Mixed QCD-EW corrections for Higgs leptonic decay via $HW^+W^-$ vertex

We consider the two-loop corrections to the $HW^+W^-$ vertex at order $\alpha\alpha_s$. We construct a canonical basis for the two-loop integrals using the Baikov representation and the intersection theory. By solving the $\epsilon$-form differential equations, we obtain fully analytic expressions for the master integrals in terms of multiple polylogarithms, which allow fast and accurate numeric evaluation for arbitrary configurations of external momenta. We apply our analytic results to the decay process $H \to \nu_e e W$, and study both the integrated and differential decay rates. Our results can also be applied to the Higgs production process via $W$ boson fusion.


Introduction
After the discovery of the Higgs boson in 2012, one of the primary goals of current and future high energy collider experiments is the precision measurements of various Higgs properties. These include the Large Hadron Collider (LHC) and its future upgrade High-Luminosity LHC (HL-LHC) [1], as well as future electron-positron colliders such as the Circular Electron Position Collider (CEPC) [2], the Future Circular Collider (FCC-ee) [3] and the International Linear Collider (ILC) [4]. An important property of the Higgs boson is its coupling to weak gauge bosons, namely, the Z 0 boson and the W ± boson. Their precision knowledge is crucial for verifying the spontaneous electroweak symmetry breaking mechanism and for exploring the nature of the Higgs sector. If the Higgs sector is more complicated than the simplest version assumed by the Standard Model (SM), e.g., if the Higgs boson is composite, these gauge couplings may differ from the SM value by a tiny amount. The precision measurement of these couplings will then help us to detect such small effects resulting from physics at higher scales.
With the current LHC data, the HZZ coupling and the HW + W − coupling have been measured with an uncertainty of 11% and 15%, respectively [5]. At HL-LHC, the precision is estimated to be about 1.5% and 1.7%, respectively [6]. At future electron-position colliders, the precisions for these couplings can be greatly improved. For the HZZ coupling, the accuracy is estimated to be about 0.25% at the CEPC with a center-of-mass energy √ s = 250 GeV [2], and about 0.17% at the FCC-ee with √ s = 365 GeV [3]; while for the HW + W − coupling, the precision is estimated to be about 1.4% at the CEPC and about 0.43% at the FCC-ee. To get the most out of the future high precision data, it is necessary to provide accurate theoretical predictions for production and decay processes sensitive to these couplings. Those related to the HZZ vertex have been extensively studied in the literature [7][8][9][10][11][12][13][14][15][16][17][18].
In this paper, we concentrate on the HW + W − vertex, and study its higher order perturbative corrections within the SM.
To setup the context, we will consider the decay process H → W * W → ν e eW in this paper, although our results can also be applied to the production process e + e − → Hν eνe . The leading order (LO) results [19][20][21] and the next-to-leading order (NLO) electroweak (EW) corrections [7] for this decay process have been known many years ago. Beyond the NLO, contributions enhanced by the large top quark mass have been considered at two-loops [11,13] and three-loops [12], albeit the assumption m H > 2m W . For a light Higgs boson, Ref. [14] reconsidered these calculations, and also gave predictions for the differential decay rates.
In this paper, we provide exact analytic results for the next-to-next-to-leading order (NNLO) mixed QCD-EW corrections to the HW + W − vertex at order αα s . The analytic expressions are valid for arbitrary values of the external momenta. Hence it is applicable to cases where both W bosons are off-shell, and also to cases where the large top mass approximation may not hold. We apply these analytic results to the Higgs decay process, and provide numeric predictions for both integrated and differential decay rates. These results will be useful for extracting the HW + W − coupling from future high precision experimental data. This paper is organized as follows. In Section 2 we briefly introduce our notations and the LO and NLO results. In Section 3 we present the analytic calculation of the mixed QCD-EW corrections for the HW + W − vertex, with generic momenta configurations. In Section 4, we apply our analytic results to the H → ν e eW decay process, and obtain numeric predictions. The summary and outlook come in Section 5, and we leave some lengthy expressions to the Appendix.

