Higgs-pair production via gluon fusion at hadron colliders: NLO QCD corrections

Higgs-pair production via gluon fusion is the dominant production mechanism of Higgs-boson pairs at hadron colliders. In this work, we present details of our numerical determination of the full next-to-leading-order (NLO) QCD corrections to the leading top-quark loops. Since gluon fusion is a loop-induced process at leading order, the NLO calculation requires the calculation of massive two-loop diagrams with up to four different mass/energy scales involved. With the current methods, this can only be done numerically, if no approximations are used. We discuss the setup and details of our numerical integration. This will be followed by a phenomenological analysis of the NLO corrections and their impact on the total cross section and the invariant Higgs-pair mass distribution. The last part of our work will be devoted to the determination of the residual theoretical uncertainties with special emphasis on the uncertainties originating from the scheme and scale dependence of the (virtual) top mass. The impact of the trilinear Higgs-coupling variation on the total cross section will be discussed.


Introduction
Since the discovery of a scalar resonance [1,2] with a mass of 125.09 ± 0.24 GeV [3] that is compatible with the Standard Model (SM) Higgs boson [4][5][6][7][8][9], the detailed study of the properties of this particle has been a high priority of the analyses at the Large Hadron Collider (LHC). Theoretical uncertainties are a limiting factor for the accuracies reachable at the LHC. This restriction can partly be compensated by increasing the diversity of processes involving the Higgs boson and a broader spectrum of Higgs couplings probed at the LHC. In order to test the nature of the Higgs boson, its self-interactions are of particular interest. It will be the first step towards an experimental reconstruction of the Higgs potential. This plays a crucial role as the origin of electroweak symmetry breaking within the SM. The initial processes that provide a direct sensitivity to the Higgs selfcouplings are Higgs-pair production processes. They involve the trilinear Higgs coupling at leading order (LO) [10][11][12][13][14]. These processes are complementary to indirect effects induced JHEP04(2020)181 by the Higgs self-interactions in radiative corrections to electroweak observables and single-Higgs processes [15,16] that are plagued by unknown interference effects with other kinds of New Physics.
The Higgs self-interactions are uniquely described by the SM Higgs potential where λ defines the self-interaction strength of the SM Higgs field. In unitary gauge, the Higgs doublet φ is given by with v ≈ 246 GeV denoting the vacuum expectation value (vev) and H is the physical Higgs field. In the SM, the self-interaction strength is given in terms of the Higgs mass M H by λ = M 2 H /v 2 . Expanding the Higgs field around its vev, the Higgs self-interactions, including the corresponding permutations, are uniquely determined as where λ H 3 (λ H 4 ) denotes the trilinear (quartic) Higgs self-coupling. While the quartic Higgs coupling λ H 4 cannot be probed directly at the LHC, due to the tiny size of the triple-Higgs production cross section [17][18][19][20][21], 1 the trilinear Higgs coupling can be accessed directly in Higgs-pair production. Higgs-boson pairs are dominantly produced in the loop-induced gluon-fusion mechanism gg → HH that is mediated by top-quark loops supplemented by a per-cent-level contribution of bottom-quark loops, see figure 1. There are destructively interfering box and triangle diagrams at LO with the latter involving the trilinear Higgs coupling [10,11]. The box diagrams provide the dominant contributions to the cross section. A rough estimate of the dependence of the cross section on the size of the trilinear coupling is given by the approximate relation ∆σ/σ ∼ −∆λ H 3 /λ H 3 in the vicinity of the SM value of λ H 3 . Therefore, in order to determine the trilinear coupling, the theoretical uncertainties of the corresponding cross section need to be small. Thus, the inclusion of higher-order corrections is mandatory. The QCD corrections are fully known up to next-to-leading order (NLO) [25][26][27] and at next-tonext-to-leading order (NNLO) in the limit of heavy top quarks [28][29][30]. While the NLO corrections are large, the NNLO contributions are of more moderate size. Very recently, the next-to-next-to-next-to-leading order (N 3 LO) QCD corrections have been computed in the limit of heavy top quarks resulting in a small further modification of the cross section [31][32][33]. This calculation uses the N 3 LO corrections to the effective Higgs and Higgs-pair couplings to gluons in the heavy-top limit (HTL) [34]. The higher-order QCD corrections increase the total LO cross section by about a factor of two. Recently, the full NLO results have been matched to parton showers [35,36] and the full NNLO results in the limit of heavy top quarks have been merged with the NLO mass effects and supplemented by the additional top-mass effects in the double-real corrections [37]. • λ Figure 1. Generic diagrams contributing to Higgs-boson pair production via gluon fusion. The contribution of the trilinear Higgs coupling is marked in red.
The goal of this paper is to present in detail the calculation of ref. [27] of the full NLO corrections to Higgs pair production in gluon fusion. We rely on a direct numerical integration of the Feynman diagrams, without any tensor reduction. We extend the results presented in ref. [27] and study not only the LHC at center-of-mass energies of 13 and 14 TeV, but also present numbers for a potential high-energy upgrade of the LHC (HE-LHC) at 27 TeV [38] and for a provisional 100 TeV proton collider within the Future-Circular-Collider (FCC) project [39,40]. Special emphasis will be given to the study of the theoretical uncertainties affecting the results and in particular the scale and scheme uncertainty related to the top-quark mass. We will also study the variation of the trilinear Higgs coupling and show that the NLO mass effects shift the minimum of the total cross section as a function of λ H 3 . They vary substantially over the range of λ H 3 values.
The paper is organized as follows. We present the notation of our calculation in section 2 and discuss the results at LO. In section 3 we move to the NLO QCD corrections. We discuss the details of the calculation of the virtual corrections in section 3.1. We describe the derivation of the real corrections in section 3.2. Our numerical analysis is performed in section 4. Finally, the conclusions are given in section 5.

