Gluon fusion into Higgs pairs at NLO QCD and the top mass scheme

We present the calculation of the full next-to-leading order (NLO) QCD corrections to Higgs boson pair production via gluon fusion at the LHC, including the exact top-mass dependence in the two-loop virtual and one-loop real corrections. This is the first independent cross-check of the NLO QCD corrections presented in the literature before. Our calculation relies on numerical integrations of Feynman integrals, stabilised with integration-by-parts and a Richardson extrapolation to the narrow width approximation. We present results for the total cross section as well as for the invariant Higgs-pair-mass distribution at the LHC, including for the first time a study of the uncertainty due to the scheme and scale choice for the top mass in the loops.


Introduction
Since the discovery of a Higgs boson at the LHC [1,2], an enormous amount of data has been collected to analyse the properties of this newly discovered particle. Up to now the Standard Model (SM) Higgs boson hypothesis [3][4][5][6] is the most favoured one, even if there is still some room for new physics effects in the Higgs-coupling measurements. The experimental determination of the Higgs boson self-couplings, one of the most important measurements in the Higgs sector and a major goal of the high-luminosity upgrade of the LHC and of future high-energy colliders, is still yet to be a e-mail: julien.baglio@uni-tuebingen.de b e-mail: francisco.campanario@ific.uv.es c e-mail: seraina.glaus@kit.edu d e-mail: milada.muehlleitner@kit.edu e e-mail: michael.spira@psi.ch f e-mail: juraj.streicher@uni-tuebingen.de performed and would give access to the scalar potential itself which is at the origin of the electroweak symmetry breaking mechanism. In order to obtain the Higgs self-couplings and in particular the triple Higgs coupling λ H 3 , the pair production of Higgs bosons needs to be considered [7,8]. For the energies reachable at the LHC the measurement of the quartic Higgs self-coupling will be hopeless, since the triple Higgsproduction cross section is too small [9][10][11].
The dominant production channel for Higgs pair production at hadron colliders is gluon fusion. It is mediated by heavy quark-loops at leading order (LO), in two different topologies: triangle diagrams containing the triple Higgs coupling, and box diagrams as an irreducible "background" for the extraction of λ H 3 [12][13][14][15]. The NLO QCD corrections have been calculated in the heavy top-quark limit (HTL) some time ago [7], a framework in which the leading term of a heavy top mass expansion is obtained so that the NLO calculation reduces to one-loop corrections to effective H Hg and H Hgg couplings. The K -factor is found to be as sizeable as the corresponding K -factor for single Higgs production, K 2. The next-to-next-to-leading order (NNLO) QCD corrections in the same HTL approximation have been computed in Refs. [16][17][18] for the total cross section and found to be of the order of +20%, and in Ref. [19] for the differential distributions. Soft gluon resummation at nextto-next-to-leading logarithmic (NNLL) accuracy has been performed [20,21], and the 3-loop matching coefficient has been derived in Refs. [18,22].
There has been an enormous effort in the past few years to reach the full NLO QCD accuracy in Higgs pair production via gluon fusion including finite top mass m t effects. After the calculation of m t -effects in the real radiation [23,24] leading to a −10% reduction of the cross section, a 1/m t expansion up to the order O(1/m 12 t ) lead to an estimate of the mass effects at the order of ±10% for the total cross section at NLO QCD [25,26]. The first calculation of the full NLO QCD corrections has been performed in [27,28] applying numerical methods based on sector decomposition and contour deformation to master integrals, and has shown that the m t -effects in the virtual corrections are of the order of ∼ −5% for the total cross section, but can amount to ∼ −25% in the tail of the invariant Higgs-pair-mass (m H H ) distribution. This fixed-order calculation has been matched to parton shower programs later [29,30] and combined with the NNLO QCD corrections in the HTL in Ref. [31]. Up to now no other independent calculation of the full NLO QCD corrections has been available, but several approximations have been able to partially reproduce the mass effects in the virtual corrections [32,33]. Analytic results in the high-energy limit are also available in Ref. [34]. This letter presents the first independent calculation of the top-mass effects at NLO QCD for Higgs boson pair production via gluon fusion at the LHC since the original work of Refs. [27,28]. Our method is based on the direct numerical integration of the Feynman integrals of the full diagrams, using integration-by-parts to improve the stability beyond the virtual thresholds, defined by intermediate gluon pairs and top-quark pairs, and a Richardson extrapolation [35] to obtain the final numbers in the narrow-width approximation of the virtual top quarks. We will present our results for the total cross section and for the invariant Higgs-pair-mass distribution at a center-of-mass (cm) energy of 14 TeV. We will perform for the first time an analysis of the uncertainties due to the scheme and scale chosen for the top mass.
The paper is organised as follows. In Sect. 2 we will describe the technical details of our calculation. In Sect. 3 we will present our numerical results. We present the renormalisation and factorisation scale uncertainties in Sect. 4 and derive the uncertainty due to the top mass in Sect. 5. In Sect. 6 we will close with our conclusions.

