QCD corrections to the Golden decay channel of the Higgs boson

Future colliders aim to provide highly precise experimental measurements of the properties of the Higgs boson. In order to benefit from these precision machines, theoretical errors in the Higgs sector observables have to match at least the experimental uncertainties. The theoretical uncertainties in the Higgs sector observables can be reduced by including missing higher-order terms in their perturbative calculations. In this direction, we compute the mixed QCD-electroweak corrections at ${\mathcal O}(\alpha \alpha_s)$ to the Higgs decay into four charged leptons by considering the golden decay channel, $ H \to e^+e^-\mu^+\mu^-$. Due to color conservation, these corrections receive contributions only from the two-loop virtual diagrams. In the complex mass scheme, we find that the mixed QCD-electroweak corrections to the partial decay width, relative to the leading order predictions, are positive and about $0.27\%$ for $\alpha_s(M_Z)$. Relative to the next-to-leading order electroweak corrections, the mixed QCD-electroweak corrections are found to be approximately $18\%$ for $\alpha_s(M_Z)$. With respect to the leading order, we observe a flat effect of the mixed QCD-electroweak corrections on the invariant mass distribution of the lepton pairs with fixed QCD coupling. The $\phi$ distribution, due to the mixed QCD-electroweak corrections, follows a $(1-\cos \phi)$ dependence.