Leading-order cross section
At LO, Higgs-boson pair production via gluon fusion is mediated by the generic diagrams of figure 1, including all permutations of the external lines. There are triangle and box diagrams with the former involving the trilinear Higgs coupling through an s-channel Higgs exchange. The LO matrix element of g(q 1 )g(q 2 ) → H(p 1 )H(p 2 ) can be cast into the form with Q = m HH denoting the invariant Higgs-pair mass. Here a, b denote the color indices of the initial gluons, 1/2 their polarization vectors, Γ H the total Higgs width, 2 G F the Fermi JHEP04(2020)181 constant and α s (µ R ) the strong coupling at the renormalization scale µ R . Since in this work we neglect the small bottom-quark contribution, the LO function of the triangle-diagram contribution is given by the top-quark contribution, with τ t = 4m 2 t /Q 2 and the basic function where m t denotes the top mass, while the more involved analytical expressions for F and G can be found in ref. [11]. In the HTL, the LO form factors approach the values There are two tensor structures contributing which correspond to the total angular-momentum states with S z = 0 and 2, where p T is the transverse momentum of each of the final-state Higgs bosons. Working in n = 4 − 2 dimensions, the following projectors on the two form factors can be constructed, Using these projectors, the explicit results of the two form factors F 1,2 can be obtained in a straightforward manner. The analytical expressions can be found in refs. [10,11]. Working out the polarization and color sums of the matrix element of eq. (2.1), the LO partonic cross sectionσ LO is given bŷ with the integration boundarieŝ

JHEP04(2020)181
where the symmetry factor 1/2 for the identical Higgs bosons in the final state is taken into account. The LO hadronic cross section σ LO can then be derived by a convolution with the parton densities with the gluon luminosity, given in terms of the gluon densities g(x, µ F ), at the factorization scale µ F and the integration boundary τ 0 = 4M 2 H /s, where s denotes the hadronic center-of-mass (c.m.) energy squared. The differential cross section with respect to the invariant squared Higgs-pair mass Q 2 can be obtained as As can be expected from single Higgs-boson production via gluon fusion (see [41][42][43][44][45]), the NLO QCD corrections to these LO expressions will be large.

Next-to-leading-order corrections
The NLO QCD corrections to Higgs-pair production via gluon fusion have been computed in the HTL, a long time ago [12]. The NLO result for the gluon-fusion cross section can be generically expressed as [12] σ NLO (pp → HH + X) = σ LO + ∆σ virt + ∆σ gg + ∆σ gq + ∆σ qq , withσ LO (Q 2 ) denoting the partonic cross section at LO and the strong coupling α s (µ R ) is evaluated at the renormalization scale µ R . The objects dL ij /dτ (i, j = g, q,q) denote the JHEP04(2020)181 parton-parton luminosities, defined analogously to dL gg /dτ of eq. (2.11), using the quark densities q(x, µ F ), at the factorization scale µ F and P ij (z) (i, j = g, q,q) are the specific Altarelli-Parisi splitting functions [46].
The quark-mass dependence is in general encoded in the LO cross sectionσ LO (Q 2 ) and the terms C virt (Q 2 ), d ij (Q 2 , z) for the virtual and real corrections, respectively. These expressions can easily be converted into the differential cross section with respect to Q 2 , while the differential cross section at LO is given in eq. (2.12).
Within the HTL, the Higgs coupling to gluons can be described by an effective Lagrangian [42,[47][48][49][50] involving the Wilson coefficients (L t = log µ 2 R /m 2 t ) [12,30,34,[51][52][53][54][55]] that are known up to N 4 LO [34,53,54]. Since the top quark is integrated out, the number of active flavours has been chosen as N F = 5. If these effective Higgs couplings to gluons in the calculation of the NLO QCD corrections are used, the calculation of these is simplified to a one-loop calculation for the virtual corrections and a tree-level one for the matrix elements of the real corrections. The terms C virt (Q 2 ) and d ij (Q 2 , z), for the virtual and . Typical two-loop triangle (left), one-particle reducible (middle) and box (right) diagrams contributing to Higgs-pair production via gluon fusion at NLO. real corrections, approach in the HTL the simple expressions whereŝ,t,û (ŝ = Q 2 at LO and for the virtual corrections) denote the partonic Mandelstam variables and C is the contribution of the one-particle reducible diagrams, see figure 2. At NLO QCD, the full mass dependence of the LO partonic cross section has been taken into account, while keeping the virtual corrections C virt and the real corrections d ij in the HTL ("Born-improved" approach) [12]. This yields a reasonable approximation for smaller invariant Higgs-pair masses and approximates the full NLO result of the total cross section within about 15% [25][26][27]. The NLO QCD corrections in the HTL increase the cross section by 80 − 90% [12]. Within the Born-improved HTL, the NNLO QCD corrections have been obtained in refs. [28][29][30] increasing the total cross section by a moderate amount of 20 − 30% [29]. Beyond these NNLO QCD corrections, the soft-gluon resummation (threshold resummation) has been performed at next-to-next-to-leading logarithmic (NNLL) accuracy for the total cross section and invariant mass distribution, modifying the total cross section further by a small amount if the central scales are chosen as µ R = µ F = Q/2 [56,57]. Very recently, the N 3 LO QCD corrections have been computed in the Born-improved HTL resulting in a small modification of the cross section beyond NNLO [31][32][33][34]. These N 3 LO QCD corrections in the HTL have been merged with the full top-mass effects of the NLO calculation [33].
The calculations in the HTL have been improved by several steps including mass effects partially at NLO. The full mass effects in the real correction terms d ij have been included by means of the full one-loop real matrix elements for gg → HHg, gq → HHq, qq → HHg. This improvement reduces the Born-improved HTL prediction for the total cross section by about 10% [58,59] and is called the "FTapprox" approximation. The calculation of the full real matrix elements has been performed by using the MG5 aMC@NLO framework [60,61].

JHEP04(2020)181
Another improvement has been achieved by an asymptotic large-top-mass expansion of the full NLO corrections at the level of the integral [62] and the integrand [63]. This indicated sizable mass effects in the virtual two-loop corrections alone. In addition, the large top-mass expansion has been extended to the virtual NNLO QCD corrections resulting in 5% mass effects estimated on top of the NLO result [63]. The large-top-mass expansion of the NLO QCD corrections has been used to perform a conformal mapping of the expansion parameter and to apply Padé approximants. In this way, an approximation of the full calculation has been achieved for Q values up to about 700 GeV [64]. Another approximation builds on an expansion in terms of a variable that dominantly corresponds to the transverse momentum of the Higgs bosons. The results of this approach show good agreement with the full calculation for Q values up to about 900 GeV [65]. Analytical results are also available in the large-Q limit [66]. The latter have recently been combined with the numerical results of refs. [25,26] for the full QCD corrections [67]. In the following, we will discuss the details of our NLO calculation.

Virtual corrections
Typical diagrams of the two-loop virtual corrections are shown in figure 2. They can be arranged in three different classes: (a) triangle, (b) one-particle-reducible and (c) box diagrams. 3 They contribute to the coefficient C virt (Q 2 ) of eq. (3.1), where ∆F , ∆F and ∆G denote the virtual corrections to the corresponding LO form factors. While ∆F involves only virtual corrections to the triangle diagram, ∆F and ∆G acquire contributions from the one-particle-reducible and box diagrams.

Triangle diagrams
The generic 2-loop triangle diagrams contributing to the virtual coefficient C virt (Q 2 ) are shown in figure 3. They only contribute to the spin-0 form factor F 1 of eq. (2.1) and can be parametrized as the correction ∆F to the form factor F , where C virt (Q 2 ) denotes the complex virtual coefficient relative to the LO form factor F of the amplitude. This virtual coefficient is related to the single-Higgs case so that the relative QCD corrections can be simply obtained from the known (complex) virtual JHEP04(2020)181 (3.9) In the HTL, this virtual coefficient (before renormalization) approaches the expression with the 't Hooft scale µ 0 , where the (infinitesimal) regulator¯ defines the proper analytical continuation of this expression. This result has to be followed by the renormalization of the strong coupling α s and the top mass m t that will be discussed in section 3.1.4. In addition, we have subtracted the HTL to obtain the pure top-mass effects at NLO (relative to the massive LO expression F ) to ensure that in the end the results of the program Hpair [68] can be added back. This last step will be discussed in section 3.1.4, too.
3.1.2 One-particle-reducible diagrams The one-particle-reducible contribution is depicted in figure 2 (middle diagram), where a second diagram with the initial gluons interchanged has to be added. These will constitute JHEP04(2020)181 thet-andû-channel parts where the second is related to the first just by the interchangê t ↔û [see C of eq. (3.6)]. The analytical expression of the coefficient c 1 can be related to the top contribution of the process H → Zγ [69,70]. The basic building block will be the one-loop contribution of the Higgs coupling to an on-shell and an off-shell gluon that is described, after translating all couplings and masses, by the "effective" Feynman rule, where the functions I 1,2 are defined as [71] with τ = 4m 2 t /m 2 H , λ = 4m 2 t /q 2 2 and the basic functions and f (τ ) defined in eq. (2.3). Implementing this building block for the two top loops of the one-particle-reducible diagrams, one arrives at the final coefficient c 1 of eq. (3.6), with λt = 4m 2 t /t (and λû = 4m 2 t /û for thet ↔û interchanged contribution accordingly). This expression, inserted in the coefficient C of eq. (3.6), determines the contribution of the one-particle-reducible diagrams analytically and agrees with the previous calculation of ref. [72]. In the HTL, this coefficient approaches the value c 1 → 2/9 in accordance with eq. (3.6). We have subtracted the HTL with c 1 = 2/9 from the coefficient C in order to account for the NLO top-mass effects only so that eventually the results of the program Hpair [68] can be added back. While the total effect of the one-particle-reducible contributions on the total cross section ranges below the per-cent level, the finite mass effects at NLO contribute less than one per mille.
Reference [73] has proposed an approximation of this one-particle-reducible contribution in terms of the triangle form factor of two on-shell external gluons,  Comparison of the approximation of ref. [73] (blue) for the one-particle-reducible contributions and the HTL (red), both normalized to the full analytical expression. The singularity at about 720 GeV is due to a sign change of the exact expression.
2) evaluated at half the invariant Higgs-pair mass Q/2 instead of Q], where the function F can be found in eq. (2.2). Since ref. [73] works in the HTL, the contribution of the second form factor F 2 vanishes, i.e. G → 0, and the approximation V 2 eff /2 is in fact treated as an approximation for the coefficient c 1 of the exact expression of C as given in eq. (3.6). 5 Thus, the approximate expression involving the coefficient c 1 has to be compared to the corresponding expression involving the exact coefficient c 1 of eq. (3.13). This comparison is presented, normalized to the exact expression, in figure 4 and shows that the approximation of ref. [73] is not better than the HTL.

