Efficient computation of two-loop amplitudes for Higgs boson pair production

We present a precise and efficient computation of the two-loop amplitudes entering the Higgs boson pair production process via gluon fusion. Our approach is based on the small-Higgs-mass expansion while keeping the full dependence on the top quark mass and other kinematic invariants. We compare our results to the up-to-date predictions based on a combination of sector decomposition and high-energy expansion. We find that our method provides precision numeric predictions in the entire phase space, while at the same time is highly efficient as the computation can be easily performed on a normal desktop or laptop computer. Our method is valuable for practical phenomenological studies of the Higgs boson pair production process, and can also be applied to other similar processes.


INTRODUCTION
Higgs boson pair production via gluon fusion is one of the most important scattering processes yet to be discovered at the Large Hadron Collider (LHC). It provides information on the Higgs trilinear self-coupling (which is an important parameter in the Higgs potential, relevant for electroweak symmetry breaking and vacuum stability), or more generically, on the Wilson coefficients of the relevant effective operators encoding physics at higher scales beyond the standard model (SM) [1][2][3][4]. Hence, precision theoretical predictions for this process (both within the SM and within the effective field theory of beyond-the-SM physics) are highly important for a future extraction of these crucial information from experimental measurements of the differential cross sections.
At the leading order (LO) in the SM couplings, the process goes mainly through a top quark loop, with a small contribution from a bottom quark loop. The relevant LO Feynman diagrams can be categorized into two types: the triangle diagrams which are sensitive to the Higgs trilinear self-coupling, and the box diagrams which are sensitive to the top quark Yukawa coupling. In the SM, these two types of diagrams numerically cancel against each other, giving a rather small cross section [5,6]. New physics beyond the SM will modify these couplings and may spoil this accidental cancellation. In that case the production rate will be significantly enhanced. An early observation of this process at the LHC is then a clear evidence of new physics lying not too far above the electroweak scale.
Due to the importance of the Higgs boson pair production process, its precision theoretical prediction has been an active research area in the past few years. The LO results have been available since about 30 years ago [5,6]. The calculation of the next-to-leading order (NLO) corrections in quantum chromodynamics was quite challenging, due to the appearance of difficult twoloop Feynman integrals whose analytic expressions are still unknown so far. Because of that, people have at-tempted different methods to approximately calculate these two-loop integrals. The most straightforward way is to compute the integrals numerically using the method of sector decomposition [7][8][9]. This approach is rather resource-consuming and has only become possible recently with the help of modern computer clusters [10][11][12][13]. Besides the automated computation using sector decomposition, Ref. [14] calculated these integrals by carefully manipulating the Feynman-parameterized integrals supplemented by numeric integration and extrapolation.
Apart from the purely numerical method, it is also possible to approximate the integrals by working with a series expansion in terms of one (or more) small parameter(s) in certain kinematic limits. The simplest one is the 1/m t expansion [15][16][17][18][19][20][21], whose leading power is dubbed the Higgs effective field theory (HEFT) in the literature. People have also studied the small transverse momentum (p T ) expansion [22] and the high-energy expansion [23][24][25]. All these expansions work in a specific kinematic region, and cannot be directly applied to the whole phase space. Recently, the high-energy expansion has been combined with the results of sector decomposition via the Padé approximation, aiming at a good description of the kinematic distribution from low to high energies [26]. Later on we will compare our results to the last one in more detail.
In [27], some of the authors of this Letter proposed a novel method to calculate loop integrals with a heavy top quark loop and lighter external particles such as the Higgs boson and the vector bosons in the SM. In this Letter, we apply this method to calculate the complete NLO virtual corrections to the Higgs boson pair production process via gluon fusion. We show that our method provides accurate numeric predictions for the two-loop amplitudes in the entire physical phase space. At the same time, the numeric performance of our method is rather efficient, and can be accomplished on a desktop computer or even a laptop. The precision and efficiency of our approach makes it convenient for phenomenological studies. Our method can also be implemented for other 2 → 2 processes with a heavy quark loop, such as gg → ZZ, gg → ZH and gg → Hj.