Introduction
With the groundbreaking discovery of the Higgs boson by the CMS and ATLAS [1,2] collaborations at the LHC in 2012, particle physics has entered the realm of precision studies.The Higgs boson is one of its kind in the Standard Model (SM); therefore, precision measurements in the Higgs sector provide an opportunity to look for the physics beyond the Standard Model (BSM).All the present and future colliders, such as the Large Hadron Collider (LHC) [3], High-Luminosity LHC (HL-LHC) [4], Future Circular Collider (FCCee) [5], Circular Electron Positron Collider (CEPC) [6,7], and the International Linear Collider (ILC) [8] aim to explore the uncharted territory of the fundamental interactions by measuring various Higgs properties with higher statistics.On the theory side, highly precise predictions for the Higgs production and decay channels are needed for a fairer comparison with future experimental data.
Among the five prominent decay modes of the Higgs, a rare but important one is the decay of the Higgs boson into four charged leptons, also known as the "Golden decay channel".This decay channel played a significant role during the Higgs discovery in 2012 as it provided a particularly clean signature around the Higgs mass (∼ 125 GeV) in the invariant mass spectrum of the final state leptons.Furthermore, kinematic distributions of the final state leptons for this decay mode not only allow for precision mass measurements of the Higgs boson but also serve as a powerful tool to study its spin and CP properties [9][10][11][12][13].The data collected in the off-shell production and decay of the Higgs to four leptons via the Z-boson pairs can constrain its total decay width [14][15][16].Thus, improved predictions for the golden decay mode are paramount to having a better understanding of the Higgs properties.In this direction, several works have been reported in the literature.
The exact one-loop QED corrections of O(α) to the Higgs decay into four leptons with the off-shell Z-bosons have been evaluated in [17,18].The complete one-loop electroweak corrections for the leptonic, semi-leptonic, and hadronic final states, and the one-loop QCD corrections for the semi-leptonic and hadronic final states to the decay H → W W/ZZ → 4f have already been evaluated in [19,20] and are encoded in a Monte Carlo (MC) code Prophecy4f [21,22].The O(α) corrections, reported for the case of four charged leptons in the final state, are of the order of 2-4% for moderate Higgs masses (M H ≤ 200 GeV) and increase with the growing Higgs mass, reaching up to 14%.In addition to that, the one-loop electroweak and QCD corrections to the Higgs decay into four fermions in the context of a simple extension of the SM have also been studied and are implemented in the code Prophecy4f [23][24][25].The next-to-leading order (NLO) electroweak corrections to the Higgs decay into charged leptonic final states H → Z ( * ) Z ( * ) → 4ℓ with 4ℓ = 4e, 4µ, 2e2µ matched with QED Parton Shower (PS), have also been calculated, for which the results are available in a public event generator, Hto4l [26].
In the present work, we compute the QCD corrections to the decay H → e + e − µ + µ − on top of the electroweak corrections that mainly receive contributions from the two-loop diagrams, appearing at O(αα s ) in the perturbation theory.These mixed QCD-electroweak corrections are expected to be small because of the two-loop effect, as compared to the NLO electroweak corrections.However, these are essential to provide precise predictions for the Higgs sector observables at the LHC and future colliders and to test the validity of the perturbative QFT calculations.Our motive is to quantify these corrections, simulate decay events, and provide improved numerical predictions for the partial decay width of H → e + e − µ + µ − , and relevant kinematic distributions with an accuracy of O(αα s ).The twoloop diagrams contributing to the amplitude at O(αα s ) are very similar to those appearing in the processes e + e − → ZH and e + e − → µ + µ − H, for which the mixed QCD-electroweak corrections of O(αα s ) have been evaluated in [27][28][29].In this work, numerical calculation of the two-loop amplitude is performed systematically using our in-house codes, and finally, to provide improved predictions for the partial decay width, phase-space integration over the final-state leptons is performed by interfacing our codes with the publicly available code Hto4l [26].
The rest of the paper is organized as follows: In section 2, we classify the Feynman diagrams that contribute to O(αα s ).The organization of the matrix elements in terms of form-factors and their divergence structure is discussed in section 3.In section 4, the UV renormalization of the two-loop matrix elements, along with the opted renormalization scheme, is described.The numerical implementation of the two-loop matrix elements for event generation and the checks performed on them are given in section 5 followed by our numerical results in section 6.Finally, we draw conclusions from our work in section 7.
2 QCD correction to H → e + e − µ + µ − In the SM, the leading order (LO) contribution to the on-shell decay of the Higgs boson into four charged leptons (H → e + e − µ + µ − ), mediated by a pair of Z-bosons, comes from a tree-level diagram, shown in Fig. 1.Due to energy conservation, at least one of the two Z-bosons has to be off-shell, depicted by Z * .In this work, we consider a more general case by treating both the Z-bosons off-shell.We choose the following momenta assignments for the particles in the decay.
where momentum conservation requires q = (p 1 + p 2 ), p 1 = (p 3 + p 4 ) and p 2 = (p 5 + p 6 ).In our calculation, we are neglecting the masses of the final-state leptons.Various scalar products of interest are given by, where m H is the Higgs mass.As the decay width is proportional to the squared amplitude, the amplitude for H → Z ( * ) Z ( * ) → e + e − µ + µ − in the perturbative expansion up to twoloop order can be written as, Here M 0 , M 1 , and M 2 are the LO, one-loop, and two-loop amplitudes.In the SM, M 1 receives contribution only from the electroweak (EW) sector, as the particles involved at the LO are color neutral.However, at the two-loop level, both the EW and QCD sectors can contribute.Therefore, there can be contributions of O(α 2 ) and O(αα s ) at the twoloop level, but one can neglect the contributions of O(α 2 ) due to the smallness of the EW coupling α in comparison to the strong coupling α s .Thus, we only consider the mixed QCD-electroweak corrections of O(αα s ) at the two-loop level and focus on the evaluation of M ααs 2 .
In going beyond the LO, O(α) amplitude receives contributions from several one-loop diagrams mediated by the weak bosons and quarks.However, only quarks are susceptible to couple with the gluons and take part in the QCD corrections; therefore, at O(αα s ), only the diagrams with quark loop along with the gluon dressing will contribute.The contribution to the amplitude at O(αα s ) can be divided into three categories as follows: where δM ααs HV 1 V 2 consists of corrections coming from the HV 1 V 2 vertex, δM ααs S.E.contains corrections due to the self-energy insertions on the vector-boson legs, and δM ααs Zℓ l appears due to O(αα s ) counter-term for the Zℓ l vertex, respectively.These contributions are described below.