Partonic leading order cross section
At LO, the gluon fusion process is mediated by heavy quark loops. There are triangle diagrams involving the triple Higgs coupling and box diagrams with two Yukawa couplings. As the Yukawa coupling is proportional to the mass of the quarks in the loop, we only consider the top-quark contribution. Following the conventions of Ref. [7], the cross section can be decomposed into form factors after applying two tensor projectors on the matrix elements, leading to the following LO expression for the partonic cross sectionσ (gg → H H),σ where G F = 1.1663787 × 10 −5 GeV −2 is the Fermi constant, α s (μ 2 R ) is the strong coupling constant evaluated at the renormalisation scale μ R , and the Mandelstam variablesŝ andt are given bŷ with the scattering angle θ in the partonic c.m. system and where m H is the Higgs boson mass. The integrations limits are given bŷ The GeV being the vacuum expectation value of the Higgs field, and the form factors reduce to F = −F = 2/3 and G = 0 in the HTL approximation. The full m t -dependence can be found in Refs. [13,15].

Hadronic cross section
The NLO QCD corrections include the two-loop virtual corrections to the triangle and box diagrams, the virtual oneparticle-reducible double-triangle diagrams, and the oneloop 2 → 3 real corrections. All these contributions are convolved with the parton distributions functions (PDFs) f i evaluated at the factorisation scale μ F , that are included in the parton luminosities dL i j /dτ , with τ = Q 2 /s, s being the hadronic cm energy. The hadronic cross section can be cast into the form with [7] for i j = gg, q,q qg, and q qq, z = Q 2 /τ s, and τ 0 = 4m 2 H /s. We include five external massless quark flavours. The coefficients C virt of the virtual and C i j of the real corrections in the HTL have been obtained in Ref. [7] and are given by where C ∞ denotes the contribution of the one-particle reducible diagrams (see Fig. 1b) in the HTL with the transverse momentum and P gq (z) are the corresponding Altarelli-Parisi splitting functions [36], given by with N F = 5 in our calculation. For the LO cross sectionσ LO (Q 2 ) the full quark-mass dependence is taken into account at the integrand-level. These expressions can easily be converted into the differential cross section with respect to the invariant Higgs-pair mass,