Notations and lower-order results
In this section we setup our notations and present the LO results as well as the NLO EW corrections for the partial decay width. In our calculation we will neglect masses of all light fermions except that of the top quark. At LO we assign the momenta as The scalar products of the momenta lead to The LO partial decay width is then written as where the 3-body phase-space integration measure is given by In general the 3-body phase-space integral can be performed numerically with a Monte Carlo program, which will be useful when we compute differential decay widths. For the fully integrated decay width at LO, the phase space integral can be calculated analytically and we find where s W = sin θ W and r = m H /(2m W ). The NLO EW corrections to this decay process involve closed fermion loops, exchanges of electroweak gauge bosons and Higgs boson, as well as real photon emissions. For these corrections we invoke the program MadGraph5 aMC@NLO [22] which can automatically compute NLO QCD and EW corrections to standard model processes. In the rest of the paper, we present our calculation of the NNLO mixed QCD-EW corrections at order αα s .

The HW + W − vertex at two loops
In this section, we present the calculation of the two-loop form factors entering the mixed QCD-EW corrections for the HW + W − vertex. This serves as a main new result of this paper. A typical Feynman diagram is depicted in Fig. 1. The most generic Lorentz structure of the HW + W − vertex can be written as: Among the scalar coefficient functions T i in the above formula, only T 1 and T 4 contribute to the partial decay width. At NNLO, the coefficient functions can be written as linear combinations of scalar two-loop integrals of the form where d = 4 − 2 in dimensional regularization, k 1 and k 2 are loop momenta, and (IBP) reduction with the program packages FIRE6 [23] and LiteRed [24], and crosscheck using Kira [25]. We find 38 master integrals in total. We note that some of the sub-topologies have been studied in [26], with 31 master integrals. 1 The additional master integrals we get are: The last integral in the above does not appear in the amplitude, and is included because we perform the reduction from the top topology with 7 propagators. The other six integrals are required in the calculation of the two-loop amplitude. They can be related to T 11 , T 25 , T 27 -T 30 of Figure 4 in [26] by an exchange of external momenta p 1 ↔ p 2 (p 3 ↔ q in our notation).

Canonical basis of the two-loop integrals
We would like to find linear combinations of integrals as in Eq. (3.2) which satisfy the canonical-form differential equations. Namely, the canonical basis f (u, v, w; ) of the master integrals satisfy differential equations of the form where A i are constant matrices independent of kinematic variables and the dimensional regulator. The "letters" α i ≡ α i (u, v, w) are algebraic functions of the kinematic variables u, v and w.
To find the canonical basis f (u, v, w; ), we use the method of [39] supplemented by suitable linear transformations to further simplify the differential equations. The method of [39] utilizes the Baikov representation [40] and look for d log-form integrals which serve as candidates for the canonical basis. The Baikov representation for integrals of the form Eq. (3.2) can be written as where z i ≡ D i are the propagator denominators which become integration variables in the Baikov representation, and the vector z is the collection of z i . The prefactor N and the function u(z) are given by The integration domain C in Eq. (3.5) has the property that u(z) vanishes on its boundary ∂C.
For illustration purposes, we consider the {1, 0, 1, 1, 1, 0, 1} sub-topology, which turns out to be the most complicated one. This includes all integrals of the form Eq. (3.5) with the powers a 1 , a 3 , a 4 , a 5 and a 7 being positive, and the powers a 2 and a 6 being non-positive. For the construction of d log-form integrals, it is reasonable to first consider the cases where either a 2 or a 6 is fixed to zero. As an example, we will fix a 2 = 0. It is then possible to perform the integration over z 2 in Eq. (3.5), resulting in where z is given by z with z 2 removed. The new prefactor and the polynomials are given by The integrals in this sub-topology with a 2 = 0 can then be expressed as which is equivalent to a loop-by-loop Baikov representation.
We now want to search for rational functionsφ(z ) such that the integrand u 2 (z )φ(z ) takes the d-dimensional generalized d log-form [39]. We perform the construction in the variable-by-variable approach in the order z 6 , z 1 , z 5 , z 3 , z 4 , z 7 , and find 4 candidates for the canonical Feynman integrals in this sub-topology: These candidates can be converted to Feynman integrals of the form (3.5) using intersection theory [41][42][43][44][45][46][47] or via IBP relations. The results are given by Where we have introduced the following two square roots: Applying the same procedure to other (sub)-topologies, we are able to construct 38 d log-form integrals in the Baikov representation, and convert them to Feynman integrals using either traditional IBP method or the intersection theory. We have checked that they indeed satisfy the canonical-form differential equations (3.4). Further (simple) linear transforms can be applied to simplify the matrices A i . The final form of the canonical basis f (u, v, w; ) in terms of linear combinations of Feynman integrals is given in Appendix A.

