Higgs boson pair production at NNLO with top quark mass effects

We consider QCD radiative corrections to Higgs boson pair production through gluon fusion in proton collisions. We combine the exact next-to-leading order (NLO) contribution, which features two-loop virtual amplitudes with the full dependence on the top quark mass $M_t$, with the next-to-next-to-leading order (NNLO) corrections computed in the large-$M_t$ approximation. The latter are improved with different reweighting techniques in order to account for finite-$M_t$ effects beyond NLO. Our reference NNLO result is obtained by combining one-loop double-real corrections with full $M_t$ dependence with suitably reweighted real--virtual and double-virtual contributions evaluated in the large-$M_t$ approximation. We present predictions for inclusive cross sections in $pp$ collisions at $\sqrt{s}$=13, 14, 27 and 100TeV and we discuss their uncertainties due to missing $M_t$ effects. Our approximated NNLO corrections increase the NLO result by an amount ranging from +12% at $\sqrt{s}$=13TeV to +7% at $\sqrt{s}$=100TeV, and the residual uncertainty from missing $M_t$ effects is estimated to be at the few percent level. Our calculation is fully differential in the Higgs boson pair and the associated jet activity: we also present predictions for various differential distributions at $\sqrt{s}$=14 and 100TeV. Our results represent the most advanced perturbative prediction available to date for this process.


Introduction
One of the primary goals of the LHC programme in the next decades is the detailed study of Higgs boson properties. In particular, the high-luminosity upgrade of the LHC is expected to provide direct constraints on the Higgs boson trilinear coupling from Higgs boson pair production [1, 2], which may reveal whether the Higgs potential is indeed Standard Model-like. A detailed theoretical understanding of Higgs boson pair production processes is thus mandatory.
Considering the magnitude of the total Higgs boson pair production cross sections at √ s = 14 TeV [3,4], the most promising process to constrain the Higgs trilinear coupling is pair production via gluon fusion. Due to the smallness of the corresponding production cross sections, it has been recently suggested to additionally harness complementary information on the trilinear Higgs coupling from higher-order contributions to single Higgs boson production [5][6][7][8][9][10] or electroweak precision observables [11,12].
For the gg → hh production channel, the leading order (LO) calculation was performed some time ago in Refs. [13][14][15]. The next-to-leading-order (NLO) corrections with full top quark mass (M t ) dependence, involving two-loop diagrams with several mass scales, became available only recently [16,17], and have been supplemented by soft-gluon resummation at small transverse momenta of the Higgs boson pair [18] and parton shower effects [19,20].
In the M t → ∞ limit, also called Higgs Effective Field Theory (HEFT) approximation, point-like effective couplings of gluons to Higgs bosons arise. In this limit, the NLO corrections were first calculated in Ref. [21] and rescaled by a factor B FT /B HEFT , where B FT denotes the LO one-loop matrix element squared in the full theory. This procedure is often called "Bornimproved HEFT" approximation.
In Refs. [4,22] an approximation for Higgs boson pair production at NLO, labelled "FT approx ", was introduced, in which the real radiation matrix elements contain the full top quark mass dependence, while the virtual part is calculated at NLO in the HEFT approximation and rescaled at the event level by the re-weighting factor B FT /B HEFT . At the inlusive cross section level this approximation suggests at the LHC a correction with respect to the "Born-improved HEFT" approximation of about −10%, close to the corresponding correction of −14% later obtained in the full NLO calculation [16,17].
The next-to-next-to-leading-order (NNLO) QCD corrections in the HEFT approximation have been computed in Refs. [23][24][25][26], where Ref. [26] provides fully differential results. The NNLO HEFT results for the total cross section have been supplemented by an expansion in 1/M 2 t in Ref. [27]. Approximations for the top-quark mass dependence of the two-loop amplitudes in the NLO calculation have been studied in Ref. [28] via a Padé ansatz. Soft gluon resummation has been performed at NLO+NNLL in Ref. [29] and at NNLO+NNLL in Ref. [30]. The NNLO+NNLL HEFT results lead to K-factors of about 1.2 relative to the Born-improved NLO HEFT result.
the NNLO+NNLL HEFT results [30], corrected by a factor δ t accounting for top quark mass effects, extracted from Ref. [16]. However, this procedure is somewhat ad hoc, and not viable to study kinematical distributions. In order to account for the NNLO K-factor in the HEFT calculation as well as for the correct description of the tt threshold and the high-energy tails of the distributions, where the top quark loops are resolved, a first attempt to combine the two calculations has been made in Ref. [17], where the full NLO result for a particular distribution was reweighted by the NNLO K-factor obtained from Ref. [26] on a bin-by-bin basis. However, this procedure, called "NLO-improved NNLO" has its drawbacks, as it needs to be repeated for each observable (and binning) under consideration.
The aim of this paper is to study alternative methods to combine the two results, i.e. to incorporate top quark mass effects in the calculation of the production of Higgs boson pairs at NNLO. One of the studied approximations comprises exact top-quark mass dependence up to NLO and also exact top quark mass dependence in the double-real emission contributions to the NNLO cross section at differential level. The results of this approximation can be regarded as the most advanced prediction currently available for Higgs boson pair production in gluon fusion.
This work is organized as follows: In Section 2 we describe the technical details of our calculation, and present the different approximations we will consider to incorporate mass effects in the NNLO contribution. In Section 3 we present our numerical predictions, both for the total cross section and differential distributions. Finally, in Section 4 we summarise our results.

