An Analytical Method for the NLO QCD Corrections to Double-Higgs Production

We propose a new method to calculate analytically higher-order perturbative corrections and we apply it to the calculation of the two-loop virtual corrections to Higgs pair production through gluon fusion. The method is based on the expansion of the amplitudes in terms of a small Higgs transverse momentum. This approach gives a very good approximation (better than per-mille) of the partonic cross section in the center of mass energy region $\sqrt{\hat{s}} \lesssim 750$ GeV, where $\sim95\%$ of the total hadronic cross section is concentrated. The presented method is general and can be applied in a straightforward way to the computation of virtual higher-order corrections to other $2\to2$ processes, representing an improvement with respect to calculations based on heavy mass expansions.


INTRODUCTION
The experimental exploration of the properties of the Higgs boson is one of the major targets of the Large Hadron Collider (LHC). However the self-couplings of the Higgs boson, which in the Standard Model (SM) are fully determined in terms of the mass of the Higgs boson and the Fermi constant, have not been probed yet. While the quartic Higgs self-coupling is not directly accessible at the LHC [1,2], the trilinear self-coupling might be measurable from Higgs pair production processes [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].
Therefore, a precise prediction of the gluon fusion channel is essential to determine the Higgs trilinear selfcoupling and constrain new physics. At leading order (LO) the gluon fusion process is known since the 80s [29]. At next-to-leading order (NLO) this process is fully known only numerically [30,31], while analytical results are available in the heavy top mass (m t ) limit [32][33][34] and partially in the light m t limit [35]. In [36] a method was proposed for obtaining an analytical result combing large top mass expansion and a threshold expansion by means of Padé approximants.
We propose a new approach for the analytical calcu-lation of the virtual NLO corrections to the Higgs pair production through gluon fusion. The method is based on the expansion of the amplitudes around a small Higgs transverse momentum p T and Higgs mass m h . After properly expanding, the resulting amplitudes are functions of only m t and √ŝ and can be calculated analytically without resorting to further expansions. With this method we are able to correctly describe the Higgs pair production in the region √ŝ 750 GeV, nicely complementing the present literature. It must be also noted that, due to the shape of the gluon pdfs, this region represents 95% of the total hadronic cross section.
Our approach has the virtue of covering larger regions of the phase space with respect to approaches based on heavy mass expansions or high energy expansions and can be easily implemented to other 2 → 2 processes.
In this letter we describe the basics of the method, and the main results of our calculation, while we will reserve a more detailed discussion of the computation to future works.

NOTATION AND DEFINITIONS
In this section we introduce the notation we will use in the rest of the paper and define a set of kinematical variables. The amplitude g µ a (p 1 )g ν b (p 2 ) → H(p 3 )H(p 4 ) can be written as where G µ is the Fermi constant, α s (µ R ) is the strong coupling defined at the renormalization scale µ R and T F = 1/2 is the normalization factor for the fundamental representation of SU (N c ). In eq. (1) A µν 1,2 are the orthogonal projectors onto the spin-0 and spin-2 states, respectively, while the corresponding form factors We defined A µν 1 and A µν 2 as with p T the transverse momentum of the Higgs particle, that can be expressed in terms of the Mandelstam variables as The Born cross section, then, is For our purpose, it is particularly convenient to introduce the prime Mandelstam variables: for which s + t + u = 0. In these variables the Higgs transverse momentum becomes Our ultimate goal is to make an expansion for small Since the final result is symmetrical in t ↔ u , the latter can be achieved expanding 2 for t ∼ 0, u ∼ −s . This is going to restrict F 1 and F 2 in eq. (3) to a forward kinematic, namely to be function of s/m 2 t only, reducing the computational difficulty from a three scales problem to a single scale one. To perform the expansion we need to express the momenta in terms of the parallel and transverse components w.r.t the beam 1 All momenta are assumed incoming. 2 Expanding only in t ∼ 0 would not be correct if the final result were not symmetrical in t ↔ u axis. For this purpose we define the combination of momenta It is easy to show that r 2 =t,r 2 =û, and that where r µ ⊥ =r µ ⊥ is perpendicular to p 1 and p 2 and, as expected, Finally, in this reparametrization, A µ 1,2 assume particularly simple forms: EXPANSION From eq. (6), assuming real valuedt andû, we obtain the condition that allows us to expand for p 2 T /s 1 and m 2 h /s 1. Although our program is clear, it is hindered by the fact that p T does not appear directly at the amplitude level. However it is possible to show that an expansion for r µ ∼ 0 µ is equivalent to an expansion in p 2 T ∼ 0. Using eqs. (12) and (13), and noticing that r ⊥ is purely space-like, we can exchange the expansion in p 2 T ∼ 0 with an expansion in r µ ∼ 0 µ or, equivalently, p µ 3 ∼ −p µ 1 . This observation is one of the main results of this letter, and it allows us to proceed. We can then rewrite the form factors in eq. (1) as Although eq.(15) is always valid, the expansion proposed in (16) requires an hierarchy between r 2 and m 2 t . We are going to estimate the range of validity of the small p T expansion by comparing at the LO the result obtained via eq. 16 with the exact LO result (see next section).
We conclude this section with an important remark on how correctly truncate the series in eq. (16). Since the final result should be symmetrical for p 1 ↔ p 2 and ∂F 1,2 /∂p µ 3 is a rank 1 tensor, the second term in eq. (16) r.h.s can be rewritten as r µ F 1,2 (p µ 1 + p µ 2 ) = F 1,2 (t +û), with F 1,2 a function of s and m 2 t . For similar arguments, the third term should be instead proportional to r µ r ν (g µν + . . . ) = p 2 T + . . . . It is clear, then, that to expand to the first order in p 2 T one has to expand to the second order in p µ 3 , or, more in general, an order n expansion in p 2 T needs the order 2n expansion in p µ 3 .

