Three-loop mixed QCD-electroweak corrections to Higgs boson gluon fusion

We compute the contribution of three-loop mixed QCD-electroweak corrections ($\alpha_S^2\alpha^2$) to the $gg \to H$ scattering amplitude. We employ the method of differential equations to compute the relevant integrals and express them in terms of Goncharov polylogarithms.


Introduction
It is an open question if the scalar boson discovered by ATLAS and CMS collaborations in 2012 is indeed the Higgs boson of the Standard Model. To answer this question, Higgs boson production cross sections and decay rates are measured experimentally and compared to theoretical predictions. Thanks to the renormalizability of the Standard Model, properties of the Higgs boson, including its quantum numbers and its couplings to gauge bosons and matter particles, are completely fixed. As a result, it is possible, at least in principle, to describe production and decay of the Higgs boson in the Standard Model with unlimited precision.
On the other hand, if the Standard Model is not a complete theory, New Physics at scale Λ that couples to the Higgs boson is expected to modify its couplings to matter and gauge particles by an amount δg/g ∼ O(v 2 /Λ 2 ) where v = 246 GeV is the scale of electroweak symmetry breaking. Taking Λ = 1 TeV, we find a typical modification in Higgs couplings of about 5%, which sets a precision goal for both experimental measurements and theoretical computations.
The major Higgs-boson production channel at the LHC is gluon fusion gg → H; it contributes more than 90% to the total Higgs production cross section at 13 TeV. As discussed in Refs. [1,2], almost 5% of the gluon fusion cross section comes from mixed QCD-electroweak (QCD-EW) contributions where the Higgs boson couples to electroweak gauge bosons that, later, convert to gluons through a quark loop. The remaining 95% of the gg → H cross section is generated by pure QCD interactions.
The importance of the gg → H channel for the investigation of the Higgs boson properties motivated its extensive studies in the past. Since the Higgs boson is light and the top quarks are relatively heavy, it is possible to integrate them out and work with an effective Lagrangian where the ggH coupling is point-like. The gg → H production cross section in the m t → ∞ approximation is known through N 3 LO in QCD while other, smaller contributions, are known less precisely.
At present, the theoretical uncertainty in gg → H cross section originates from (see [1,2]) the O(2%) residual scale uncertainty in pure QCD contributions, the O(1%) uncertainty caused by unknown mass effects of b and c quarks in higher orders of QCD perturbation theory, and the O(1%) uncertainty in QCD-EW contributions. The latter uncertainty is particularly peculiar since computation of QCD corrections to mixed QCD-EW contributions exists [3]. However, although the leading order QCD-EW O(α S α 2 ) contribution is known for arbitrary values of the Higgs and electroweak gauge boson masses [4][5][6][7], the next-to-leading order (NLO) QCD corrections to it have only been computed in the unphysical limit m H m W,Z [3]. It turns out [3] that in this approximation the QCD corrections to mixed QCD-EW contributions are large and, in fact, similar to NLO corrections to leading QCD contributions. Nevertheless, the magnitude of QCD corrections and the approximate nature of the calculation of Ref. [3] are major reasons for assigning a large uncertainty to the QCD-EW contribution to gg → H in Ref. [1]. To remove this theoretical uncertainty, the NLO QCD corrections to QCD-EW contribution to gg → H should be re-computed for physical values of m H , m W and m Z .
The goal of this paper is to present one of the ingredients of such a computation -the virtual QCD corrections to mixed QCD-EW contribution to gg → H amplitude. The NLO QCD-EW corrections are given by three-loop O(α 2 S α 2 ) diagrams where gluons couple to electroweak gauge bosons through quark loops, and electroweak gauge bosons fuse into Higgs boson. In principle, all quark flavors contribute to mixed QCD-EW corrections but, since the Higgs boson is light, the contribution of top quark loops is small [3]. For this reason we focus on the contribution of massless quarks in what follows.
To compute the three-loop mixed QCD-EW contribution to the gg → H amplitude, we use the method of differential equations [8][9][10] to calculate the relevant master integrals (MIs). We employ the concept of canonical basis introduced in Ref. [11], and express the MIs in terms of Goncharov's multiple polylogarithms (GPLs) [12][13][14], up to weight 6. We use the mathematical limit of small Higgs mass m H m W,Z to determine the boundary values for the solutions of the differential equations. Our final three-loop result for QCD-EW contributions to gg → H amplitude is expressed in terms of GPLs up to weight five.
The paper is organized as follows. We introduce the gluon fusion amplitude gg → H and discuss the QCD-EW contributions in Section 2. In Section 3 we describe the classes of integrals that contribute to the QCD-EW amplitude. In Section 4 we discuss differential equations for MIs and their simplification to a canonical form. We introduce the large-mass expansion in Section 5 and explain how we use it to provide boundary conditions for the solutions of the differential equations. The final result for the gg → H amplitude is given in Section 6. We conclude in Section 7. Additional information about integral topologies and MIs can be found, respectively, in Appendices A and B. We provide the analytic results for leading and next-to-leading order QCD-EW contributions to gg → H amplitude in the ancillary file.
Additional material that includes the expressions for MIs, the basis of canonical functions, the matrix of coefficients of the differential equations and the integration constants is available at https://www.ttp.kit.edu/_media/progdata/2017/ttp17-047/anc.tar.gz.