HV 1 V 2 vertex corrections
At the two-loop level, in addition to the decay of the Higgs into e + e − µ + µ − through the Z ( * ) Z ( * ) channel, we also have to consider the contributions coming from the Z ( * ) γ ( * ) , and γ ( * ) γ ( * ) channels.Thus, we consider the most general amplitude for the HV 1 V 2 vertex corrections denoted by δM ααs , where both V 1 and V 2 are taken off-shell, and write, δM ααs HV 1 V 2 = δM ααs HZZ + δM ααs HZγ + δM ααs HγZ + δM ααs Hγγ . (2.5) In the above, the HZγ and HγZ contributions are written explicitly to take care of the fact that e + e − can come from either Z or γ.As the leptonic decay of the vector bosons is not affected by the QCD corrections to the HV 1 V 2 vertex, we can decompose the amplitude for the HV 1 V 2 vertex corrections as where, M µν is the two-loop amplitude for H(q) → V 1 (p 1 )V 2 (p 2 ) decay; J µ (p 1 ) and J ν (p 2 ) are the fermionic currents corresponding to V 1 → e + e − and V 2 → µ + µ − respectively.Since coupling of the top quark with the Higgs is the largest among all the quark flavors, we neglect the contributions from diagrams with any quark other than the top quark in the loop.There are in total 48 two-loop triangle diagrams contributing to the H(q) → V 1 (p 1 )V 2 (p 2 ) decay at O(αα s ), out of which some representative diagrams are shown in Fig. 2. The diagrams are generated with the help of the QGRAF [30] package.The amplitude for each of these diagrams is organized in FORM [31,32] and is manipulated using Mathematica.
Instead of using the conventional methods, the projector technique has opted to perform the amplitude evaluation more systematically.In this technique, the amplitude for H → V 1 V 2 using the Lorentz covariance can be expressed as where A, B, C, D, E and F are scalar functions called form-factors and ϵ µνp 1 p 2 = ϵ µνρσ p 1ρ p 2σ .The form-factors can be obtained by applying suitable projectors P i µν (i = A, B, C, D, E, F ) on the amplitude M µν .These form-factors are, in general, functions of the Mandelstam variables present in the problem under consideration.The two-loop form-factors contributing to the amplitude at O(αα s ) are discussed in section 3.
with the top quark running in the loop.Diagrams with the reversed direction of the fermionic current are not shown.

Self-energy corrections
The self-energy part of the two-loop amplitude, denoted by δM ααs S.E., receives contributions from O(αα s ) corrections to the ZZ and mixed Zγ self-energies.This part of the amplitude receives contributions from the top quark as well as from the light quarks.In total, 72 self-energy diagrams contribute to the δM ααs S.E., out of which some are shown in Fig. 3.In order to calculate δM ααs S.E., we need O(αα s ) expressions of the gauge boson self-energies.Their analytical expressions are available in [33,34].u, d, c, s and b) and the top quark contribute.The light quarks are taken to be massless.

Zℓ l vertex correction
The amplitude at O(αα s ) also receives contribution from the diagrams shown in Fig. 4 involving a Zℓ l vertex counterterm.The contribution to the vertex counterterm at O(αα s ) comes from the self-energies of the vector bosons.It depends on the wave function renormalization constants of the vector bosons, the charge renormalization constant, and the renormalization constant for the weak mixing angle [34].The pure electroweak nature of the Zℓ l vertex does not allow any kind of QCD corrections of O(αα s ); therefore, the counterterm for the Zℓ l vertex is non-divergent in nature.One must be careful about the renormalization scheme in which the vertex correction is computed.