Box diagrams
The third class of two-loop contributions to the virtual corrections is given by the box diagrams. The generic box diagrams are shown in figures 19-21 in the appendix. The simultaneous exchange of the gluons and Higgs bosons has to be added to complete the set of diagrams. The only exception is diagram 44 that is already totally symmetric so that in the final end there are 93 two-loop box diagrams. The generic 47 diagrams are grouped into 6 topology classes. The first 5 topologies contain only a virtual threshold for Q 2 > 4m 2 t . The diagrams of topology 6 on the other hand develop a second threshold for Q 2 > 0, because two virtual gluon lines next to the external gluons can be cut. This implies that the form factors are complex in the entire Q 2 range. Therefore, a dedicated 5 Since V eff is symmetric with respect tot ↔û the additional factor 2 emerges from the second term in the numerator of C in eq. (3.6). Figure 5. Explicit definitions of the virtual momenta in box 39. treatment of this last topology in terms of a suitably constructed infrared subtraction term to isolate the associated infrared singularities is required.

JHEP04(2020)181
In the following, we will exemplify our method for the boxes 39 of topology 5 and 45 of topology 6. The diagrams of topologies 1-5 are treated analogously to box 39 and those of topology 6 analogously to box 45. The algebraic manipulation of the traces and projections onto the form factors have been performed with the help of the symbolic tools FORM [74,75], Reduce [76], and Mathematica [77]. Our method of Feynman parametrization and endpoint subtraction to isolate the ultraviolet singularities for the numerical integration has first been applied to the NLO two-loop QCD corrections to H → γγ, Zγ in refs. [78,79] and later to the squark-loop contributions to h, H ↔ gg, γγ within the minimal supersymmetric extension of the SM [80]. The method of the infrared subtraction as applied to topology 6 originates from numerical cross checks of the full NLO QCD corrections to single Higgs production in refs. [41,42,80,81]. The stabilization of virtual thresholds by integration by parts of the integrand has first been applied to the SUSY-QCD corrections to single Higgs production in refs. [82,83]. The basic idea behind the integration by parts is to reduce the power of the threshold-singular denominator and in this way to stabilize the numerical integration. The treatment of the thresholds in our approach is performed by replacing the squared top mass m 2 t by a complex counter part with a positive regulator¯ > 0 to ensure proper micro-causality. This defines the analytical continuation of our two-loop box integrals. In the following, the parameter¯ will be kept finite in our numerical analysis, while the narrow-width limit¯ → 0 is achieved by a Richardson extrapolation [84]. This will be discussed in more detail in the following paragraphs.
Box 39. Using the definition of real and virtual momenta as in figure 5, the contribution to the tensor A µν [see eq. (2.1)] of the virtual two-loop corrections is given by where k, q are the loop momenta that are integrated over. The Feynman parametrization is first performed for the integration over k. We provide Feynman parameters x 1 , . . . , x 4 for the first four propagators in the denominator and 1 − i x i for the last one (k 2 − m 2 t ). Performing the substitutions we arrive at a four-dimensional integral over x, y, z, r with integration boundaries from 0 to 1. To symmetrize the n-dimensional k-integration, we have to perform the shift in both the numerator and denominator. The residual (properly normalized) denominator after the k-integration is treated as a propagator for the second loop integration over q.
We attribute additional Feynman parameters x 5 , x 6 to this residual propagator and the next one [(q − q 1 ) 2 ] and 1 − x 5 − x 6 for the last one (q 2 ) in eq. (3.16). Performing the substitution 6 we again arrive at integrals over s, t from 0 to 1. This latter parametrization requires the shift in the numerator and denominator to be able to perform the loop integration over q symmetrically. After projecting on the two form factors, we finally arrive at integrals of the type with x = (x, y, z, r, s, t) and d 6 x = dx dy dz dr ds dt. H i ( x) denotes the full numerator, including regular factors of the Jacobians due to the Feynman parametrization and substitutions, and singular as well as higher powers of the dimensional regulator , and N ( x) the final denominator, The singular powers in of H i ( x) arise from powers of k 2 and q 2 in the numerators of the final integrations of the loop momenta k and q. It is important that the final denominator develops the form of 1 + O(1/m 2 t ) to ensure that no further ultraviolet nor infrared singularities arise from this part of the integrand.
The integral for ∆F i of eq. (3.21) is singular for s → 0. To separate this singularity from the integral, we perform an endpoint subtraction, where in the second term ∆F i,2 the integration over s has been performed analytically and the integration measure is given by d 5 x = dx dy dz dr dt. It should be noted that in the terms L, L 0 , L 1 the logarithms of the denominator N need to be linear in N to be consistent with the analytical continuation along the proper Riemann sheet. We have checked numerically that the first (subtracted) part ∆F i,1 is finite for each order in the dimensional regulator by introducing cuts in the integration boundaries, i.e. integrating from˜ to 1 −˜ , varying˜ down to 10 −10 and checking that the integrals become independent of˜ .
These integrals are numerically stable below the virtual tt-threshold, i.e. for Q 2 < 4m 2 t or ρ s < 4. However, above this threshold, the integrals have to be stabilized. We have achieved this stabilization by means of integration by parts with respect to the Feynman parameter z. The denominator is a quadratic polynomial in z, To simplify the integration by parts, we insert a unit factor ∆/∆ with ∆ = 4ac − b 2 in the integrand and replace ∆ in the numerator by the expression