Virtual corrections
The two-loop virtual corrections split into three different types of contributions: triangle diagrams with a virtual Higgs propagator in the s-channel, box diagrams, and one-particle reducible diagrams. Typical examples of these three classes are shown in Fig. 1. According to Eq. (9) the two-loop triangle diagrams coincide with the production of a single Higgs boson with virtuality Q 2 so that they can be taken from this NLO calculation [37][38][39][40][41]. On the other hand the one-particle reducible diagrams involving two triangle loops of the Higgs coupling to one on-shell and one off-shell gluon can be constructed from the LO (one-loop) calculation of the Higgs decay H → Z γ [42,43] where the couplings and colour factors have to be adjusted appropriately and the Z mass has to be replaced by the virtuality of the gluon in the t/u-channel. This results in the one-particle reducible contribution to the virtual corrections with τ = 4m 2 t /m 2 H , λt = 4m 2 t /t and the functions This expression has to be inserted in the C coefficient of Eq. (7) and agrees with the previous analytical result of Ref. [44].
The cumbersome part of this calculation is the computation of the two-loop box diagrams. We have performed a Feynman parametrisation of the full virtual diagrams individually without any tensor reduction to master integrals. The ultraviolet singularities are extracted by appropriate end- Typical two-loop diagrams contributing to Higgs-boson pair production via gluon fusion: a triangle diagram, b one-particle reducible diagram, c box diagram point subtractions. Special care, however, is required for the infrared divergent diagrams of the type of e.g. Fig. 1c that correspond to gluon-rescattering involving a threshold starting at Q 2 = 0, i.e. in the whole kinematical range of the process. The Feynman parametrisation can be adjusted such that the (singular) denominator is expressed as a second-order polynomial in terms of one of the Feynman parameters (here x 6 , while x 1 , . . . , x 5 correspond to the other parameters), where H (x, x 6 ) denotes the numerator related to the full spinorial structure of the diagram and where the coefficient c can be expressed as a pure product of a kinematical ratio of O(1/m 2 t ) and other Feynman parameters. Based on the fact that the relative infrared singularities are universal we constructed a subtraction term by keeping only b and c in the denominator N 0 (x, x 6 ) = bx 6 + c thus rewriting the integral of Eq. (12) as I = I 1 + I 2 , The integral I 1 can now be expanded in leading to numerically finite integrals for each expansion coefficient, while the integral I 2 can be integrated over x 6 yielding hypergeometric functions. The transformation properties of the latter (related to an inversion of the kinematical argument) allow for a clean separation of the infrared singularities. The remaining singularities in the coefficient c can be treated by end-point subtractions. The cumbersome treatment of these diagrams can be related to the appearance of two different scales, m t and Q, that control the high-scale and low-scale parts of the whole calculation in the sense of an effective theory, the HTL, where the high scale is given by the top mass and the low scale by Q. This method of constructing the infrared subtraction term has been developed within the old NLO calculation of single Higgs production [37,38] for internal numerical checks.
Since the integration overt is not finite for the individual two-loop diagrams we have introduced a cut at the integration bounds as well as a suitable substitution to stabilise the integration close to the bounds. By varying the cut around its central value we have checked that the total sum of all twoloop box diagrams is finite and independent of the particular value of this cut. For the analytical continuation of our virtual amplitudes above the thresholds we introduced a small imaginary part of the virtual top mass, m 2 t → m 2 t (1−i¯ ) in our numerical integration. However, to stabilise the numerical integration above the thresholds we had to perform integration by parts in one of the involved Feynman parameters in order to reduce the power of the corresponding denominator for each diagram individually. In this way we achieved a reliable accuracy of our numerical integrations for¯ values down to about 0.1. In order to obtain the result in the narrow-width approximation (¯ → 0) we performed a Richardson extrapolation [35] applied to the results for different values of¯ . 1 In particular, we use¯ in the range given by¯ n = 0.05 × 2 n , with n = 0, . . . , 9. In the dominant region, we use the set n = 1, . . . , 9, with the exception of the bins in the range Q ∈ [300 − 475] GeV where the complete set of values is used. Starting at Q = 700 GeV, we restrict ourselves to values in the range n = 1, . . . , 5. This allows us to obtain a series of extrapolated results up to the ninth order in the dominant region and up to the fifth order in the tails. We define a theoretical error estimate due to the Richardson extrapolation as the difference of the fifth and the fourth order extrapolated results. Moreover, this error is multiplied by a factor of two close to the top threshold, in order to be conservative. The obtained estimated Richardson extrapolated error falls below the percent level of accuracy and is added in quadrature to the statistical integrated error.
The top mass has been renormalised in the on-shell scheme and the strong coupling constant in the MS scheme with 5 active flavours, i.e. with a decoupled top quark. We have achieved a finite result for the virtual corrections by subtracting the virtual correction in the HTL so that our numerical integration yields the NLO mass effects only. The virtual corrections in the HTL have then been added back by the results of HPAIR. 2 The calculation of each two-loop box diagram has been performed independently at least twice with different Feynman parametrisations.