Details on the method and approximations
We start by presenting the different technical ingredients entering our computation, as well as the definition of the various approximate ways to include mass effects in the NNLO calculation introduced and used in this work. Finally, we also discuss the numerical stability of our predictions.

Technical ingredients
Our calculation is based on the publicly available computational framework Matrix [32], which allows the user to perform fully differential NNLO calculations for a wide class of processes at hadron colliders. For the purpose of the present work, the public version of the code has been extended, based on the calculation of Ref. [26], to include the production of a pair of Higgs bosons via gluon fusion. For the calculation of the NNLO corrections the code implements the q T -subtraction formalism [33], in which the genuine NNLO singularities, located where the transverse momentum of the Higgs boson pair, p T,hh ≡ q T , vanishes, are explicitly separated from the NLO-like singularities in the hh + jet contribution. The q T subtraction formula reads where in particular the contribution dσ hh+jet NLO can be evaluated using any available NLO subtraction procedure to handle and cancel the corresponding infrared (IR) divergencies * . The remaining q T → 0 divergence is canceled by the process-independent counterterm dσ CT NNLO . The process-dependence of the hard-collinear coefficient H hh NNLO enters only via the NNLO (HEFT) two-loop virtual corrections [23] through an appropriate subtraction procedure [37].
The difference in the square bracket of Eq. (1) is finite when q T → 0, but each of the terms exhibits a logarithmic divergence. Therefore, a technical cut, r cut , needs to be introduced on q T /Q, where the scale Q is chosen to be the invariant mass of the final-state system. More details about the r cut → 0 extrapolation are provided in Section 2.3.
At variance with the calculation of Ref. [26], which was strictly done within the HEFT, this time all the routines needed to compute the full NLO cross section as well as the different NNLO reweightings have been implemented. This includes linking the code to the NLO two-loop virtual corrections obtained via a grid interpolation [19] and to several loop-induced amplitudes provided by the OpenLoops amplitude generator [38]. Within this framework we reproduced the differential NLO results of Refs. [16,17] at the per mille level.
The grid for the NLO virtual two-loop amplitudes is based on the calculation presented in Refs. [16,17], which in turn for the calculation of the two-loop amplitudes relies on an extension of the program GoSam [39,40] to two loops [41], using also Reduze2 [42], SecDec3 [43] and the Quasi-Monte Carlo technique as described in Ref. [44] for the numerical integration. These amplitudes (for fixed values of the Higgs boson and top quark masses) are provided in a two-dimensional grid together with an interpolation framework, which allows us to evaluate them at any phase space point without having to perform the computationally costly two-loop integration. For more details, see Refs. [19,45].
All tree and one-loop amplitudes in the HEFT and also all loop-squared amplitudes in the full theory as discussed below are obtained via a process independent interface to Open-Loops [38,46,47]. For the latter this comprises loop-squared amplitudes for pp → hh + 1, 2 jets, that need to be evaluated in IR divergent unresolved limits. In particular the limit q T → 0 represents a significant challenge for the numerical stability of the hh + 2 jets amplitudes in the full theory. Thanks to the employed algorithms the numerical stability is under control, as discussed in detail in Section 2.3. A major element of this stability originates from the employed tensor integral reduction library COLLIER [48].

