Double Higgs boson production at NLO in the high-energy limit: complete analytic results

We compute the NLO virtual corrections to the partonic cross section of $gg\to HH$, in the high energy limit. Finite Higgs boson mass effects are taken into account via an expansion which is shown to converge quickly. We obtain analytic results for the next-to-leading order form factors which can be used to compute the cross section. The method used for the calculation of the (non-planar) master integrals is described in detail and explicit results are presented.


Introduction
Higgs boson pair production is a promising channel to investigate the self interaction of the Higgs boson. Although it is very challenging from the experimental point of view it is expected that after the high-luminosity upgrade of the LHC constraints on the Higgs boson tri-linear coupling will be able to be obtained. In order to determine whether or not the Higgs sector is Standard Model-like it is therefore important to have the higher order corrections to double Higgs boson production under control. A further building block towards this goal is considered in this paper by providing analytic results in the high-energy limit.
Higgs boson pairs are predominantly produced by the gluon-fusion channel and in the recent years a number of higher order corrections have been computed to gg → HH, both for the total cross section and for differential distributions. We refrain from providing a detailed review but refer to Ref. [1] where several recent results are combined and approximate next-to-next-to-leading order (NNLO) expressions are constructed.
From the technical side the main new ingredients from this paper are analytic results for the two-loop non-planar master integrals for gg → HH which, in combination with the findings of Ref. [2], allows one to obtain the next-to-leading order (NLO) amplitude for this process in the high-energy limit. This complements the NLO results obtained from the large top quark-mass expansion [3,4], from the threshold expansion [5] and from an expansion for small Higgs transverse momentum [6]. Furthermore, it provides an important cross check and eventually an alternative approach to the exact result obtained in Refs. [7,8] using a numerical approach. Recently it has been suggested to expand the gg → HH amplitude only in the Higgs boson mass but keep the dependence on the kinematic invariants and the top quark mass [9]. This also leads to simpler expressions, however, one still has to solve integrals involving three scales.
As described in more detail in Subsection 2.3 we perform an expansion in the Higgs boson mass. This means that we use the kinematics defined in Eqs. (1) and (2) when evaluating the amplitude, but before evaluating the Feynman integrals we set m H = 0 and obtain the following variables which are relevant for the computation of the integrals 1 s = 2q 1 · q 2 , t = 2q 1 · q 3 , u = 2q 2 · q 3 = −s − t .
Thus the integrals will only depend on the variables s, t and m 2 t , and when computing them we further assume that m 2 t s, |t|. It is convenient to introduce the scattering angle θ of the Higgs boson in the center-of-mass frame which leads to the following relation in terms of these variables, Due to Lorentz and gauge invariance it is possible to define two scalar matrix elements M 1 and M 2 as where a and b are adjoint colour indices and the two Lorentz structures are given by with The Feynman diagrams involving the triple Higgs boson coupling only contribute to A µν 1 and, thus, it is convenient to decompose M 1 and M 2 into "triangle" and "box" form factors with where T = 1/2 and µ is the renormalization scale. We furthermore define the expansion in α s of the form factors as and similarly for M i . Throughout this paper the strong coupling constant is defined with six active quark flavours. Note that the form factors are defined such that the one-loop colour factor T is contained in the prefactor X 0 .
The main results of this paper can be summarized as follows: • We compute all planar (see Ref. [2]) and non-planar master integrals for gg → HH in the limit m 2 t s, |t| and m H = 0.
• We obtain analytic results for the NLO form factors which are used to parametrize the process gg → HH. These results can be used to construct the partonic cross section in the high-energy limit.
• We perform an expansion in the Higgs boson mass which converges very quickly in the region in which our result is valid. Here the relevant expansion parameter is m 2 H /(2m t ) 2 ≈ 0.13. In fact, at LO very good agreement with the exact result is obtained after including only the quadratic term.
• We provide input for the Padé method suggested in Ref. [5] for the process gg → HH.
The remainder of the paper is organized as follows: in Section 2 we describe the method we used to compute the amplitude and master integrals and discuss the ultraviolet and infra-red structure of the amplitude. Additionally, we explain our approach to obtain an expansion of the amplitude in the Higgs boson mass. Afterwards, in Section 3 we discuss our results for the form factors and present both analytic and numerical results. Our conclusions are presented in Section 4. In Appendix A we define our non-planar master integrals and provide graphical representations, and in Appendix B we describe the basis change which facilitates the computation of the boundary conditions.