QCD-EW contribution to gg → H amplitude
We show the mixed QCD-electroweak corrections to the gg → H amplitude in Fig. 1. They appear for the first time at two loops. The two gluons annihilate to electroweak vector bosons through a quark loop. The vector bosons fuse into a Higgs boson in a second step. These leading order QCD-electroweak contributions are well-known [4][5][6][7]. Our goal in this paper is to compute the QCD corrections to the leading order QCD-EW contribution to gg → H amplitude.
Both W and Z bosons appear in diagrams that contribute to the QCD-EW amplitude. In the case of W -bosons, we account for quarks of the first two generations; when considering Z-bosons, we include, in addition, bottom quarks. We treat all quarks as massless. To write an expression for the gg → H amplitude we assume that the two incoming gluons g 1 and g 2 carry on-shell momenta p 1,2 , color indices c 1,2 and polarization labels λ 1,2 . The momentum of the outgoing Higgs boson is taken to be p 3 = −p 1 − p 2 , with p 2 3 = m 2 H = s. QCD gauge invariance ensures that the gg → H amplitude depends on a single form factor where ε λ 1,2 are gluon polarization vectors that satisfy the transversality conditions ε λ i ·p 1 = ε λ i ·p 2 = 0, for i = 1, 2. It is convenient to construct a projection operator that can be applied to individual Feynman diagrams to extract their contributions to the form factor F. The projection operator reads where d = 4 − 2ε is the space-time dimensionality. Using the standard expression for the sum over polarizations for each of the two gluons, consistent with the transversality conditions it is easy to show that the following equation holds We can further simplify this equation if we write the amplitude as apply the projection operator and sum over helicity labels of the two gluons, to obtain This formula can be used to extract the contribution of individual diagrams to the gauge-invariant form factor F. Examples of three loop diagrams that contribute to the NLO QCD corrections to the mixed QCD-electroweak gg → H amplitude are shown in Fig. 2. The majority of the diagrams contributes both to Z and W channels with a few exceptions that contribute only in the Z channel; some examples are shown in Fig. 3. However, all diagrams that distinguish Z and W contributions vanish either because of color algebra (Fig. 3a), or Furry's theorem (vector current contribution in case of Fig. 3b) or the fact that we consider complete generations of massless quarks (   axial current). We note that in the latter case, because of the axial anomaly, the contribution of the b-quark is not well defined without the top quark; in this paper we ignore this issue and completely discard all contributions that lead to diagrams of the type ( Fig. 3b) but we include the b-quark in all other three-loop QCD-EW diagrams that we consider in this paper. After these simplifications, there remain 47 three-loop non-vanishing diagrams for both W -bosons and the Zbosons. All the relevant contributions can be obtained by considering diagrams where a massive vector boson interacts with the quark loop through a vector current.
The [j]s denote inverse Feynman propagators that we specify in Appendix A. Since W and Z bosons never appear in the same Feynman diagram, we use a generic notation M for the mass of the vector bosons when we discuss Feynman diagrams and integrals. All Feynman integrals are analytic functions in the variables s = (p 1 + p 2 ) 2 and M 2 . Discontinuities common to all parent topologies occur at s ≥ 0 (on-shell massless intermediate states), s ≥ M 2 (production of an on-shell vector boson plus two massless quarks), and s ≥ 4M 2 (production of an on-shell pair of vector bosons).
To compute the different Feynman integrals contributing to the form factor, we identify eleven different topologies; these topologies can be written using three different sets of the inverse Feynman propagators. The topologies are described in Appendix A. We use the program Reduze2 [15] to express all Feynman integrals that contribute to the invariant form factor F through 95 master integrals (MIs). The full list of MIs can be found in Appendix B.