Approximations for top-mass effects at NNLO
In the following we present three approximations for the NNLO Higgs boson pair production cross section, which take into account finite top quark mass effects in different ways. In all cases, we always include the full NLO result when computing the NNLO prediction, and only apply the different approximations to the O(α 4 S ) contribution.
NNLO NLO-i The NLO-improved NNLO approximation (NNLO NLO-i ) has already been presented in Ref. [17]. It can be constructed based on an observable-level multiplicative approach. In this approximation, for each bin of each histogram we multiply the full NLO result by the ratio between the HEFT NNLO and NLO predictions for this bin.
NNLO B-proj A different approximation can be obtained by reweighting each NNLO event by the ratio of the full and HEFT Born squared amplitudes. We denote this procedure as Bornprojected approximation (NNLO B-proj ). Of course, in order to do so and due to the different multiplicities involved, an appropriate projection to Born-like kinematics is needed; for this purpose we make use of the q T -recoil procedure defined in Ref. [49]. Following this prescription, the momenta of the Higgs bosons remain unchanged, and the new initial-state parton momenta are obtained by absorbing the recoil due to the additional radiation. Specifically, denoting the momenta of the incoming partons by p 1 and p 2 , and the momentum of the Higgs boson pair system by q, the new momentum to be used for the LO projection k 1 (then, k 2 = q − k 1 ) is given by where and k µ 1T is a two-dimensional vector in the q T plane which needs to fulfill the condition k 1T → 0 when q T → 0, and we set k 1T = q T /2 (and therefore k 2T = q T /2). This condition guarantees that the subsequently applied reweighting does not spoil the NNLO q T -cancellation. More details about this procedure can be found in Ref. [49].
NNLO FTapprox The third approximation we consider is constructed to profit from the fact that the double-real emission contributions to the NNLO cross section require only one-loop amplitudes in the full theory (FT) and can thus be computed by using OpenLoops. Of course, the inclusion of these loop-induced amplitudes needs to be done in such a way that the dipole cancellations in the NLO hh + j calculation and the low-q T cancellation for hh at NNLO are not spoiled.
We will define our approximation by using the following procedure: working in the HEFT, for each n-loop squared amplitude that needs to be computed for a given partonic subprocess HEFT (ij → HH + X), we apply the reweighting where A Born Full stands for the lowest order (loop-induced) squared amplitude for the corresponding partonic subprocess, computed in the full theory. † We note that, contrary to what happens in the Born-projected approach, here the reweighting is defined using amplitudes that correspond to the same subprocess under consideration. Therefore, the kinematics is always preserved and there is no need to define a Born projection. Moreover, for amplitudes that are of tree-level type in the HEFT (as it is the case for the double-real emission contributions), this reweighting simply implies using the exact loop-induced amplitudes with full top mass dependence. The reweighting procedure defined by Eq. (4) agrees at NLO with the so-called FT approx introduced in Ref. [22], therefore we will use the same notation.
Given that the performance of the Born-projection and FT approximations was already studied in Ref. [17] at NLO, we directly present NNLO predictions in Section 3. We point out that, based on the ingredients entering each of the approximations, the NNLO FTapprox is expected to be the most advanced prediction for Higgs boson pair production via gluon fusion. By contrast, the NNLO B-proj is expected to be the less accurate, since it is based on a simple Born level reweighting procedure. Nevertheless, and for comparison purposes, we always present results for the three approximations described above.

