On the two-loop virtual QCD corrections to Higgs boson pair production in the Standard Model

We compute the next-to-leading order virtual QCD corrections to Higgs pair production via gluon fusion. We present analytic results for the two-loop contributions to the spin-0 and spin-2 form factors in the amplitude. The reducible contributions, given by the double-triangle diagrams, are evaluated exactly while the two-loop irreducible diagrams are evaluated by an asymptotic expansion in heavy top quark mass up to and including terms of $\mathcal{O}(1/m_t^8)$. Assuming that the finite top-quark mass effects are of similar size in the entire range of partonic energies we estimate that mass effects can reduce the hadronic cross section by at most $10\%$.


Introduction
After the discovery of the Higgs boson in Run 1 of the Large Hadron Collider (LHC) [1,2], one of the major targets of Run 2 is the experimental exploration of its properties. In Run 1, the measured Higgs boson production rate and the extracted values of the Higgs couplings to fermions and to gauge bosons have been found to be compatible with the predictions of the Standard Model (SM) within an experimental accuracy of (10 -20)% [3]. On the other hand, the self-couplings of the Higgs boson, which in the SM are determined in terms of the mass of the Higgs boson and the vacuum expectation value of the Higgs field and are thus fully predicted, have not been probed yet. They are accessible in multi-Higgs production processes [4,5] though a measurement of the quartic Higgs self-coupling lies beyond the reach of the LHC [6,7]. Instead, for the trilinear Higgs self-coupling various studies showed that it might be accessible at the LHC in Higgs pair production in bbγγ [8][9][10][11][12][13], bbττ [9,14], bbW + W − [15] and bbbb [16][17][18] final states.
Higgs pair production is not only interesting as a probe of the trilinear Higgs selfcoupling, but its rate can be significantly modified by new physics effects. For the dominant Higgs pair production mode, gluon fusion, this can, for instance, occur due to new loop contributions [19], in models with novel hhtt coupling [20][21][22] or if the Higgs boson pair is produced through the decay of a heavy new resonance. The latter two possibilities can lead to a strong increase of the cross section. First limits on such scenarios have been given in refs. [23][24][25][26][27].
A precise prediction of the gluon fusion Higgs-pair production channel is essential to constrain new physics or to determine the Higgs self-coupling. The gluon fusion process is mediated by heavy fermions via diagrams with box and triangle topologies and is hence loop-induced already at the leading order (LO). In the "triangle" contribution a single Higgs boson splits via an s-channel exchange into two Higgs bosons, thus it contains the trilinear Higgs self-coupling. The "box" contribution plays the role of an irreducible background, as it does not incorporate the trilinear Higgs self-coupling.
In the SM, the LO cross section is fully known since the late eighties [28]. However, similarly to what happens in single Higgs production, one expects the LO contribution to be subject to large radiative corrections. A computation of a 2 → 2 process at higher orders is extremely challenging. The next-to-leading order (NLO) "triangle" contribution can be borrowed from the production of a single Higgs boson [29][30][31][32], whereas a full computation of the NLO "box" form factors is at the moment not available and technically much more difficult. Higher order corrections to Higgs pair production are, however, available in the effective theory with infinite top mass, m t , or, equivalently, in the limit of vanishing external momentum, at NLO [33] and more recently also at next-to-next-to-leading order (NNLO) [34,35]. 1 Soft gluon resummation at next-to-next-to-leading-logarithmic (NNLL) accuracy has been performed in refs. [41,42]. Whereas the approximation of small external momenta was shown to work quite well for single Higgs production [29], it can be expected to be less effective for pair production, due to the larger energy scale that characterizes the latter process. The approximation can, however, be improved by factoring out the full LO cross section. The error due to the infinite top-mass limit for the part related to the real corrections has been estimated in refs. [43,44] to be roughly −10% by comparing the m t → ∞ limit result with the numerical calculation of the real corrections with full top-mass dependence. Instead, the uncertainty of the effective-theory result for the virtual corrections has been estimated in refs. [45] by the inclusion of higher orders in an expansion in small external momenta finding a positive shift with respect to the m t → ∞ result. This leads to an estimate of the uncertainties due to mass corrections at NLO, including also the the real contributions expanded in small external momenta, of ±10%, with a reduction to ±5% when the NNLO effective theory result is included [46] .
In this paper we reexamine the evaluation of the virtual NLO QCD corrections in Higgs pair production. We present an exact result for the reducible contribution given by the double-triangle diagrams, while the irreducible diagrams are evaluated via an asymptotic expansion in the top mass. Our work differs from the similar previous analyses in refs. [45,46] by the fact that we perform the asymptotic expansion up to and including terms O(1/m 8 t ) at the level of the amplitudes and not of the cross section, allowing us to derive simple analytic expressions for the spin-0 and spin-2 form factors in the amplitudes. The latter could be used in the future as a check of the result, in the relevant center-of-mass partonic energy region, when a complete calculation of the virtual corrections will be available. Furthermore, our expressions can be easily implemented in Monte-Carlo codes that compute the hadronic cross section in order to achieve a better description of the partonic center-of-mass energy region below the the 2 m t threshold.
In order to quantify the finite top-mass effects in the NLO corrections to the hadronic cross section we make two different comparisons: i) We compare the NLO cross sections computed using different orders in the top-mass expansion. ii) We compare the cross section including the O(1/m 8 t ) terms with the one computed factorizing the exact LO cross section while evaluating the NLO correction factor in the m t → ∞ limit.
The paper is organized as follows: in section 2 we give general formulae for the Higgs pair production cross section. In the next section we discuss different large-mass evaluations of the LO cross section comparing them with exact result. In section 4 we outline our method of calculation of the NLO corrections that are presented in the next section where we also discuss their numerical impact and the estimate of the error due to the mass effects in the virtual corrections. Finally, in section 6 we draw our conclusions. The paper is completed with an appendix where we present the analytic result for the expanded NLO form factors up to and including terms of O(1/m 8 t ).