CALCULATION OF THE AMPLITUDES
We consider the partonic process g µ a (p 1 ) + g ν b (p 2 ) → H(p 3 ) + H(p 4 ), where a and b are color indices. The amplitude can be written as where G F is the Fermi constant. The Mandelstam variables are defined asŝ = ( H , which satisfyŝ +t 1 +û 1 = 0. The tensor structures in Eq. (1) are given by where p T denotes the transverse momentum of the Higgs boson with respect to the beam axis and can be expressed in terms of the Mandelstam variables as H . The corresponding form factors F 1,2 are functions ofŝ,t 1 , m 2 t and m 2 H , and can be expressed as power series in the strong coupling α s : The NLO terms F (1) i receive contributions from two-loop diagrams with a top quark loop and virtual gluon exchanges, which are the main concern of this work.
The two-loop Feynman integrals contain both ultraviolet (UV) and infrared (IR) divergences. The UV divergences are removed by renormalization, for that we adopt the on-shell scheme for m t and the five-flavor MS scheme for α s . The IR divergences will be eventually cancelled by the real-emission contributions and the renormalization of the parton distribution functions (PDFs). We follow the dipole subtraction method [28] to remove the IR divergences from the two-loop amplitudes, and define the finite part of the NLO form factors as [26] where F (1),ren i denote the UV-renormalized form factors, Here and below, one should remember that the physical branch is chosen byŝ →ŝ + iδ.
In order to describe the contributions of the NLO form factors to the partonic cross section, we follow the notation of Refs. [26] and define the finite part of the two-loop virtual corrections as with The quantity V fin involves various ingredients, among which the one-loop diagrams, the two-loop one-particle reducible diagrams, and the two-loop triangle diagrams can be calculated analytically [5,6,21,[29][30][31]. The only difficult part is the two-loop one-particle irreducible box diagrams. We denote those contributions as F (1) i,box , and employ the small-Higgs-mass expansion proposed in [27] for their calculation.
We generate the two-loop virtual diagrams using FeynArts [32], and manipulate the resulting amplitudes with FeynCalc [33] and FORM [34]. The form factors are then expressed in terms of linear combinations of scalar integrals, which are functions ofŝ,t 1 , m 2 t and m 2 H . These scalar integrals are rather complicated, and cannot be evaluated analytically using existing methods. In fact, even their reduction to master integrals is highly nontrivial [10,11]. To deal with them, we observe that these integrals, and hence the form factors, can be expanded as Taylor series in m 2 H : .
(8) When acting on the loop integrals, the partial derivative operator can be exchanged with the integration, which then acts on the integrands as After expansion, we need to deal with two-loop box diagrams with massless external legs. They can be categorized into 6 integral topologies shown in Fig. 1. We use program packages FIRE6 [35,36] and LiteRed [37] to perform the integration-by-parts (IBP) reduction of these loop integrals. There are 29, 32, 37, 27, 54, 36 master integrals in topologies A-F, respectively. The first five topologies have been discussed in [27,38,39]. They can be solved using the method of canonical differential equations [40] in terms of iterated integrals [41]. Most of these iterated integrals can be expressed as generalized polylogarithms (GPLs) [42]. We evaluate these GPLs using in-house routines and the program package handyG [43]. The remaining iterated integrals are expressed as one-fold integrals over GPLs of transcendental weight 2.
Topology F involves elliptic Feynman integrals and is more complicated. We have constructed a basis of the master integrals such that the 0 part of their differential equations are as simple as possible, where is the dimensional regulator. The solutions of these differential equations involve iterated integrals over elliptic integrals of the first and second kinds, which we perform numerically. These integrals turn out to be the most time-consuming part of the numeric evaluation of the amplitudes. It is possible to speed up the computation using methods of, e.g., Refs. [44,45].

NUMERIC RESULTS AND DISCUSSIONS
In this section, we present some numeric results for the finite part of the two-loop amplitude, i.e., V fin defined in Eq. (6), computed using our method. While these results are not new, our aim is to demonstrate the precision and efficiency that can be achieved from the small-mass expansion. We will compare our results with those computed using sector decomposition in [10][11][12] (optionally supplemented by the Padé-improved high-energy expansion [23][24][25][26]). For that purpose we choose the same input parameters as in [26]: the Higgs mass m H = 125 GeV, the top quark mass m t = 173 GeV, and the renormalization scale µ = √ŝ /2. Our approach is a mixture of analytic and numeric methods. Its main advantage lies in the fact that it is much more efficient than the purely numeric integration with sector decomposition, yet provides even better precision in a major portion of the whole phase space. To see that, we show V fin as a function of √ŝ and p T in Fig. 2. In order to generate the surfaces, we plot the results of V fin at 6320 phase-space points corresponding to the grid of [26], and then simply join them without any interpolation or fitting procedure.
The upper plot in Fig. 2 uses the values of V fin taken from the grid file of [26], which were computed using sector decomposition. It can be seen that the surface has a lot of spikes and is far from smooth. This is due to the uncertainties of numeric integrations, which are most severe when √ŝ is above the 2m t threshold. The integration uncertainties can in principle be reduced with more sampling points in the (quasi-)Monte Carlo methods. However, this requires much more computational resources which prevents its viability. To cure this spiky behavior, Ref. [26] has created an interpolation code which applies a Clough-Tocher interpolator resulting in a smooth distribution. This, however, does not reduce the intrinsic uncertainties of the results.
The lower plot in Fig. 2, on the other hand, uses our results of V fin computed with the small-mass expansion up to O(m 4 H ). One can see that the double distribution is rather smooth in the whole phase-space region probed by the 6320 points. We emphasize again that the surface is generated without any interpolation or fitting procedure. This demonstrates the high precision and stability of our method, which is crucial for practical phenomenological applications. Due to the high efficiency of our method, it is straightforward to generate more phase-space points on a modern laptop computer. This helps to better describe the shape of the double distribution, which can then be easily integrated over without worrying too much about the interpolation.
The 3-D plots in Fig. 2 only give a qualitative impression of the two approaches. To make a more quantitative comparison, in Fig. 3, we present V fin as a function of √ŝ with fixed values of p T . The black curves are obtained from our results with the small-mass expansion up to O(m 4 H ), while the red ones are obtained using the interpolation code associated with [26] (with the grid file containing results from the Padé-improved high-energy expansion). The error bars of the interpolated results are calculated following the approach outlined in [26]. We note that when no error bar is attached to a given point, it does not necessarily mean that the result is accurate, but rather means that the result is not supported by nearby grid points and might not be reliable. The uncertainties of our results come mainly from the missing higher order terms at O(m 6 H ), which would not be visible on the plot. Overall, we find that our results are consistent with the interpolated ones wherever the interpolation is reliable (i.e., where there are enough grid points to support the interpolation and the associated numeric uncertainties are small). Also, our curves are smooth in all cases. These observations lend us confidence that the small-mass expansion provides a rather good approximation to the exact result across the whole phase space.
In Fig. 3, one can also see that there are discrepancies between the two results in certain phase-space regions, which we now discuss. In the small p T region (represented by the p T = 100 GeV curves), the high-energy expansion does not apply and the interpolated results are fully controlled by numeric integrations with sector decomposition. Since the numeric integration is timeconsuming, the number of grid points are not sufficient to support the whole range of √ŝ (here only those points with error bars are supported). Also, numeric integrations lead to large uncertainties at moderate and high √ŝ , which can be clearly seen in the plot. These two factors result in the wiggly behavior of the p T = 100 GeV curve for √ŝ > 700 GeV. On the other hand, our result remains smooth and is much more reliable in this region.
As p T increases, the Padé-improved high-energy expansion comes into play. This helps to make the curve smoother, albeit that the uncertainties in the region of moderate √ŝ are still large. The deviation in the tail (high √ŝ ) region is due to the lack of grid points there to support the interpolation. Note that the interpolation in the tail region can in principle be improved by adding more data from high-energy expansion. Finally, at high p T (represented by the p T = 800 GeV curves), the two results almost completely coincide with each other, although the interpolated curve still has some small wiggles.
Given the high precision of our results in the whole phase space, we now turn to discuss the efficiency of our approach compared to that based on sector decomposition. Ref. [10] quoted their performance of 4680 GPGPU hours (using NVIDIA Tesla K20X GPUs) for 665 phasespace points, giving an average of about 7 GPGPU hours per phase-space point. On the contrary, we only need about 10 seconds for one phase-space point using one CPU core on a desktop or even laptop computer. Note that a majority of the time in our approach is spent on the integration over elliptic functions, which may be further reduced using methods similar to those proposed in, e.g., [44,45]. The high efficiency of our method makes it valuable for phenomenological analyses, in particular, for studies involving different values of m t and m H .
In all previous numerics, we have truncated the smallmass expansion at O(m 4 H ). It is now time to present the convergence of the expansion to justify this choice. In Fig. 4

CONCLUSION
In this Letter, we present an efficient method to calculate the two-loop amplitudes for Higgs boson pair production via gluon fusion. Our method is based on a series expansion in terms of the Higgs boson mass m H . After expansion, the remaining integrals can be calculated with a combination of analytic and numeric methods.
The main advantage of our approach lies in two aspects: precision and efficiency. To demonstrate these, we compare our results with the most up-to-date calculation which combines the results from sector decomposition and those from Padé-improved high-energy expansion. While in principle the method of sector decomposition can lead to the exact answer for the relevant integrals, in practice it is highly challenging to reduce the uncertainty of the Monte-Carlo integration, especially in the moderate and high energy regions. On the other hand, the integration uncertainty in our approach is negligible, and the uncertainty due to the expansion is undercontrol and can be easily reduced by including higher power terms in m 2 H . This leads to precision numeric results across the whole phase-space, as can be seen by the smooth behavior of kinematic distributions and the excellent agreement between our predictions and those from other approaches. Our method also gives rather good numeric performance, such that the computation can be carried out by a normal desktop or laptop computer. The precision and efficiency of our method is highly valuable for practical phenomenological analyses.
Our calculation can be easily extended to the cases with dimension-6 effective operators describing new physics beyond the SM, as the required two-loop integrals are the same. Our method can also be applied to Higgs boson pair production in a concrete new physics model, where new heavy particles such as the top quark partners may enter the loop. Finally, our method can be implemented for other 2 → 2 processes with a heavy quark loop, such as gg → ZZ, gg → ZH and gg → Hj, which we leave for future investigations.