JHEP04(2020)181
Then the following manipulation can be performed, and analogously for integrals involving additional powers of log N factors in the numerator of the integrand. The progress achieved with these integrations by parts is that the maximal power of the denominator in the new integral is reduced by one compared to the original integral. One could perform additional integrations by parts with respect to another Feynman parameter. However, we did not investigate this further, since the stability we achieved at this point has been sufficient for the numerical integrations for the top loops. 7 After performing the integrations by parts, the integral is stable for regulators¯ [see eq. (3.15)] down to 0.05 for the relevant Higgs mass, top mass and Q 2 range. Since this is still apart from the plateau of the narrow-width limit, we performed a Richardson extrapolation [84] from finite values of¯ down to zero. Richardson extrapolation is possible since the¯ -dependence of the integral is polynomial for small values of¯ . The basic principle behind this extrapolation method is very simple: let a function f (¯ ) behave for small¯ as If we know f (¯ ) for two different values¯ and t¯ , we can construct the new function This function shows a better convergence towards the value at¯ = 0, Our integrals I(¯ ) behave for small values of¯ as so that the first new extrapolation function in our case is given by Using an additional value of¯ , this method can be repeated iteratively for the new function obtained by applying eq. (3.28), Figure 6. Explicit definitions of the virtual momenta in box 45. In this way, the estimated error is reduced by each additional iteration. We have used this method for a set of¯ separated by factors of t = 2. Then, we obtain the following extrapolation polynomials, and so on. We have used extrapolation polynomials up to R 9 (¯ ). To determine the extrapolation error, we have chosen different sets of¯ values and derived the spread of the extrapolated values appropriately (see section 4 for more details).
Box 45. Based on the distribution of the loop and external momenta of figure 6, the contribution to the two-loop matrix element is given by Following the same procedure as for box 39 for the Feynman parametrization, we have first performed the parametrization of the k-integration following the ordering of the denominator of eq. (3.34). The shift in the loop momentum k and the corresponding substitutions of the Feynman parameters are given by

JHEP04(2020)181
Performing the second loop integration over q with the residual (normalized) denominator of the k integration as the first propagator of the q integration, attributing the additional Feynman parameters x 4 , x 5 , x 6 to the remaining propagators in eq. (3.34) and applying the substitutions 8 we arrive at the final expressions for the shift of q and the denominator that contribute to the two form factors, and the final integrals of the two form factors (i = 1, 2) can be cast into the form where H i ( x) contains all additional regular Feynman-parameter factors from Jacobians and the normalization of the denominator of the first loop-integration over k. It develops a singular Laurent-expansion in . The final denominator exhibits the basic form of r + O(1/m 2 t ), so that the additional singular behavior is entirely controlled by the limit of small r. Since the denominator is of the form with a, c = O(1/m 2 t ) and b = 1 + O(1/m 2 t ) and the infrared singularities are universal (relative to the LO expressions) the coefficient a does not contribute to the infrared singularity structure, because a is subleading relative to b in the limit r → 0. Thus, we can construct infrared subtraction terms that turn the contributions to the form factors into