Two-loop form-factors and divergences
The form-factors for the HV 1 V 2 vertex corrections given in Eq. 2.7 can be obtained by applying the projectors given in Appendix A on the bare two-loop amplitude of all the contributing diagrams shown in Fig. 2. In order to calculate the trace involving the γ 5 matrix, the prescription given in [35] is used.Furry's theorem forbids the appearance of a single γ 5 in the trace over a closed fermionic loop due to charge invariance, which gives C = 0 on adding the contributions from all the triangle diagrams together.Furthermore, as the current conservation is associated with massless leptons in the final state, only the form-factors A and B contribute at the squared amplitude level.These form-factors can be written as linear combinations of scalar two-loop integrals of the type Here, k 1 and k 2 are the loop momenta, d denotes the space-time dimensions, γ E is the Euler-Mascheroni constant, µ represents an arbitrary scale introduced to maintain the canonical dimension of the integral, {ν i } = {ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 ν 7 } are powers of the inverse propagators P i , and ν = 7 i=1 ν i .The inverse propagators P i are given by The set of two-loop integrals in these form-factors is reduced to a minimal set of integrals called master integrals (MIs) using the integration-by-parts (IBP) [36,37] and Lorentz-invariance (LI) [38] identities with the programs LiteRed [39,40] combined with Mint [41] and FIRE [42][43][44].We find a total of 41 master integrals after the IBP reduction.The choice of master integrals is not unique.The basis set ⃗ I of 41 master integrals we obtain is listed below.The involvement of two-loop integrals makes the evaluation of these form-factors highly non-trivial.The master integrals can be evaluated numerically with the help of publicly available codes such as pySecDec [45] and AMFlow [46] etc.The analytical results for these 41 master integrals in the canonical basis are now also available in terms of iterated integrals [47].In our work, for an efficient evaluation, keeping the required accuracy of the MIs in mind, we use an in-house code for the numerical evaluation of all the two-loop master integrals involved, using the sector-decomposition method given in [48][49][50].The results obtained for the MIs using the numerical integration are in good agreement with the analytical results.
Due to the presence of loop integrals and massless particles in the loop, the twoloop form-factors develop both the Ultraviolet (UV) and Infrared (IR) divergences.We regularize these divergences in dimensional regularization by taking d = 4 − 2ϵ.After regularization, the divergences are encoded in the two-loop master integrals as poles in ϵ, 1 ϵ 4 being the highest order of pole that can appear.The 1 ϵ 3 and 1 ϵ 4 poles are exclusively due to the IR singularities, while 1 ϵ and 1 ϵ 2 poles can be due to both IR and UV singularities.According to the KLN (Kinoshita-Lee-Nauenberg) theorem [51][52][53][54], the IR singularities eventually get cancelled against the real emission Feynman diagrams to give IR safe observables.
In our case, the possible real emission Feynman diagrams involve the emission of a gluon from the closed quark loop, as shown in Fig. 5. Due to a closed fermionic loop, the amplitudes of these diagrams are proportional to the trace over T a , the generator of the SU (3) gauge group.Since tr(T a ) = 0, diagrams for the real corrections give zero.Moreover, these real emission diagrams can contribute only at the amplitude-squared level, which is an O(α 2 α s ) effect with respect to the LO.Therefore, the real emission diagrams do not contribute at O(αα s ).The absence of real corrections thus demands the cancellation of virtual IR divergences among the contributing diagrams, and the two-loop amplitude is expected to be free from any IR divergences in our case.Since the form-factors A and B for the HV 1 V 2 vertex are independent, at the two-loop, they should be separately free from 1 ϵ 3 and 1 ϵ 4 poles.This fact will provide one of the important checks on our calculation.Furthermore, since the form-factor B is zero at the tree-level, and the first non-zero contribution to it arises at the one-loop, we expect that the two-loop form-factor B does not have 1 ϵ 2 UV pole dependence.This can serve as another consistency check on our calculation.