Numerical stability
Before presenting our quantitative predictions, we briefly discuss the numerical stability of our results. From the computational point of view, the most challenging of the three approaches to incorporate mass effects at NNLO is the NNLO FTapprox procedure, as it involves loop-induced double-real contributions in the full theory. In particular the dominant gg → hhgg amplitude comprises computationally very challenging six-point loop integrals with internal masses. In fact, these contributions have to be evaluated in the numerically intricate NNLO unresolved limits and to the best of our knowledge, the present calculation is the first application of a six-point one-loop amplitude integrated over its IR divergent unresolved limits in an NNLO calculation.
Thanks to the numerical stability of the applied algorithms in OpenLoops together with Collier, the bulk of the phase-space points remains stable in double precision when approaching q T → 0, even close to the dipole singularity, i.e. in the NNLO double-unresolved limits. On average the runtime per phase space point for the gg → hhgg amplitude is ∼ 1 sec. In principle OpenLoops provides a rescue system, such that remaining numerically unstable phase-space points can be reevaluated in higher numerical precision based on reduction with † Strictly speaking, the reweighting is applied to the finite part of the loop amplitudes. However, at one-loop level this procedure reproduces the loop structure of the full theory. CutTools [50]. However, the runtime of the loop-induced gg → hhgg amplitude in Open-Loops is significantly increased when CutTools is used in quadruple precision (to the level of ∼ 10 minutes per phase-space point), rendering the quadruple precision stability system prohibitive for this amplitude for practical purposes ‡ . Therefore, we restrict the evaluation to double precision and replace potentially unstable phase-space points close to the dipole singularities, quantified by α L-i = (p i · p j /ŝ) min , where the minimum among all potential emitter parton combinations i and j is taken, with an approximation: Below a technical cut α L-i, cut we switch from the (loop-induced) double-real amplitude in the FT to the (tree-level) double-real amplitudes in the HEFT, reweighted at LO. This approach could in principle introduce a bias in the NLO hh+jet cross section, thereby hampering the low-q T cancellation of the NNLO computation. We have checked that this is not the case, as detailed in the following.
For the predictions presented in Section 3 we use α L-i, cut = 10 −4 and we varied this parameter in the range 10 −3 to 10 −5 , finding independence of all results. In Fig. 1 we illustrate the resulting dependence of the NNLO FTapprox total cross section on the q T -subtraction cut, r cut , for √ s = 14 TeV. Due to the previously discussed stability challenges, we considered values of r cut between 1% and 3.5%, which are larger than the ones typically used in previous q T -subtraction calculations (compared for instance with the default values in the public Matrix release [32]). Nevertheless our results present a good stability, with effects that are below 0.2% in the whole q T /Q range under study, validating this choice. The r cut → 0 extrapolation is performed using a linear least χ 2 fit. The fit is repeated varying the upper bound of the interval (in this case starting from a minimum of 25 points, which corresponds to an upper bound of r cut = 1.6%, and up to r cut = 3.5%). Then, the result with the lowest χ 2 /degrees-of-freedom value is taken as the best fit, and the rest is used to estimate the extrapolation uncertainty. In the case shown in Fig. 1 the extrapolation uncertainty for r cut → 0, indicated with the dotted lines, is ±0.14%.
A further uncertainty arises due to the numerical evaluation of the two-loop integrals with full top-quark mass dependence in the virtual corrections of the NLO contribution. The error ‡ Here we want to note that these stability issues will be strongly mitigated in the future based on the new OpenLoops on-the-fly reduction method introduced in Ref. [47].
of the numerical integration of the amplitudes is propagated to the total cross section using Monte Carlo methods, varying the amplitude level results according to the corresponding error estimates. This leads to changes of the NLO cross section at the per mille level. Furthermore, we have checked that, within this uncertainty, results based on the grid for the virtual amplitude are consistent with the ones directly obtained from the amplitude results calculated in Refs. [16,17]. We want to point out that the uncertainties can be somewhat larger in differential results, in particular in the tails of p T and invariant-mass distributions.
This discussion shows that the uncertainties due to the q T -subtraction method and the numerical evaluation of the NLO virtual contribution and grid interpolation are clearly under sufficient control.

Results
In this section we present our numerical predictions for inclusive and differential cross sections for Higgs boson pair production in pp collisions. We consider centre-of-mass energies of 13, 14, 27 and 100 TeV. For the sake of brevity, differential distributions are presented only for 14 TeV and 100 TeV. We use the values M h = 125 GeV and M t = 173 GeV for the Higgs boson and top quark masses, respectively. We do not consider bottom quark loops, whose contribution at LO is below 1%. We also neglect top quark width effects, which at LO are at the level of 2% for the total cross section [22]. We use the PDF4LHC15 sets [51][52][53][54][55][56] of parton distribution functions (PDFs), with parton densities and α S evaluated at each corresponding perturbative order (i.e., we use the (k + 1)-loop running α S at N k LO, with k = 1, 2). As renormalization and factorization scales, we use the central value µ 0 = M hh /2, and we obtain scale uncertainties via the usual 7-point scale variation.