Double Higgs Production via gluon fusion
In this section we summarize some general results on the Higgs boson pair production via the gluon fusion mechanism in proton-proton collisions, pp → HH. The hadronic cross section for the process p + p → H + H + X at center-of-mass energy √ s, can be written as: where M 2 HH is the invariant mass of the two Higgs system, τ = M 2 HH /s, µ F is the factorization scale, f a (x, µ 2 F ), the parton density of the colliding proton for the parton of type a, (a = g, q,q) andσ ab is the cross section for the partonic subprocess ab → H + H + X at the center-of-mass energyŝ = x 1 x 2 s. The partonic cross section can be written in terms of the LO cross section σ (0) as: where, up to NLO terms, with µ R denoting the renormalization scale. The LO contribution is given by the gluongluon (gg) channel only, i.e. G (0) The amplitude for g µ a (p 1 )g ν b (p 2 ) → H(p 3 )H(p 4 ) can be written as: where T F is the matrix normalization factor for the fundamental representation of SU (N c ) (T F = 1/2) and the form factors F 1 , F 2 are functions, besides of m 2 t , of the partonic Mandelstam variablesŝ In eq. (5) the orthogonal projectors A 1 and A 2 onto the spin-0 and spin-2 states, respectively, in n d = 4 − 2 dimension and normalized to 2 read A µν with p T the transverse momentum of the Higgs particle that can be expressed in terms of the Mandelstam variables as The spin-2 state receives contributions only from box topologies (see fig. 1) while in the spin-0 case both box and triangle diagrams contribute such that F 1 takes the form where F ∆ (F ) is the contribution of the triangle (box) diagrams.
The Born cross section is written as The one-loop form factors F 1 1 , F 1 2 are fully known analytically [28,47] and their values in the limit of vanishing external momentum can be obtained via a low energy theorem (LET) calculation [48][49][50] giving F 1 ,LET The NLO terms include, besides the gg channel, also the one-loop induced processes gq → qHH and qq → gHH. The gg-channel contribution, involving two-loop virtual corrections to gg → HH and one-loop real corrections from gg → HHg, can be written as where In eq. (12), C A = N c (N c being the number of colors), β 0 = (11 C A − 2 N f )/6 (N f being the number of active flavors) is the one-loop β-function of the strong coupling in the SM, R gg is the contribution of the real corrections, P gg is the LO Altarelli-Parisi splitting function and The first line of eq. (12) displays the two-loop virtual contribution regularized by the infrared singular part of the real-emission cross section. In eq. (13) the terms represents the contribution of the two-loop double-triangle diagrams with a t/u-channel gluon exchange ( fig. 2b) to the spin-0 (spin-2) part of the amplitude. In the limit of vanishing external momenta the double-triangle diagrams can be expressed in terms of The second line in eq. (12) contains the non-singular contribution from the real gluon emission in the gluon-fusion process. The function R gg is obtained from one-loop diagrams where quarks circulate in the loop, and in the limit of vanishing external momenta it becomes R gg → −11(1 − z) 3 /(6z). The contributions of the gq → qHH and qq → gHH channels are given by: where The functions R qq and R qg in (17) are obtained from one-loop quark diagrams, and in the limit of vanishing external momenta become