UV renormalization and the complex-mass scheme
In order to remove the UV divergences from the matrix elements contributing at O(αα s ), the standard on-shell renormalization scheme is used to renormalize all the fields and masses involved.The renormalization process involves evaluating the required counterterm (CT) diagrams and then adding the CT amplitude back to the UV divergent amplitude to obtain finite results after the cancellation of all the divergences.In total, there are 48 one-loop triangle, 96 one-loop self-energy, and 8 tree-level CT diagrams contributing to the amplitude at O(αα s ).The representative CT diagrams are shown in Figs. 6, 7 and 8.
As shown in Fig. 6, the triangle counterterm diagrams mainly involve the V t t, Ht t vertex counterterms and counterterms for the top-quark mass and wave function.However, the renormalization of the quark wave function is related to the vertex renormalization.As a result, the vertex counterterms involving quarks completely cancel out the contributions from the quark wave function counterterms.It is not surprising, as in our calculation, the total amplitude does not depend on the quark wave function.Thus, the complete cancellation of the quark wave function counterterm provides another important check on our calculation.Therefore, at the one-loop level, we only need to evaluate diagrams with the top-quark mass counterterm insertions.On the other hand, for the evaluation of the self-energy and tree-level counterterm diagrams, one needs O(αα s ) expressions for the renormalization constants δZ e , δZ ZZ , δZ H , δM 2 Z , δM 2 W , δZ γZ and δZ Zγ , which we have deduced from the self-energies of the Higgs and gauge bosons given in [33,55].
The diagrams for the process under consideration involve unstable Z-bosons in the propagators.Therefore, one needs to introduce the finite Z-width in the propagators to ensure the stability of the perturbative calculation at the Z-pole.However, this incorporation of the finite Z-width can lead to several problems, such as the violation of gauge invariance due to the mixing of different perturbative orders [56].These problems are tackled via the adoption of the "complex-mass scheme" (CMS) [57][58][59][60].In this scheme, we analytically continue the masses of the weak gauge bosons to the complex plane as    where, V = W, Z and Γ V is the corresponding constant decay width.This scheme also makes the weak mixing angle complex as cos . Therefore, the CMS scheme is just a generalization of the on-shell scheme with complex renormalization constants.There is no need to introduce a finite width for the top quark running in the loop, as the Higgs mass is taken to be 125 GeV, which is below the t t threshold.
After combining the CT amplitude with the UV divergent amplitude in the complex mass scheme, we get the finite results for the amplitudes of the HV 1 V 2 vertex and the self-energy corrections.

Numerical implementation and checks
We combine the UV finite amplitudes for the HV 1 V 2 vertex and self-energy corrections with the fermionic currents to get the two-loop amplitude for H → e + e − µ + µ − .Adding the finite Zℓ l contribution to it, we get the matrix element M ααs 2 given in Eq. 2.4.In the perturbative expansion of the amplitude-squared up to two-loop, this matrix element interferes with the LO amplitude as ). (5.1) The interference term is organized in terms of the two-loop form-factors using the symbolic manipulation program FORM, and a FORTRAN output is obtained for the numerical evaluation.
To obtain the partial decay width, we need to perform the phase-space integration over the final state leptons.This is done using the publicly available Hto4l [26] code, which is a Monte Carlo code that generates events for the process H → 4ℓ (ℓ = e, µ).To achieve a good accuracy on the observables of interest, we need to sample a huge number of phase-space points.However, calculating the form-factors for every phase point is very time-consuming.To manage this issue, we have prepared a two-dimensional grid for the form-factors A and B using pre-specified phase-space points (p 2 1 and p 2 2 values to be more precise) by numerically evaluating 41 MIs for input parameters given in section 6.The grid is prepared with an accuracy of O(10 −3 ).The grid is then used to estimate the formfactors at random phase-space points with the help of a linear interpolation code developed in-house.It is important to note that any change in the input parameter set would require a new grid.We interface the squared matrix elements, the grid of the form-factors, and the interpolation code with the Hto4l code to perform the phase-space integration.This allows us to obtain the partial decay width and kinematical distributions for the final state leptons at O(αα s ).In order to prove the reliability of our implementation, we have performed the following checks: 1.To good numerical accuracy, we find that the 1 ϵ 4 and 1 ϵ 3 poles cancel in both the form-factors A and B. In the form-factor B, the 1 ϵ 2 pole also vanishes.The UV poles in the form-factors A and B cancel after adding the CTs, and the result does not depend on the choice of the scale µ in the dimensional regularization.These checks have been performed for several phase-space points.
2. By taking the gluon propagator in the R ξ gauge, we find that the two-loop formfactors, and consequently the two-loop amplitude, are gauge-parameter independent.Additionally, the Ward identity for the HV 1 V 2 vertex demands A+p 1 .p 2 B+p 2 1 D = 0, and we have verified this relation for the two-loop amplitude at O(αα s ).
3. As mentioned earlier, the two-loop diagrams for H → e + e − µ + µ − are closely related to the ones appearing in the production process e + e − → ZH.In Ref. [28], analytical expressions for the contributing form-factors are given up to order m 0 t after series expanding them in1 mt .In order to check the accuracy of the grid prepared for the form-factors, we produced the grid for e + e − → ZH, taking a very large value of the top-quark mass (m t ).Further, we matched the numerical values of the form-factors from the grid with those given in [28].We found an excellent agreement between the two for different values of the center-of-mass energies.
4. The correctness of our numerical implementation is checked via reproducing the results for the mixed QCD-electroweak corrections for the e + e − → ZH process given in Ref. [27] in the G µ and α(0) schemes.We performed this check by implementing our calculation in MadGraph [61], and we found that the calculated corrections matched the available results in both the schemes with a relative error of less than 1%.