Non-Planar Master Integrals
Details on the calculation of the NLO amplitude gg → HH and in particular on the reduction to master integrals can be found in Ref. [2]. An algorithm is provided which minimizes the number of families and yields 10 one-loop and 161 two-loop master integrals. At one-loop order all integrals are planar. At two-loop order we obtain 131 planar integrals, which are discussed in detail in [2], and 30 non-planar master integrals. The computation of the latter, which is based on differential equations, is described in the following. A detailed description of the computation of the boundary conditions can be found in Ref. [10].
Graphical representations of the non-planar master integrals can be found in Appendix A, see Fig. 10. Note that the 30 non-planar master integrals can be divided into two sets; 16 integrals for which actual calculation is needed, and 14 integrals which can be obtained with the help of crossing relations. Among the 16 integrals there are 9 seven-line and 7 six-line master integrals (cf. Fig. 10). We have computed all 30 integrals directly, however, and use the crossing relations as a cross check.
The main idea to obtain the high-energy expansion is the same as for the planar integrals; for each integral we make an ansatz which reflects the expected functional form of the expansion. This ansatz is inserted into the differential equation obtained by differentiating the master integrals with respect to m t . It is a new feature of the non-planar integrals that the ansatz requires both odd and even powers in m t (see, e.g., Ref. [11]) whereas for the planar integrals just even powers were sufficient. Note that due to the structure of the differential equations w.r.t. m t the even-and odd-power ansatz terms decouple and can be treated independently.
For the computation of the planar master integrals in Ref. [2] we followed two approaches.
In the first we computed the boundary integrals in the limit m t → 0 for a fixed values of s and t and used differential equations in t to reconstruct the t-dependence (still in the limit m t → 0). The differential equations in m t were then used to construct the expansion terms in the high-energy limit. In the second approach t-dependent boundary conditions were computed using asymptotic expansion and Mellin-Barnes techniques. For the non-planar master integrals we follow only this second approach, which can be used largely without modification. There are a few peculiarities, however, mainly connected to the presence of additional regions in the asymptotic expansion. This requires an extension of the method, which is described in detail in Ref. [10]. We note that this method has many algorithmic elements, which are certainly more generally applicable beyond the computation of the amplitude described in this paper.
For the computation of the non-planar master integrals (at least for those with seven lines) it is crucial to choose a basis in which the master integrals do not contain poles in their prefactor in the amplitude. This guarantees that only the constant ( 0 ) terms of the master integrals are required, which contain objects with transcendental weight of at most four. We obtain such a basis by replacing dotted propagators, which are present in the original basis chosen by FIRE [12], with numerator scalar products, see Appendix B for details. Our choice of basis for the 4×4 and 5×5 coupled blocks are given in Appendix A.
An important cross check of our results is provided by the explicit expressions from Ref. [11] where NLO corrections to Higgs plus jet were considered in the high-energy limit. Unfortunately, it is not possible to simply take over the results from [11] since our amplitude has single poles in in the master integral coefficients if we use their integral basis. This means that we would require O( ) terms of these master integrals, which are not known. We nonetheless compare our results to those of [11], to the orders possible, and they agree. Note that the results of [11] are given in terms of kinematics where t > 0, s < 0, u < 0, so require analytic continuation to our physical kinematics. We have also successfully compared our "triangle" master integrals to Ref. [13]. All of our non-planar results could additionally be cross checked numerically using both FIESTA [14] and pySecDec [15]. Analytic results for the master integrals can be found in the ancillary file to this paper [16].
In order to illustrate the structure of our results we present the explicit expression for the pole part of G 51 (1, 1, 1, 1, 1, 1, 1, 0, −2) (see Appendix A for the definition of this integral) in the limit m t → 0. We include the first and second terms of the small-m t expansion, and set s = 1. The s dependence can easily be restored by making the replacements t → t/s, m t → m t / √ s and multiplying by an overall factor of (−µ 2 /s) /s to fix the mass dimension of the integral. Our result reads where H a (x) denote Harmonic Polylogarithms as defined in [17]. Note that here one observes that the leading term is proportional to 1/m t ; as explained above, these odd powers of m t are particular to the non-planar master integrals and do not appear in the planar results of Ref. [2].
For illustration we show in Fig. 1 the real and imaginary part of the 0 term of G 59 (1, 1, 1, 1, 1, 1, 1, 0, 0) as a function of √ s for θ = π/2. We include successively higher orders in the m t expansion, which improves the agreement with the exact result shown as dots (pySecDec) and crosses (FIESTA). We want to stress that the odd m t terms are numerically significant and are needed to reach the agreement. It is, furthermore, interesting to mention that after including an odd expansion term the agreement gets worse and improves only after adding also the next even m t term. Thus, if the combination of the m 2n−1 t and m 2n t terms are considered a steady improvement is observed. We have obtained similar plots for all 30 non-planar master integrals.

