Two-loop electroweak corrections to Higgs-gluon couplings to higher orders in the dimensional regularization parameter

We compute the two-loop electroweak correction to the production of the Higgs boson in gluon fusion to higher orders in the dimensional-regularization parameter $\epsilon = (d-4)/2$. We employ the method of differential equations to compute the relevant integrals and express them in terms of Goncharov polylogarithms. Our result provides one of the necessary inputs for the computation of mixed three-loop QCD-electroweak corrections to $gg \to H$.


Introduction
The recent discovery of the Higgs boson and the non-observation of any New Physics at the LHC establishes the validity of the Standard Model as the low-energy effective theory of Nature. At the same time, the apparent inability of the Standard Model to explain several experimental facts makes the need for physics beyond the Standard Model (BSM) as strong as ever. Searching for clues about BSM physics is in the focus of contemporary particle physics. The Higgs boson is bound to play an important role in this endeavor. Indeed, the Higgs mechanism in the Standard Model is very simplistic and rather ad hoc. At the same time, there are many extensions of the Standard Model where the Higgs boson is the only particle that is sensitive to rich physics beyond it. More generally, if the Higgs boson is responsible for generating masses not only of the Standard Model but also of some BSM particles, which appears to be necessary for protecting the Higgs boson mass from large radiative corrections, these new particles may affect the couplings of the Higgs boson to gauge bosons and fermions through radiative corrections. A percent modification of the Higgs couplings is a generic consequence of physics beyond the Standard Model at the energy scale of about 1 TeV. Therefore, measurement of the Higgs bosons couplings to Standard Model particles to this level of precision is an important goal of the LHC physics program.
The major production channel of Higgs bosons at the LHC is gluon fusion. The recent computation of the three-loop QCD corrections to σ gg→H significantly reduces theoretical uncertainty in the predicted cross section. According to Ref. [1], the theory uncertainty in σ gg→H is close to 5% and the uncertainties related to imprecise knowledge of parton distribution functions and the strong coupling constant are close to 4%. The theoretical uncertainty has several sources such as the residual scale dependence of the three-loop QCD result, imperfect knowledge of the bottom quark contribution to gg → H and the mixed three-loop QCD-electroweak corrections which are known in the unphysical limit m Z,W m H [2]. Each of these sources contributes similar amount to the final uncertainty which implies that a better understanding of all of them is required for reducing the uncertainty to ∼ 1-2%.
In this paper we focus on the computation of the two-loop electroweak correction to the production of the Higgs boson in gluon fusion. This contribution arises because gluons couple to electroweak vector bosons W and Z through a quark loop; a subsequent fusion of the electroweak bosons to the Higgs boson gives rise to electroweak-mediated ggH coupling. The quark loop receives contributions from both light and heavy quarks but the relatively small mass of the Higgs boson leads to a strong dominance of the light quark contributions. 1 The electroweak contributions to ggH have been evaluated analytically at leading (two-loop) order in Refs. [3][4][5][6]. Since the QCD corrections increase the leading, top-quark mediated, contribution to gg → H by almost a factor two, it is essential to understand if a similar enhancement is present in case of the electroweak corrections to gg → H . To clarify this issue, we need a computation of QCD corrections to the electroweak contribution to ggH . However, since the electroweak contribution starts at two loops, calculation of NLO QCD corrections requires dealing with three-loop diagrams with massive internal lines. Given the complexity of the required computation, one can try to simplify it by considering different kinematic limits: mixed QCD-electroweak corrections in the unphysical limit of a vanishingly small Higgs boson mass m Z,W m H were calculated in Ref. [2]. However, recent progress in the theoretical understanding of QCD effects in gg → H and continuous developments in the technology of multi-loop computations make it worthwhile and interesting to attempt an exact computation of the NLO QCD corrections to the electroweak contribution to ggH . In this paper, we make an important step in this direction by setting up a modern calculational framework for this problem that employs canonical bases for master integrals and differential equations, and computing the two-loop electroweak contribution to ggH to higher orders in the dimensional regularization parameter ε = (4 − d)/2. The knowledge of the two-loop amplitude to higher orders in ε is necessary for subtracting infrared and collinear singularities from the electroweak contributions to the gg → Hg inelastic process or, alternatively, for extracting the relevant finite remainder, defined by the Catani formula [7], from the three-loop mixed QCD-electroweak contribution to ggH amplitude.
Specifically, we derive the two-loop electroweak correction to gg → H through O(ε 2 ) and show that only GPLs up to weight five appear in this amplitude. We perform our calculation using the method of differential equations [8][9][10], augmented by the choice of a canonical basis of master integrals, introduced in Ref. [11]. 2 A canonical basis of master integrals is presented and the master integrals are calculated in terms of Goncharov's multiple polylogarithms (GPLs) [13][14][15]. In order to fix analytically all boundary conditions we make extensive use of the large mass expansion the PSLQ algorithm. This allows us to derive the expansion for the master integrals in the dimensional regularization parameter ε through weight six. From our calculation we can easily reproduce the O(ε 0 ) results which have been known for a long time in the literature [5].
The paper is organized as follows. We introduce the notation and discuss the structure of the scattering amplitude gg → H in Section 2. We describe the master integrals and the differential equations in Section 3. We explain how the boundary conditions can be fixed using the large mass expansion and outline the analytic continuation of GPLs, required to obtain results in physical kinematics, in Section 4. The gg → H finite part of the amplitude is given in Section 5. We provide constants of integration for the master integrals up to weight six in Appendix A. The explicit expressions for the master integrals up to this weight, and the gg → H amplitude through O ε 2 are available in the ancillary file.