Inclusive cross sections
In Table 1 we present results for the total cross sections at NLO and NNLO in the various approximations. At NLO we report the exact result, including the full M t dependence, and also the FT approx result. By comparing the two NLO predictions, we see that the FT approximation overestimates the exact NLO result by 4% (6%) at 14 (100) TeV. At NNLO the largest prediction is obtained in the NNLO B-proj approximation, resulting in an increase with respect to the exact NLO result of about 20% at 14 TeV. For this collider energy, the increase within the NNLO NLO-i approach (which is computed based on the M hh distribution) is smaller, being about 18%. Finally, the NNLO FTapprox prediction is the lowest one, with a 12% increase with respect to the NLO cross section at 14 TeV. For all the considered approximations and collider energies the scale uncertainties are significantly reduced when including the O(α 4 S ) NNLO corrections. This reduction is largest for the NNLO B−proj and NNLO FTapprox approximations § . For instance at 14 TeV, the total scale uncertainty is reduced from about ±13% at NLO to +2% − 5% at NNLO FTapprox , i.e. by about a factor of three. This reduction of the scale uncertainties is stronger as we increase the collider energy, being close to a factor of five at 100 TeV.
As is well known, scale uncertainties can only provide a lower limit on the true perturbative uncertainties. In particular, from Table 1 we see that the difference between the NNLO and NLO central predictions is always larger than the NNLO scale uncertainties (although within the NLO uncertainty bands). In any case, the strong reduction of scale uncertainties, together with the moderate impact of NNLO corrections, suggests a significant improvement in the perturbative convergence as we move from NLO to NNLO.
It is also worth mentioning that the three approximations have a different behaviour with √ s. For instance at 100 TeV, the increase with respect to the NLO prediction for the NNLO B-proj and NNLO NLO-i approaches is 23% and 17%, respectively, values that are close to the ones for 14 TeV (20% and 18%, respectively). By contrast, the NNLO FTapprox result increases the NLO prediction by 7% at 100 TeV, i.e. the correction is smaller by almost a factor of two than at 14 TeV (12%), which also means a larger separation with respect to the other two NNLO approximations. The smaller size of the NNLO corrections in the FT approx at higher energies is also consistent with the observed reduction of scale uncertainties.
As was mentioned already in Section 2.2, the NNLO FTapprox result is expected to be the most accurate one among the approximations studied in this work, and therefore it is considered to be our best prediction. In order to estimate the remaining uncertainty associated with finite top quark mass effects at NNLO, we start by considering the accuracy of the FT approx approximation at NLO. At 14 TeV the NLO FT approx result (see Table 1) overestimates the full NLO total cross section by only about 4%, or equivalently by about 11% of the pure O(α 3 S ) contribution. If we assume that FT approx performs analogously at one order higher, we obtain a ±11% uncertainty on the O(α 4 S ) contribution ¶ . Given that the relative weight of the O(α 4 S ) contributions to the total NNLO cross section is definitely smaller than the weight of the O(α 3 S ) contributions to the NLO cross section, we obtain a significantly smaller overall uncertainty, in this case of ±1.2%. In order to be conservative, we can increase this estimate by a factor of two. The relative difference between the FT approx and the full NLO result slightly increases with the collider energy. However, at the same time the relative size of the O(α 4 S ) correction decreases. The NNLO uncertainty obtained with this procedure ranges from ±2.3% at 13 TeV to ±3.1% at 100 TeV.
We can repeat the above procedure to estimate the uncertainty of the NNLO B−proj approximation, which displays the largest differences with respect to the NNLO FTapprox result. Similarly to what we do for FT approx , we can assign an uncertainty to the NNLO B−proj result by relying on the accuracy of the same approximation at NLO, and conservatively multiplying by a factor of two. The ensuing uncertainties range from ±14% at √ s = 13 TeV to ±36% at √ s = 100 TeV. We find that the NNLO FTapprox prediction (always evaluated at µ R = µ F = µ 0 ) is fully contained in the NNLO B−proj uncertainty band. Actually, there is a large overlap between the two approximations, which includes in all the cases the central value of the NNLO FTapprox , even when the conservative factor of two is not included. This can be regarded as a non-trivial consistency check for our procedure. We may be tempted to conclude our discussion by adopting the above procedure for the uncertainty estimate of our NNLO FTapprox result.
However, we have already pointed out that, as √ s increases, the difference between the NNLO FTapprox and the other approximations increases. In particular, the difference between the NNLO FTapprox result and our "next-to-best" NNLO prediction, NNLO NLO−i , is 5.2% at √ s = 13 TeV, and it becomes 9.2% at √ s = 100 TeV. The significant increase of this difference with the collider energy suggests us a more conservative approach. Our final estimate for the finite top quark mass uncertainty of our NNLO FTapprox result is defined as half the difference between the NNLO FTapprox and the NNLO NLO−i approximations, and is reported in Table 1 for the different values of √ s. At √ s = 13 and 14 TeV these uncertainties are ±2.6% and ±2.7%, and thus very similar to the ones obtained with the method discussed above. At √ s = 100 TeV, however, the uncertainty increases to ±4.6%, which appears to be more conservative than the ±3.1% obtained with the previous procedure.