Real corrections
We are left with the calculation of the coefficients C i j , or more specifically the calculation of the finite m t -effects in the real corrections, Δσ mass involves the full mass-dependent LO contribution as exemplified in Eqs. (6,9). The HTL expressions for the coefficients C i j are given in Eq. (7) [7] and can be calculated using the program HPAIR.
To obtain the mass effects, we use the fact that the infrared divergences in the real corrections are universal and are the same in the full calculation and in the HTL approximation. At any given phase-space point we can subtract the HTL result from the full calculation, obtaining an infrared-finite result which is exactly the remainder due to the mass effects in the full real corrections, In order to calculate the LO matrix elements we need to map the full 2 → 3 phase-space onto an on-shell LO 2 → 2 subspace, leading to a transformation of the 4-momenta p i → p i . For this purpose, we use the mapping of Ref. [45] for the case of initial-state emitter and initial-state spectator to build our local subtraction term. The HTL matrix elements are calculated analytically. In order to achieve a good numerical stability we have implemented a technical collinear cut in the scattering angle integration of the extra parton in the final state, | cos θ | < 1 − δ, with δ = 10 −4 . It has been checked that this value does not affect the physical results by varying δ around our nominal choice. The full one-loop matrix elements contain triangle, box, and pentagon diagrams. They are generated using Feyn-Arts [46] and FormCalc [47], while the scalar oneloop integrals are numerically calculated using the library 2 The program can be downloaded at http://tiger.web.psi.ch/hpair/. COLLIER 1.2 [48] interfaced to the output of FormCalc with an in-house library. We have performed two independent calculations with different parametrisations of the 2 → 3 phase-space and different versions of FeynArts and FormCalc. We have also cross-checked the mass effects of the real corrections against Refs. [23,24,27,28] and we have obtained mutual agreement.

Numerical results
We present our numerical results at the LHC for a c.m. energy of √ s = 14 TeV. We use m H = 125 GeV and m t = 172.5 GeV. We have performed the calculation using two NLO PDF sets, MMHT2014 [49] and PDF4LHC15 [50] as implemented in the LHAPDF-6 library [51]. Our central scale choice is μ R = μ F = μ 0 = m H H /2, and α s (M 2 Z ) is set according to the PDF set chosen, with a NLO running in the five-flavour scheme. We recall that Q and m H H refer to the same physical quantity. We have performed the whole calculation twice independently using also different parametrisations of the virtual Feynman integrals and of the real phase-space. Both independent calculations agree with each other within the numerical errors. Note that our calculation has been performed in the narrow-width approximation for the top quark. The finite-width effect amounts to ∼ −2% at LO for the total cross section, corresponding to a maximum deviation of ∼ −4% at the tt threshold in the invariant Higgs-pair-mass distribution [24].
We have calculated a grid of Q-values from Q = 250 GeV to Q = 1500 GeV, so that we obtain the invariant Higgspair-mass distribution depicted in Fig. 2. We compare our full results, displayed in red, to three different approximations: The HTL results according to Ref. [7] are displayed in blue; the HTL plus the mass effects in the real corrections only are shown in yellow; the HTL plus the mass effects in the virtual corrections only are presented in green. The error bars due to the total numerical errors are also given for the green and red histograms. They are determined from the individual integration errors and the errors due to the Richardson extrapolation -added in quadrature. The numerical errors are negligible for the other predictions (HTL and HTL plus the mass effects in the real corrections only). The red band indicates the renormalisation and factorisation scale uncertainties of our prediction for the full NLO QCD predictions, see next section for more details. In the case of the MMHT2014 PDF set we also compare to the LO results. We reproduce the results of Refs. [23,24] for the mass effects in the real corrections, which are mildly varying from m H H = 400 GeV to m H H = 1500 GeV, and are of the order of −10%. The global K -factor displayed on the left-hand-side of Fig. 2 [27,28]. The mass effects are negative in accordance with the expected restoration of partial-wave unitarity in the high-energy limit.
From the differential distribution we calculate the total cross section using a numerical integration over Q with the trapezoidal method supplemented by a Richardson extrapolation [35] for Q > 300 GeV, while for Q < 300 GeV we have used the extension of Boole's rule to six nodes [52]. We obtain the following results for the NLO cross sections (the numbers in parenthesis indicate the numerical errors, i.e. the quadratic sum of the statistical and the Richardson extrapolation error), Comparing Eqs. (16) and (17), we derive −15% top-mass effects at NLO on the total cross section. We are in mutual agreement with the results of Refs. [27,28] within the numerical uncertainties and taking into account the small difference in the choice of the top mass value.