Feynman diagrams and master integrals
We consider the electroweak correction to the gg → H amplitude mediated by a light-quark loop. The relevant contributions are shown in Fig. 1. The fermionic lines represent up, down, strange and charm quarks, that are taken to be massless. 3 The incoming gluons g 1 and g 2 are on-shell and carry momenta p 1 and p 2 , with color indices c 1 and c 2 and polarizations ε λ 1 (p 1 ) and ε λ 2 (p 2 ), respectively. The momentum of the Higgs boson is taken to be p 3 = p 1 + p 2 with p 2 3 = m 2 H = s. Thanks to gauge-invariance and parity constraints, the gg → H scattering amplitude is expressed in terms of a single form factor It is possible to extract the form factor F by contracting M c 1 c 2 λ 1 λ 2 with the projection operator The form factor F is a linear combination of integrals which depend on the scalar products between loop and external momenta, and on the scalar products of loop momenta between themselves. All the integrals in F are obtained starting from the two topologies shown in Fig. 2. At variance with Feynman diagrams in Fig. 1, when we consider topologies and master integrals, we use wavy (solid) lines to denote massless (massive) propagators, respectively. We take all momenta to be incoming, i.e.
The planar and non-planar integrals are parametrized as (2.5) where (2.7) In both cases, the propagator [7] is auxiliary; it is only needed for the parametrization of tensor integrals with (otherwise) irreducible numerators. Both planar and non-planar integrals are analytic functions in the complex plane of the variable s with the cut along the real axis starting at s = 0. This discontinuity corresponds to massless intermediate states in Feynman diagrams. At s ≥ 4M 2 , it also becomes possible to produce pairs of vector bosons on the mass shell; this leads to additional contributions to the discontinuities of I P,NP . We use the program Reduze2 [16] to express all integrals that appear in the evaluation of gg → H amplitude through master integrals (MIs). We also use the integration-by-parts reduction identities to derive the differential equations in s and M 2 satisfied by the master integrals.

