1 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 \(b\bar{b} \gamma \gamma \) [8–13], \(b\bar{b} \tau \bar{\tau }\) [9, 14], \(b\bar{b}W^+ W^-\) [15] and \(b\bar{b}b\bar{b}\) [16–18] final states and in Higgs-pair production in association with two jets [19, 20] and a \(t\bar{t}\) pair [21].

Higgs-pair production is not only interesting as a probe of the trilinear Higgs self-coupling, 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 [22], in models with novel \(hh\bar{t}t\) coupling [23–25] 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. [26–30].

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 1980s [31]. 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 \rightarrow 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  [32–35], 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 [36] and more recently also at next-to-next-to-leading order (NNLO) [37, 38].Footnote 1 Soft gluon resummation at next-to-next-to-leading-logarithmic (NNLL) accuracy has been performed in Refs. [44, 45]. Whereas the approximation of small external momenta was shown to work quite well for single Higgs production [32], 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. [46, 47] to be roughly \(-10~\%\) by comparing the \(m_t\rightarrow \infty \) 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. [48] by the inclusion of higher orders in an expansion in small external momenta finding a positive shift with respect to the \(m_t\rightarrow \infty \) result. This leads to an estimate of the uncertainties due to mass corrections at NLO, including also the real contributions expanded in small external momenta, of \(\pm 10~\%\), with a reduction to \(\pm 5~\%\) when the NNLO effective-theory result is included [49].

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 similar previous analyses in Refs. [48, 49] by the fact that we perform the asymptotic expansion up to and including terms \(\mathcal{O} (1/m_t^8)\) 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 easily be 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 \(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 \(\mathcal{O} (1/m_t^8)\) terms with the one computed factorizing the exact LO cross section while evaluating the NLO correction factor in the \(m_t\rightarrow \infty \) limit.

The paper is organized as follows: in Sect. 2 we give general formulas 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 Sect. 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 Sect. 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 \(\mathcal{O} (1/m_t^8)\).

2 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 \rightarrow H H\). The hadronic cross section for the process \(p + p \rightarrow H+ H+X\) at center-of-mass energy \(\sqrt{s}\) can be written as:

$$\begin{aligned} M_{HH}^2 \frac{\mathrm{d}\,\sigma }{\mathrm{d}\, M_{HH}^2}= & {} \sum _{a,b}\int _0^1 \mathrm{d}x_1 \mathrm{d}x_2 \,\,f_{a}(x_1,\mu _F^2)\, f_{b}(x_2,\mu _F^2)\nonumber \\&\times \int _0^1 \mathrm{d}z~ \delta \left( z-\frac{\tau }{x_1 x_2} \right) M_{HH}^2 \frac{\mathrm{d}\,\hat{\sigma }_{ab}}{\mathrm{d}\, M_{HH}^2}, \end{aligned}$$
(1)

where \(M_{HH}^2\) is the invariant mass of the two Higgs system, \(\tau = M_{HH}^2/s\), \(\mu _F\) is the factorization scale, \(f_{a}(x,\mu _F^2)\), the parton density of the colliding proton for the parton of type \(a, \,(a = g,q,\bar{q})\) and \(\hat{\sigma }_{ab}\) is the cross section for the partonic subprocess \( ab \rightarrow H+ H +X\) at the center-of-mass energy \(\hat{s}=x_1 x_2 s\). The partonic cross section can be written in terms of the LO cross section \(\sigma ^{(0)}\) as:

$$\begin{aligned} M_{HH}^2 \frac{\mathrm{d}\,\hat{\sigma }_{ab}}{\mathrm{d}\, M_{HH}^2}= \sigma ^{(0)}(z \hat{s})\,z \, G_{ab}(z), \end{aligned}$$
(2)

where, up to NLO terms,

$$\begin{aligned} G_{a b}(z) = G_{a b}^{(0)}(z) + \frac{\alpha _s (\mu _R)}{\pi } \, G_{a b}^{(1)}(z) \, \end{aligned}$$
(3)

with \(\mu _R\) denoting the renormalization scale. The LO contribution is given by the gluon–gluon (gg) channel only, i.e.

$$\begin{aligned} G_{a b}^{(0)}(z) = \delta (1-z) \,\delta _{ag}\, \delta _{bg} \, . \end{aligned}$$
(4)

The amplitude for \(g_a^\mu (p_1) g_b^\nu (p_2) \rightarrow H(p_3) H(p_4)\) can be written

$$\begin{aligned} A^{\mu \nu } = \frac{G_\mu }{\sqrt{2}} \frac{\alpha _s (\mu _R)}{2 \pi } \delta _{ab}\, T_F\, \hat{s}\left[ A_1^{\mu \nu } \,F_1 + A_2^{\mu \nu }\, F_2 \right] \end{aligned}$$
(5)

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_t^2\), of the partonic Mandelstam variables

$$\begin{aligned} \hat{s} = (p_1 + p_2)^2,\quad \hat{t} = (p_1 - p_3)^2,\quad \hat{u} = (p_2 - p_3)^2. \end{aligned}$$
(6)

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\,\epsilon \) dimension and normalized to 2 read