JHEP04(2020)181
Numerically, we have tested that the subtracted integral G 1 (after expansion in the dimensional regulator ) is finite for each coefficient of the expansion in individually by integrating the Feynman-parameter integrals from˜ to 1 −˜ with˜ varied down to 10 −10 . The second integral G 2 can be integrated over the Feynman parameter r analytically giving rise to hypergeometric functions, with d 5 x = dx dy dz ds dt. Since this integral is singular for c → 0, we have to invert the last argument of the hypergeometric function. Using the transformation relation the special property 2 F 1 (a, 0; c; z) = 1 (3.43) and suitable end-point subtractions of the residual singular integrals analogous to box 39, we arrive at the final decomposition of the initial Feynman-parameter integral
Box 45 contains a second threshold for Q 2 > 0 so that even below the tt-threshold, integrations by parts are required to stabilize the integrand numerically. These integrations by parts are performed for the Feynman parameter r in the contributions S 1,2 along the same lines as for box 39, while the integrals S 3−6 are stable without integrations by parts.

Renormalization
The strong coupling α s has been renormalized in the MS scheme with the top quark decoupled, i.e. the renormalization constant is given by where the LO form factors F i have to be used in n dimensions, i.e. including higher orders in the dimensional regulator . For our default prediction, we have renormalized the top mass on-shell so that the renormalization constant is given by The explicit contribution of the mass counterterm can either be obtained by calculating the corresponding counterterm diagrams or, in much more elegant manner, by differentiating the LO form factors with respect to the top mass, where we followed the second option. For the renormalization of the top mass in terms of the MS mass, a counterterm

JHEP04(2020)181
has to be used with the LO and NLO expressions of the form factors expressed in terms of the MS top mass m t (µ t ). For the evaluation of the MS top mass, we use the N 3 LO relation between the pole and MS mass [85][86][87][88], with K 2 ≈ 10.9 and K 3 ≈ 107.11. The scale dependence of the MS mass is treated at N 3 LL, with the coefficient function [89,90] c Since we are interested in the finite top-mass effects on top of the LO ones, we have subtracted in addition the Born-improved HTL of the virtual corrections involving the full top-mass dependence at LO [12]. This yields the additional subtraction term (3.55) After adding this subtraction term, the result of Hpair can simply be added back to the NLO top-mass effects obtained in this way for the virtual corrections. Thus, the total counterterm plus HTL-subtraction is given by (3.56) The addition of this term results in an infrared and ultraviolet finite result for the virtual corrections as we have explicitly checked numerically. It should be noted that we have defined this total subtraction term with the imaginary part¯ for the top mass to be consistent with our treatment of the two-loop diagrams. For the two-loop triangle diagrams, this total subtraction term is included in the narrow-width approximation according to the known result for the single-Higgs case.

Differential cross section
The final numerical integrations have been performed by Vegas [91] for the differential cross sections dσ/dQ 2 of eq. (3.3), i.e. the integration overt is included. Each individual box diagram is divergent int at the lower and upper bound of thet-integration in general.
To stabilize thet-integration, we have performed a suitable substitution to smoothen the integrand,t JHEP04(2020)181 where the integration boundariest ± are given in eq. (2.9). By means of this substitution, we can rewrite the integration over t 1 generically as 9 where f (t 1 ,û 1 ) denotes the corresponding virtual matrix element with the (singular) denominatort 1û1 −ŝM 2 H extracted and the integration boundaries read where we have introduced a cut˜ for the upper and lower bound of thet 1 -integration (after rewriting this into an integral from 0 to 1 and replacing these integration boundaries by˜ and 1 −˜ ). We have checked that the total sum of all box diagrams becomes independent of this cut by varying˜ down to 10 −10 , i.e. that the total sum is again finite. 10

Real corrections
We are left with the evaluation of the real contributions to complete the picture of the NLO QCD corrections. As we are interested in the calculation of the top-mass effects on top of the HTL calculation that is provided by Hpair, we use the universality of the infrared divergent pieces to subtract the Born-improved HTL contributions dσ HTL ij in such a way that our integration of the real contributions d∆σ mass ij = dσ ij − dσ HTL ij is finite. We construct a local subtraction term for the partonic channels dσ ij , where p k denote the four-momenta from the full 2 → 3 phase-space andp k stand for the mapping of the momenta p k on a 2 → 2 sub-phase-space. As the results in the HTL limit are given in the Born-improved approximation in which the pure HTL is rescaled with the full LO matrix elements, we need to map the full 2 → 3 phase-space onto a projected 2 → 2 phase-space to construct the subtraction term involving this rescaling to the full LO contribution dσ LO . The mapping is done by using the transformation formulae for initial-state emitter and initial-state spectator in the construction of dipole subtraction terms, i.e. using eqs. (5.137-5.139) of ref. [92]. The (mapped) momenta of the initial-state partons are p 1/2 JHEP04(2020)181 (p 1/2 ), the (mapped) momenta of the final-state Higgs bosons are p 3/4 (p 3/4 ), and the momentum of the radiated parton is p 5 . For the initial-state partons, we use the following mapping,p In order to transform the Higgs momenta, we introduce the variables K andK, allowing us to definep The HTL matrix elements are calculated analytically. We introduce the partonic centerof-mass energyŝ, and the Mandelstam variablest = (p 1 − p 5 ) 2 andû = (p 2 − p 5 ) 2 . The invariant squared Higgs-pair mass is Q 2 =ŝ +t +û. The real spin-and colour-averaged matrix elements are (3.65) The full one-loop matrix elements have been generated with FeynArts [93] and FormCalc [94]. They contain triangle, box, and pentagons diagrams. Generic diagrams for the contribution gg → HHg are given in figure 8, generic diagrams for the contributions qg → HHq and qq → HHg are displayed in figure 9. The numerical evaluation of the scalar integrals [95] as well as the tensor reduction has been performed using the techniques developed in refs. [96][97][98][99] and implemented in the library Collier 1.2 [100]. The latter has been interfaced to the analytic expressions generated by FormCalc with an in-house routine. In order to improve our numerical stability, we have implemented a technical collinear cut in the phase-space parametrization. The integration of the scattering angle θ of the radiated parton in the c.m. system is restricted to the range |cos θ| < 1 − δ with δ = 10 −4 . We have checked that our results are stable against a variation of δ from 10 −4 to 10 −6 and therefore they are not affected by our choice for this technical cut. We have cross-checked the final mass-effects of the real corrections against the results presented in the literature [25,26,58,59] and we have obtained agreement.