Large mass evaluation of the LO cross section
Even though the one-loop form factors F 1 1 , F 1 2 are fully known analytically [28,47], we will give here approximate results in order to inspect the validity range of the applied  ) and (20)) normalized to the real part of exact F 1 1 form factor.
approximations. This will later on allow us to apply the same approximations to the NLO cross section, where the full form factors are yet unknown. We discuss the large top-mass-expansion evaluation of the LO cross section. We start by reporting the expressions that we obtained via a Taylor expansion forŝ,t,û, m 2 The evaluation of the LO cross section using for F 1 and F 2 the values obtained via the LET calculation, i.e. the leading term in the large top-mass expansion in eqs. (19)(20)(21), gives a poor approximation of the exact result. Furthermore, the validity of this approximation is quite sensitive to the hadronic center-of-mass energy and to the choice of the renormalization and factorization scales [51]. This is at variant with the case of single Higgs production where the LET result gives a quite accurate estimate of the cross section. Indeed, the LET result is expected to be reliable in the region of partonic energies below the √ŝ < 2 m t threshold. In Higgs pair production, also the region above the 2 m t threshold contributes significantly to the hadronic cross section up to √ŝ ∼ 600 − 700 GeV. In this latter region the vanishing external momenta condition is obviously not satisfied and therefore the result obtained in this approximation is unreliable.
The inclusion of more terms in a large top-mass expansion of the form factors does not improve the evaluation of the LO cross section [51,52]. The reason is easily understood looking at the plots in Fig. 3. They are obtained evaluating the F 1 form factor with p T randomly generated but distributed as for the the integration of the full LO cross section. The spread in the points for equal √ŝ is induced by the difference in the value of p T for fixed √ŝ . The LET result 2 ( fig. 3a) approximates relatively well the exact result for F 1 1 in the region √ŝ 2 m t but it fails in describing the region √ŝ > 2 m t when √ŝ 450 GeV. The sum of the first five terms in the large top-mass expansion of F 1 1 ( fig. 3b) reproduces quite well the exact results when √ŝ 400 GeV while the region √ŝ > 400 GeV is described very badly, worse than in the LET case. Similar considerations apply to F 1 2 . We remark that the evaluation of F 1 and F 2 via a large mass expansion has a range of validity up to the 2 m t threshold. Describing the region above this threshold via the LET results means to replace the exact form factors by constant values. Instead using the sum of few terms in the large mass expansion means to replace F 1 and F 2 by a powerlike combination ofŝ/m 2 t that has a wrong behavior whenŝ grows. As a consequence, the partonic cross section in eq. (11) grows, for large values of the partonic center-of-mass energy, asŝ in the former case, while asŝ n+1 /m 2n t in the latter case with n the order of the expansion. Although in both cases the behavior of the partonic cross section in the region √ŝ > 2 m t is not described correctly, it is evident that in this region the cross section is much better (or less worse) approximated by its LET value than by including additional terms in the large mass expansion. As a further remark, we recall that the full form factors develop an imaginary part above √ŝ > 2 m t which cannot be described by an expansion in small external momenta. This imaginary part is however smaller than the real part up to √ŝ ≈ 450 GeV.
In Fig. 4 we present the partonic cross section as a function of √ŝ . The exact cross section (solid black line), σ (0) ex , is compared with the approximated ones (dashed colored lines), σ (0) app,n , obtained using for the form factors the expansions in eqs. (19)(20)(21) to the order n. The figure tells us that the validity of an estimate of the hadronic cross section from eq. (1) based on the use of σ √ŝ are going to contribute more to the hadronic cross section, so that the LET approximation is going to grow in size and therefore become either closer to the full cross section or overestimating it. For instance for √ s = 100 TeV the LET result overestimates the full cross section by a factor ∼ 2.2. Figure 4 indicates that an estimate of the LO hadronic cross section obtained employing  the large-mass expanded results for F 1 and F 2 in the entire range of partonic energies is not going to be realistic. An alternative estimate, based on the use of the maximal approximate information available and on simplicity, can be obtained by evaluating F 1 and F 2 via a large mass expansion only up to a cut √ŝ c in the partonic center-of-mass energy while above √ŝ c , where we do not trust any more the expansion, setting them to their LET values 3 . This can be considered an improvement with respect to an evaluation based only on the LET result because we are describing better the region √ŝ < 2 m t . In Table 1 we report the values of the LO hadronic cross section computed employing different orders in the expansion of F 1 and F 2 from eqs. (19 -21) in the region below √ŝ Note that the bottom quark loops contribute with less than 1%. One can see from the first column in the table that the use of the large mass expansion in the entire range of partonic energies gives rise to a non convergent result. The table also shows that the if √ŝ c is taken around 400 GeV the LO cross section obtained in this way is closer to the exact result than the one that is obtained using the LET results (last column of the table).