Ultraviolet and Infrared Divergences
The bare two-loop expressions for the form factors are both ultraviolet and infrared divergent. We take care of the ultraviolet poles by renormalizing the top quark mass in the on-shell scheme and the strong coupling constant in the MS scheme using standard one-loop counterterms. Since we consider the high-energy region we renormalize α s with six active quark flavours.
Note that after top quark-mass renormalization the C F colour factor of the two-loop form factors are finite. However, there are still infrared divergences in the C A colour factor. We have checked that they agree with the poles predicted in Ref. [18]. We thus construct the (infrared finite) soft-virtual corrections as where F (1),IR is one of the ultraviolet-renormalized, but still infrared divergent, form factors introduced in Eq. (8). K (1) g can be found in Ref. [18]. For the normalization introduced in Eq. (10) it is given by where γ E is Euler's constant. Note that since infrared and ultraviolet divergences are regulated with the same parameter and since scaleless integrals are set to zero, the poles in the terms proportional to n f from Eq. (13) cancel against the counterterm contribution induced by the α s renormalization. However, finite terms proportional to log(µ 2 /(−s−iδ)) and the LO result remain. We thus cast F (1) in the form with β 0 = 11C A /12 − T n f /3. Only F (1),C F and F (1),C A contain new information and thus only these will be discussed in the following. Note that F (1),C F and F (1),C A are independent of µ.

Expansion in m H
In Ref. [2] the calculation has been performed for a massless Higgs boson which constitutes a good approximation since the relevant expansion parameter m 2 H /(2m t ) 2 ≈ 0.13 is sufficiently small. In the present calculation we incorporate finite Higgs mass effects via an expansion in m 2 H /m 2 t . For our process the dependence on the Higgs boson mass is analytic, i.e., there are no log(m H ) terms in the limit m H → 0 since the Higgs boson couples only to the massive top quark. It is thus possible to perform a simple Taylor expansion (in contrast to a more involved asymptotic expansion) which we have implemented as follows: • We generate the amplitude using the kinematics for a finite Higgs boson mass as given in Eq. (2). In particular, we use m H = 0 in the projectors onto the individual tensor structures and express the amplitude as a linear combination of scalar integrals, which depend ons,t, m t and m H .
• Next, the pre-factors of the scalar integrals are expanded about m 2 H = 0. Expressions for the Taylor expansion of the scalar integrals themselves are constructed using LiteRed's [19,20] derivative function Dinv.
• At this point the amplitude is expressed as a linear combination of scalar integrals which only depend on s, t and m t ; m H only appears in their prefactors. All scalar integrals can be mapped to one of the families defined in Ref. [2]. We can thus use the same procedure to obtain the reduction tables with the help of FIRE 5.2 [12] and FIRE 5.7. 2 Note, however, that the number of scalar integrals is significantly increased; at two-loop order one has about 25,000 scalar integrals to reduce to master integrals, for the m 0 H contribution. A further 70,000 integrals were reduced in order to produce differential equations for the master integrals. For the m 2 H and m 4 H contributions, one must reduce an additional 123,000 and then 457,000 integrals respectively.
At one-loop order we performed an expansion up to O(m 4 H ). We show below that the contribution from the m 4 H terms is very small in the kinematic region where the small-m t expansion is valid (see the discussion regarding Fig. 4). For this reason, at two loops we consider only the m 2 H terms of the expansion, and do not perform the computationally expensive reduction of the above-mentioned additional 457,000 scalar integrals to masters.

Numerical Results for the Form Factors
In the following we discuss the √ s dependence of the form factors at one-and two-loop order. If not stated otherwise we use m t = 173 GeV and m H = 0 or m H = 125 GeV for the top quark and Higgs boson masses, respectively.