$$\begin{aligned} A_1^{\mu \nu }= & {} \sqrt{\frac{2}{n_d-2}} \left[ g^{\mu \nu }-\frac{p_1^{\nu }\,p_2^{\mu }}{\left( p_1\cdot p_2\right) } \right] \end{aligned}$$
(7)
$$\begin{aligned} A_2^{\mu \nu }= & {} \sqrt{\frac{n_d-2}{2(n_d-3)}} \left\{ \frac{n_d-4}{n_d-2} \left[ g^{\mu \nu }-\frac{p_1^{\nu }\,p_2^{\mu }}{\left( p_1\cdot p_2\right) } \right] +g^{\mu \nu } \frac{ p_3^2\, p_1^{\nu }\, p^{\mu }_2 - 2\left( p_3\cdot p_2\right) p_1^{\nu }\,p_3^{\mu } - 2\left( p_3\cdot p_1\right) p_3^{\nu }\,p_2^{\mu } + 2 \left( p_1\cdot p_2\right) p_3^{\mu }\,p_3^{\nu }}{p_{\scriptscriptstyle T}^2\left( p_1\cdot p_2\right) } \right\} \nonumber \\ \end{aligned}$$
(8)

with \(p_{\scriptscriptstyle T}\) the transverse momentum of the Higgs particle that can be expressed in terms of the Mandelstam variables as

$$\begin{aligned} p_{\scriptscriptstyle T}^2 = \frac{\hat{t}\hat{u} - m_{\scriptscriptstyle H}^4}{\hat{s}}~. \end{aligned}$$
(9)
Fig. 1
figure 1

Generic Feynman diagrams for box and triangle topologies for Higgs-pair production

Full size image

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

$$\begin{aligned} F_1 = F^{ }_\Delta \frac{3 m_{\scriptscriptstyle H}^2}{\hat{s} -m_{\scriptscriptstyle H}^2 } + F_\Box \end{aligned}$$
(10)

where \(F^{}_\Delta (F^{}_\Box )\) is the contribution of the triangle (box) diagrams.

The Born cross section is written as

$$\begin{aligned} \sigma ^{(0)}(\hat{s}) = \frac{G^2_\mu \alpha _s^2 (\mu _R)}{512\, (2 \, \pi )^3}\int _{\hat{t}_{-}}^{\hat{t}_{+}} \mathrm{d} \hat{t} \left\{ \left| T_F \, F^{1\ell }_1 (\hat{s}) \right| ^2 + \left| T_F\, F^{1\ell }_2 (\hat{s})\right| ^2 \right\} \nonumber \\ \end{aligned}$$
(11)

with \(\hat{t}_{\pm } = -\hat{s}/2 (1-2\,m_{\scriptscriptstyle H}^2/\hat{s} \mp \sqrt{1-4\,m_{\scriptscriptstyle H}^2/\hat{s}})\). The one-loop form factors \( F_1^{1\ell },\, F_2^{1\ell }\) are fully known analytically [31, 50] and their values in the limit of vanishing external momentum can be obtained via a low energy theorem (LET) calculation [51–53] giving \(F^{1\ell ,\mathrm{LET}}_\Delta =- F^{1\ell ,\mathrm{LET}}_\Box = 4/3, \, F^{1\ell ,\mathrm{LET}}_2 =0\), which correspond to the effective theory \(m_t\rightarrow \infty \) result.

The NLO terms include, besides the gg channel, also the one-loop induced processes \(gq \rightarrow q H H\) and \(q \bar{q} \rightarrow g H H\). The gg-channel contribution, involving two-loop virtual corrections to \(g g \rightarrow H H\) and one-loop real corrections from \( gg \rightarrow H H g\), can be written as

$$\begin{aligned}&G_{g g}^{(1)}(z) = \delta (1-z) \left[ C_A \, \frac{~\pi ^2}{3} \,+ \beta _0 \, \ln \left( \frac{\mu _R^2}{\mu _F^2} \right) + \mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\right] \nonumber \\&\quad + P_{gg} (z)\,\ln \left( \frac{\hat{s}}{\mu _F^2}\right) + C_A\, \frac{4}{z} \,(1-z+z^2)^2 \,\mathcal{D}_1(z) + C_A\, \mathcal{R}_{gg}, \end{aligned}$$
(12)

where

$$\begin{aligned}&\!\!\mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\nonumber \\&= \,\frac{\int _{\hat{t}_-}^{\hat{t}_{+}} \mathrm{d}\hat{t}\,\left[ \left( T_F\, F_1^{1\ell }\right) ^* \,T_F \left( \,F_1^{2\ell } {+} F_1^{2\Delta } \right) ~{+}~ \left( T_F\,F_2^{1\ell }\right) ^* \,T_F\left( F_2^{2\ell } {+} F_2^{2\Delta } \right) ~ \right] }{\int _{\hat{t}_-}^{\hat{t}_{+}} \mathrm{d}\hat{t}\, \left( \left| T_F \,F_1^{1\ell }\right| ^2 {+} \left| T_F\,F_2^{1\ell }\right| ^2\right) }\nonumber \\&\quad +\,\mathrm{h.c.}. \end{aligned}$$
(13)