Outline of calculation
An exact analytic evaluation of the two-loop QCD corrections to the F 1 and F 2 form factors is presently not available. Exact expressions for F 2∆ 1 and F 2∆ 2 can be derived given the structure of the double-triangle diagrams ( fig. 2b) that allows to express the result in terms of products of one-loop Passarino-Veltman functions [57]. An exact analytic result for F 2l ∆ can be obtained by adapting the corresponding calculation in single-Higgs production [29][30][31][32]. Instead the exact analytic evaluations of F 2 and F 2 2 seem, at the moment, beyond our computational ability. However, it seems feasible to obtain an approximate evaluation of latter form factors using the method of asymptotic expansions [58,59]. Two different kind of expansions must be employed according to the region of partonic energy one is considering: for √ŝ 2 m t a large mass expansion in the top mass has to be performed while in the complementary region ( √ŝ 2 m t ) a large momentum expansion is required 5 . Here we provide a first step in the evaluation of the O(α s ) corrections to F 1 and F 2 via asymptotic expansions addressing the large mass case.
The large top-mass expansions of the two-loop diagrams contributing to F 2 and F 2 2 is performed using the strategy described in ref. [60] that we briefly recall here. The relevant diagrams are generated with the help of FeynArts [61], and contracted with the projector A µν 1 (A µν 2 ) to extract the F 1 (F 2 ) contribution. Then they are separated in two classes: i) those that can be evaluated via an ordinary Taylor expansion in powers ofŝ/m 2 t ,t/m 2 t andû/m 2 t ; ii) the diagrams that require an asymptotic expansion, i.e. those that when Taylor-expanded in the external momenta exhibit an infrared (IR) divergent behavior.
Class-i diagrams require the evaluation of the generic integral v(j 1 , . . . , j 9 , m 1 , m 2 , where any exponent j 1 − j 9 is either 0 or a positive integer and the propagator masses, m 1 − m 3 , are either m t or 0. The integral (23) can be reduced to vacuum integrals, i.e. v(0, . . . , 0, j 7 , j 8 , j 9 , m 1 , m 2 , m 3 ), using the tensor reduction formula presented in ref. [62]. The two-loop vacuum integrals obtained from the reduction can be evaluated using the results of ref. [63].
The two-loop diagrams that belong to class ii) are those containing either two triplegluon vertices ( fig. 2c) or one four-gluon vertex (fig. 2d). These diagrams become more and more IR divergent when the gluon propagators are Taylor-expanded with respect to the external momenta. The evaluation of the class ii) diagrams is obtained by supplementing the Taylor-expanded result by the exact computation of their IR divergent contribution. 6 The IR-divergent part of any diagram is constructed by a repeated application of the identity in eq. (3.1) of ref. [60] controlled by the power counting in the IR-divergent terms. The outcome of the procedure is that the IR-divergent part of any diagram is expressed in terms of products of one-loop integrals with numerators that contain terms of the form (k i ·q j ) m (k 1 ·k 2 ) n (i = 1, 2 , j = 1, 2, 3) where m, n are generic powers. Finally, the Passarino-Veltman reduction method is applied to eliminate the numerators and express the result in terms of the known one-loop scalar integrals [57].