COMPUTATION AND RESULTS
We generated the relevant amplitudes for the virtual NLO corrections to gg → HH with FeynArts [37]. The amplitudes were contracted with the two orthogonal projectors eq. (4) and (5), using FeynCalc [38], and reduced to a combination of scalar integrals. The integrals were then Taylor expanded, as described in the previous section. Subsequently, the resulting integrals were reduced in terms of a basis of Master Integrals (MI) using FIRE [39] and LiteRed [40]. All of the MIs, of which nearly the totality can be expressed in terms of multiple polylogarithms, were already known in the literature [41][42][43][44][45][46][47]. However, we evaluated them again directly in the phase space region of interest. We cross-checked our results using SecDec [48]. The details of the calculation presented here will be the topic of a second paper on this argument [49], while in this letter we will focus on the final result.
In order to show that our method correctly describes the partonic cross section for √ŝ < 750 GeV, we will start applying it to the LO.
In fig. 1 we report our calculation for the partonic cross section using eq. (16). As discussed in the introduction, while the heavy m t expansion describes well only the range √ŝ < 2m t , with our method we are able to correctly describe a wider range. It is also interesting to note that an expansion up to order O(p 4 T ) is already sufficient to describe the complete result with enough precision.
The range of validity of the small p T expansion can be estimated comparing the partonic cross section calculated with our method with the one from the full LO calculation.
In table I, we show where σ full is the cross section calculated without expansion, and σ exp is the one calculated in this letter. The FIG. 1: Partonic cross section of gg → HH as a function of the partonic center of mass energy. The black continuous line is the full result [29]. The dotted lines represent two orders of approximation in the heavy mt limit. The dashed lines are the result of the small p 2 T approximation presented in this letter.  ∼ 750 GeV. Moreover, in the region of interest, the approximation rapidly improves as one considers higher order in the expansion in p 2 T and m 2 h . The range of validity of our formulas is complementary to the one present in the literature, and represents 95% of the total hadronic cross section.
This behaviour is confirmed (and even improved) in the comparison with the full numerical result at NLO. It is well known that the NLO virtual corrections are IR divergent and these divergences cancel against the ones that come from real corrections [32][33][34]. Following [34], we cancel the IR divergences by adding the counterterm 1/(2 2 )F LO 1,2 ( 2 )(ŝ) − , where F LO 1,2 ( ) are the LO form factors with the inclusion of the O( , 2 ) terms. In fig. 2 we compare our result to the numerical results from [30], at the partonic level, using the grid and the interpolation function for the finite part of the virtual corrections V f in provided in [50]. As can be inferred from the figure, our expansion perfectly agrees with the full result when the first correction in p T and m h is included. It can clearly be seen that our lines smooth out the error on the full  result stemming from the interpolation.

CONCLUSION
In this letter we have proposed a novel approach for the analytical computation of the NLO virtual corrections to Higgs pair production through gluon fusion. This method, based on a expansion for small p 2 T , allows us to describe accurately the regionŝ 750 GeV that until now has been explored only numerically. In particular we showed that a few terms in the expansion already reproduce the full LO within 10 −3 , in the region of interest. At NLO we find excellent agreement already at O(p 2 T + m 2 h ) comparing to the full result of [30]. We remark that this method is general and can be useful for the analytic computation of radiative corrections to other fundamental processes for the physics programme of the LHC.