Results
Our numerical results will be presented for the invariant Higgs We define an estimate of the theoretical error due to the Richardson extrapolation as the difference of the extrapolated results at fifth and fourth order. In addition, we multiply this error by a factor of two close to the virtual tt threshold in order to be conservative. The total estimated Richardson-extrapolation error ranges below the per-cent level and is added in quadrature to the statistical integration error. Since we have subtracted the (Born-improved) HTL consistently from the virtual and real corrections, we are left with the pure top-mass effects at NLO that are infrared and ultraviolet finite individually after renormalization. This part has then been added to the results of Hpair [68] to derive the full NLO cross section. The final invariant Higgs-pairmass distributions are displayed in figures 10-12 for the three c.m. energies, 14, 27, 100 TeV. The blue curves show the Born-improved result in the HTL of ref. [12] as implemented in Hpair [68], the yellow ones the Born-improved HTL result plus the mass effects of the real corrections, the green curves the Born-improved HTL result plus the mass effects of the virtual corrections and the red curves the full NLO results. The plots on the left side of JHEP04(2020)181 each figure have been obtained by using MMHT2014 PDFs [101] and the ones on the right with PDF4LHC PDFs [102]. The lower panel on the left shows the consistently defined Kfactors K = dσ NLO /dσ LO . The lower panel on the right shows the ratio of the differential NLO cross section to the one obtained in the Born-improved HTL.
While the Born-improved HTL provides a reasonable approximation for Q-values close to threshold, the real corrections add a negative mass effect of about −10% for √ s = 14 TeV (yellow curves) that is approximately uniform in the entire Q range. The (negative) mass effects of the virtual corrections (green curves), however, become large at large values of Q reaching a level of more than 20% for Q beyond about 1 TeV. While the relative mass effects of the virtual corrections at NLO are independent of the collider energy (see the right plots showing the ratios to the HTL in the lower panels) in agreement with eq. While (as for the ratios) the full NLO K-factors shown in the left plots are close to the Born-improved HTL (blue curves) at Q values close to the production threshold, they deviate significantly at larger values of Q due to the additional NLO top-mass effects that decrease the total size of the NLO QCD corrections compared to the HTL as expected from unitarity arguments.
To estimate the theoretical uncertainties, we have varied the renormalization and factorization scales for each bin in Q by a factor of 2 up and down around the central scale µ R = µ F = Q/2 and derived the envelope of a 7-point variation, i.e. excluding points where the renormalization and factorization scales differ by more than a factor of two. The residual uncertainties are shown by the red band around the full NLO results (red curves) in figures 10-12. They range at the level of 10-15% in total as can be inferred from the explicit numbers for √ s = 14 TeV (using PDF4LHC PDFs), dσ NLO dQ Q=300 GeV = 0.02978(7) +15.3% −13.0% fb/GeV, dσ NLO dQ Q=400 GeV = 0.1609(4) +14.4% −12.8% fb/GeV, dσ NLO dQ Q=600 GeV = 0.03204(9) +10.9% −11.5% fb/GeV, dσ NLO dQ Q=1200 GeV = 0.000435(4) +7.1% −10.6% fb/GeV . We have analyzed the structure of the NLO QCD corrections in more detail by comparing the K-factor with the one of the triangle diagrams alone, i.e. with the K-factor of single-Higgs production with mass M H = Q, in all individual approximations. This will determine the amount of universal NLO top-mass effects, common in the triangle and box diagrams. We define the ratio of the NLO triangle-diagram K-factor to the one including all diagrams as K-fac /K-fac. This is shown, as a function of Q = m HH , in figure 13 (left). It is visible that the triangle-diagram K-factor provides an acceptable approximation to      the full NLO K-factor only for Q values below about 500-600 GeV if maximal deviations of about 15% are allowed (red histogram). The break down into the different mass effects of the virtual (green histogram) and real (yellow histogram) corrections singles out the origin of non-universal mass effects in the virtual corrections, while the non-universal mass effects beyond the single-Higgs case of the real corrections are limited to less than about 5% (apart from the virtual tt-threshold region). In comparison to the contribution of the triangle diagrams alone, we also present the ratio of the K-factor obtained by including only the continuum diagrams (box diagrams of the virtual corrections and all box and pentagon diagrams of the real corrections without trilinear Higgs couplings) to the full K-factor in figure 13 (right). The different curves show the results for the various approximations, i.e. the blue curves for the Born-improved HTL, the yellow ones with the inclusion of the NLO mass effects of the real corrections, the green curves with only the virtual NLO mass effects and the red curves the full NLO results. The right figure shows that the full NLO K-factor (red curve) is well-described (within 5%) by the one for the continuum diagrams alone which coincides with the observation that the continuum diagrams play a significant role for small values of Q (where the K-factor does not deviate much from the single-Higgs case) and are dominant for large Q. This result shows that the K-factor cannot be approximated well by the one of single-Higgs production for large values of Q due to the large mass effects of the virtual corrections.

Total cross section
The total cross section has been obtained from the invariant Higgs-pair mass distribution by means of a numerical integration of the bins in Q with the trapezoidal method for Q > 300 GeV. For a reliable result, we used a Richardson extrapolation [84] in terms of the bin size in Q also for this step. For Q < 300 GeV, we have adopted the extension of Boole's rule to six nodes [104]. We obtain the following values for the total cross section at various c.m. energies, Comparing the results of eqs. (4.2) and (4.3), we observe a reduction of the total cross section by about 15% due to the top-mass effects at NLO and a reduction of the scale uncertainty. These numbers, as well as the differential distributions presented in section 4.1, agree with the results of refs. [25,26]. 13 It should be noted that a comparison of the full JHEP04(2020)181 virtual corrections with the analytical large top-mass expansion presented in ref. [63] was performed in refs. [25,26] and shows a convergence to the full result below the tt-threshold, as expected.

