Two-loop massless QCD corrections to the $g+g \rightarrow H+H$ four-point amplitude

We compute the two-loop massless QCD corrections to the four-point amplitude $g+g \rightarrow H+H$ resulting from effective operator insertions that describe the interaction of a Higgs boson with gluons in the infinite top quark mass limit. This amplitude is an essential ingredient to the third-order QCD corrections to Higgs boson pair production. We have implemented our results in a numerical code that can be used for further phenomenological studies.


Introduction
The discovery of the Higgs boson at the Large Hadron Collider [1,2] is an important milestone in particle physics. It puts the Standard Model (SM) in a firm position to describe the dynamics of all the known elementary particles. Of course, there are several shortcomings in the SM which lead physicists to explore physics beyond the SM. There have been tremendous efforts to construct models that address these shortcomings and at the same time demonstrate rich phenomenology that can be explored at present and future colliders. All these culminated into dedicated experimental searches for hints of new physics which in turn constrain the parameters of beyond the SM scenarios [3].
By measuring the mass of the Higgs boson, one can predict the trilinear self-coupling in the Higgs sector of the SM. This is a crucial parameter that describes the shape of the Higgs potential. In order to better understand the Higgs sector and the nature of the electroweak symmetry breaking mechanism, it is important to measure this self-coupling independently. At hadron colliders, one of the potential channels that can probe this self-coupling is the production of a pair of Higgs bosons [4][5][6][7]. The dominant production channel in the SM is through the loop-induced gluon fusion subprocess [8,9]. At leading order (LO), this process involves two mechanisms, with the scattering amplitude for one of the them depending on the trilinear Higgs boson coupling. Since both mechanisms are loop-induced through heavy quarks and there is destructive interference between their respective amplitudes, the SM production cross section at LHC energies is only few tens of a femtobarn. In addition, a large and irreducible background [10][11][12][13][14][15] makes its detection an experimentally demanding task. Double Higgs boson production can receive substantial contributions from physics processes beyond the SM, and there are already several detailed studies indicating scenarios for a substantial increase in its production rate (see [16] and the references therein).
Theoretically, it is a challenging task to compute higher order QCD effects when taking into account the exact top quark mass dependence, since the Born-level contribution appears only at one loop. The first computation of next-to-leading order (NLO) QCD corrections was performed in the infinite top quark mass limit in [4]. In this limit, the top quark is integrated out, resulting in a field theory that contains effective operators coupling the Higgs field to the gluon field. These early results were then improved upon by considering various NLO contributions from finite top quark mass effects [17][18][19][20][21][22]. Recently, the full NLO corrections with exact top quark mass dependence could be completed [23,24], owing to technical progress in the numerical evaluation of two-loop integrals and amplitudes with internal masses. At next-to-next-to-leading order (NNLO) level, results are available only in the heavy top limit. The prediction at NNLO level in the soft plus virtual (SV) approximation can be found in [25], the leading top quark mass corrections were then included in [26], while in [27] the impact of the remaining hard contributions were studied. The relevant Wilson coefficients at NNLO were obtained in [28]. For the fully differential results at NNLO level, see [29][30][31]. By using a re-weighting approach, these fixed-order NNLO results for infinite top quark mass can be combined with the exact NLO top quark mass dependence to quantify [32] the top quark mass effects at NNLO. Effects of threshold resummation at next-to-next-to-leading logarithm (NNLL) level using soft collinear effective theory were obtained in [31,33].
Going beyond NNLO level in QCD is a challenging task owing to the technical difficulties involved in computing the loop integrals for the virtual subprocesses and the phase space integrals when there are real emissions. In this article we make a first step towards computing the third-order correction to the production of a pair of Higgs bosons in the gluon initiated channels. In particular we compute virtual amplitudes for the subprocess g + g → H + H, resulting from two effective operator insertions, at the two-loop level. The paper is structured as follows. In Section 2, we introduce the notation, describe the effective field theory that results in the limit of an infinite top quark mass, and discuss the different purely virtual contributions to Higgs boson pair production up to next-to-next-tonext-to-leading order (N 3 LO). Section 3 describes in detail the calculation of the two-loop amplitude for g + g → H + H, and the numerical evaluation of the results is discussed in 4. We conclude with an outlook on future applications in Section 5.