Virtual corrections to gg → HH
In this section we give the analytical results for the double-triangle form factors F 2∆ 1 and F 2∆ 2 and the two-loop form factors F 2 1 and F 2 2 and discuss their numerical impact.

Analytic results for the two-loop form factors
We present, for the first time, the exact computation of the double-triangle diagrams, i.e. keeping the full dependence on the quark masses. The top contribution to the form factors can be expressed in terms of one-loop integrals so that, defining we find for F 2∆ 1 and F 2∆ 2 in eq. (13) The finite parts of the scalar one-loop integrals appearing in eq. (24) are given by with β x = 1 − 4 m 2 /x and β y defined in analogy. The bottom contribution can be obtained thorough the substitution m t → m b in eq. (24).
The two-loop form factors F 2 ∆ , F 2 and F 2 2 can be written as where F 2 i,C F is directly obtained from the two-loop virtual diagrams and depends upon the renormalized top-mass parameter employed that we choose to be the on-shell mass. The term F 2 i,C A represents the IR regularized results after subtraction of the IR poles, i.e.
where F virt i,C A is the contribution of the two-loop virtual diagrams and δF i,C A the counterterm required to make it finite that reads where F 1 i (ŝ, ) is the one-loop result including the O( , 2 ) terms. Employing the method described in sect. 4 we obtained the large top-mass expansion of the two-loop spin-0 and spin-2 form factors up to and including terms O(1/m 8 t ). The results are presented in appendix 1. As in the one-loop case the form factors are expressed in terms ofŝ, p 2 T , m 2 H , and m 2 t . The computation was performed first using orthogonal projectors in n d = 4 − 2 dimension (see eqs. (7,8)) and then using orthogonal projectors in n d = 4 dimension. We found that, after the addition of the counterterm pieces from eq. (31), the two results are identical. We checked that, once the IR counterterm is chosen as in eq. (31), F 2 ∆,C A reproduces the result for the triangle form factor that can be obtained directly adapting the known results on single Higgs production (cfr. eq. (A2) with ref. [32]).
Finally, we want to comment on the comparison of our results with the ones of refs. [45,46]. These references deal with the large top-mass evaluation of the NLO cross section while we concentrated only on the virtual corrections. We use a different method for the asymptotic expansion compared to refs. [45,46]. Instead of adding subgraphs and co-subgraphs we followed ref. [60], where we only add the IR divergent parts, evaluated fully, to the diagrams that exhibit IR divergencies. We work at the level of amplitudes while in ref. [45] the total cross section is computed by deriving the imaginary part of the gg → gg amplitude and connecting it via the optical theorem to the total cross section. 7 In refs. [45,46] the phase space integrals are computed analytically which requires an expansion in δ = 1 − 4 m 2 H /ŝ including the s-channel Higgs propagator in the triangle contributions. The result is then expressed as an expansion both in ρ = m 2 H /m 2 t and δ. Our approach is instead to compute the phase space integrals numerically via Monte Carlo methods. Performing a Monte Carlo integration of the phase space integrals will allow us in future to include expansions in other regimes or exact results, once available, in a straightforward way. We remark that the integration overt in eq. (13) of the expanded form factors can be easily done analytically as the latter are given as power series int. A precise numerical comparison with refs. [45,46] cannot be performed, since we did not compute the real radiation part of the gg amplitude.