Uncertainties originating from the top-mass definition
An uncertainty that has been neglected or underestimated often previously is the intrinsic uncertainty due to the scheme and scale choice of the virtual top mass. This does not play a large role for single on-shell Higgs-boson production via gluon fusion, gg → H, since the Higgs mass is small and thus the HTL works well, i.e. top-mass effects are suppressed. This uncertainty, however, plays a significant role for the larger values of Q in Higgs-pair production. Top-mass effects are already sizeable at LO, but the NLO corrections add additional relevant top-mass dependences on top of the LO result as we have discussed in the previous subsection. The top mass is a scheme and scale dependent quantity so that the related uncertainties need to be estimated for a reliable determination of the total theoretical uncertainties. For this analysis, we have evaluated the differential cross section for the top mass defined in the on-shell scheme (default) and in the MS-scheme at the scale µ t , i.e. adjusting the counterterms and input parameters to the choices m t (m t ) and m t (µ t ) with µ t in the range between Q/4 and Q according to section 3.1.4. 14 Since the scale dependence on µ t is a monotonously falling function, we evaluated the differential cross section for four choices of the top mass, m t , m t (m t ), m t (Q/4) and m t (Q), for each bin in Q.
For the three c.m. energies of 14, 27 and 100 TeV the differential cross sections are presented in figures 14, 15 as a function of Q = m HH for the various definitions of the top mass. The lower panels exhibit the ratios of the differential cross sections to the ones in terms of the top pole mass (OS scheme). It is clearly visible that the scale and scheme dependence of the top mass induces sizeable variations of the NLO Higgs-pair production cross section and thus contributes to the theoretical uncertainties. For small Q values, the size pattern of the differential cross section due to the different scale and scheme choices is varying. For large values of Q, the maximum is always given by the on-shell scheme and the minimum in terms of the MS-top mass m t (Q) with sizeable differences to the on-shell scheme. Adopting the related uncertainties as the envelope of the cross sections for our four choices, we arrive at the following uncertainties of the differential cross section for a c.m. energy √ s = 14 TeV, Since these uncertainties are given relative to the on-shell results, the upper uncertainty vanishes for Q ≥ 400 GeV, because the on-shell results provide the maximal values. These uncertainties turn out to be significant and at a similar level as the usual renormalization and factorization scale uncertainties. Thus, they constitute an additional contribution to the total theoretical uncertainties that has to be taken into account. The uncertainties due to the top-mass scheme and scale are about a factor of two smaller than at LO,  using PDF4LHC PDFs. A further reduction of these uncertainties can only be achieved by the determination or reliable estimate of the full mass effects at NNLO. Since these uncertainties are sizeable, one may wonder why this has not been observed already for single-Higgs boson production gg → H. The measured value of the Higgs mass M H = 125 GeV is small compared to the top mass so that for single on-shell Higgs production we are close to the HTL, i.e. finite top-mass effects are small and thus the related uncertainties, too. However, going to larger virtualities Q for off-shell Higgs production JHEP04(2020)181 that have been obtained with PDF4LHC PDFs as in the Higgs-pair case. On the other hand, the uncertainties for Q = 125 GeV confirm that they are small for on-shell Higgs production via gluon fusion (already at LO) in agreement with the analysis of the LHC Higgs Cross Section Working Group [105,106]. A relevant issue is the theoretical background of the different scale choices for the top mass. For small values of Q, the matrix element will be closer to the HTL such that the NLO corrections get closer to the HTL calculation. The HTL on the other hand can be treated by starting from the effective Lagrangian of eq. (3.4) which is the residual effective coupling of Higgs bosons to gluons after integrating out the top quark. Thus, the corresponding Wilson coefficients C 1 and C 2 are determined by matching the full SM with the top quark to the effective theory without the top quark. The matching scale is naturally given by the top mass. Performing the proper matching at the scale of the top mass, i.e. using either the top pole mass or the top MS mass at the scale of the top mass itself leads to non-logarithmic (in the top mass) matching contributions [see eq. (3.5) for µ R = m t ] also for higher powers in 1/m 2 t , i.e. higher-dimensional operators contributing to the gluonic Higgs couplings at the subleading level. This implies that the top mass is the preferred scale choice for small values of Q. This is confirmed by the heavy top expansion of the form factors of refs. [62,63,66].
At large Q values, on the other hand, we can use the results for the high-energy expansion of ref. [66]. In the regime of large Q, the triangle-diagram contributions are suppressed by the s-channel Higgs propagator so that the box diagrams provide the dominant contributions. In our normalization, the explicit results of the virtual box-form factors in the JHEP04(2020)181 high-energy limit (Q m t , M H ) in terms of the top pole mass m t are given by 16 where G 1,2 (ŝ,t) denote explicit and lengthy functions of the kinematical variablesŝ andt that do not depend on the top mass [66]. The logarithms L ts , L 1ts are defined as Transforming the top pole mass m t into the MS mass m t (µ t ), we arrive at the LO expressions for F 1/2,LO with m t replaced by m t (µ t ) and the appropriately transformed NLO coefficients To minimize the logarithms of µ t , a dynamical scale of the order of √ŝ = Q has to be chosen, but not the top mass. A coefficient κ in front of the dynamical scale choice µ t = κQ is still arbitrary (but should not be large) since additional finite parts of the functions G 1,2 (ŝ,t) may be absorbed in the scale choice. Thus, the dynamical scale Q can be identified as the preferred central scale choice of the Yukawa couplings for large Q values.
The uncertainties originating from the scheme and scale dependence of the top mass can be reduced by calculating the NNLO mass effects. Such a three-loop calculation is beyond everything that has been performed so far with current methods, but for Q values close to threshold a large-mass expansion at NNLO could be used to reach an approximate estimate of the finite top-mass effects at NNLO. As a first step, partial results of the NNLO top-mass effects are known in the soft+virtual approximation [63]. For Q values around the JHEP04(2020)181 virtual tt threshold Q ∼ 2m t , non-relativistic Green's functions could be used that allow the introduction of higher-order corrections to the QCD potential [107][108][109][110][111]. This may lead to an improved description of the threshold region. However, for the triangle diagrams, the threshold behaviour is determined by P -wave contributions, since the tt-ground state appears as a CP-odd configuration that does not mix with the virtual CP-even threshold state of the triangle diagrams. For the box diagrams, the dominant S-wave contributions have to be considered. Moreover, it is unclear how large the impact of top-mass effects of the remainder beyond the non-relativistic Green's functions will be. Finally, for the high-energy tail, the approximate calculation of ref. [66] could be extended to NNLO.