In Eq. (12), \(C_A =N_c\) (\(N_c\) being the number of colors), \(\beta _0 = (11\, C_A - 2\, N_f)/6 \) (\(N_f\) being the number of active flavors) is the one-loop \(\beta \)-function of the strong coupling in the SM, \(\mathcal{R}_{gg}\) is the contribution of the real corrections, \(P_{gg}\) is the LO Altarelli–Parisi splitting function

$$\begin{aligned} P_{gg} (z) =2\, C_A\,\left[ \mathcal{D}_0(z) +\frac{1}{z} -2 + z(1-z) \right] , \end{aligned}$$
(14)

and

$$\begin{aligned} \mathcal{D}_i (z) = \left[ \frac{\ln ^i (1-z)}{1-z} \right] _+ . \end{aligned}$$
(15)

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 \(F_{1}^{2\ell }\) and \(F_{2}^{2\ell }\) contain the contribution of irreducible two-loop diagrams (see Fig. 2a, c, d) and in the limit of vanishing external momenta they read \(F^{2\ell }_\Delta = -F^{2\ell }_\Box = - C_F + 5/3 \,C_A, \, F^{2\ell }_2 = 0\) [36] with \(C_F = (N_c^2-1)/(2\,N_c)\). The term \(F_{1}^{2\Delta }\, (F_{2}^{2\Delta })\) 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 \(F_{\Delta }^{1\ell ,\mathrm{LET}}\) as

$$\begin{aligned}&F_1^{2\Delta }\,\rightarrow \, \frac{1}{2} T_F\, \left( F^{1\ell , \mathrm{LET}}_\Delta \right) ^2\quad \text {and}\nonumber \\&F_2^{2\Delta }\,\rightarrow \, - \frac{1}{2} T_F\,\frac{p_{\scriptscriptstyle T}^2}{2 \hat{t}\hat{u}}(\hat{s}-2\,m_{\scriptscriptstyle H}^2) \, \left( F^{1\ell , \mathrm{LET}}_\Delta \right) ^2~. \end{aligned}$$
(16)
Fig. 2
figure 2

Sample of Feynman diagrams for the virtual two-loop corrections to Higgs-pair production via gluon fusion

Full size image

The second line in Eq. (12) contains the non-singular contribution from the real gluon emission in the gluon-fusion process. The function \(\mathcal{R}_{gg}\) is obtained from one-loop diagrams where quarks circulate in the loop, and in the limit of vanishing external momenta it becomes \(\mathcal{R}_{gg} \rightarrow -11 (1-z)^3/(6 z)\). The contributions of the \(gq \rightarrow q H H\) and \(q \bar{q} \rightarrow g H H\) channels are given by

$$\begin{aligned} G_{q \bar{q}}^{(1)}(z)= & {} \mathcal{R}_{q \bar{q}},\nonumber \\ G_{q g}^{(1)}(z)= & {} P_{gq}(z) \left[ \ln (1-z) + \frac{1}{2} \ln \left( \frac{\hat{s}}{\mu _F^2}\right) \right] + \mathcal{R}_{qg}, \end{aligned}$$
(17)

where

$$\begin{aligned} P_{gq} (z)= C_F \,\frac{1 + (1-z)^2}{z}. \end{aligned}$$
(18)

The functions \(\mathcal{R}_{q \bar{q}}\) and \(\mathcal{R}_{q g}\) in (17) are obtained from one-loop quark diagrams, and in the limit of vanishing external momenta become \(\mathcal{R}_{q \bar{q}} \rightarrow 32 \,(1-z)^3/(27 z)\),  \(\mathcal{R}_{q g} \rightarrow 2\,z/3 - (1-z)^2/z\).

3 Large-mass evaluation of the LO cross section

Even though the one-loop form factors \( F_1^{1\ell },\, F_2^{1\ell }\) are fully known analytically [31, 50], we will give here approximate results in order to inspect the validity range of the applied 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 \(\hat{s}, \hat{t}, \hat{u}, m_{\scriptscriptstyle H}^2 \ll m_t^2\) up to and including \(\mathcal{O}(1/m_t^8)\) terms