Solution to the canonical differential equations
The differential equations (3.4) involve 18 independent letters as following: where we have introduced the following polynomials: The solution to the canonical form differential equations (3.4) can be formally written as a path-ordered exponential integrated from the boundary point We choose the boundary point to be x 0 = (0, 0, 0), which corresponds to zero external momenta. The values of the mater integrals at the boundary can be easily obtained and are given by The path-ordered exponential (3.16) can be expanded as power series in . At each order in we need to evaluate iterated integrals in terms of generalized polylogarithms (GPLs) [48]. For that it is convenient to perform a change of variables such that the integration kernels are rational functions without square roots. There are two ways to achieve that, which lead to the same final results. The first is to perform variable changes on u, v and w, such that both R 1 and R 2 become rational functions of the new variables. This can be done using the method of [49,50]. The second way, which we'll describe in more details, amounts to choosing a specific integration path in the (u, v, w) space, such that the integration reduces to a single-variable problem.
Observing that R 2 is homogeneous in u, v and w, it is reasonable to pick a straight path from (0, 0, 0) to (u, v, w) parameterized by a variable t as (tu, tv, tw). During the integration over t, R 2 simply becomes t multiplied by a constant factor: R 2 (tu, tv, tw) = tR 2 (u, v, w) . (3.18) On the other hand, we now need to perform a variable change on t such that R 1 (tu) becomes a rational function. This can be simply done with When t goes from 0 to 1, x goes from 0 to 2(R 1 (u) − u). This determines the integration domain over x.
After the variable change, the iterated integrals at order n have the generic with G (n−1) (x) being functions at order n−1 . These integrals can be naturally evaluated in terms of GPLs, with the branch points x i given by 2v , 2w , The resulting GPLs can be evaluated numerically using program libraries GiNaC [51][52][53] and handyG [54]. We have checked that our results for the master integrals are in good agreement with the numeric results of pySecDec [55].
4 Application to the decay process 4

.1 Analytic results at NNLO
In this section, we apply the two-loop integrals calculated in the previous section to the decay process H → eν e W at order αα s . The decay process receives three kinds of contributions at this order: the corrections to the HW + W − vertex, the corrections to the propagator of the off-shell W ± field, and the corrections to the eν e W vertex (arising from counter-terms). In practice these can all be described by their contributions to the coefficient functions in Eq. (3.1) (with the third kind of contributions only affecting T 1 ). The squared amplitude at this order can then be written as To calculate the coefficient functions, we generate the relevant Feynman diagrams using FeynArts [56]. The resulting amplitudes are further manipulated with Mathematica. In addition to the Feynman integrals calculated in the previous section, we also need to calculate a couple of extra integrals arising in the W ± propagator and in the counter-terms. We have checked these results against those in the literature [57,58].
We renormalize the fields and the masses in the on-shell scheme, with pole masses m W , m Z , m H and m t as input parameters. The weak mixing angle is defined by the on-shell relation c W ≡ cos(θ W ) = m W /m Z , whose renormalization is given by For the renormalization of the fine structure constant α, we choose two different schemes: the α(m Z )-scheme with α(m Z ) = α(0)/(1−∆α(m Z )) as the input parameter, where ∆α(m Z ) encodes the contributions to the renormalization of α from light fermion loops (including leptons and 5 flavors of light quarks); and the G µ -scheme with the Fermi-constant G µ determined from the muon decay data as the input parameter. In the G µ -scheme, the coupling is given by