Numerical results
In this section, we will present the numerical results for our calculation obtained by its implementation in the Hto4l code.We work in the G µ scheme and use the following set of input parameters, In the G µ scheme, α is given by which is scale-independent in the case of on-shell renormalization.From the above list of input parameters, the transformation of the on-shell values of the masses and widths (M OS V , Γ OS V ) of the W and Z bosons to the corresponding values in the CMS, denoted by M V and Γ V , is provided in Appendix B. It is worth mentioning that the difference between using (M OS V , Γ OS V ) or (M V , Γ V ) is hardly visible in the numerical results presented in this section.The results for the mixed QCD-electroweak corrections are obtained using both fixed and running α s (Q).We take Q = M Z to obtain results with a fixed value of the QCD coupling.Since the mixed QCD-electroweak corrections are leading order in α s , we employ the one-loop result for the running of α s (Q).For the choice of the running scale, we have used p T of the lepton and the invariant mass of the ℓ + ℓ − pair (M ℓ + ℓ − ).
We find that the mixed QCD-electroweak correction to the partial decay width is around 0.27% of the LO contribution for α s at Q = M Z 1 .With the running coupling, it becomes 0.30% for Q = M ℓ + ℓ − and 0.35% for Q = p T (ℓ − ).To draw a comparison, we note that the inclusive two-loop QCD corrections in H → Zγ decay have been found around 0.22% of the LO [62][63][64].The two-loop QCD corrections in the H → γγ decay lie in the range of 1-2% for the intermediate Higgs mass below the t t threshold [65].With respect to the NLO EW correction, the mixed QCD-electroweak correction amounts to 18% for the fixed and 21% (24%) for the running QCD coupling with It is well known that the higher-order corrections are sensitive to the kinematics of the events.In Fig. 9, we investigate the impact of the two-loop corrections with respect to the LO predictions and NLO EW corrections on the invariant mass distribution of the final state lepton pair.We define δ i = Γ two-loop /Γ i (i = LO, NLO) to indicate the relative correction with respect to the LO contribution and NLO EW correction.For the fixed α s , the mixed QCD-electroweak corrections relative to the LO are roughly the same in all the bins and are of the order of 0.27%, as seen in Fig. 9(a).This suggests that the nature of the events for the LO and two-loop corrections is kinematically similar.The two-loop corrections for the running α s differ in each bin and are higher in the lower mass bins.It can go beyond 0.35% in the lower bins, depending on the choice of the running scale.For the invariant mass above the Z pole, the corrections for the running α s with Q = M ℓ + ℓ − become slightly smaller than the corrections for the fixed α s .This behavior is dictated by the one-loop running of α s .
In Fig. 9(b), we plot the two-loop corrections with respect to the NLO EW corrections in the upper panel.Since the electroweak corrections are also sensitive to the kinematics of the events, we note that in certain bins, between 20 GeV and 40 GeV, δ NLO becomes quite large independent of the choice of the scale Q.The relative effect of the NLO EW corrections with respect to the LO contributions in each bin, as shown in the lower panel, can be referred to understand the features of the distribution in the upper panel of Fig. 9(b).Significantly large values of δ NLO in certain bins are simply due to the fact that the NLO EW corrections are negligible in those bins.This observation can be useful for the bin-wise analysis of the data.The two-loop corrections appear flat in the bins between 40 GeV and 80 GeV.However, a closer look reveals that δ NLO decreases in higher mass bins due to the larger NLO EW corrections in those bins.We have shown this in Fig. 9(c).The inverse behavior of the distributions in the upper and lower panels can be attributed to the kinematic similarity between the two-loop events and the LO events noted earlier.
Apart from the invariant mass distribution, angular distributions are also effective in studying the Higgs properties.Therefore, we study the effect of the mixed QCD-electroweak corrections in the ϕ distribution, which is one of the most sensitive observables for the BSM studies.It is defined as the angle between the decay planes of the intermediate Z-bosons in the rest frame of the Higgs boson.This angle ϕ is the main observable for spin-parity assignment of the Higgs boson [66][67][68][69][70][71][72].We have plotted the ϕ dependence of the relative corrections δ i in Fig. 10.
In Fig. 10(a), we see that the mixed QCD-electroweak corrections relative to the LO do not exhibit a flat behavior in the ϕ distribution for the fixed α s .This is in contrast to what we see in the case of the invariant mass distribution in Fig. 9(a).We observe a (1 − cos ϕ) dependence in the shape of the ϕ-distribution due to the mixed QCD-electroweak corrections.Also, the shape is independent of the scale choice for α s .The stated dependence is different from the LO behavior, which follows a cos 2 ϕ dependence [72,73].The difference can be attributed to the change in the effective HZZ coupling due to the two-loop corrections.The two-loop corrections with respect to the LO are insignificant at the edges for α s , respectively.The dotted straight lines mark the results at the inclusive level.
and large in the central region.The relative correction δ LO peaks at ϕ = π.It is 0.49% for the fixed and 0.56% (0.64%) for the running QCD coupling with Compared to the fixed α s , relative corrections are higher for the running α s across all the bins.
)) for α s , respectively.The dotted straight lines mark the results at the inclusive level.
In the upper panel of Fig. 10(b), we have shown the relative effect of the two-loop corrections with respect to the NLO EW corrections.In the lower panel of the figure, the NLO EW corrections with respect to the LO are displayed.In the mid-region of the distribution, where the NLO EW corrections are negligible, more prominent peaks for the two-loop corrections can be seen.The numerical values of these peaks should not be taken very seriously.It is just that the two-loop corrections are more relevant in the bins with peaks, and they should be taken into account in the bin-wise analysis of the data aimed at BSM searches.δ NLO for ϕ looks flat near the edges.However, it is indeed not the case, as shown in Fig. 10(c) for the bins between 0 and π 2 .In this range, δ NLO rises as dΓ NLO /dΓ LO goes down with increasing ϕ.Similarly, in the region beyond 250 • , as dΓ NLO /dΓ LO increases, δ NLO decreases (not shown explicitly) with an increasing value of ϕ.These features are independent of the scale choice for α s .