Factorisation and renormalisation scale dependence
The last result is in mutual agreement with Refs. [27,28] within the numerical uncertainties and taking into account the small difference in the choice of the top mass value.
using PDF4LHC parton densities. The top-quark scheme uncertainty is significant over the whole range of m H H . Note that a similar result has been observed in single Higgs production for large Higgs masses which correspond to our triangle diagram involving the triple Higgs coupling. Furthermore, this scheme uncertainty is reduced by roughly a factor of two from LO to NLO. The prediction involving the top pole mass, that we take as our central prediction, is the maximal prediction for high m H H values. The uncertainties induced by the top-mass scheme and scale choice on the total cross section at NLO will be given in a forthcoming publication [53].

Conclusions
We have presented the calculation of the full NLO QCD corrections to Higgs-boson pair production via gluon fusion for the top-loop contributions. This has been performed by numerical integrations of the involved virtual two-loop corrections to the four-point functions, while the results of the single-Higgs case have been translated to the three-point contributions that involve the trilinear Higgs self-coupling. The one-particle reducible contributions that appear for the first time at NLO have been inferred from the explicit analytical one-loop results for H → Z γ , where the Z -boson mass plays the role of the virtuality of the gluon in the dressed Hgg * vertex. In order to isolate the ultraviolet, infrared and collinear divergences, we have performed appropriate endpoint subtractions at the integrand level and described the explicit construction of infrared subtraction terms that allow for a clean separation of the infrared singularities from the regular rest. The real corrections have been obtained by generating the full matrix elements with automatic tools. We have constructed the infrared and collinear subtraction term as the heavy-top limit of the real matrix elements involving the fully massive LO sub-matrix element. Adding back the full results in the heavy-top limit completed the full real corrections. The final results we have obtained agree with previous calculations for the individual finite parts of the real and virtual corrections. We find finite NLO mass effects that are up to −30% for large invariant Higgs-pair masses, while the total NLO top-mass effects modify the total cross section by about −15%.
We have studied the theoretical uncertainties related to variations of the renormalisation and factorisation scales and have found agreement with the previously known results finding uncertainties at the level of 10 − 15%. A novel outcome of our calculation is the additional uncertainty induced by the scheme and scale dependence of the top mass that can be significant, amounting to +6%/ − 34% at m H H = 300 GeV and +0%/ − 35% at m H H = 1200 GeV. The induced uncertainty on the total cross section will be given in a forthcoming publication [53].
In the future we plan to extend our calculation to beyondthe-SM models as e.g. the 2HDM or MSSM.