General considerations
To calculate the MIs, we employ the differential equations method. The MIs depend on two dimensionfull variables s and M 2 and it is convenient to trade them for one dimensionfull and one dimensionless variable. We take s to be the dimensionfull variable, extract the mass dimensions of the MIs and write The dimensionless variable y reads Values of the variable y for different relations between M 2 an s are shown in Table 1, together with prescriptions for analytic continuation.

Variable Prescription Euclidean region Minkowski region
Below threshold Above threshold To study the non-trivial dependence of the MIs on y, we differentiate the vector of rescaled MIs J(ε, y) with respect to y, express the resulting integrals again through MIs and obtain a closed system of differential equations dJ(ε, y) dy = A(ε, y)J(ε, y).
To compute the MIs we need to solve this system of differential equations. To do so, we first determine a canonical basis of MIs. With the integrals in the canonical basis, the system of differential equations assumes the ε-homogeneous form [11] dF(ε, y) dy = εA(y)F(ε, y). The necessary and sufficient conditions for the existence of a canonical basis are not known, but it is expected that for all cases where solutions can be written in terms of Chen iterated integrals, a canonical basis exists. Procedures to determine canonical bases were discussed in Ref. [11,[16][17][18][19][20][21].
Once the canonical basis is found and integrals are normalized in such a way that the finite ε-limit exists lim ε→0 F(ε, y) = F 0 , it becomes straightforward to obtain the solution to the system of differential equations as series expansion in the dimensional regularization parameter ε. An important property of the canonical form of differential equations in Eq. (4.4) is that the matrix A(y) contains only a finite number of simple poles in the variable y [16] dF(ε, y) where B k are y-independent matrices. Equations of the type shown in Eq.(4.5) are called fuchsian equations.
The solution of Eq.(4.5) is obtained upon iterative integration Thanks to the fuchsian form of the matrix A, the nested integrations in Eq.(4.6) can be performed in terms of multiple polylogarithms [22,23], best known in the physics literature as Goncharov's polylogarithms (GPLs) [12][13][14]24]. GPLs are defined recursively using the following equations where m w = (m w , m w−1 ). The weight w of a GPL is the length of the vector m w . GPLs evaluated at rational points are expressed by constants of the same weight. In Section 5, we use this property to efficiently fix the integration constants F

The system of differential equations
We are now in position to discuss the system of differential equations and the steps that are required to transform it into a canonical form. To achieve that, we used a combination of different methods that we now explain. The matrix A(ε, y) that appears on the right hand side of the differential equation has a blockdiagonal form. This form implies that differential equations for simpler MIs close and that simpler integrals appear as inhomogeneous terms in differential equations for more complex ones. This form of the differential equations suggests that one should start with bringing simplest topologies to a canonical form and then moving on to more complex topologies, step by step.
For MIs with four, five or six propagators coupled together in at most 3 × 3 blocks, we first adjusted powers of propagators and introduced ε-dependent prefactors, as described in Ref. [17,18]. It is well understood how to do this to reach canonical form for simple integrals, such as bubbles and triangles where the guiding principle is to ensure that integrals are free from ultraviolet divergences. After that, the candidate integrals are multiplied by a polynomial in ε of the form ε a 1 (c 1 + c 2 ε) a 2 , to bring the system into a linear-ε form We then integrate out the ε-independent part of the differential equation and obtain the canonical form We stress that it is not always possible to integrate out A 0 (y) in the presence of off-diagonal terms and that, even if successful, this procedure is not guaranteed to keep the Fuchsian form of the system of differential equations. Nevertheless, it provides a practical way to achieve canonical form for some differential equations relevant for our problem. For larger sets of coupled MIs or for MIs with a larger number of denominators, the differential equations become too complex to reach the linear-ε form by a clever choice of propagator powers and prefactors. For these differential equations, we employ the algebraic algorithm presented in [19] and, in particular, its computer-algebraic implementation Fuchsia [25]. We also use techniques described in Ref. [20] and their implementation in the Mathematica package CANONICA [26]. Despite the power and versatility of these algorithms, their blind application to our system of differential equations results either in a failure of the procedure to reach a canonical fuchsian form, or in an incorrect redefinition of the MIs that have already been fixed. In the second case, the new system shows the correct differential structure, but the constants F (n) 0 are, in general, not of an uniform weight n. Again, by carefully selecting the MIs and choosing their prefactors, it is possible to overcome these problems.
A useful procedure to select candidate integrals for the canonical basis is to inspect their generalized unitarity cuts, as explained in [16]. The idea is that if one replaces a propagator in an integral by a delta-function 1/(p 2 − m 2 ) → δ(p 2 − m 2 ), the differential equation that this integral satisfies does not change except that all integrals on the right hand side of the differential equation where this propagator is not present have to be set to zero. Thanks to this observation, by cutting a MI in different ways we can inspect different subsets of the differential equation that it satisfies. By replacing all propagators of a given integral with the δ-functions (the maximal cut), we obtain the differential equation which involves only MIs of the highest topology since all other integrals drop out [27]. The analysis of the cuts can be simplified by use of the so-called Baikov representation [28,29]. The upshot of these discussions is that if all the cut integrals are of the d log-iterated form , a 1 , . . . , a n ) d log R n−1 (x, a 1 , . . . , a n−1 ) . . . d log R 1 (x, a 1 ), (4.12) where C(x) is a function of the kinematic invariants x and R(x, a 1 , . . . , a k ) are rational functions of x and of the integration variables a i , the integral I is a valid candidate for the canonical basis. For our purposes, we require integrals with large number of propagators; they satisfy complex differential equations with up to four MIs coupled. For this reasons, the study of all the different cuts is prohibitive and in most of the cases we limit ourselves to the maximal cut. As explained above, the latter only allows us to inspect the part of the differential equation that corresponds to the highest topology. Once the d log-form is achieved for the maximal cut of all the coupled highest-level MIs, the homogeneous part of the equations is ensured to be canonical. Of course, integrals of lower topologies that appear in the differential equations still do not have a a canonical form, but the study of the maximal cut has proven to be sufficient to provide a starting points for a successful application of either Fuchsia or CANONICA.
Performing manipulations that we just described, we were able to transform the system of differential equations into a canonical Fuchsian form. We write it as (4.13) We have chosen arguments of the logarithms to ensure that the logarithms are real-valued for 0 < y < 1 (s < 0). It follows from this form that integration kernels are identical to what was found in the calculation of MIs for leading order QCD-electroweak amplitude [7]. They read where the function f (r; y) is a short-hand notation for a particular linear combination of two integration kernels f (r; y) = 2y − 1 that make the alphabet complex-valued. Note that this expansion is only needed for numerical evaluations, while the integration of the differential equation can be performed using the quadratic form that we associate with the alphabet symbol r [30,31]. Indeed, thanks to the linearity of both the differential equations and the definition of the GPLs, it is immediately possible to express the obtained generalized GPL in terms of the canonical ones