Variation of the cross section with λ H 3
Higgs-pair production at the LHC is directly sensitive to the trilinear Higgs coupling. The dependence of the total and differential cross sections on the trilinear coupling λ H 3 is modified by the NLO QCD corrections and in particular by the finite mass effects at LO and NLO. Finite top-mass effects result in a non-vanishing matrix element at threshold, while in the HTL the matrix element of eq. (2.1) vanishes exactly [10,11,112], where we have used that the second form factor G vanishes in the HTL [see eq. (2.4)]. The cancellation is induced by the destructive interference between the triangle and box diagrams at LO. This property is modified by finite subleading O(1/m 2 t ) terms but explains why the matrix element itself is suppressed at the production threshold. As a function of λ H 3 , the cross section develops a minimum at λ H 3 -values around 2.4 times the SM-value in the Born-improved HTL [12,14] since the phase-space integration adds contributions from above the production threshold. The NLO QCD corrections will shift the minimum of the cross section as a function of λ H 3 and finite top-mass effects play a prominent role in the amount of these cancellations. For the determination of the trilinear coupling, the variation of the cross section with λ H 3 is of interest. As mentioned in the introduction, the total cross section behaves approximately as ∆σ/σ ∼ −∆λ H 3 /λ H 3 for λ H 3 close to the SM value.
In the following, we will analyze the NLO results, where only the trilinear coupling has been varied. In general, however, several coupling modifications contribute to the Higgspair production cross section. This could be treated consistently by extending the SM Lagrangian by all contributing dimension-6 operators as has been studied in ref. [113] in the HTL at NLO and in ref. [73] at NNLO. Recently the HTL analysis has been extended to the inclusion of finite top-mass effects at NLO [114]. However, we will neglect all dimension-6 operators but the one modifying the Higgs self-interactions. A proper and consistent effective model of this type has been discussed in ref. [16] that adds higherdimension operators to the scalar Higgs sector only. Thus, a sole variation of the Higgs self-interactions could be realized within Higgs portal models with additional heavy scalar states that couple only to the SM-like Higgs field and are integrated out.    account when determining the value of λ H 3 from the experimental data at the HL-LHC. This agrees with the findings of ref. [114]. The NLO mass effects on the variation of the total cross section with λ H 3 become larger with rising c.m. energy of the hadron collider.
In figure 18, we display the consistently defined K-factors K = σ NLO /σ LO as a function of λ H 3 in units of the SM coupling. The full curves show the NLO K-factors including the NLO top-mass effects for various c.m. energies. The dotted curves exhibit the corresponding K-factors in the Born-improved HTL as computed in refs. [12,113]. The impact of the NLO mass effects on the K-factors ranges at the level of 10-15% for negative λ H 3 values, where the triangle and box diagrams interfere constructively. For positive values of λ H 3 (destructive interference), the size and sign of the NLO mass effects is changing considerably as can be inferred from the comparison to the dotted curves. The full Kfactors develop a larger dependence on λ H 3 than the Born-improved HTL due to the NLO top-mass effects. This confirms the findings of ref. [114]. The NLO top-mass effects of the total cross section increase with rising collider energy in general except for the regions of destructive interference between the triangle and box diagrams (positive λ H 3 ).
The full NLO cross section as a function of λ H 3 can be parametrized as It should be noted that the final numerical errors of the cross sections as shown in figures 16-18 are smaller than the ones emerging from using the coefficients of eqs. (4.14), (4.15) since the combinations of each bin in Q before integration reduces them.

Conclusions
In this work, we have discussed the full QCD corrections to Higgs-pair production at NLO. We have explained the details of our numerical approach to solve the multi-scale two-loop integrals involving ultraviolet and infrared singularities. The ultraviolet singularities could be extracted from the finite parts by suitable end-point subtractions, while the infrared singularities have been isolated by means of dedicated subtraction terms. The ultraviolet singularities have been absorbed by the proper renormalization of the strong coupling and the top mass, while the infrared ones cancel against the one-loop real corrections involving an additional gluon or quark in the final state of the Higgs-boson pair. We have performed the evaluation of the virtual corrections diagram by diagram without tensor reduction.
The emerging integrals develop thresholds if the virtual tt-threshold is crossed, but also at small virtualities due to the presence of purely gluonic intermediate states. The numerical stabilization of the virtual two-loop integrals has been achieved through integrations by parts of the integrands such that the power of the threshold-singular denominators is reduced. The narrow-width limit of the virtual top quarks has been obtained by a Richardson extrapolation of the results for different sizes of an auxiliarly introduced width parameter. This has allowed a numerical integration of the virtual two-loop corrections with an accuracy of less than one per cent.
The matrix elements for the real corrections have been generated with FeynArts and FormCalc and integrated using the library Collier. The collinear region of the phase-space integration has been regularized numerically by a technical cut.
We have subtracted the Born-improved HTL from the virtual and real corrections individually so that we have been left with the pure NLO top-mass effects beyond the Born-improved HTL that is implemented in the public tool Hpair. Thus, the final NLO results have been obtained by adding back the numbers from Hpair.
The final results have been analyzed in detail for the differential cross section in the invariant Higgs-pair mass and the total cross section. Finite top-mass effects beyond the Born-improved HTL decrease the total cross section by about 15% at the LHC. However, the negative mass effects are larger for the differential cross section reaching a level of −30% or −40% for large invariant Higgs-pair masses. This implies that the inclusion of the NLO top-mass effects is crucial for a reliable analysis at the LHC and future proton colliders. We have discussed the usual renormalization and factorization scale uncertainties that are in agreement with previous calculations. However, we have identified an additional scale and scheme uncertainty due to the virtual top mass. This uncertainty reaches a level of 15% for the total cross section but can be larger (up to 35%) for the differential cross section. Based on the heavy-top and high-energy expansions, we have discussed the preferred scale choices of the running top mass and identified a large dynamical scale as the proper choice for large invariant Higgs-pair masses. This additional uncertainty has to be combined with the usual renormalization and factorization scale uncertainties. Since the (relative) scheme and scale uncertainties originating from the top mass only mildly depend on the renormalization and factorization scale choice, the addition of this uncertainty may lead to JHEP04(2020)181 about a linear addition to the other uncertainties, if the total uncertainty is defined as the envelope. This, however, has to be analyzed in more detail which is left for future work.
We have investigated the total cross section as a function of the trilinear coupling varied from its SM value. We have found significant NLO mass effects beyond the Born-improved HTL that result in a shift of the minimum of the cross section at various present and future c.m. energies of the hadron colliders. While the main effect of shifting the minimum originates from the NLO top-mass effects of the real corrections, the more symmetric virtual mass effects mainly affect the size of the total cross section as a function of λ H 3 . The full K-factors develop a larger dependence on λ H 3 than those of the Born-improved HTL due to the NLO top-mass effects.