Differential distributions
In this section we present predictions for differential Higgs boson pair production at 14 TeV and 100 TeV. We consider the following kinematical distributions: the invariant mass (M hh , Fig. 2) and rapidity (y hh , Fig. 3) of the Higgs boson pair, the transverse momenta of the Higgs boson pair and the leading jet (p T,hh and p T,j1 , Figs. 4 and 5), the transverse momenta of the ¶ We point out that in order to obtain the pure O(α 4 S ) corrections, we have subtracted the lower order contributions computed with NNLO parton distributions and strong coupling. The corresponding numbers are a few percent lower than the ones given in Table 1 for the NLO results. harder and the softer Higgs boson (p T,h1 and p T,h2 , Figs. 6 and 7), and the azimuthal separation between the two Higgs bosons (∆φ hh , Fig. 8). For the sake of clarity, we only show the scale uncertainty bands corresponding to the NLO and NNLO FTapprox predictions.
We start our discussion from the invariant-mass distribution of the Higgs boson pair, reported in Fig. 2. We observe that the NNLO B-proj and NNLO NLO-i approximations predict a similar shape, with very small corrections at threshold, an approximately constant K-factor for larger invariant masses, and only a small difference in the normalization between them, which increases in the 100 TeV case. The NNLO FTapprox , on the other hand, presents a different shape, in particular with larger corrections for lower invariant masses, a minimum in the size of the corrections close to the region where the maximum of the distribution is located, and a slow increase towards the tail. The different behavior of the NNLO FTapprox in the region close to threshold is more evident at 100 TeV, where the increase is about 30% in the first bin. Naively we could expect that if this region is dominated by soft parton(s) recoiling against the Higgs bosons, the Born projection and FT approx should provide similar results. We have investigated the origin of this difference, and we find that in the region M hh ∼ 2M h the cross section is actually dominated by events with relatively hard radiation recoiling against the Higgs boson pair (for example, at √ s = 100 TeV, the average transverse momentum of the Higgs boson pair in the first M hh bin is p T,hh ∼ 100 GeV at NLO). In this region the exact loop amplitudes behave rather differently as compared to the amplitudes evaluated in the HEFT: As the production threshold is approached, they go to zero faster than in the mass-dependent case, thus explaining the differences we find. Within the NNLO FTapprox , the corrections to the M hh spectrum range between 10% and 20% at 14 TeV. The scale uncertainty is substantially reduced in the NNLO FTapprox , and this reduction is particularly strong for large invariant masses. As observed at the inclusive level, the NNLO FTapprox corrections are smaller at 100 TeV (except only for the first bin) and the difference with respect to the other approximations is larger.
Next we move to the rapidity distribution of the Higgs boson pair, reported in Fig. 3. The NNLO results are similar for all three approximations. This is not unexpected as the shape of the rapidity distribution is mainly driven by the PDFs. Besides the obvious difference in the normalization, the largest effect in the shape of the NNLO NLO−i distribution is observed in the central region, which is particularly evident in the 100 TeV case. Again we observe a clear reduction of scale uncertainties over the whole range under study. More significant differences between the three approximations are obtained in the p T,hh distribution, reported in Fig. 4. The NNLO B-proj approximation predicts huge corrections for large transverse momentum, the result being almost an order of magnitude larger than the NLO prediction and the other approximations for p T,hh ∼ 500 GeV. This behavior is hardly surprising since already at NLO the Born-projected result deviates from the exact NLO prediction in this way [17]. In fact, given that the p T,hh distribution is not defined at LO, the NNLO B-proj corrections cannot inherit any information about the (full) lowest-order prediction for this distribution. This is of course not the case for the other two approximations, which in fact make an almost identical prediction at large p T,hh , with large corrections that can be well above 50%, and sizable uncertainties at the level of 30%-40%, reflecting the NLO-nature of this observable. At lower transverse momenta, however, the NNLO NLO−i and NNLO FTapprox deviate from each other, and the latter approaches the NNLO B−proj prediction. Once again, the different behavior of these approximations is more pronounced in the 100 TeV distribution, for which the central NNLO NLO−i curve lies outside the NNLO FTapprox uncertainty band below p T,hh ∼ 200 GeV. Of course, in order to obtain reliable results in the low-p T,hh region, the corresponding logarithmically enhanced contributions need to be properly resummed to all orders in the strong coupling constant.
The transverse momentum distribution of the leading jet p T,j1 , reported in Fig. 5, has similar features as the p T,hh distribution. Again we observe the unphysical excess predicted by the NNLO B−proj approximation, which can be understood using the same arguments as presented for the p T,hh distribution, and the agreement between NNLO B−proj and NNLO FTapprox at low p T,j1 . The difference between the NNLO NLO−i and NNLO FTapprox results is more pronounced here, with the FT approx predicting a softer spectrum for this observable, and small corrections that are almost always contained in the NLO scale uncertainty band.
The transverse-momentum distributions of the harder and the softer Higgs boson are reported in Figs. 6 and 7, respectively.
As can be expected from the p T,hh spectrum, the NNLO B-proj result for p T,h1 features very large corrections as p T,h1 increases. The effect, however, is less severe than the one observed in p T,hh because the p T,h1 observable is already well defined at LO. The NNLO NLO-i curve is overall in good agreement with the NNLO FTapprox prediction: It shows moderate corrections with respect to the NLO result which increase as p T,h1 increases, while the scale uncertainties are about ±15%. At very small p T,h1 the higher-order corrections become perturbatively unstable as the available phase space for the real radiation is severely restricted in this regime yielding large logarithms that should be resummed in order to get a reliable prediction, see also the discussion in Section 3.4 of Ref. [19]. For the transverse momentum of the softer Higgs boson, p T,h2 , the NNLO effect is rather uniform in all three approximations, especially at 14 TeV. The NNLO FTapprox predicts small corrections of order 10%, while the other two approximations show larger corrections with a similar shape. In the tail of the distribution the scale uncertainty at NNLO is larger than at NLO, most likely due to an accidentally small size of the NLO scale variation (in fact, in this region the NLO corrections almost vanish).
Finally, the distribution in the azimuthal angle between the two Higgs bosons, ∆φ hh , is shown in Fig. 8. At LO we have ∆φ hh = π, due to the back-to-back production of the two Higgs bosons at Born level. Real contributions allow ∆φ hh to be smaller than π, and again we  observe that the NNLO B-proj approximation predicts larger corrections in the region dominated by hard radiation compared to the other two results, which again are in good agreement with each other in that region, whereas they start to deviate for larger angles. For values of ∆φ hh close to π, this observable receives large corrections from soft-gluon emission, and the corresponding large logarithms should be resummed in order to get a reliable prediction.
We conclude this section by adding a few comments on the finite-M t uncertainties at NNLO for the various differential distributions. The analysis that was performed for the total cross section cannot be easily extended to differential distributions. On one hand, any accidental agreement between the FT approx and the full result at NLO in a given phase-space region would likely lead to an underestimation of the top quark mass effects; on the other hand, the regions in which the NLO corrections are very small due to cancellations between different contributions can present very large relative differences in the O(α 3 S ) contribution of the NLO FTapprox and NLO results, thus leading to artificially large uncertainties at NNLO. In addition, there are observables that are by definition reproduced in an exact way by the FT approx at NLO (in our case p T,hh , p T,j1 and ∆φ hh ), and the uncertainty estimate procedure that we defined for the inclusive case is therefore not applicable. Despite these facts, and based on the performance of the FT approx at NLO [17] as well as on the observed differences between our NNLO approximations, we can try to assess the order of magnitude of the expected missing M t effects for the distributions presented above.
In the Higgs boson pair invariant-mass distribution, for values of M hh below 500 GeV the level of accuracy of the FT approx at NLO is similar to the inclusive case, and therefore the M t uncertainty at NNLO is expected to be of a comparable size. In the tail of the distribution, however, the quality of the FT approx decreases (see Fig. 5 of Ref. [17]), and we thus expect the finite top quark mass effects to be of O(10%) in this region.
The shape of the rapidity distribution of the Higgs boson pair is correctly described by the FT approx at NLO (see Fig. 8 of Ref. [17]), and the difference to the full result is only the overall normalization. Based on this, the estimated top quark mass uncertainty for the NNLO FTapprox result is constant in the whole y hh range and of the same size as for the inclusive cross section.
The transverse momentum of the harder Higgs boson is very well described at NLO by the FT approx (see Fig. 7 of Ref. [17]), being always within the NLO scale uncertainty band. This fact, together with the close agreement between the NNLO FTapprox and NNLO NLO−i predictions, suggests that the missing top quark mass effects at NNLO are probably of moderate size. The same holds true for the transverse-momentum distribution of the softer Higgs boson, except for the tail where at NLO the FT approx overestimates the full NLO corrections, which in fact almost vanish in this region.
The remaining distributions, which are either not defined or trivial at LO, are by definition reproduced in an exact way by the FT approx at NLO, and this makes the estimate of the missing top quark mass effects at NNLO more difficult. In this case, a possible approach can be to use the difference between the NNLO FTapprox and NNLO NLO−i prediction as an estimate of the uncertainty (as discussed before, the NNLO B−proj prediction is not expected to be reliable in the regions dominated by hard real radiation, where it largely deviates from the other two approximations). This procedure would imply relatively low top quark mass uncertainties for the p T,hh and ∆φ hh distributions, except for the low p T,hh and the ∆φ hh ∼ π regions, typically below the size of the scale uncertainties, and larger uncertainties for the leading-jet transverse momentum, for which the difference between the two approximations is larger.