Differential equations
We denote a vector of master integrals by I, a set of kinematic variables by x ∈ {s, M 2 }, and write the differential equations as It was conjectured in Ref. [11] that in many physically relevant cases a canonical basis of master integrals I exists with the property that the right hand side of the differential equation has a simple, factorized dependence on the regularization parameter ε. While the statement has not been rigorously proved, it is expected to be true at least for those cases that can be expressed in terms of Chen iterated integrals. The differential equations for the canonical basis assume the following form so that the iterative construction of I as series in ε becomes straightforward. General criteria to find candidate canonical integrals are given in Ref. [17] and, under certain conditions for ordinary differential equations, in Ref. [18]. We do not use this last algorithm in this paper; instead, we begin by constructing canonical bases for the simplest integrals in the set and gradually move to more complex ones, as described extensively in [19]. Since the original matrices A i are relatively sparse, this approach turns out to be quite practical for finding the canonical basis.
It is convenient to choose as independent variables the center of mass energy squared s and the dimensionless ratio ω = −M 2 /s. Since the dependence of any master integral on s follows uniquely from its mass dimension, we write master integrals as where a is an integer determined by the canonical mass dimension of the integral. The non-trivial information is contained in the functions I i (ω), which are dimensionless quantities. By choosing these functions to be appropriately re-scaled versions of the master integrals found by Reduze2 we can cast the system of differential equations for I(ω) in the following form The matrices A 0,1 are rational functions of ω, and have a block-triangular structure.
To construct a systematic expansion of master integrals in ε, it is convenient to change basis of master integrals and transform the system of differential equations into a canonical form. This requires A 0 to be removed. We can do that in a symbolic form by first solving the matrix differential equation where P ω is the path-ordering operator defined as By defining a new set of master integrals it is easy to see that F satisfies differential equations in the canonical form The non-trivial part of this procedure is to find the matrix Ŝ A 0 . A systematic way to do that, based on the Magnus exponentiation was suggested in Ref. [20]. Instead, we do that iteratively, using the sparse nature of the matrix A 0 and considering different blocks of A 0 separately. As an illustration, consider two integrals from the list of master integrals, I 2 and I 3 . Neglecting the matrix A 1 , we find that they satisfy the system of coupled differential equations d dω (3.10) Integrating this equation, we find where C 1,2 are the two integration constants. Since the above solution satisfies the system of differential equations for arbitrary C 1 , C 2 , the matrix Ŝ 0 satisfies the original differential equation and, therefore, can be taken to be a part of Ŝ A 0 . Finally, we compute F =Ŝ −1 0 I and find The system of differential equations for the integrals F 2,3 is then guaranteed to be in the canonical form. We find d dω We apply this procedure block by block, to the block-triangular matrix A 0 + εA 1 and obtain the canonical system of differential equations that we write in the following form dF(ω) = ε dB(ω)F(ω). (3.14) The canonical basis of the master integrals F(ω) reads ⎛ For the matrix B(ω) we obtain 4 The ω-independent matrix coefficients in the above equation are given by   4 The signs of the arguments of the logarithms are chosen to ensure that for positive ω the logarithms are real.
It is convenient to remove the square roots from the Eqs. (3.14), (3.15) by changing variables ω → y where 5 The differential equations (3.14) take the following form dF (y) = ε dC(y)F(y), (3.19) where the matrix C reads C(y) = C 1 log y + C 2 log(1 − y) + C 3 log(1 + y) + C 4 log(1 − y + y 2 ). (3.20) The y-independent matrices C 1,...,4 are It is straightforward to write the solution of the system of differential equations (3.19) as Taylor series in ε (3.23) Table 1 Different kinematic regions in s, ω and y.

Variable
Euclidean Minkowski, below threshold Above threshold 0 are integration constants that can not be fixed from the differential equations. Given the iterative structure of the solution, it can be written as a linear combination of the socalled Goncharov's polylogarithms (GPLs), also known as hyperlogarithms, defined as [13][14][15]21] , (3.24) where m w indicates the vector (m w , m w−1 ). The functions f (a; τ ) represent the integration kernels; for our system of differential equations they span the following set The last term in the set is quadratic in the integration variable; it is possible to re-write it in a usual linear form at the expense of introducing complex-valued poles. This last step is essential for numerical evaluation of the GPLs but it is not required to integrate the system of differential equations since, thanks to the linearity of the differential equations and the definition of the GPLs, This implies that we can perform the analytic integration using the symbol r and then, for the numerical evaluation of the final result, switch to r ± using Eq. (3.27). 6 In Table 1 the relations among s, ω and y, as well as the different kinematic regions, are summarized.