Numerical results
We discuss now the numerical impact of the corrections we computed. The numerical results are obtained with a private version of the code HPAIR [53] where we implemented our results. The inputs in the code are the same as in table 1, but using the NLO value for the strong coupling, α N LO s (m Z ) = 0.12018. We start analyzing the NLO contribution due to the double-triangle diagrams, i.e. σ (2∆) = σ (0) C NLO with F 2 1 = F 2 2 = 0 (see eq. (13)). In order to quantify the impact of the inclusion of the finite mass effects we plot, in fig. 5, σ (2∆) computed in two ways: i) exactly (solid line), i.e. with all the form factors evaluated in full mass dependence, namely we use for F 2∆ 1 (F 2∆ 2 ) eq. (25) (eq. (26)). ii) Computing F 2∆ 1 and F 2∆ 2 using their LET approximation as given in eq. (16), while employing the exact expressions for F 1 1 and F 1 2 (dashed line). The figure shows that the inclusion of the finite top mass effects changes the double-triangle contribution to the partonic cross section by ∼ 20 − 30%. We remark that the double-triangle contribution to the hadronic cross section is actually always very small. Indeed, it amounts to ∼ −0.18 fb while the NLO cross section is around 40 fb.
We turn now to discuss the contribution in C NLO due to F 2 1 and F 2 2 in eq. (13). We expect our results for F 2 1 and F 2 2 to be quite accurate for √ŝ 400 GeV, in analogy with the LO case as shown in fig. 3. This allows us to evaluate the contribution induced by the mass effects in the virtual part of the NLO corrections by computing σ (0) C NLO at various order in the large mass expansion.
In table 2 we report the contribution of     [28,47], while F 2 ,n i the expression for the two-loop form factor we derived (eqs. (A1-A6)) to the relevant order n. Comparing the LET row with the 1/m 8 t one, we find that the mass effects induce a relative variation with respect to the m t → ∞ result up to ∼ 40%.
Based on the experience gained in single Higgs production one expects that the factorization of the exact LO cross section can improve the m t → ∞ determination of the hadronic cross section. Applying the same procedure to a large mass expansion determination amounts to evaluate σ (0) C NLO employing for σ (0) the exact LO cross section while evaluating C NLO at the same order of approximation both in the numerator and in the denominator. The contribution to the hadronic cross section of σ 0 C NLO computed factorizing the exact LO cross section is presented in table 3. Looking at tables 2 and 3 we notice that the factorization of the exact LO cross section in the m t → ∞ result has the tendency to overestimate the NLO cross section as approximated by the 1/m 8 t rows in the tables. Both tables show in the first two columns a good convergence with respect to the order of the expansion, with the exception of the 1/m 2 t term. As expected, for M c HH > 2 m t the convergence starts to downgrade.
Our analysis cannot say anything about the region of partonic energies √ŝ 400 GeV where our results cannot be trusted. Concerning the hadronic cross section we can only make a guess assuming that in both regions, √ŝ 400 GeV and √ŝ 400 GeV, the variation induced by mass effects in σ 0 C NLO will be of a similar size and behaviour so that compensations between the two energy regions are not going to happen. Comparing the first row in table 3 with the last one in table 2 we find a relative variation ∼ 20%. We notice that the contribution of σ 0 C NLO to the NLO cross section is about 10% of the total, that the contributions we did not discuss, i.e. R gg , R qq and R qg , when evaluated in the limit of vanishing external momenta contribute to the total NLO cross section by ∼ 2% and that, according to the analysis in refs. [43,44], the finite mass effects reduce the size of the real contributions with respect their LET estimate. Considering a maximal case we expect that mass effects are going to reduce the m t → ∞ value of the NLO cross section by less than 10%. This size of variation is indeed found if one compares the NLO cross section evaluated in the m t → ∞ limit with the LO term factorized, σ N LO LET = 40.00 fb, with the NLO cross section computed as in table 1 with √ŝ c = 400 GeV, the cut in the partonic energy that at LO gives a result close to the exact LO value. The latter cross section, that is computed evaluating F 2 i in the region below √ŝ c using the 1/m 8 t order in the expansion while above √ŝ c employing the LET values, amounts to σ N LÔ sc = 37.86 fb. Finally we comment on larger hadronic center-of-mass energies. At LO, the LET result, e.g. at √ s = 100 TeV, approximates the true one worse than at √ s = 14 TeV. Even though a large center-of-mass energy gives a stronger weight to the region where the approximation of large top mass is not valid, we can expect that our conclusion on the uncertainty on the hadronic cross section due to mass effects is not going to change significantly, since the parts of the NLO cross section that are actually mass dependent are small.