Conclusions and outlook
In this paper, we have computed the mixed QCD-electroweak corrections to the partial decay width of the H → e + e − µ + µ − channel.It is one of the most crucial decay channels to study the Higgs properties at the LHC and for new physics searches in the Higgs sector.In-house codes are developed to systematically compute the contributing two-loop matrix elements at O(αα s ) using the projector technique.The bare two-loop matrix elements are found to be free from IR divergences but contain UV divergences, which are regularized using dimensional regularization and eliminated by the on-shell counter-term renormalization procedure.The whole calculation is implemented in the publicly available event generator Hto4l code to perform the phase-space integration over the final state leptons and to obtain the improved predictions for the partial decay width.In the G µ scheme, the mixed QCD-electroweak correction to the partial decay width for the Higgs mass of 125 GeV is found to be 0.27% (18%) of the LO prediction (NLO EW correction) for the fixed QCD coupling.With the running coupling, the corrections become larger.The corrections can be significantly larger for the Higgs mass above the t t threshold, which is relevant to various BSM scenarios.
For the invariant mass distribution of the final state lepton pairs, the two-loop corrections with respect to the LO can cross 0.35% in the lower bins depending on the scale choice.For the angular distributions, which are crucial in measuring the spin-CP properties of the Higgs boson, the corrections with respect to the LO are of the order 0.6% in the bins around ϕ = 180 • when using the running α s .We also note that there are kinematic and angular bins in which the mixed QCD-electroweak corrections dominate the NLO EW corrections.These results may be helpful in probing new physics in the Higgs sector.Needless to say, our computational framework also allows predictions for H → γγ, γZ, γℓ + ℓ − , Zℓ + ℓ − decays.In addition to that, we can also use the ingredients calculated in this paper to predict O(αα s ) corrections for the H → ℓ + ℓ − ℓ + ℓ − (ℓ = e, µ).