Boundary conditions and the large-mass expansion
To fully solve the system of linear differential equations Eq. (4.13), we need to determine the integration constants F (n) 0 . The easiest way to fix them is to compute the MIs at some kinematic point y 0 using alternative methods and then compare the result of the computation with the solution of the differential equations. As we already explained in Ref. [7], it is convenient to compute the required integrals at y = 1, corresponding to s/M 2 = 0. This kinematic condition can be studied using the large-mass expansion procedure [32] to independently compute the MIs in this limit. In this Section we explain how to apply the large-mass expansion procedure to obtain boundary values for the MIs.
A typical MI that we need to compute depends on two different scales: the energy of the external gluons p 1 ∼ p 2 ∼ √ s, and the mass of the electroweak vector bosons M . We consider the hypothetical limit M 2 s. To provide a non-vanishing contribution in dimensional regularization, the loop momentum k that flows through any sub-set of internal lines has to scale either as k ∼ √ s or as k ∼ M . All possible combinations of scalings for the three loop momenta are allowed and have to be considered, provided that momentum conservation is satisfied at each vertex. 4 Once a valid momentum scaling is identified for a particular Feynman integral, the integrand is Taylorexpanded in all small momenta (both external and loop ones) and then integrated. This procedure is performed for all possible momenta assignments and the results are summed over to obtain the large-mass expansion of the original integral.
As an example of the large-mass expansion procedure, we consider MI I 7 where a dot on a line implies that the corresponding propagator is raised to second power. Among all allowed choices of large and small momenta only two give non-zero contributions. Indeed, it is possible to have large momenta flowing in the two-loop self-energy sub-graph and small momentum flowing in the "outer" loop, or to have all the three loop momenta to be large. It is clear that lines with large momenta "decouple" from adjacent propagators: they become tadpoles that multiply a diagram obtained by shrinking the corresponding "large momentum" lines to a point. For the two possible kinematic regions described above, the large mass expansion reads (bold internal lines in the original integral indicate the flow of large momenta) We note that in the first case the result is exact since, upon expansion, the tadpole produces a term that cancels one of the propagators of the massless bubble. We obtain the large-mass expansion of the integral (5.1) by summing over the two possible momenta flows described above (5.4) All integrals that appear on the right hand side of this equation can easily be evaluated.
Clearly, the number of different momenta flows that need to be considered increases as one moves from relatively simple to more complex integrals. Nevertheless, the complexity remains manageable. Indeed, once the large mass expansion is performed, the most complex integrals are three-loop tadpoles and two-loop massless triangles that are well-known, see e.g. Refs. [33,34].
The large-mass expansions of the MIs computed to relevant order in s/M 2 expansion, is compared with the y → 1 limit of the solutions of the differential equations; requiring that the two agree, we obtain all the integration constants F (n) 0 , n ≤ 6, c.f. Eq.(4.6). Since analytic expressions for GPLs with letters e ±iπ/3 , evaluated at y = 1, are unknown at present, we follow a numerical approach to determine the constants. The fact that the integrals of the canonical basis are of uniform weight is crucial, since in this case only a small variety of constants can appear at each order of their ε-expansion. In particular, F must be proportional to π 2 , ζ(3) and π 4 , respectively, F 0 should be given by a linear combination of π 2 ζ(3) and ζ(5), and F (6) 0 by a linear combination of π 6 and ζ 2 (3). We evaluate with high precision the canonical integrals written in terms of GPLs at y = 1 and match their numerical values to rational linear combinations of analytic constant using PSLQ algorithm, to at least 750 digits. As a final check, all MIs have been compared for at least two different values of s and M 2 to the numerical results obtained using the programs SecDec [35] and pySecDec [36]. In all cases agreement was found. 6 The final result for the gg → H amplitude As explained in Section 2, the mixed QCD-electroweak contributions to gg → H amplitude are parametrized in terms of a single form factor. The form factor receives contributions from diagrams with W and Z bosons; such diagrams appear for the first time at two loops. We account for quarks of the first two generations and include b-quark contributions in diagrams with Z-boson. 5 All quarks are taken to be massless. The CKM matrix is approximated by the identity matrix.
Computation of three-loop contribution to the form factor yields divergent results. The divergences are of both ultraviolet and infrared origin. The ultraviolet divergences are removed by ultraviolet renormalization; the infrared divergences remain but will, eventually, get canceled by the real emission contributions when infrared safe observables are computed.
To remove the ultraviolet divergences, we only need to renormalize the strong coupling constant; this is so because for massless quarks both vector and axial currents are conserved and, for this reason, weak couplings do not need an ultraviolet renormalization. To renormalize the strong coupling constant, we use the relation between bare and MS coupling constants that reads As we already mentioned, we discard certain anomalous-type diagrams also for the third generation, see Section 2.
where γ E is the Euler-Mascheroni constant and β 0 is the first coefficient of the QCD β-function. It reads where C A = N c = 3, T F = 1/2 and N f = 5 is the number of massless fermions that contribute to the renormalization of the QCD coupling constant. We write the UV-renormalized form factor as where α is the QED fine structure constant, ϑ W is the Weinberg angle and v = m W sin ϑ W / √ πα is the Higgs field vacuum expectation value. In addition, The two terms in Eq.(6.3) denote contributions of diagrams with W and Z bosons. The function A can be written as an expansion in the strong coupling constant It is interesting to remark that individual integrals that contribute to gg → H amplitude are strongly divergent; the strongest divergence is O(ε −6 ). When the integrals are combined to produce a form factor, all divergences stronger than O(ε −2 ) cancel out. The ultraviolet renormalization described above removes some of the 1/ε poles but infra-red singularities that start at 1/ε 2 still remain. The general structure of these singularities in QCD is described by Catani's formula [37]. This formula, applied to the form factor that describes mixed QCD-electroweak corrections to gg → H amplitude reads where Note that the leading order amplitude in Eq.(6.6) is needed through order O(ε 2 ); it was computed in Ref. [7]. An analytic expression for the finite part of the NLO amplitude A fin NLO is given in the ancillary file.
For future reference, we give numerical values of the two functions A LO and A fin NLO , for physical values of the Higgs boson and vector boson masses. Taking m H = 125.09 GeV, m W = 80.385 GeV, m Z = 91.1876 GeV (see [2]), N c = 3, N f = 5 and µ = m H , we find  H > m 2 Z > m 2 W , so that the first two cuts s = 0 and s = m 2 V should both contribute to the imaginary parts in Eq. (6.8). However, comparing absolute values of real and imaginary parts in Eq. (6.8), the difference in the relative importance of imaginary parts at NLO compared to LO is striking. The reason for this can be traced back to the fact that, contrary to expectations, the LO amplitude is actually real in the region 0 < s < m 2 V , and that the discontinuity at s = 0 contributes only to the imaginary part of NLO amplitude. This is easy to understand. Indeed, the imaginary part of the LO amplitude in this region is obtained by cutting through the fermion loop and schematically corresponds to the process gg → qq | qq → H , see Fig. 4a. Clearly, since the Higgs boson cannot couple to massless fermions, this contribution must vanish. We emphasize that this is a feature of the amplitude that does not hold at the level of individual MIs. The same is not true at NLO, where more cuts contribute to the discontinuity that starts at s = 0. In particular, a cut through a gluon loop provides a non-vanishing contribution to the imaginary part of the amplitude for 0 < s < m 2 V from the re-scattering process gg → gg | gg → H, see

Conclusions
We computed the three-loop virtual contributions to next-to-leading order mixed QCD-electroweak corrections to Higgs boson production amplitude in the annihilation of two gluons gg → H. The analytic result for the amplitude is obtained for arbitrary relation between the Higgs boson and the electroweak gauge boson masses. We computed the required integrals using the method of differential equations and obtained the boundary constants using the large-mass expansion valid in the limit m H m W,Z . The finite part of the three-loop amplitude is written in terms of Goncharov polylogarithms and is easy to evaluate numerically.
To understand the impact of the mixed QCD-electroweak corrections on the Higgs boson production cross section, the virtual corrections computed in this paper will have to be combined with the real emission contributions of the type gg → Hg, qg → Hq where, again, the Higgs boson couples to electroweak vector bosons. The computation of the relevant real emission amplitude is non-trivial as it involves two-loop box diagrams with both internal and external masses. Nevertheless, it is conceivable that, with the current computational technology, analytic results for these contributions can be obtained.

B Master integrals
The form factor F can be expressed in terms of 95 MIs, depicted in Figs. 8, 9 and 10. A circled number besides the graph indicates that the corresponding denominator features in the numerator.
The MIs with an asterisk ( * ) do not contribute to the NLO form factor, but appear in the differential equations.