Conclusions
In this paper we computed the virtual NLO QCD corrections in Higgs pair production. The double-triangle contribution was computed exactly while the spin-0 and spin-2 twoloop form factors in the amplitude were computed via an asymptotic expansion in the top mass up to and including terms O(1/m 8 t ). Analytic results are presented for both contributions. Before this work F 2 1,2 and F 2∆ 1,2 were known only in the m t → ∞ limit [33]. Our results allow a more precise evaluation of the NLO cross section for partonic energies up to √ŝ 400 GeV. This energy region is not the one contributing most to the hadronic cross section, however, its investigation enabled us to quantify the difference between the NLO result obtained in the m t → ∞ limit and the true one, where the top mass is kept finite. Although we did not discuss the large mass evaluation of the real contributions R gg , R qq and R qg , their size, as estimated from their LET values, is quite small so that even a 100% error on these terms will not make a large difference in the hadronic cross section. Under the assumption that in both energy regions, √ŝ 400 GeV and √ŝ 400 GeV, the finite top-mass effects are of similar size and behavior, we estimated that the true NLO result is going to be smaller than the one obtained in the LET limit by less than 10%. We remark that while our results for √ŝ 400 GeV are solid, our estimate of the hadronic cross section, as any other based on results obtained via a large top-mass expansion, should be understood just as a guess.
Our analysis differs in several points from previous works in the literature [45,46]. The main differences are: i) we performed the asymptotic expansion at the level of form factors and not of the cross section as in refs. [45,46]. ii) We did not discuss the real contributions as instead was done in those works. Point i) allowed us to compute the virtual NLO contribution as in eq. (32) without making use of the factorization of the exact LO cross section neither at the partonic level for the total cross section as in ref. [45] nor at the level of differential factorization, i.e. before the integration over the Higgs pair invariant mass, as in ref. [46]. The factorization of the LO cross section is known to work fine in single Higgs production where the exact NLO result is known [29][30][31][32], however there is no proof that the same happens also in double Higgs production. From the comparison of table 2 and table 3 in section 5 it seems that the differential factorization, that is expected to lead to a better result than the other possibility since it gives rise to a better-behaved integrand [46], when the LET result is employed tends to overestimate the result. Although a detail comparison of our results with those of refs. [45,46] is not possible, we notice that our results in table 2 and 3 exhibit the same behavior with respect to the order of the expansion of the soft-virtual cross section of ref. [46].
Finally, we would like to point out that our work should be seen as one of the first steps towards a complete calculation of the two-loop virtual corrections in Higgs pair production. A complete calculation of the NLO corrections, requires to address, besides the real contributions that were studied in refs. [43,44], the computation of the virtual corrections in the energy region √ŝ 400 GeV. These corrections are very difficult to compute but can be attacked either via a large momentum expansion calculation or by numerical methods.

Two-loop form factors in the large top-mass expansion
In this appendix we provide the two-loop form factors appearing in eq. (29) expanded in the small external momenta up to O(1/m 8 t ). The triangle form factors read (A2) The spin-0 box form factors are given by