Figure 2 :
Figure 2: Representative two-loop triangle diagrams contributing to the bare amplitude ofH → V 1 V 2 (V 1 , V 2 = Z, γ)with the top quark running in the loop.Diagrams with the reversed direction of the fermionic current are not shown.

Figure 3 :
Figure 3: Representative diagrams contributing to δM ααs S.E. with a quark running in the loop.In these diagrams, five light quark flavors(u, d, c, s and b) and the top quark contribute.The light quarks are taken to be massless.

Figure 4 :
Figure 4: Tree level diagrams involving O(αα s ) Zℓ l vertex counterterm denoted by a crossed circle.

Figure 5 :
Figure 5: Representative diagram for real corrections to the amplitude for H → ZZ * .

Figure 6 :
Figure 6: Representative one-loop triangle CT diagrams.The counterterm vertex proportional to α s is denoted by a crossed circle.

Figure 7 :
Figure 7: Representative one-loop self-energy CT diagrams.The counterterm vertex proportional to α s is denoted by a crossed circle.The diagrams with counterterm insertions on the lower leg are not shown.

Figure 8 :
Figure 8: Tree level CT diagrams.The counterterm vertex proportional to αα s is denoted by a crossed circle.

-Figure 9 :
Figure9: Effect of the mixed QCD-electroweak corrections on the invariant mass distribution of the final state lepton pair ℓ + ℓ − .In (a), δ LO is the ratio of the mixed QCDelectroweak correction and the LO contribution.In the upper panels of (b) and (c), δ NLO is the ratio of the mixed QCD-electroweak correction and the NLO EW correction.The quantity dΓ NLO /dΓ LO in the lower panels of (b) and (c) gives the ratio of the NLO EW correction and the LO contribution.For clarity, plot (c) displays the information of plot (b) in the region between 40 GeV and 80 GeV.The distributions in blue and red (black) color are the results obtained using the fixed scale Q = M Z and the running scaleQ = M ℓ + ℓ − (p T (ℓ − ))for α s , respectively.The dotted straight lines mark the results at the inclusive level.

-Figure 10 : 2 .
Figure 10: Effect of the mixed QCD-electroweak corrections on the distribution for the angle ϕ between the decay planes of the intermediate Z-bosons.In (a), δ LO is the ratio of the mixed QCD-electroweak correction and the LO contribution.In the upper panels of (b) and (c), δ NLO is the ratio of the mixed QCD-electroweak correction and the NLO EW correction.The quantity dΓ NLO /dΓ LO in the lower panels of (b) and (c) gives the ratio of the NLO EW correction and the LO contribution.For clarity, plot (c) displays the information of plot (b) in the region between 0 and π 2 .The distributions in blue and red (black) color are the results obtained using the fixed scale Q = M Z and the running scaleQ = M ℓ + ℓ − (p T (ℓ −)) for α s , respectively.The dotted straight lines mark the results at the inclusive level.