One-Loop Form Factors
In Figs The triangle form factor (Fig. 2) is approximated very well by the asymptotic results. The solid and dashed curves lie on top of each other for the entire √ s region above the threshold 2m t . In Fig. 3 one observes that for F and √ s ≈ 500 GeV for the real and imaginary parts, respectively. Below these energies the curves diverge from each other. In general, one can trust the expansions in the regions where successive approximation orders agree with each other. Due to the very marginal improvement of the m 16 t approximation relative to the m 14 t approximation, we expect that computing higher order terms of the expansion will not improve the approximation, and that the small-m t expansion has a finite radius of convergence.
In Fig. 4 we consider the m H dependence of the partonic cross section for θ = π/2. Since this quantity is non-zero for the whole √ s range we can consider the ratio of our approximations to the exact result, evaluated for m H = 125 GeV. For √ s = 1000 GeV one observes that the m 0 H approximation (purple dashed curve) reproduces the m H = 0 exact curve well, and that these curves deviate from the m H = 125 GeV exact curve by about 2%. Including m 2 H terms in the approximation is sufficient to describe the m H = 125 GeV exact curve very well. Including also m 4 H terms provides a very small correction. Based on this observation, we compute m 2 H contributions to NLO quantities but not contributions proportional to m 4 H . We want to remark that the numerical values for Fig. 4 have been  √ s for θ = π/2. The notation is the same as in Fig. 2.
obtained by using the relation and performing a consistent expansion in m H . In this way we obtain the form factors as a function of s, θ and m H .

Two-Loop Form Factors
For simplicity we set m H = 0 in the following discussion of the two-loop corrections. The two-loop form factors F (1),C F and F (1),C A are shown in Figs. 5, 6 and 7, where approximations including terms up to m 14 t and m 16 t are shown. For the triangle form factor (Fig. 5) the approximations can be compared to the exact result from Ref. [21] and, as at one-loop order, good agreement is found down to √ s ≈ 2m t . For the box form factors no exact results are available. For the C F contribution we observe a similar behaviour as at one-loop order; the two highest expansion terms agree down to √ s ≈ 800 GeV and √ s ≈ 500 GeV for real and imaginary parts respectively, and diverge for smaller √ s values. For the C A contribution the convergence properties for real and imaginary part are reversed; we find agreement of the highest expansion terms down to values √ s ≈ 750 GeV and √ s ≈ 800 GeV for the real and imaginary parts respectively.
In order to illustrate the size of the m 2 H terms we show in Tab. 1 for two values of √ s the relative corrections for the real part of the NLO box form factors 3 as compared to the m H = 0 result. One observes corrections up to a few percent, in agreement with the one-loop results discussed in Fig. 4.  as a function of the partonic center-of-mass energy for θ = π/2. The same notation as in Fig. 2 is adopted. In these plots the exact curves are not known.  as a function of the partonic center-of-mass energy for θ = π/2. The same notation as in Fig. 2 is adopted. In these plots the exact curves are not known.

θ Dependence of the Form Factors
In the previous subsection we have chosen θ = π/2 where t = −s/2, i.e., the absolute value of t is maximal and our approximation works best. In Fig. 8 we show the "box1" form factors as a function θ with 0 ≤ θ ≤ π/2. Symmetric results are obtained for π/2 ≤ θ ≤ π.
The form factors F are shown in the three columns and the rows correspond to three different choices of √ s: 800 GeV, 1000 GeV and 1500 GeV. We show both the real and imaginary part for expansion depths m 14 t and m 16 t and assume that our approximation is good if the two curves agree. At one-loop order we can compare to the exact result.
In the case of F the θ range where our approximation works well is significantly smaller for √ s = 800 GeV.
However, for √ s = 1000 GeV and √ s = 1500 GeV similar results are obtained as for F (0) box1 and F (1),C F box1 . Fig. 9 shows analogous results to Fig. 8 for the "box2" for factors. We observe very similar convergence properties.

Conclusions
We consider Higgs boson pair production in gluon fusion at NLO and compute analytic results in the high-energy limit where the squared top quark mass is much smaller than s and |t|. We compute analytic results in this limit for all non-planar master integrals, which complement the results for the planar integrals, already presented in Ref. [2]. Analytic expressions for the master integrals are provided in an ancillary file to this paper Our expressions allow for a fast numerical evaluation of the form factors and thus provide an alternative to the exact, numerically expensive calculation of Ref. [8] in the high-energy region of the phase-space. It is in particular tempting to combine our results with other approximations [3][4][5][6] to cover the full phase space. Such investigations are the subject of ongoing research. with