Summary
In this work we considered Higgs boson pair production through gluon fusion in proton collisions. We presented new QCD predictions for inclusive and differential cross sections, which include the full NLO contribution and also account for finite top quark mass effects at NNLO. Our best prediction, denoted NNLO FTapprox , retains the full top quark mass dependence in the double-real emission amplitudes, while the remaining real-virtual and two-loop virtual HEFT amplitudes are treated via a suitable reweighting for the corresponding subprocesses with a given final-state multiplicity. This approximation represents the most advanced prediction available to date for this process.
The numerical results we obtained for the NNLO FTapprox are quantitatively different from the results obtained in previous combinations. In particular, as far as the total cross section is concerned, the corrections turn out to be smaller than previous estimates, increasing the NLO result by about 12% at 13 TeV and 7% at 100 TeV. The reduction of the scale uncertainties is significant, by about a factor of three for LHC energies. Given that our NNLO FTapprox prediction includes top quark mass effects in an approximated way, it is important to assess the corresponding uncertainty. We carefully examined the performance of our approximations at both the inclusive and differential levels. The uncertainty on our reference inclusive NNLO FTapprox prediction is estimated to be about ±2.7% at 14 TeV, increasing with the collider energy to reach ±4.6% at 100 TeV.
Regarding differential distributions, in most of the cases we can observe clear qualitative differences with respect to the bin-by-bin reweighting procedure introduced in Ref. [17], in the shape and/or the normalization. For some of the distributions, however, specifically the tails of the p T,hh and p T,h1 spectra, both approximations are in very good agreement. We discussed an estimate of the uncertainty associated with top quark mass effects at NNLO at the differential level, and we found that in most of the cases its magnitude is comparable to the size of the scale uncertainties, except for the tails of some distributions where the uncertainty from missing M t effects can be dominant.