Boundary conditions and analytic continuation
Linear differential equations allow us to restore the dependencies of master integrals on ω up to a single constant. The constant can be fixed by computing the required integral at any point ω = ω 0 and comparing the result with the solution of the differential equation.
It turns out to be convenient to determine the boundary conditions by computing the master integrals in ω → ∞ or y → 1 limit. Since ω = −M 2 /s, this corresponds to M 2 → ∞ at fixed s; in this limit the integrals can be easily computed using the so-called large-mass expansion procedure [24]. The large-mass expansion procedure can be formulated as follows: consider two different scalings for each of the loop momenta k 1,2 ∼ M and k 1,2 ∼ √ s and, for a chosen scaling, systematically expand the integrand in Taylor series in all small variables. The set of small variables will differ from scaling to scaling; for example, if k 1 ∼ k 2 ∼ M, one Taylor expands in external momenta but if the scaling k 1 ∼ M, k 2 ∼ √ s is considered, one Taylor expands in the external momenta and k 2 . It is easy to see that, to leading order in s/M 2 , most of the master integrals are expressed in terms of two-loop tadpole integrals and, in some cases, in terms of products of one-loop three-and two-point integrals and one-loop tadpole integrals.
As an illustration, consider the non-planar master integral We are interested in determining its behavior in the y → 1 limit. By applying the large mass expansion, we find that the non-planar integral scales as in the large-M limit. This implies that the large-M limit of the F 9 non-planar master integral reads where we used relations among s, ω and y variables summarized in Table 1. F 9 vanishes as (1 − y) 3 as y goes to 1, and this information is sufficient to fix the constant of integration for this master integral. A similar analysis reveals that only three master integrals F 1 , F 2 , F 10 possess non-vanishing y → 1 limit. These limits are To fix the boundary conditions, we need to evaluate the GPLs of the form G(m, y), where m is composed of the elements of the set and the limit y → 1 is taken where possible. For weights higher than three, not all the Goncharov polylogarithms at y = 1 with the arguments from Eq. (4.7) can be analytically expressed in terms of canonical irrational numbers such as π and ζ(n). Nevertheless, we expect that the boundary constants are linear combinations of these irrational numbers; to find them we follow a numerical approach. We use our solution of the differential equations and the boundary conditions discussed above to find the integration constants numerically with high precision. 7 We then fit the resulting numerical value to a linear combination of π 2 and ζ(n) of a well-defined weight. For example, at weight two we only have π 2 , at weight three ζ(3), at weight four π 4 , at weight five ζ(5) and π 2 ζ(3) and at weight six π 6 and ζ(3) 2 . For each of the master integrals, we have achieved the matching of the numerical and the analytic results to at least 750 digits. Our final remark concerns the analytic continuation of the master integrals F(y). So far, we have studied them in the Euclidean region but we need them in the region where s = m 2 H > 0 and, yet, s < 4M 2 . The correct analytic continuation is achieved by replacing s → s + i0 at fixed M 2 . It is easy to see that this implies ω → ω + i0 and y → y + i0.
All the master integrals evaluated here have been compared for at least four different values of s/M 2 , both in the Euclidean and Minkowski region, to the numerical results obtained with the program Secdec [25]. In all cases we found excellent agreement.

Form factor for gg → H
The gg → H amplitude is described by a single form factor, as explained in Section 2. This form factor receives contributions from loops with W and Z bosons. The form factor is finite in four dimensions (ε → 0) and can be written as: We take the CKM matrix to be an identity matrix. The contributions of W bosons is computed in Eq. (5.1) taking into account first and second generations. The contribution of the Z boson is calculated for five massless quarks (u, d, s, c and b).
The expression for F in Eq. (5.1) has been compared with a previous calculation in Ref. [4] and agreement was found. The terms A 1 (y) and A 2 (y) are new. They can be found in the ancillary file provided with this paper.

Conclusions
We have presented a calculation of the mixed two-loop QCD-electroweak corrections mediated by massless quarks to the production of the Higgs boson in gluon fusion. We extended the known result for these corrections to two higher orders in the dimensional regularization parameter ε. This is one of the ingredients required for the computation of the NLO mixed QCD-electroweak corrections to gg → H amplitudes. We employed the method of differential equations, determined a canonical basis of master integrals and expressed all the relevant functions in terms of Goncharov polylogarithms. Finally, we used a mixed numerical and analytical approach, based on the PSLQ algorithm, in order to fix all necessary boundary conditions. This establishes the necessary framework to successfully address the calculation of the missing threeloop virtual contributions, whose calculation is ongoing.
In this appendix we present the boundary conditions for the master integrals defined in Eq. (3.23). The weights 0, 1 and 2 were determined analytically. For weights 4, 5 and 6 the results were obtained by fitting numerical results to an analytic Ansatz to at least 750 digits.