$$\begin{aligned} F^{1\ell }_\Delta (\hat{s})= & {} \frac{4}{3} + \frac{7}{90} \frac{\hat{s}}{m_t^2} + \frac{1}{126} \frac{\hat{s}^2}{m_t^4} + \frac{13}{12600} \frac{\hat{s}^3}{m_t^6} + \frac{8}{51975} \frac{\hat{s}^4}{m_t^8}, \end{aligned}$$
(19)
$$\begin{aligned} F^{1\ell }_\Box (\hat{s})= & {} - \frac{4}{3} - \frac{7}{15} \frac{m_{\scriptscriptstyle H}^2}{m_t^2} - \frac{45 \,m_{\scriptscriptstyle H}^4 - 14\, m_{\scriptscriptstyle H}^2 \hat{s} + 6 \hat{s}^2}{315 \,m_t^4} + \frac{13}{630} \frac{p_{\scriptscriptstyle T}^2\,\hat{s} }{m_t^4} \nonumber \\&-\frac{780 \,m_{\scriptscriptstyle H}^6 -620\, m_{\scriptscriptstyle H}^4 \,\hat{s} + 355\, m_{\scriptscriptstyle H}^2\, \hat{s}^2- 16 \, \hat{s}^3}{18900\, m_t^6} \nonumber \\&- \frac{p_{\scriptscriptstyle T}^2 ( 11 \,\hat{s}^2 -36 \,m_{\scriptscriptstyle H}^2\, \hat{s})}{1890 \, m_t^6} \nonumber \\&- \frac{ 2400\, m_{\scriptscriptstyle H}^8 - 3480\, m_{\scriptscriptstyle H}^6 \,\hat{s} + 2955\, m_{\scriptscriptstyle H}^4 \,\hat{s}^2 - 704 \, m_{\scriptscriptstyle H}^2 \,\hat{s}^3 + 120\, \hat{s}^4}{207900 \, m_t^8} \nonumber \\&+ \frac{p_{\scriptscriptstyle T}^2 \hat{s} (114\, m_{\scriptscriptstyle H}^4 - 85\, m_{\scriptscriptstyle H}^2\, \hat{s} + 16 \, \hat{s}^2 - 8\, p_{\scriptscriptstyle T}^2 \,\hat{s})}{10395 \,m_t^8}, \end{aligned}$$
(20)
$$\begin{aligned} F^{1\ell }_2 (\hat{s})= & {} \frac{p_{\scriptscriptstyle T}^2}{m_t^2} \left\{ -\frac{11}{45} - \frac{ 62\, m_{\scriptscriptstyle H}^2 - 5 \, \hat{s}}{630 \,m_t^2} - \frac{400 \, m_{\scriptscriptstyle H}^4 - 156\, m_{\scriptscriptstyle H}^2 \, \hat{s} + 49 \, \hat{s}^2}{ 12600 \,m_t^4}\right. \nonumber \\&\left. + \frac{103}{18900} \frac{ p_{\scriptscriptstyle T}^2 \,\hat{s} }{ \,m_t^4} - \frac{980 \, m_{\scriptscriptstyle H}^6 - 867\, m_{\scriptscriptstyle H}^4 \, \hat{s} + 469\, m_{\scriptscriptstyle H}^2 \, \hat{s}^2 - 34 \, \hat{s}^3 }{103950 \,m_t^6}\right. \nonumber \\&\left. + \frac{p_{\scriptscriptstyle T}^2\,\hat{s}( 24\, m_{\scriptscriptstyle H}^2 \, -7 \,\hat{s}) }{4950 \,m_t^6} \right\} . \end{aligned}$$
(21)
Fig. 3
figure 3

a LET result for \(F_1^{1\ell }\) normalized to the real part of the exact \(F_1^{1\ell }\) form factor. b The sum of first five terms of the large top-mass expansion of \(F_1^{1\ell }\) (Eqs. (19) and (20)) normalized to the real part of exact \(F_1^{1\ell }\) form factor