Numeric results
We now present the numeric predictions from our calculations.  80.379 GeV, G µ = 1.1663787 × 10 −5 GeV −2 , α(m Z ) = 1/128.929 and α s (m Z ) = 0.1179 [59]. The default renormalization scale for α s is chosen as µ = m H . In Table 1 and Table 2, we show the LO, NLO EW, and NNLO QCD-EW predictions to the partial decay widths in the α(m Z )-scheme and G µ -scheme, respectively. We find that the mixed QCD-EW corrections amount to about 1% of the LO prediction in the α(m Z )-scheme, and about 0.2% of the LO prediction in the G µ -scheme. The smallness of the contributions in the G µ -scheme can be partly attributed to the fact that corrections to the W propagator and to the eν e W vertex are absorbed into the coupling constant. Nevertheless, we find that the NNLO results in the two schemes are rather close to each other, showing the stability of our prediction against scheme changes.
We now proceed to study the differential decay rates. In Fig. 2 we depict the M eν distribution at various orders in perturbation theory, in the two renormalization schemes mentioned above. This distribution gives us information about the virtuality of the off-shell W boson. Here we observe similar behaviors as for the integrated decay rate. The NNLO corrections are much smaller in the G µ scheme, as shown in the right plot in Fig. 2, where the green and blue curves almost completely overlap with each other. Since the neutrino is not observable (although could be reconstructed), it is also interesting to study the invariant mass of the visible part of the decay products. In Fig. 3 we show the M eW distribution, again in the two schemes at different orders. This distribution can be measured at the LHC and the future Higgs factories, with the W boson decaying hadronically. Our results then provide the high precision theoretical predictions to be compared with the experimental data.

Summary and Outlook
In this paper, we have studied a class of two-loop triangle integrals entering the O(αα s ) corrections to the HW + W − vertex. We have constructed a canonical basis consisting of 38 master integrals using the Baikov representation and intersection theory. We have derived the -form differential equations for the master integrals, and are able to find fully analytic solutions in terms of GPLs. The fully analytic form allows fast and accurate numeric evaluation for any combination of external momenta and internal masses, which is important for phenomenological studies. We apply our results to the H → ν e eW decay process. For the integrated decay width, we find that the NNLO mixed QCD-EW corrections can reach 1% of the LO result in the α(m Z ) scheme. The size of the corrections are reduced to about 0.2% in the G µ scheme. This seems to indicate that the perturbative convergence in the G µ scheme is better. However, this needs to be confirmed by the behavior of the purely electroweak two-loop corrections of order α 2 . We have also studied the differential decay rates and drawn similar conclusions.
For the decay process, the complete O(αα s ) corrections are not expected to be quite different in size compared to those calculated in the heavy top limit [11][12][13][14]. However, the fully analytic form of our results enable us to study cases where the virtuality of W ± bosons exceeds the top quark mass, where the heavy top approximation is no longer valid. This is the case for, e.g., the W boson fusion process W + * W − * → H. The O(αα s ) corrections for this production process will be presented in a forthcoming article.
where u, v, w, R 1 and R 2 are defined by Eq. (3.3) and Eq. (3.13), respectively, and the functions g(u, v, w) and h(u, v, w) are given in Eq. (3.15).