{'fraction_non_alphanumeric': 0.03233684210526316, 'fraction_numerical': 0.011915789473684211, 'mean_word_length': 4.232650363516193, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mixed scalar-pseudoscalar Higgs boson production through next-to-next-to-leading order at the LHC

We study the production of a mixed scalar-pseudoscalar Higgs boson in gluon fusion at the LHC, through next-to-next-to-leading order (NNLO) in QCD. We obtain fully differential results, including the decay of the Higgs boson to two charged lepton pairs. We discuss the impact of the interference between the scalar and pseudoscalar states. We also show differential distributions for several kinematic variables whose shape is sensitive to the parity of the Higgs boson, and assess the impact of the NNLO QCD corrections on these shapes.


Introduction
Comprehensive studies of the properties of the Higgs boson will be a focus of the particle physics community in the forthcoming years. The large dataset that will be recorded by the Large Hadron Collider (LHC) during Run 3 and in its high-luminosity phase will allow precise probes of the quantum numbers and couplings of the Higgs boson. One of the most interesting properties of the Higgs boson is its parity. Initial measurements indicate that the 125 GeV Higgs boson is a scalar state J P = 0 + , while the pseudoscalar state J P = 0 − has been ruled out [1][2][3][4]. On the other hand, the possibility that the Higgs boson is an admixture of scalar and pseudoscalar states has not been excluded yet by measurements, although constraints on the parameter space for such a mixed state exist [5,6].
Such mixing between scalar and pseudoscalar states implies CP violation in the Higgs sector. This would be a signal for physics Beyond the Standard Model (BSM), and might explain the origin of CP violation within the Standard Model (SM). Moreover it could shed light on the physics of the early universe, given the large enough amount of CP violation that is required for baryogenesis. However, we also stress that, irrespective of explicit BSM scenarios, determining the behavior of the Higgs boson under parity is important as a matter of principle, in order to build a complete picture of this new particle.
From a phenomenological point of view, the possibility of producing mixed scalarpseudoscalar Higgs states through gluon fusion has been considered in Refs. [7][8][9]. These references studied a number of angular observables whose shapes are sensitive to the parity of the Higgs boson. However, Ref. [7] considers the leading-order (LO) production, while Refs. [8,9] include the next-to-leading order (NLO) QCD corrections and parton shower effects. It is well known that Higgs production is subject to large perturbative QCD corrections and that results at LO and NLO are not always reliable. Indeed, the results for the scalar Higgs boson are known at next-to-next-to-next-to-leading order (N3LO) accuracy [10], and these indicate that the perturbative expansion in α s only begins to converge at next-to-next-to-leading order (NNLO). The NNLO corrections for both scalar and pseudoscalar Higgs boson production have been known for some time [11][12][13], and are implemented in the public codes SuSHi [14,15] and its extension SuSHiMi [16]. These results, however, are limited in two ways: first, they do not include the possibility of an admixture between scalar and pseudoscalar Higgs states, including possible interference effects. Second, they only include the decays of the Higgs boson as overall branching ratios. As a result, they cannot be used to investigate the potential to determine the parity of the Higgs boson using information from its decay products, nor can they be used to compute fiducial cross sections defined by cuts on the decay products of the Higgs boson.
In this paper, we aim to bridge this gap by presenting the first NNLO QCD-accurate fully differential predictions of mixed scalar-pseudoscalar Higgs boson production, including the Higgs decay into two charged lepton pairs. In particular, we consider distributions in the angles ∆φ l 1l2 , Φ, and cos θ 1 , which describe the geometry of the decay of a spin zero particle into two charged lepton pairs and are known to be sensitive to its parity [7], and we assess the impact of the NNLO QCD corrections on these observables.
The remainder of the paper is organized as follows. In Section 2, we briefly summarize the technical details involved in this calculation, before presenting results in Section 3, and concluding in Section 4.

Technical Details
In order to describe the production of an arbitrary mix of scalar (0 + ) and pseudoscalar (0 − ) Higgs states, we make use of the Higgs Characterization model introduced in Ref. [8]. The Lagrangian describing the interaction between a spin zero particle and two heavy fermions is where ψ f is a fermionic field of flavor f , X H/A is the Higgs field, g Hf f and g Af f are the scalar and pseudoscalar Higgs couplings to the fermions respectively, and we have used the notation c α = cos(α) and s α = sin(α). The mixing between the scalar and pseudoscalar Higgs states is therefore controlled solely by the parameter α, with the pure scalar and pseudoscalar states given by c α = 1 and c α = 0, respectively. It was shown in Ref. [8] that this Lagrangian can be used to build an effective Lagrangian for the interaction of a spin-0 mixed scalar-pseudoscalar state with vector bosons below a cutoff scale Λ. In this paper, we will focus on the production of the Higgs boson through gluon fusion and its subsequent decay into two charged lepton pairs. The relevant terms in the effective Lagrangian are where the field strength tensors are defined as and the dual tensor isṼ The factors κ HV V and κ AV V (V = g, Z, γ) allow modifications of the (dimensionful) couplings g HV V and g AV V . We comment on the values of these couplings in the following section. We now discuss a technical detail concerning the renormalization of the pseudoscalar amplitudes. Neglecting the Yukawa interaction between the pseudoscalar Higgs boson and the light quarks, the interaction between the 0 − state and light QCD particles is mediated by a top quark loop only. In the heavy top quark limit, this loop can be integrated out, leading to the two effective operators [17], where we denote bare operators with O B i . The first operator describes the interaction of the pseudoscalar Higgs boson with gluons, and is present in Eq. (2.2). The second operator describes the loop-induced interaction with light quarks, and first appears at NNLO in QCD. These operators mix under renormalization as where O R i is the renormalized operator. Note that Z 21 = 0, and Z 22 is determined in such a way as to preserve the absence of higher order corrections [18] to the axial anomaly, to all orders in perturbation theory.
In the limit of massless light quarks, contributions from the squared operator O B This contribution has to be added to the two-loop amplitudes for pseudoscalar Higgs boson production [12,19,20], which we did by expressing it in terms of the leading order contribution of O R 1 2 using Eq. (2.7).
We will use Eq. (2.2) and Eq. (2.5) to describe the process pp → X H/A → e − e + µ − µ + to NNLO in QCD. We will consider the Higgs boson to be produced onshell, so that the production and decay processes factorize. We briefly discuss our implementation of these two processes below.
The production is governed by the g Hgg and g Agg terms in Eq. (2.2). To compute the NNLO corrections, we need the double-real, real-virtual and double-virtual amplitudes for the production of a scalar and a pseudoscalar Higgs boson. We take these from Refs. [19][20][21][22][23][24][25], and use the nested soft-collinear subtraction scheme [26][27][28][29][30] to extract and remove the infrared singularities associated with these contributions. We have checked our results for scalar and pseudoscalar Higgs boson production through NNLO against the program SuSHi [14, 15] and found full agreement.

Results
We now present numerical results for the production of mixed scalar-pseudoscalar Higgs states. We consider the production of a Higgs boson of mass m H = 125 GeV at the LHC operating at 13 TeV center-of-mass energy. We use a factorization and renormalization scale µ = m H /2 throughout the paper, and estimate the scale uncertainty by varying this scale by a factor of 2 in either direction [10,32]. We use NNPDF3.0 NNLO parton distribution functions [33] for all results, i.e. we compute LO and NLO cross sections with NNLO PDF's. The Higgs vacuum expectation value is taken to be We choose the mass of the Z boson to be m Z = 91.1876 GeV and its width to be Γ Z = 2.4952 GeV. We use the weak coupling the relevant interaction terms are G µν G µν X H/A and G µνG µν X H/A , which are controlled by five parameters: the dimensionful couplings g Hgg and g Agg , the dimensionless parameters κ Hgg and κ Agg which allow the modifications of the couplings, and the scalar-pseudoscalar mixing parameter c α . The dimensionful couplings have the values We set κ Hgg = κ Agg = 1, and present results for three representative values of the mixing parameter, c α = {1, 0, 1/2}, which correspond to a pure scalar, pure pseudoscalar, and an equal scalar-pseudoscalar admixture, respectively. The cross sections at LO, NLO and NNLO in QCD are shown in Table 1. As is well known for Higgs boson production, the NLO and NNLO corrections are large, with NLO and NNLO k-factors of approximately 2.3 and 1.25, respectively, while the scale uncertainties at LO and NLO underestimate the missing higher order corrections. We also note that the cross sections in the scalar case are smaller than those in the pseudoscalar case by a factor of about 0.44, due to the coupling of gluons to a scalar Higgs boson being suppressed by a factor g 2 Hgg /g 2 Agg = 4/9 relative to the pseudoscalar Higgs boson. We note that at all three orders, the result for c α = 1/2 is the arithmetic average of the results for c α = 0 and c α = 1. This implies that there are no interference contributions ∼ c α s α . This is immediately obvious at LO and NLO as the relevant matrix elements do not include such interference terms. However, both the double-real and the real-virtual matrix elements which enter the NNLO calculation contain such terms ∼ c α s α , which only vanish upon integration over the phase space. We have confirmed this observation for the case of LO gluon fusion Higgs production in association with two jets (which corresponds to the fully resolved double-real contributions to the NNLO corrections), both using our own code and using MadGraph [8,31]. We therefore conclude that, if one considers the production of a Higgs boson and neglects its decay, the results up to NNLO for an arbitrary value of c α may be obtained by simply rescaling the scalar and pseudoscalar results We now turn to the case of the Higgs boson decaying into two charged lepton pairs pp → X H/A → Z/γ * Z/γ * → e − e + µ − µ + , with the Higgs boson being onshell. We impose  minimal kinematic cuts on the final state leptons, inspired by a recent ATLAS analysis [6]. We require all leptons to have transverse momentum p T,l > 15 GeV and pseudorapidity |η l | < 2.5. Moreover, we require the invariant mass of each lepton pair to be in a window around the Z mass peak, 50 GeV < m l − l + < 106 GeV. This last cut implies that the contribution of the offshell photons is negligible; for simplicity, we set κ Hγγ = κ Aγγ = κ HZγ = κ AZγ = 0 in Eq. (2.2). Moreover, since we are interested in CP -violation in the Higgs sector, we will set κ H∂Z = κ H∂γ = 0 as these derivative terms do not have a pseudoscalar counterpart. Therefore, we will only consider the terms in Eq. (2.2) which are governed by κ HZZ and κ AZZ , together with the SM term with coupling g HZZ = 2m 2 Z /v, and the production terms with couplings g Hgg and g Agg which we have already discussed. We then set κ SM = κ HZZ = 1, κ AZZ = 1, Λ = 1 TeV, and consider the three benchmark scenarios with values of c α = {0, 1, 1/2}. We also consider the case c α = 0.6 together with κ AZZ = 20 in order to illustrate the sensitivity of the shape information of certain angular observables to the parity of the Higgs boson. All other choices of parameters, scales, and the pdf set are the same as for the undecayed Higgs boson, described above.
We show the fiducial cross sections for this setup in Table 2. We first note that, in contrast to the results presented in Table 1, the cross sections for the pure pseudoscalar case are smaller than those for the pure scalar case by three orders of magnitude. This can be understood by looking at Eq. (2.2). The (scalar) SM interaction between the Higgs boson and the Z boson pair has a coupling given by g HZZ = 2m 2 Z /v, as mentioned previously. The pseudoscalar interaction Z µνZ µν X H/A leads to a factor f ({p})/Λ in the decay amplitude, where f ({p}) is a kinematic factor with dimension of mass-squared. The value of f ({p}) is generally smaller than m 2 Z , and moreover Λ > v, leading to the pseudoscalar decay to Z bosons being suppressed relative to the SM scalar decay by several orders of magnitude.
From Table 2 one can also see that the fiducial cross section for c α = 1/2 is no longer given by an arithmetic average of the fiducial cross sections for c α = 1 and c α = 0. This is clear from the fact that the degree of scalar-pseudoscalar mixing is controlled by c α both in the production as well as in the decay. This implies that terms ∼ c α s α do appearmost notably, from the combination of the pseudoscalar interaction in production and SM interaction in the decay. This means that, in general, a simple reweighting formula like  Figure 1: Normalized distribution at NNLO accuracy and differential k-factors dk NLO and dk NNLO for the angle Φ. See text for further details.
Eq. (3.2) cannot be used anymore due to the interplay between the production and decay of the Higgs boson. We note that the scale uncertainties and the impact of the NLO and NNLO corrections are similar to those for the undecayed results (see Table 1). Moreover, both the scale uncertainties and the effects of the NLO and NNLO corrections are the same for all four values of c α in Table 2. This, together with the fact that for the SM Higgs, the N3LO corrections lie within the NNLO scale uncertainty bands [10,32], lead us to conclude that NNLO is the first order at which the results for any value of c α are reliable.
It is clear that, for this choice of parameters, the cross sections provide enough information to discriminate between the pure scalar and pure pseudoscalar scenarios. If we compare the results for the c α = 1 and c α = 1/2 cases, we see that they are compatible within the scale uncertainties at LO and NLO. The NNLO corrections, however, lead to reduced scale uncertainties, and the results for these two values of c α are no longer compatible at this order. This emphasizes the need for higher order corrections in determining the properties of the Higgs boson. On the other hand, the results for c α = 1 and c α = 0.6 (with κ AZZ = 20) are compatible within the scale uncertainties at LO, NLO and NNLO, meaning that one cannot differentiate between these two cases based on the rates alone, and additional information from the shape of kinematic distributions is required.
We now show differential distributions for three observables Φ, cos θ 1 and ∆φ l 1l2 . The first two observables were proposed in Ref. [7], where they have been shown to be particularly sensitive to the spin and parity of the Higgs boson. The observable Φ is the azimuthal angle between the planes constructed by the respective decay products of the two Z bosons in the rest frame of the Higgs boson, while cos θ 1 is the polar angle of the decay products of the first Z boson in its own rest frame. Identifying the Z boson that decays to electrons as Z 1 and the one that decays into muons as Z 2 , the angle Φ is defined as [7,34] where q i is the three-momentum of Z i , and q i1 (2) is the three-momentum of the lepton (antilepton) resulting from the decay of Z i . Here, all three-momenta are defined in the rest frame of the Higgs boson. The observable cos θ 1 is defined as [7,34] cos where q 2 and q 11 are now defined in the rest frame of a scattering axis whose direction is that of Z 1 in the Higgs rest frame. We begin by showing the Φ distribution in Fig 1. The upper pane shows the NNLO results for c α = 1.0, c α = 0.6 (with κ AZZ = 20) and c α = 0.0, normalized to their respective NNLO cross sections. The distribution obtained for the scale µ = m H /2 is depicted as the thicker central line, while the band around it is the envelope from varying the scale by factors of two and 1/2 around this central value. As expected from Ref. [7], there is a notable shape difference between the three values of c α . The value Φ = 0 corresponds to a maximum for the scalar distribution and a minimum for the pseudoscalar, while the scalar has minima and the pseudoscalar has maxima at Φ = ±π/2. The minima and maxima of the c α = 0.6 distribution are shifted relative to the pure scalar and pure pseudoscalar cases, giving the c α = 0.6 case a distinct shape. This shift originates from the interference between the 0 + and 0 − production and decay contributions. We recall from Table 2 that one cannot tell the c α = 1 and c α = 0.6 scenarios apart based on overall rates alone, as the cross sections for these values of c α lie within each others' uncertainty bands, even once NNLO corrections are taken into account. The shape of the Φ distribution may provide a means to distinguish between these two scenarios. Of course, the choice of parameters c α = 0.6 and κ AZZ = 20 is somewhat contrived, but it does illustrate the importance of shape information in determining the parity of the Higgs boson in a large EFT parameter space.
It follows from this discussion that it is important to have reliable predictions for the shapes of distributions, meaning that the impact of the NNLO corrections needs to be known. In the lower two panes of Fig. 1, we show the differential NLO and NNLO kfactors, defined as In order to simplify the plot, the k-factors are only shown for the central scale choice. We observe that the value of Φ has a mild effect on the NLO k-factor, which peaks at Φ = ±π/2 and is smallest at Φ = 0, and therefore tends to reduce the difference between the scalar and pseudoscalar distributions. This behavior of the NLO k-factor appears to be an effect of the real radiation in conjunction with the cuts; without the cuts, we observe the k-factor to be perfectly flat, as expected. We conclude that the real emission moves the final state inside or outside of the fiducial volume defined by the kinematic cuts, in a way which is similar for the scalar and the pseudoscalar cases. The NNLO k-factor is quite flat, and amounts to a simple rescaling of the NLO results by a factor of approximately 1.25, implying that the additional radiation present at NNLO does not dramatically change the acceptance rates for this fiducial volume. We also note that the NLO k-factor is slightly larger for the pseudoscalar Higgs than for the scalar one, while the NNLO k-factors are almost identical for all three values of c α . We now turn to the cos θ 1 distribution, shown in Fig. 2. As for the Φ distribution, the shape of this distribution is significantly different for the pure scalar and pure pseudoscalar scenarios, with the former having a maximum and the latter a minimum at cos θ 1 = 0. The shape difference between the pure scalar and the c α = 0.6 distributions is much milder, implying that this observable is less sensitive to the parity of the Higgs than Φ is, given our setup. The NLO corrections have a mild dependence on the value of cos θ 1 and appear to have a slightly larger impact at low values of this angle. The NNLO k-factor is flat, and the k-factors at both NLO and NNLO are the same for all values of c α .
Finally, we show the distribution in the opening angle ∆φ l 1l2 in the lab frame between the e − and the µ + leptons in Fig. 3. Unlike the Φ and cos θ 1 observables considered previously, this observable can be measured even if the final state lepton configurationand thus the Higgs boson rest frame -cannot be fully reconstructed. As such, it is also an  interesting proxy for the W + W − decay channel of the Higgs boson. We see a noticeable shape difference between the pure scalar and pure pseudoscalar cases, but the difference between the c α = 0.6 and the pure scalar or pure pseudoscalar cases is much milder and is covered by the scale uncertainty bands of the distributions. Therefore, as expected, this observable has a lower sensitivity to the parity of the Higgs boson than either Φ or cos θ 1 .
Looking at the NLO k-factor, we see that the NLO corrections enhance the distribution at small angles. This is due to the additional radiated parton, and the effect is made more pronounced by the kinematic cuts that we impose. On the other hand, the k-factor at NNLO is relatively flat, implying that the presence of a second radiated parton has less of an impact, as we saw for the Φ and cos θ 1 distributions. Again, the value of c α does not seem to affect the differential k-factors at NLO or NNLO.

Conclusions
We have presented the first fully differential results for the production of a mixed scalarpseudoscalar Higgs boson X H/A through gluon fusion to NNLO accuracy in QCD. We made use of an effective Lagrangian to parametrize the coupling of the mixed state to gluons as well as its decay into Z bosons. In particular, the mixing between the scalar and pseudoscalar Higgs states is controlled by a single parameter c α . This allows us to make precise predictions for a generic observable at the LHC for an arbitrary admixture of scalar and pseudoscalar Higgs states.
For the production of a stable X H/A boson we observe that the cross section for an arbitrary mixing angle can be obtained by an appropriate reweighting of the cross sections for scalar and pseudoscalar Higgs production. This is true even at NNLO in QCD, where interference effects between the scalar and pseudoscalar amplitudes occur, meaning that such interference vanishes upon integration over the full phase space.
Furthermore, we considered the subsequent decay X H/A → ZZ * → e − e + µ − µ + for onshell intermediate X H/A . We observe the fiducial cross section for the pure pseudoscalar Higgs boson to be smaller than that for the pure scalar one by several orders of magnitude, as a result of the pseudoscalar decay amplitudes being suppressed relative to the scalar ones. This implies that the scalar and pseudoscalar bosons may be distinguished through the rates alone. On the other hand, for certain values of the mixing parameter c α and of the coupling of the Z bosons to the Higgs, the cross sections for the pure scalar and the mixed state are comparable. In these situations, angular observables are known to provide additional discriminating power. We considered three differential distributions, and found that for the setup considered here, the angle Φ showed the most noticeable sensitivity to the parity of the Higgs. We observed that, while the NLO corrections showed some dependence on the observable for these three distributions, the k-factors for the NNLO corrections were relatively flat. Furthermore, the corrections are largely independent of the value of c α , implying that the differential corrections to pure scalar production are a good approximation for the differential corrections for any value of c α . However, one should be cautious in drawing this conclusion, as the situation may change if different kinematic cuts are applied, or if different values of the parameters in the effective Lagrangian were chosen.