Full size image

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)–(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 [54]. 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 \(\sqrt{\hat{s}} < 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 \(\sqrt{\hat{s}} \sim \) 600–700 GeV. In the 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 [54, 55]. The reason is easily understood looking at the plots in Fig. 3. They are obtained evaluating the \(F_1\) form factor with \(p_{\scriptscriptstyle T}\) randomly generated but distributed as for the integration of the full LO cross section. The spread in the points for equal \(\sqrt{\hat{s}}\) is induced by the difference in the value of \(p_{\scriptscriptstyle T}\) for fixed \(\sqrt{\hat{s}}\). The LET resultFootnote 2 (Fig. 3a) approximates relatively well the exact result for \(F_1^{1\ell }\) in the region \(\sqrt{\hat{s}} \lesssim 2\, m_t\) but it fails in describing the region \(\sqrt{\hat{s}} > 2 \, m_t\) when \(\sqrt{\hat{s}} \gtrsim 450\) GeV. The sum of the first five terms in the large top-mass expansion of \(F_1^{1\ell }\) (Fig. 3b) reproduces quite well the exact results when \(\sqrt{\hat{s}} \lesssim 400\) GeV while the region \(\sqrt{\hat{s}} > 400\) GeV is described very badly, worse than in the LET case. Similar considerations apply to \(F_2^{1\ell }\).

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 \(\hat{s}/m_t^2\) that has a wrong behavior when \(\hat{s}\) grows. As a consequence, the partonic cross section in Eq. (11) grows, for large values of the partonic center-of-mass energy, as \(\hat{s}\) in the former case, while as \(\hat{s}^{n+1}/m_t^{2 n}\) 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 \(\sqrt{\hat{s}} > 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 \(\sqrt{\hat{s}}> 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 \(\sqrt{\hat{s}}\approx 450 \text { GeV}\).

Fig. 4
figure 4

Leading order partonic cross section as a function of the partonic center-of-mass energy. The solid line corresponds to the exact result, the dashed ones to the results obtained using different terms in the large top-mass expansion

Full size image

In Fig. 4 we present the partonic cross section as a function of \(\sqrt{\hat{s}}\). The exact cross section (solid black line), \(\sigma ^{(0)}_{\mathrm{ex}}\), is compared with the approximated ones (dashed colored lines), \(\sigma ^{(0)}_{\mathrm{app},n}\), obtained using for the form factors the expansions in Eqs. (19)–(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 \(\sigma ^{(0)}_{\mathrm{app},n}\) depends on the relative weights in the hadronic integral of the regions where \(\sigma ^{(0)}_{\mathrm{app},n} < \sigma ^{(0)}_{ex}\) vs. \(\sigma ^{(0)}_{\mathrm{app},n} > \sigma ^{(0)}_{\mathrm{ex}}\) and how these two regions can compensate each other. With the increase in the hadronic energy, regions with larger \(\sqrt{\hat{s}}\) 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 \(\sqrt{s}=100 \text { TeV}\) the LET result overestimates the full cross section by a factor \(\sim 2.2\).

Table 1 Values in fb of the LO cross section computed using the large-mass expansion results of \(F_1\) and \(F_2\) of Eqs. (19)–(21) for partonic energies up to \(\sqrt{\hat{s}_c} \) (in GeV) while for partonic energies greater than \( \sqrt{\hat{s}_c} \) approximating \(F_1\) and \(F_2\) with their LET values
Full size table

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 \(\sqrt{\hat{s}_c}\) in the partonic center-of-mass energy while above \(\sqrt{\hat{s}_c}\), where we do not trust any more the expansion, setting them to their LET values.Footnote 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 \(\sqrt{\hat{s}} < 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 \(\sqrt{\hat{s}_c}\) while above it the LET values are used. The values for the cross section are obtained using a modified version of the code HPAIR [56], with \(\sqrt{s} = 14\) TeV, \(m_t=173.2\) GeV, \(m_{\scriptscriptstyle H}=125\) GeV and employing the parton distribution functions (pdf) MSTW08 [57–59]. The \(\alpha _s\) value is taken as the default in the pdf set, namely \(\alpha _s^{LO}(m_{\scriptscriptstyle Z}) = 0.13939\). The renormalization and factorization scales have been set to \(\mu _R=\mu _F=M_{HH}/2\) as suggested by the NNLL threshold expansion performed in Ref. [45]. The numbers in the table should be compared with the exact LO resultFootnote 4 that, including also the bottom contribution, reads

$$\begin{aligned} \sigma _{\mathrm{LO}}^{\mathrm{full}}=23.38 \text { fb}. \end{aligned}$$
(22)

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 \(\sqrt{\hat{s}_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).

4 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_{1}^{2 \Delta }\) and \(F_{2}^{2 \Delta } \) can be derived given the structure of the double-triangle diagrams (Fig. 2b) that allows one to express the result in terms of products of one-loop Passarino–Veltman functions [60]. An exact analytic result for \(F^{2l}_\Delta \) can be obtained by adapting the corresponding calculation in single-Higgs production [32–35]. Instead the exact analytic evaluations of \(F^{2\ell }_\Box \) and \(F^{2\ell }_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 [61, 62]. Two different kind of expansions must be employed according to the region of partonic energy one is considering: for \(\sqrt{\hat{s}} \lesssim 2 \,m_t\) a large-mass expansion in the top mass has to be performed while in the complementary region (\(\sqrt{\hat{s}} \gtrsim 2 \,m_t\)) a large momentum expansion is required.Footnote 5 Here we provide a first step in the evaluation of the \(\mathcal {O}(\alpha _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\ell }_\Box \) and \(F^{2\ell }_2\) is performed using the strategy described in Ref. [63] that we briefly recall here. The relevant diagrams are generated with the help of FeynArts [64], and contracted with the projector \(A_1^{\mu \nu } \, (A_2^{\mu \nu })\) 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 \(\hat{s}/m_t^2, \, \hat{t}/m_t^2\) and \(\hat{u}/m_t^2\); (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

$$\begin{aligned}&v(j_1,\dots ,j_9,m_1,m_2,m_3) = \int \mathrm{d}^4k_1 \,\mathrm{d}^4 k_2\nonumber \\&\quad \times \,\frac{ (k_1. p_1)^{j_1} (k_1. p_2)^{j_2} (k_1. p_3)^{j_3} (k_2. p_1)^{j_4} (k_2. p_2)^{j_5} (k_2. p_3)^{j_6} }{ (k_1^2 - m_1^2)^{j_7} (k_2^2 - m_2^2)^{j_8}((k_1+k_2)^2 - m_3^2)^{j_{9}}}\nonumber \\ \end{aligned}$$
(23)

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,\ldots ,0,j_7,j_8,j_9,m_1,m_2,m_3)\), using the tensor reduction formula presented in Ref. [65]. The two-loop vacuum integrals obtained from the reduction can be evaluated using the results of Ref. [66].

The two-loop diagrams that belong to class (ii) are those containing either two triple-gluon 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.Footnote 6 The IR-divergent part of any diagram is constructed by a repeated application of the identity in Eq. (3.1) of Ref. [63] 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}\cdot q_{j})^m \, (k_1 \cdot 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 [60].

5 Virtual corrections to \(gg\rightarrow HH\)

In this section we give the analytical results for the double-triangle form factors \(F^{2\Delta }_1\) and \(F^{2\Delta }_2\) and the two-loop form factors \(F_{1}^{2\ell }\) and \(F_{2}^{2\ell }\) and discuss their numerical impact.

5.1 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

$$\begin{aligned} F^{2 \Delta }(x)= & {} \frac{8\, m_t^4}{(m_{\scriptscriptstyle H}^2-x)^2} \bigg [1+ \frac{x}{(m_{\scriptscriptstyle H}^2-x)}\nonumber \\&\times \left( B_0(m_{\scriptscriptstyle H}^2, m_t^2, m_t^2) - B_0(x, m_t^2, m_t^2)\right) \nonumber \\&+\frac{1}{2} \left( 4\, m_t^2-m_{\scriptscriptstyle H}^2+x \right) \, C_0(0, x, m_{\scriptscriptstyle H}^2, m_t^2, m_t^2, m_t^2) \bigg ]^2,\nonumber \\ \end{aligned}$$
(24)

we find for \(F_{1}^{2\Delta }\) and \(F_{2}^{2\Delta }\) in Eq. (13)

$$\begin{aligned} F_1^{2\Delta }= & {} F^{2\Delta }(\hat{t}) + F^{2\Delta }(\hat{u})~, \end{aligned}$$
(25)
$$\begin{aligned} F_2^{2\Delta }= & {} \frac{p_{\scriptscriptstyle T}^2 }{\hat{t}} F^{2\Delta }(\hat{t})+ \frac{p_{\scriptscriptstyle T}^2 }{\hat{u}} F^{2\Delta }(\hat{u})~. \end{aligned}$$
(26)

The finite parts of the scalar one-loop integrals appearing in Eq. (24) are given by

$$\begin{aligned}&B_0(x, m^2, m^2)= 2 +\beta _x \log \frac{\beta _x-1}{\beta _x+1}-\log \frac{m^2}{\mu _R^2}, \end{aligned}$$
(27)
$$\begin{aligned}&C_0(0,x,y,m^2, m^2, m^2)= \frac{1}{2(x-y)}\nonumber \\&\quad \times \left( \log ^2\frac{\beta _x+1}{\beta _x-1}- \log ^2\frac{\beta _y+1}{\beta _y-1}\right) , \end{aligned}$$
(28)

with \(\beta _x=\sqrt{1-4\, m^2/x}\) and \(\beta _y\) defined in analogy. The bottom contribution can be obtained thorough the substitution \(m_t\rightarrow m_b\) in Eq. (24).

The two-loop form factors \(F^{2\ell }_\Delta , \, F^{2\ell }_\Box \) and \(F^{2\ell }_2\) can be written as

$$\begin{aligned} F_i^{2\ell } (\hat{s}) = C_F\, F_{i,C_{F}}^{2\ell }(\hat{s}) + C_A\, F_{i,C_{A}}^{2\ell } (\hat{s})\quad (i= \Delta ,\Box ,2) \end{aligned}$$
(29)

where \( F_{i,C_{F}}^{2\ell }\) 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_{i,C_{A}}^{2\ell }\) represents the IR regularized results after subtraction of the IR poles, i.e.

$$\begin{aligned} F_{i,C_{A}}^{2\ell } (\hat{s}) = F_{i,C_{A}}^{\mathrm{virt}} (\hat{s}) + \delta F_{i,C_{A}} (\hat{s}) \end{aligned}$$
(30)

where \(F_{i,C_{A}}^{\mathrm{virt}}\) is the contribution of the two-loop virtual diagrams and \( \delta F_{i,C_{A}}\) the counterterm required to make it finite that reads

$$\begin{aligned} \delta F_{i,C_{A}} (\hat{s}) = \frac{1}{2 \epsilon ^2} F_{i}^{1\ell } (\hat{s},\epsilon ) (\hat{s})^{-\epsilon } \end{aligned}$$
(31)

where \(F_{i}^{1\ell } (\hat{s},\epsilon )\) is the one-loop result including the \(\mathcal{O}(\epsilon , \epsilon ^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 \(\mathcal{O}(1/m_t^8)\). The results are presented in Appendix 1. As in the one-loop case the form factors are expressed in terms of \(\hat{s}, \, p_{\scriptscriptstyle T}^2,\, m_{\scriptscriptstyle H}^2,\) and \(m_t^2\). The computation was performed first using orthogonal projectors in \(n_d=4 -2\,\epsilon \) 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_{\Delta ,C_{A}}^{2\ell }\) reproduces the result for the triangle form factor, which can be obtained directly adapting the known results on single Higgs production (cf. Eq. (A2) with Ref. [35]).

Finally, we want to comment on the comparison of our results with the ones of Refs. [48, 49]. 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. [48, 49]. Instead of adding subgraphs and co-subgraphs we followed Ref. [63], 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. [48] the total cross section is computed by deriving the imaginary part of the \(gg\rightarrow gg\) amplitude and connecting it via the optical theorem to the total cross section.Footnote 7 In Refs. [48, 49] the phase space integrals are computed analytically which requires an expansion in \(\delta = 1- 4\,m_{\scriptscriptstyle H}^2/\hat{s}\) including the s-channel Higgs propagator in the triangle contributions. The result is then expressed as an expansion both in \(\rho =m_{\scriptscriptstyle H}^2/m_t^2\) and \(\delta \). 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 the future to include expansions in other regimes or exact results, once available, in a straightforward way. We remark that the integration over \(\hat{t}\) in Eq. (13) of the expanded form factors can easily be done analytically as the latter are given as power series in \(\hat{t}\). A precise numerical comparison with Refs. [48, 49] cannot be performed, since we did not compute the real radiation part of the gg amplitude.

5.2 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 [56] 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, \(\alpha _s^{\mathrm{NLO}} (m_Z) = 0.12018\).

We start analyzing the NLO contribution due to the double-triangle diagrams, i.e. \(\sigma ^{(2 \Delta )} =\sigma ^{(0)} \mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\) with \( F_1^{2\ell }=F_2^{2\ell }=0\) (see Eq. (13)). In order to quantify the impact of the inclusion of the finite mass effects we plot, in Fig. 5, \(\sigma ^{(2 \Delta )}\) 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_1^{2\Delta }\,(F_2^{2\Delta })\) Eq. (25) (Eq. (26)). (ii) Computing \(F_1^{2\Delta }\) and \(F_2^{2\Delta }\) using their LET approximation as given in Eq. (16), while employing the exact expressions for \( F_1^{1\ell }\) and \(F_2^{1\ell }\) (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 \(\sim \)20–30 %. We remark that the double-triangle contribution to the hadronic cross section is actually always very small. Indeed, it amounts to \(\sim -0.18\) fb while the NLO cross section is around 40 fb.

We turn now to discuss the contribution in \(\mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\) due to \( F_1^{2\ell }\) and \(F_2^{2\ell }\) in Eq. (13). We expect our results for \( F_1^{2\ell }\) and \(F_2^{2\ell }\) to be quite accurate for \(\sqrt{\hat{s}} \lesssim 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 \(\sigma ^{(0)}\, \mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\) at various order in the large-mass expansion.

In Table 2 we report the contribution of

$$\begin{aligned} \sigma ^{(0)}\, \mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}= & {} \frac{G_{\mu }^2\alpha _s^2(\mu _R)}{512\, (2\pi )^3} \frac{\alpha _s(\mu _R)}{\pi }\nonumber \\&\times \int _{t_-}^{t_+} 2\text {Re} \left( T_F\,F_{1}^{1\ell , \mathrm{full}}(T_F\,F\,_1^{2\ell , n})^{*} \right. \nonumber \\&\left. + T_F\,F_{2}^{1\ell , \mathrm{full}}(T_F\,F_2^{2\ell , n})^* \right) \end{aligned}$$
(32)

to the hadronic cross section for few values of an upper cut on the invariant mass of the two Higgs system,Footnote 8 \(M_{HH}^c\), and various orders in the expansion. In Eq. (32) \(F_{i}^{1\ell , \mathrm{full}},\, (i=1,2)\) indicates the exact expression of the one-loop form factor [31, 50], while \(F_i^{2\ell , n}\) is for 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_t^8\) one, we find that the mass effects induce a relative variation with respect to the \(m_t\rightarrow \infty \) result up to \(\sim \)40 %.

Fig. 5
figure 5

Double-triangle contribution to the partonic cross section as a function of the partonic center-of-mass energy. The solid line represents the exact result using Eqs. (25) and (26), while the dashed one the result obtained in the LET approximation using Eq. (16)

Full size image
Table 2 Contribution (in fb) of \(\sigma ^{(0)}\, \mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\) as defined in Eq. (32) to the hadronic cross section for few values of an upper cut on the invariant mass of the two Higgs system (in GeV)
Full size table
Table 3 As Table 2 but with \(\sigma ^{0}\mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\), computed factorizing the LO cross section (see text)
Full size table

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\rightarrow \infty \) determination of the hadronic cross section. Applying the same procedure to a large-mass expansion determination amounts to evaluate \(\sigma ^{(0)}\, \mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\) employing for \(\sigma ^{(0)}\) the exact LO cross section while evaluating \(\mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\) at the same order of approximation both in the numerator and in the denominator. The contribution to the hadronic cross section of \(\sigma ^{0} \mathcal{C}_{\scriptscriptstyle \mathrm{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\rightarrow \infty \) result has the tendency to overestimate the NLO cross section as approximated by the \(1/m_t^8\) 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_t^2\) term. As expected, for \(M_{HH}^c > 2 \, m_t\) the convergence starts to downgrade.

Our analysis cannot say anything about the region of partonic energies \(\sqrt{\hat{s}} \gtrsim 400\) GeV where our results cannot be trusted. Concerning the hadronic cross section we can only make a guess assuming that in both regions, \(\sqrt{\hat{s}} \lesssim 400\) GeV and \(\sqrt{\hat{s}} \gtrsim 400\) GeV, the variation induced by mass effects in \(\sigma ^{0}\mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\) will be of a similar size and behavior 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 \(\sim 20~\%\). We notice that the contribution of \(\sigma ^{0}\mathcal{C}_{\scriptscriptstyle \mathrm{NLO}}\) to the NLO cross section is about \(10~\%\) of the total, that the contributions we did not discuss, i.e. \(\mathcal{R}_{gg}, \, \mathcal{R}_{q \bar{q}}\) and \(\mathcal{R}_{qg}\), when evaluated in the limit of vanishing external momenta contribute to the total NLO cross section by \(\sim 2~\%\), and that, according to the analysis in Refs. [46, 47], 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\rightarrow \infty \) 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\rightarrow \infty \) limit with the LO term factorized, \(\sigma ^{\mathrm{NLO}}_{\mathrm{LET}} = 40.00\) fb, with the NLO cross section computed as in Table 1 with \(\sqrt{\hat{s}_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, which is computed evaluating \(F_i^{2\ell }\) in the region below \(\sqrt{\hat{s}_c}\) using the \(1/m_t^8\) order in the expansion, while above \(\sqrt{\hat{s}_c}\) employing the LET values, amounts to \(\sigma ^{\mathrm{NLO}}_{\hat{s}_c} = 37.86\) fb.

Finally, we comment on larger hadronic center-of-mass energies. At LO, the LET result, e.g. at \(\sqrt{s}=100\) TeV, approximates the true one worse than at \(\sqrt{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.

6 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 two-loop form factors in the amplitude were computed via an asymptotic expansion in the top mass up to and including terms \(\mathcal{O} (1/m_t^8)\). Analytic results are presented for both contributions. Before this work \(F_{1,2}^{2\ell }\) and \(F_{1,2}^{2\Delta }\) were known only in the \(m_t\rightarrow \infty \) limit [36].

Our results allow for a more precise evaluation of the NLO cross section for partonic energies up to \(\sqrt{\hat{s}} \simeq 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\rightarrow \infty \) 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 \(\mathcal{R}_{gg}, \, \mathcal{R}_{q \bar{q}}\), and \(\mathcal{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, \(\sqrt{\hat{s}} \lesssim 400\) GeV and \(\sqrt{\hat{s}} \gtrsim 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 \(\sqrt{\hat{s}} \lesssim 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 [48, 49]. 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. [48, 49]. (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. [48] nor at the level of differential factorization, i.e. before the integration over the Higgs pair invariant mass, as in Ref. [49]. The factorization of the LO cross section is well known to work fine in single Higgs production where the exact NLO result is known [32–35], however, there is no proof that the same happens also in double Higgs production. From the comparison of Tables 2 and 3 in Sect. 5 it seems that the differential factorization, which is expected to lead to a better result than the other possibility since it gives rise to a better-behaved integrand [49], when the LET result is employed tends to overestimate the result. Although a detailed comparison of our results with those of Refs. [48, 49] is not possible, we notice that our results in Tables 2 and 3 exhibit the same behavior with respect to the order of the expansion of the soft-virtual cross section of Ref. [49].

Finally, we would like to point out that our work should be seen as one of the first steps toward a complete calculation of the two-loop virtual corrections in Higgs-pair production that can be easily and efficiently implemented into a Monte Carlo generator. A complete calculation of the NLO corrections, requires one to address, besides the real contributions that were studied in Refs. [46, 47], the computation of the virtual corrections in the energy region \(\sqrt{\hat{s}} \gtrsim 400 \) GeV. These corrections are very difficult to compute but can be attacked either via a large momentum expansion calculation or by numerical methods.

Note added After we submitted the preprint version of this paper, Ref. [67] appeared on the web. In this paper the authors compute the NLO QCD corrections with a full top-mass dependence by numerical methods. According to the analysis of Ref. [67] the top-mass corrections, i.e. the difference between the exact and the LET results, are slightly larger than \(10~\%\), although with the same sign that we predicted. Concerning the total cross section evaluated including \(1/m_t^8\) terms we find, using the same inputs of Ref. [67], that the NLO cross section computed factorizing the LO term, i.e. as in Table 3 with no cut on \(M_{HH}\), amounts to \(\sigma ^{\mathrm{NLO}} = 39.31\) fb. Instead, if we evaluate the NLO cross section as in Table 1 with \(\sqrt{\hat{s}}_c = 350\, (400)\) GeV, i.e. evaluating \(F_i^{2\ell }\) in the region below \(\sqrt{\hat{s}_c}\) using the \(1/m_t^8\) order in the expansion while above \(\sqrt{\hat{s}_c}\) employing the LET values, we find \(\sigma ^{\mathrm{NLO}} = 36.03 \,(36.35)\) fb. The difference between this latter result and the number quoted in Ref. [67], \(\sigma ^{\mathrm{NLO}} =32.80\) fb, is within \(10~\%\), while in the former case, using the LO factorization, is somewhat larger. We notice that given the huge amount of computer time needed to produce the result of Ref. [67], the approach pursued in this paper of looking for approximate analytic expressions that can describe well the exact result, qualifies to be probably the only one that can be followed in order to construct a Monte Carlo generator at the NLO level.