Higgs effective field theory
We compute the relevant amplitudes in an effective theory where the top quark degrees of freedom are integrated out. The effective Lagrangian that describes the coupling of one and two Higgs bosons to gluons is given by where G µν denotes the gluon field strength tensor, φ, the Higgs boson and v = 246 GeV is the vacuum expectation value of the Higgs field. Note that we have taken only those terms in the L ef f into account that are relevant for the production of two Higgs bosons in a gluongluon initiated process. The constants C H and C HH are the Wilson coefficients [28,[34][35][36][37][38] determined by matching the effective theory to the full theory and they can be expanded in powers of the renormalized strong coupling constant a s = g 2 and where n f is the number of light flavors, m t is the M S top quark mass at scale µ R and N = 3 is fixed for QCD.

Kinematics
Consider the production of a pair of Higgs bosons in the gluon fusion subprocess, where p 1 and p 2 are the momenta of the incoming gluons, and p 3 and p 4 the momenta for the outgoing Higgs bosons, respectively. The Mandelstam variables for the above process are given by They satisfy s + t + u = 2m 2 h where m h is the mass of the Higgs boson. In the following, we describe the computation of the one-and two-loop QCD corrections to the amplitude given in Eq. (2.4). We find that it is convenient to express this amplitude in terms of the dimensionless variables x, y and z

Tensors and projectors
Using gauge invariance, the amplitude can be decomposed in terms of two second rank Lorentz tensors T µν i with i = 1, 2 as follows [8]: where the tensors are given by with p 2 T = (tu − m 4 h )/s. In color space, the amplitude is diagonal in the indices (a, b) of the incoming gluons. The scalar functions M i can be obtained from M µν ab by using appropriate projectors as follows where the projectors in d dimensions are given by, 2.4 Diagrams to O(a 4 s ) When considering higher order massless QCD corrections to the g + g → H + H amplitudes in the effective theory, we encounter two topologically distinct classes of subprocesses we call Class-A and Class-B hereafter. We perform an expansion in a s to include all the contributing diagrams.
• Class-A, see Fig. 1, contains diagrams where two Higgs bosons couple to each other and to gluons. They either couple to the gluons directly through a C HH Wilson coefficient (left-hand column of Fig. 1), or through a Higgs boson propagator and the C H Wilson coefficient (right-hand column of Fig. 1). The latter diagrams are linearly proportional to the triple Higgs coupling λ.
• Class-B, see Fig. 2, contains diagrams where Higgs bosons couple to two gluons through the effective vertices proportional to C H , but do not couple to each other.
Both Wilson coefficients C H and C HH start at order a s . Consequently, to LO in a s only Class-A diagrams contribute. Beyond LO, that is from order a 2 s onwards, the class-A diagrams are only of form factor type and the results for class-A to a 4 s can be readily obtained from the three loop form factor [39,40] that appears in purely virtual contributions to single Higgs boson production. The class-B diagrams start contributing from order a 2 s , with results only available up to order a 3 s [4]. In the following, we will complete the a 4 s contributions to the g + g → H + H amplitude, by computing the class-B diagrams to this order, which amount to their two-loop corrections.
In general, the scalar amplitudes M i can be written as a sum of amplitudes resulting from the two classes A and B Since the M A i are proportional to the Higgs boson form factor, they can be expressed as (2.14) The amplitude M A 2 is identically to zero to all orders in perturbation theory due to the choice of the tensorial basis. The form factors F (j) (d) for j = 1, 2, 3 are known in the literature [39,40].
In this article, the amplitudes of class-B are presented up to two loop level in perturbative QCD. At each order, the amplitude contains a pair of vertices resulting from the first term of the effective Lagrangian L ef f and hence will be proportional to the square of the Wilson coefficient C H (a s ), expanded to the desired accuracy in a s . Beyond leading order, the one-and two-loop diagrams are not only ultraviolet (UV) divergent but also infrared (IR) divergent resulting from soft and collinear regions of the loop momenta. We use dimensional regularization to treat both UV and IR divergences and all the divergences show up as poles in , where the space time dimension is d = 4 − 2 .

Ultraviolet renormalization and operator mixing
The bare strong coupling constant in the regularized theory is denoted byâ s which is related to its renormalized counter-part bŷ where S = exp [(ln 4π − γ) ] with γ ≈ 0.5772... the Euler-Mascheroni constant. The beta function coefficients β 0 and β 1 are given by for the SU(N) color factors we have Besides coupling constant renormalisation, the amplitudes also require the renormalisation of the effective operators in the effective Lagrangian, Eq. (2.1). Both composite operators that appear in our one-and two-loop amplitudes can develop UV divergences and thus have to undergo renormalisation, as derived in detail in [41]. In particular, a new renormalisation constant Z L 11 is needed in a counter term proportional to G µν G µν φφ to renormalize the additional UV divergence resulting from amplitudes involving two G µν G µν φ type operators starting from 2-loop order in class-B amplitudes. If we denote the amplitudes computed in the bare theory byM B i , then the relation between these bare amplitudes and the UV renormalized ones is given by The overall renormalisation constant [42][43][44] is given by where and Z L 11 is given by [41], The UV renormalized amplitude M B i can be expanded in powers of a s up to the two-loop level as follows: where, In summary, the UV divergences that appear at the one-and two-loop level can be removed using coupling constant renormalisation through Z and the overall operator and the contact renormalisation constants, Z O and Z L 11 respectively.

Infrared factorization
The resulting UV finite amplitudes will contain divergences of infrared origin, which remain as poles in the dimensional regularization parameter . These will cancel when combined with the real emission processes to compute observables. While these divergences disappear in the physical observables, the amplitudes beyond leading order demonstrate a very rich universal structure in the IR region. Catani [45] predicted IR divergences for n-point two-loop amplitudes in terms of certain universal IR anomalous dimensions, exploiting the iterative structure of the IR singular parts in any UV renormalized amplitudes in QCD. These could be related [46] to the factorization and resummation properties of QCD amplitudes, and were subsequently generalized to higher loop order [47,48]. Following [45], we obtain g ( ) are the IR singularity operators given by T F n f , It is known that the terms that become finite or vanish as goes to zero, i.e., O( α ), α ≥ 0 in the subtraction operators I (1) g and I (2) g are arbitrary and they define the scheme in which these IR divergences are subtracted to obtain IR finite parts of amplitudes, M B,(j),f in i . These scheme-dependent terms in the finite part of virtual contributions will cancel against those coming from the soft gluon emission subprocesses at the observable level. The only scheme dependence that will be left in a physical subprocess coefficient function is then due to the subtraction of collinear initial state divergences through mass factorization, parametrized by a factorization scale µ F .

Calculation of the Amplitude
For the amplitudes of class-B, we needed to consider only those diagrams which involve a pair of vertices resulting from the first term of the effective Lagrangian and hence all the amplitudes are proportional to C 2 H . These Feynman diagrams up to two-loop level were obtained with help of the package QGRAF [49]. There are 2 diagrams at tree level, 37 at one loop and 865 at two-loop order in perturbation theory. The output from QGRAF was then used for further algebraic manipulations involving traces of Dirac matrices, contraction of Lorentz and color indices, using two independent sets of in-house routines based on a symbolic package FORM [50]. The entire manipulations were performed in d = 4 − 2 dimensions and most of the algebraic simplifications were done at this stage. We used the Feynman gauge throughout and hence allowed ghost particles in the loops. External ghosts are not required due to the transversal nature of the tensorial projectors Eq. (2.11).
At this stage, we obtain a large number of Feynman integrals with different sets of propagators and each containing scalar products of the independent external and internal momenta. Using the REDUZE2 package [51], we can identify the momentum shifts that are required to express each diagram in terms of a standard set of propagators (called auxiliary topology). The auxiliary topologies in the two-loop corrections to the class-B process are identical to those in equal-mass on-shell vector boson pair production at this loop order. They are described in [52] and were used to compute the two-loop corrections to qq → V V in [53,54]. They were subsequently extended towards non-equal gauge boson masses [55][56][57][58].
It is well known that the resulting Feynman integrals are not all independent and hence they can be expressed in terms of fewer scalar integrals, called Master Integrals (MIs) by using integration-by-parts (IBP) identities [59,60]. Further simplifications can be done by exploiting the Lorentz invariance of the integrands, resulting in Lorentz invariance (LI) identities [61]. These identities can be solved systematically using lexicographic ordering (Laporta algorithm, [62]) to express any Feynman integral in terms of master integrals. These are implemented in several specialized computer algebra packages, for example AIR [63], FIRE [64], REDUZE2 [51,65] and LiteRed [66], to perform suitable integral reductions such that one ends up with only MIs. We performed two independent reductions of the integrals in the two-loop class-B amplitude, one based on the Mathematica based package LiteRed [66] and the other based on REDUZE2 [51]. Counting kinematical crossings as independent integrals, we can express the one-loop amplitude in terms of 10 master integrals, while the two-loop amplitude contains 149 master integrals. These master integrals are two-loop four-point functions with internal massless propagators and two massive external legs of equal mass. They were computed analytically as Laurent series expansion in in [52,67].
These MIs were then expressed in terms of generalized harmonic polylogarithms. An alternative functional basis can be obtained in terms of logarithms, polylogarithms Li n≤4 and the multiple polylogarithm Li 2,2 by matching the original expression at the symbol level [52]. We use the master integrals in this latter representation. Substituting the MIs from [52,67], we obtain the bare amplitudesM    (right panel). The behavior of the amplitudes close to the production threshold, x = 0, is shown in the insets. We see that the finite parts of the two-loop amplitude shows stable behavior, and they display a non-trivial dependence on the process kinematics.
In the numerical evaluations, the large rational coefficients of the classical polylogarithms can introduce numerical instabilities in case we do not demand high enough precision. In particular, there are large cancellations between the numerator and denominator of rational functions. Therefore, we evaluate the polylogarithms at double, and the rational coefficients at even higher precision.

Discussion and Conclusions
The two-loop massless corrections to the g + g → H + H amplitude derived above complete the set of purely virtual amplitudes required for the prediction of the N 3 LO corrections to Higgs boson pair production in gluon fusion, in the infinite top quark mass limit. All other amplitudes relevant at this order are either (class-A) known already from the calculation of inclusive gluon fusion Higgs boson production at N 3 LO [68,69] or (class-B) amount to one-loop and tree-level amplitudes that can be computed using automated tools. The combination of these amplitudes into a fully differential N 3 LO calculation of Higgs boson pair production does still require substantial advances in the techniques for handling in-frared singular real radiation configurations at this order, with first steps being taken most recently [70,71].
More imminent applications of the newly derived results to Higgs boson pair production are the computation of fixed order soft-virtual corrections to the total cross section or of the hard matching coefficients in the resummation of corrections at low pair transverse momentum.
In this paper, we have computed all virtual amplitudes that contribute to the production of a pair of Higgs bosons from the gluon-gluon initiated partonic processes at order a 4 s . The calculation is performed in an effective field theory where the top quark is integrated out, and all other quarks are massless. The exact calculation of top quark mass effects is currently out of reach at this order, but reweighting procedures allow to reliably quantify these effects [32]. We deal with two classes of amplitudes separately, named class-A (one effective operator insertion) and class-B (two effective operator insertions). The amplitudes of class-A can be related to the gluon form factor which is already known up to three loop order while amplitudes of class-B were known previously up to one loop. Our explicit computation of the two-loop corrections to the class-B amplitudes now completes the perturbative expansion of the g + g → H + H amplitude to order a 4 s . We observe that the pole structure of the amplitude is in agreement with predictions from infrared factorization, and provide (as ancillary files with the arXiv submission of this article) a numerical code to evaluate its finite remainder piece. The newly derived amplitudes open up opportunities for a new level of precision phenomenology predictions in Higgs boson pair production.