The gluon-fusion production of Higgs boson pair: N$^3$LO QCD corrections and top-quark mass effects

The Higgs boson pair production via gluon fusion at high-energy hadron colliders, such as the LHC, is vital in deciphering the Higgs potential and in pinning down the electroweak symmetry breaking mechanism. We carry out the next-to-next-to-next-to-leading order (N$^3$LO) QCD calculations in the infinite top-quark mass limit and present predictions for both the inclusive and differential cross sections, albeit the differential distributions other than the invariant mass distribution of the Higgs boson pair are approximated at N$^3$LO. Such corrections are indispensable in stabilising the perturbative expansion of the cross section in the strong coupling $\alpha_s$. At the inclusive level, the scale uncertainties are reduced by a factor of four compared with the next-to-next-to-leading order (NNLO) results. Given that the inclusion of the top-quark mass effects is essential for the phenomenological applications, we use several schemes to incorporate the N$^3$LO results in the infinite top-quark mass limit and the next-to-leading order (NLO) results with full top-quark mass dependence, and present theoretical predictions for the (differential) cross sections in the proton-proton collisions at the centre-of-mass energies $\sqrt{s}=13,14,27$ and $100$ TeV. Our results provide one of the most precise theoretical inputs for the analyses of the Higgs boson pair events.


Introduction
In view of the null results in the beyond the Standard Model (BSM) searches so far at colliders, it seems that a realistic way of looking for new physics in the future is to precisely study the nature of the Higgs sector. Any small deviation with respect to the Standard Model (SM) predictions would indicate the signal of new physics. In particular, the electroweak symmetry breaking mechanism remains to be understood. It can be deciphered by specifying the form of the Higgs potential. In the SM, such a potential is determined by two SU(2) L ×U(1) Y gauge invariant renormalisable operators constructed from a single Higgs SU(2) L doublet H = H + , H 0 T , i.e.
This Higgs potential has a well-known shape of a "Mexican hat". The spontaneous symmetry breaking happens after the Higgs field captures a non-vanishing vacuum expectation value v, which is related to the Fermi constant G F via v = √ 2G F −1/2 = µ 2 /λ SM 1/2 . The quantum fluctuation of the real scalar field around the minimum value of the potential V (H) at H 0 = 0, v/ √ 2 T represents a physical Higgs boson h. The Higgs boson mass at tree-level is given by m 2 h = 2µ 2 , and the Higgs self-interactions become One can see that the Higgs potential in the SM is fully determined by G F and m h , whose values have been measured precisely [1]. Therefore, independent measurements on the Higgs trilinear and quartic couplings are very important to test the SM Higgs sector. In fact, several UV-complete new physics models predict modifications of the Higgs potential and the Higgs trilinear coupling λ hhh [2,3]. Some of them (see e.g. refs. [4][5][6]) can possess very different λ hhh value from the SM expectation λ SM hhh = λ SM = m 2 h /2v 2 but still have SM-compatible Higgs interactions with the massive gauge bosons and fermions. The measurement of the Higgs self couplings seems the only way to understand the dynamics of electroweak symmetry spontaneously breaking.
The Higgs trilinear coupling can be either directly probed via the Higgs boson pair production or indirectly constrained by using the loop effects in the precision observables (e.g. the single Higgs boson signal strengths at the LHC [7][8][9][10][11][12] or at an e + e − collider [13], the electroweak oblique parameters [14], or the W boson mass and the effective sine [15]). The existing direct measurements of the Higgs pair cross sections at the LHC only loosely bound λ hhh [16,17] due to the low statistics. The current best constraint −5 < λ hhh /λ SM hhh < 12 at 95% confidence level (CL) is from the ATLAS collaboration with 36.1 fb −1 Run-2 data [16]. The situation will be largely improved at the phase of the HL-LHC with 3000 fb −1 integrated luminosity [18]. Meanwhile, novel analysis techniques (e.g. new kinematic variables [19] or machine learning [20]) have been proposed to expedite the discovery. In addition, the envisaged future hadron colliders, like the FCC-hh, are expected to be the ultimate precision machines for determining λ hhh [21], strongly gaining from both the 20 times bigger cross section and the higher integrated luminosity.
Although this process is mainly limited by the low statistics at the moment, the continuous measurements at the LHC are still quite valuable, because even the loose bounds can already exclude some new physics models or corner the parameter space, which predicts the enhanced yields of pp → hh (see e.g. [22]). The indirect constraints on λ hhh from the single-Higgs data have been set in the range λ hhh /λ SM hhh ∈ [−3.2, 11.9] at 95% CL with 79.8 fb −1 Run-2 data by ATLAS [23]. These constraints are already comparable with the direct ones and impact the final bounds with the combinations of the direct and indirect measurements [24]. As opposed to the direct bounds, the improvements of the indirect bounds are limited by the systematics and thus will be harder at the HL-LHC. Nevertheless, these indirect approaches feature different systematics than direct measurements and can be thought as independent cross checks. On the other hand, the extraction of the quartic Higgs self-coupling from the triple Higgs production is much difficult (though not hopeless) at hadron colliders, because the corresponding cross sections are three orders of magnitude smaller than the double-Higgs production [25]. Similar to the single Higgs hadroproduction, the dominant di-Higgs production channel at a high-energy hadron collider is via the gluon-gluon fusion (ggF) [3,26,27]. Other channels are at least one order of magnitude lower in their yields. Due to the absence of the tree-level interactions between the Higgs boson and gluons in the SM, the leading order (LO) cross section σ(gg → hh) was computed from one-loop amplitude squared [28][29][30], where two representative LO Feynman diagrams can be seen in the first row of figure 1. Further improvements of the fixed-order perturbative calculations without any approximation are quite challenging. The full next-to-leading order (NLO) QCD calculations involving complicated two-loop Feynman integrals were carried out only recently [31][32][33][34] thanks to the new advances of the numerical approaches [35][36][37]. The NLO results were complemented with soft-gluon resummation [38] or parton-shower (PS) effects [39][40][41]. The ggF NLO predictions are plagued with the large theoretical uncertainties from the scale variations [31] and the top-quark mass scheme dependence [33]. Moreover, at NLO+PS, some differential distributions (e.g. the distribution at large transverse momentum of the Higgs boson pair) differ significantly by adopting different matching schemes [39] or shower scales [40].
Instead of starting from the loop-induced process, one can also carry out the heavy top-quark mass m t expansions in the amplitudes. We refer to the leading expansion term in 1/m 2 t as the infinite top-quark mass limit m t → +∞. In such an approximation, the two Higgs bosons can be generated by the two gluon scatterings at tree level (see the second row of figure 1), which makes the higher-order perturbative calculations more feasible.
The first NLO computation in the m t → +∞ limit was performed two decades ago [42]. Next-to-next-to-leading order (NNLO) was also available [43][44][45][46], and recently we have presented the first next-to-next-to-next-to-leading order (N 3 LO) calculation [47]. Besides these fixed-order results, the soft-gluon resummation effects are also considered in refs. [48][49][50]. In spite of the success of improving the perturbative accuracy in the cross section calculations, it is widely acknowledged that the m t → +∞ approximation is insufficient for the phenomenological applications. Many theoretical efforts have been devoted to investigate the finite m t corrections to this approximation [25,27,[51][52][53][54][55]. Moreover, there are also many well-motivated attempts to evaluate the involved two-loop gg → hh amplitudes in the analytic forms by taking other approximations (e.g. in the small top-quark mass [56,57], the small Higgs transverse momentum [58] and the small Higgs mass [59] limits).
The primary goal of this paper is to extend our previous N 3 LO results in ref. [47] and to include the top-quark mass effects for the phenomenological applications. The remaining context is organised as follows. In section 2, after the description of our method, we provide the validation of our calculations as well as the extensive numerical results in the infinite top-quark mass limit. We take into account the finite m t effects at N 3 LO based on the NLO QCD results with full m t dependence in section 3. The conclusion is drawn in section 4. Additional results and some technical details can be found in the appendices. The hard functions, in particular the new one-loop analytic expressions, are shown in appendix A. An NLO model and the R 2 Feynman rules are described in appendix B. The renormalisation scale dependence in the N 3 LO results is discussed in appendix C. Finally, appendix D collects the additional plots for the differential distributions.
2 N 3 LO corrections in the infinite top-quark mass limit

Effective Lagrangian and Wilson coefficients
The interactions between the Higgs bosons and gluons are mainly generated by top-quark loops, where two LO Feynman diagrams are shown in the first row of figure 1. The effective Lagrangian in the infinite top-quark mass limit is obtained through integrating out the top-quark loop contribution (see the second row of figure 1). For the Higgs boson pair production, the relevant effective Lagrangian can be written as where α s is the strong coupling and G a µν is the gluon field strength tensor. On the right hand side of the second equation, O(h k , k ≥ 3) means that we have ignored terms involving more than two Higgs bosons in the effective Lagrangian. The Wilson coefficients δ and η, or equivalently C h and C hh , comprise the QCD radiative corrections of the top-quark loops. C h and C hh can be easily derived in terms of δ and η as These Wilson coefficients can be perturbatively expanded in a series of α s , Their four-loop analytic expressions are already known in the literature [43,45,[60][61][62][63][64][65][66][67][68][69]. In our N 3 LO calculations, the results up to three loops are needed. They are given in the on-shell top-quark mass scheme by [68] and where L t = ln(µ 2 R /m 2 t ), m t is the top-quark pole mass, µ R is the renormalisation scale and n f is the number of the light-quark flavours. The expressions of C h and C hh are

Breakdown in three channels
The ggF Higgs boson pair production in the infinite top-quark mass limit with the effective Lagrangian defined in eq.(2.1) can be divided into three channels according to the number of effective vertices at the squared amplitude level. Three representative Born cut-diagrams are shown in figure 2. There are two (class-a), three (class-b) and four (class-c) effective vertices insertions respectively. In other words, the double-Higgs (differential) cross section can be decomposed into Because there are at least one α s power in the Wilson coefficients C h and C hh , their Born cross sections contribute to different α s orders, which are summarised in table 1. The lowest orders of class-a, -b and -c are O(α 2 s ), O(α 3 s ) and O(α 4 s ) respectively, which means that they contribute to LO, NLO and NNLO parts of the Higgs boson pair cross section. For the purpose of N 3 LO calculations in the present paper, we need to calculate N 3 LO, NNLO and NLO corrections to the class-a, -b and -c part, respectively.  Table 1. The perturbative orders in α s for different classes at the amplitude squared level. We call the O(α 3 s ) contribution in class-b as the LO in this class though it is an NLO correction to the cross section of Higgs pair production. The same rule applies to the class-c part.

The class-a part
We have two approaches to compute NNLO (i.e. up to O(α 4 s )) cross section in the class-a part. The first one is that we can perform a fully-differential NNLO calculation based on the q T -subtraction method, which was originally proposed in ref. [70]. 1 In this paper, we will use the q T -subtraction method in the framework of the soft-collinear effective theory (SCET) [92][93][94][95][96]. In this approach, the class-a (differential) cross section can be further divided into where p hh T is the transverse momentum of the Higgs pair system, i.e. q T = p hh T . The first (second) term on the right-hand side of eq.(2.8) is imposed the kinematic cut p hh T < p veto is computed with the aid of the transverse-momentum resummation formalism in SCET. The cross section of this part is factorised as a convolution 1 With qT -subtraction method, tremendous works have been done at the NNLO accuracy [46,[70][71][72][73][74][75][76][77][78][79][80][81][82][83][84][85][86][87][88].
Through solving the renormalisation equations up to N 3 LO, the small qT cross section has also been studied at N 3 LO for certain processes [89][90][91] with constant terms missing at three loops in the collinear sector.
of the hard, beam and soft functions dσ a hh dp hh where we have ignored the power-suppressed terms O . Such a factorisation formalism holds when p hh T is sufficiently smaller than the hard scale Q, which is derived by studying the IR behaviour of QCD. The transverse-momentum dependent (TMD) gluon beam function B g is universal in the sense that it is independent of the process but only relies on the species of the initial state (i.e. gluon). The soft function S is also the same for all processes only involving colourless final states with the gluon-gluon initial state. The calculations of the TMD beam and soft functions can be carried out with a rapidity regulator, while the physical results are independent of the choice of such a regulator. The two-loop analytic results for these TMD beam and soft functions can be found in [97][98][99][100][101][102], and the N 3 LO results have been obtained very recently [103,104]. On the other hand, the hard function H a is process dependent. The detailed discussions about the hard functions can be found in appendix A.
Due to the non-vanishing transverse momentum of the Higgs pair system p hh T > p veto T , only the events with additional jets will be maintained in the second piece of eq.(2.8) dσ a hh p hh T >p veto T . In our case, the NNLO computation of class-a requires us to calculate the NLO corrections to a Higgs pair plus a jet with two effective vertices insertions. Such a task can be carried out by using the automated simulation framework MadGraph5_aMC@NLO (MG5_aMC henceforth) 2 [106] with an NLO Universal FeynRules Output (UFO) model [112] based on the SM Lagrangian and the effective Lagrangian eq.(2.1). The details about the model, in particular the analytic expressions of the rational R 2 terms, can be found in appendix B. Due to the different α s orders in the three Born classes, we need the recent development [113] in MG5_aMC that is capable of handling mixed-order scenarios.
Within the q T -subtraction approach, the independence on p veto T in the finite cross section should always be guaranteed when p veto T is approaching zero. We have explicitly checked this in the NNLO class-a cross section σ a,NNLO hh shown in figure 3. Alternatively, the class-a cross section can be related to the single Higgs production cross section, because they share exactly the same topology in the infinite top-quark mass limit. In the di-Higgs case, the class-a part can be viewed as the production of an off-shell Higgs boson from ggF and its decay into two on-shell Higgs bosons. The off-shell Higgs boson has an invariant mass of the final-state Higgs boson pair m hh . The explicit relation is  where the function f h→hh accounts for the phase space factor mapping from the single Higgs production to the Higgs pair production, and σ h denotes the cross section for the single Higgs boson production. The replacement m h → m hh in eq.(2.10) means that the cross section is calculated with the Higgs mass m hh .
In the first parentheses of the right-hand side of eq.(2.10), C hh C h accounts for the Wilson coefficient difference in figure 1c, while the second term takes into account the propagator of the off-shell Higgs and the Higgs self-coupling in figure 1d. Such a method has already been used in the previous NNLO calculation of the ggF di-Higgs production in ref. [44].
We have compared the results with the above two independent approaches for NNLO class-a cross sections shown in the left panel of figure 4. The calculation with q T -subtraction matches the result by using eq.(2.10) and iHixs2 within the Monte Carlo integration errors when p veto T ≤ 16 GeV. Thus, we have validated eq.(2.10). After inclusion of class-b and class- hh ), we can compare our two calculations with the previous NNLO di-Higgs calculation in ref. [46]. As shown in the right panel of figure 4, we have obtained perfect agreement with ref. [46].
In order to compute the N 3 LO class-a cross section, we need to know the N 3 LO cross section of σ h . Since σ h is only known inclusively (i.e. total cross section) at N 3 LO, we only perform the exact N 3 LO calculations for the total inclusive cross sections and the invariant mass distributions of the class-a part. In the present paper, we will use the public code iHixs2 [114]

The class-b part
In order to achieve the N 3 LO accuracy for the di-Higgs cross sections in the infinite topquark mass limit, we have to calculate the NNLO QCD corrections to the class-b part. The NNLO cross sections for the class-b part were computed with the q T -subtraction method similarly as described in the previous section, i.e. the differential cross section is decomposed into The two pieces dσ b hh p hh T <p veto T and dσ b hh p hh T >p veto T in eq.(2.12) can be computed using the method described above in the class-a part. Therefore, we will refrain ourselves from describing them again except the hard function H b in the following equation The explicit expression of H b is shown in appendix A. The p veto T independence of dσ b hh after summing the two pieces is explicitly verified in figure 5. As opposed to NNLO cross section of the class-a part, we do not have a second independent cross section calculation for this part. The NLO cross section σ b,NLO hh however can be easily checked with MG5_aMC as shown in figure 6. The perfect agreement below permille level is achieved when p veto T ≤ 8 GeV. At N 3 LO, the renormalisation scale cancellation is guaranteed only when combining the class-a and class-b parts, which will be detailed in appendix C. It can serve as another powerful check to the NNLO class-b cross section. The class-a (differential) cross section    can be decomposed into

The class-c part
We only need the NLO QCD corrections to the class-c part in order to give N 3 LO di-Higgs cross sections. The computations can be achieved with the full-fledged NLO techniques. We have compared the NLO cross sections for the class-c part between the q T -subtraction approach and the automated calculation by MG5_aMC in figure 8. The perfect agreement is found when p veto T ≤ 6 GeV. We have summarised the independent calculations we have performed with different approaches for the three classes contributing to various orders in table 2.

Calculational setup
In our numerical calculations, we take v = 246.2 GeV and the Higgs boson mass m h = 125 GeV. The top-quark pole mass, which enters only into the Wilson coefficients, is m t = 173

NLO
NNLO GeV. Unless it is explicitly specified, the trilinear Higgs coupling λ hhh is taken to be the SM value. We use the PDF4LHC15_nnlo_30 PDF [115][116][117][118] available in the programme LHAPDF6 [119], and the associated α s . The default central scale is chosen to be the invariant mass of the Higgs boson pair divided by 2, i.e. µ 0 = m hh /2, and the scale uncertainty is evaluated through the 9-point variation of the factorisation scale µ F and the renormalisation scale µ R in the form of In the parts of ultilising the q T -subtraction method, we will use p veto We have verified that the uncertainties due to the missing power-suppressed terms of are well below the Monte-Carlo integration errors.

Inclusive total cross sections
We present the inclusive total cross sections from LO to N 3 LO at different centre-of-mass energies √ s = 13, 14, 27, 100 TeV in table 3, where the scale uncertainties are also shown.
These particular energies are either the LHC energies (13 and 14 TeV) or the nominated energies for the future hadron colliders [18,21]. The cross sections from √ s = 7 TeV to √ s = 100 TeV are also displayed in the left panel of figure 9, where the bands represent the scale uncertainties. Similarly to the case of single Higgs production, the QCD corrections in the di-Higgs process are very prominent. The NLO QCD corrections increase the LO cross sections by 87% (85%) at √ s = 13 (100) TeV. The NNLO QCD corrections improve the NLO cross sections further by 18% (16%), reducing the scale uncertainties by a factor of two to three to be below 8%. The N 3 LO QCD corrections enhance the NNLO cross section by 3.0% (2.7%). The cross sections lie well within the scale uncertainty bands of the NNLO results, and the N 3 LO scale uncertainties are less than 3% and 2% at 13 and 100 TeV respectively. In addition, the PDF parameterisation uncertainties are almost independent of the QCD corrections. Their relative sizes amount to ±3.3%, ±3.1%, ±2.2% and ±1.4% with respect to the central values at 13, 14, 27 and 100 TeV, overwhelming the remaining N 3 LO scale uncertainties. We have also shown the contribution from three different classes separately in the right panel of figure 9, where the class-b contribution has been multiplied by a factor of -1 in order to make it visible in the frame. There is a strong hierarchy among the three classes. Typically, the class-b part is only a few percent of the class-a, while the class-c is a few percent of the class-b. Such a behaviour can be understood from the effective Lagrangian eq.(2.1). One more effective vertex in the squared amplitude results in one more factor of αs 3π ∼ 1% suppression instead of the usual α s suppression.  Table 3. The inclusive total cross sections (in unit of fb) of Higgs boson pair production in the infinite top-quark mass limit at different centre-of-mass energies √ s from LO to N 3 LO. The quoted relative uncertainties are from the 9-point scale variations. The errors due to the numerical Monte Carlo integration are well below 1 .
It was proposed in ref. [120] to use the ratios of cross sections with the same final state between different centre-of-mass energies to perform precision studies (e.g. determining PDFs) and to improve the BSM sensitivities 3 . The success of such a programme relies on the large cancellations of theoretical systematic uncertainties in the ratios. In particular, the usually dominant scale uncertainties in the cross sections can be significantly reduced by fully correlating the renormalisation and factorisation scales between numerators and denominators. Such a reasonable working assumption, however, should be carefully checked when higher-order calculations become available. With the N 3 LO calculations we have done, we can readily check such a hypothesis in the double-Higgs process. In figure 10, we have plotted the cross section ratios in six different √ s pairs from LO to N 3 LO. The scale correlation assumption in the cross section ratios is indeed verified in this process. Apart from the dependence on the collision energy, it is also very interesting to know how total cross sections vary when λ hhh deviates from the SM value. At four different centre-of-mass energies √ s = 13, 14, 27, 100 TeV, we have varied κ λ = λ hhh /λ SM hhh from −4 to 8 in figure 11. The largest deconstruction between the λ hhh -independent amplitude (e.g. from figure 1c) and the λ hhh -dependent amplitude (e.g. from figure 1d) occurs when κ λ is close to 2. The N 3 LO corrections only marginally distort the NNLO predictions around κ λ = 2. This can be understood because the QCD radiative corrections to the above two kinds of different amplitudes are not very different due to the same Lorentz structure shared between figure 1c and figure 1d.

Invariant mass distributions
Besides the total cross sections, we are also able to calculate the exact N 3 LO results for the invariant mass m hh distributions, which are shown in figure 12 with the 4 different energies √ s = 13, 14, 27, 100 TeV. Due to the larger phase space, the m hh spectrum becomes harder when √ s increases. The inclusion of the N 3 LO QCD corrections dramatically stabilises the perturbative calculations of the invariant mass differential distributions. The N 3 LO corrections only marginally change the shapes, and the N 3 LO results, which have very small scale uncertainties, are completely enclosed within the NNLO uncertainty bands. Such a feature consolidates that the perturbative expansions of the invariant mass differential cross sections are converging in α s up to the fourth order.
It is also very interesting to investigate how the invariant mass distribution changes with respect to the value of κ λ = λ hhh /λ SM hhh . We have shown the LO to N 3 LO distributions with κ λ = −1 (upper left), 3 (upper right) and 5 (lower left) in figure 13. In addition, the comparison of N 3 LO m hh distributions with four values κ λ = −1, 1, 3, 5 is given in the lower right panel of figure 13. The differential distribution dramatically changes when κ λ varies. This feature can be understood qualitatively by looking at eq.(2.10). σ h | m h →m hh decreases monotonically when increasing m hh , which explains the behaviour in the large invariant mass regime. At small m hh (i.e. m hh → 2m h ), the distribution is governed by . Given the phase space boundary m hh ≥ 2m h , the second prefactor is a monotonically decreasing and then monotonically increases when m hh > √ 1 + 3κ λ m h . This explains the fact that the suppression at threshold m hh → 2m h is more dramatic in the SM case κ λ = 1 than others. On the other hand, when m hh approaches √ 1 + 3κ λ m h (395 GeV for κ λ = 3 and 500 GeV for κ λ = 5), a cancellation , which results in the dip structures in figure 13 for the κ λ > 1 cases. These interesting features can be definitely used in the BSM searches via the di-Higgs final states [20].

Other differential distributions
In order to carry out N 3 LO calculations for other differential distribution, we have to take some approximations, because the fully-differential N 3 LO corrections to single Higgs production are still unknown. Therefore, at the moment, we have to approximate the N 3 LO  class-a corrections for other differential cross sections. As we already mentioned in section 2.3.2, the class-a differential cross sections can be divided into two pieces given in eq.(2.14). The second piece dσ    ferential calculations for the NNLO class-a cross sections with the q T -subtraction method. Therefore, in our paper, we can define the approximated N 3 LO (AN  The (a, 1) piece is simply multiplying a global K factor σ (a,1),N 3 LO hh σ (a,1),NNLO hh assuming no kinematic dependence. Such an assumption is more-or-less justified given the extremely flat K factor found in the rapidity distributions of the single Higgs process [122]. In contrast, the exact fully-differential predictions are achievable for other three pieces at O(α 5 s ). Our calculations can certainly be improved as long as the fully-differential N 3 LO calculation of the single-Higgs process is available.
We have shown 6 differential distributions in figure 14  Like the dijet hadroproduction case [123,124], the region of p T (h 1 ) → 0 is largely populated by IR quanta radiations, which makes fixed-order perturbative calculations problematic. In addition, the ∆φ distribution is quite special as all the LO events locate at ∆φ = π, i.e. back-to-back of the Higgs boson pair in the transverse plane. For all the ∆φ = π bins, an N k LO calculation only gives the N k−1 LO accuracy. On the contrast, the N k LO accuracy can be achieved by a complete N k LO calculation in the end point ∆φ = π. The region is however sensitive to the soft gluon emissions, which yields large logarithms to spoil the fixed-order perturbative calculations. Such a feature can be deduced from the fact that the scale uncertainty bands do not shrink from LO to AN 3 LO. The pathological behaviour should be cured after performing the soft-gluon resummation in the region.

Top-quark mass approximations at N 3 LO
It is well known that the top-quark mass effects are important in the Higgs boson pair production. Therefore, any relevant phenomenology studies should take into account these effects. However, the direct improvements of perturbative calculations with full top-quark mass dependence are technically very challenging because the lowest order is already loopinduced. The state-of-the-art calculation without performing 1/m 2 t expansion is NLO in α s . A standard way to improve the perturbative calculations is to combine the NLO full top-quark mass calculations (denote as NLO mt ) with the higher-order infinite top-quark mass calculations. The combination of the two different calculations is not unique, and therefore relies on various approximations.  There are several approximations to combine the differential cross sections in the infinite top-quark mass limit dσ N k LO mt→∞ and those with full top-quark mass dependence dσ N l LO mt (l < k). In our case, we have k = 0, 1, 2, 3 and l = 0, 1. They are: • N k LO⊕N l LO mt : this approximation simply improves the leading m t expansion term (3.2) • N k LO⊗N l LO mt : this assumes that the QCD K factor in the infinite topquark mass limit also applies to the other top-quark mass dependent terms. It is defined as Other approximations are of course still possible (e.g. those introduced in refs. [25,54]). However, they require the knowledge of the fully-differential distributions, which is not known at N 3 LO. In particular, the "FT approximation" 5 introduced in refs. [25,54] is considered as the most advanced predictions. We leave the FT approximation at N 3 LO for a future study. Here, we decide to restrict ourselves with the above three approximations. Among them, N k LO⊗N l LO mt is expected to be the most accurate predictions, while N k LO⊕N l LO mt is the worst approximation because the finite top-quark mass effects are missing in the correction ∆σ k,l mt→∞ . In the following, we will present the results under three approximations for comparison.

Results
With the same setup as described in the section 2.4.1, the full m t -dependent NLO (differential) cross sections can be obtained by the public code [39,41] available in the Powheg-Box [125][126][127]. The scale uncertainties for each approximation in the present paper are estimated by taking the envelope of 9-point variations ξ R = µ R /µ 0 , ξ F = µ F /µ 0 with 5 In the so-called FT approximation, the matrix elements in the infinite top-quark mass limit for each partonic subprocess are improved/reweighted by the ratios of the one-loop full top-quark mass squared amplitudes over the tree-level mt → +∞ squared amplitudes. µ 0 = m hh /2, ξ R , ξ F ∈ {0.5, 1, 2}. The (differential) cross sections at each point are defined as . (3.4)

Inclusive total cross sections
The inclusive total cross sections after taking into account the top-quark mass effects are tabulated in table 4. The NLO cross section with full top-quark mass dependence (denoted by NLO mt ) is 27.56 fb at √ s = 13 TeV, 6 which is 6.8% larger than the result in the infinite top-quark mass limit (denoted by NLO) shown in table 3. However, at 100 TeV, the NLO mt cross section 7 is more than 3 times smaller than the NLO result. This indicates that the large top-quark mass approximation is not valid any more at a very high energy collider.
The remaining scale uncertainties in NLO mt cross sections are beyond 10%. Such theoretical uncertainties are expected to be reduced by including higher-order QCD corrections. We evaluated the NNLO and N 3 LO cross sections by using three approximations defined in the previous section based on the NLO mt results. The central values as well as the scale uncertainties are presented in table 4. Because the finite m t corrections in ∆σ k,1 mt→∞ , k = 2, 3 are still missing, the N k LO⊕NLO mt approximation is least accurate and even not reliable at 100 TeV, which is also implied in the shown pathological scale uncertainties. In contrast, both N k LO B−i ⊕NLO mt and N k LO⊗NLO mt approximations have partially captured the finite top mass effects in the higher-order QCD correction pieces ∆σ k,1 mt→∞ . The differences between the two different approximations can be viewed as a way of estimating the remaining 1/m 2 t uncertainties, which are around 2-3%. In particular, we take N 3 LO⊗NLO mt predictions as the state-of-the-art. The relative scale uncertainties in N k LO⊗NLO mt are identical to those in N k LO.  6 We have verified that the slightly offsets between our NLOm t results and those in ref. [54] at √ s = 13, 14 TeV can be attributed to the different PDFs. In our calculations, we always use the same NNLO PDF, while ref. [54] used a NLO PDF for the NLO calculations and a NNLO PDF in the NNLO calculations.   Table 4. The inclusive total cross sections (in unit of fb) of Higgs boson pair production at different centre-of-mass energies √ s within the considered approximations. The quoted relative uncertainties are from the 9-point scale variations.

Invariant mass distributions
is degraded to NLO accuracy when m hh becomes larger than two times of the top-quark mass where the scale cancellations are not guaranteed. For N 3 LO⊗NLO mt , because of the manner of varying ξ R , ξ F in differential cross sections eq.(3.4), their relative scale uncertainties are exactly the same as N 3 LO in section 2.4.3. Comparisons between NNLO⊗NLO mt and N 3 LO⊗NLO mt predictions are given in figure 16. Similar to what has been found at NNLO in ref. [54], the higher-order QCD corrections are quite small near the threshold region m hh 2m h . The K factors are almost constants (around 1.2) at larger m hh . A lesson from NNLO tells us that the NNLO⊗NLO mt predictions feature different shapes as the FT approximation. Therefore, it would be quite desirable to carry out the latter approximation at N 3 LO, which is however beyond the scope of the present paper.

Other differential distributions
With the approximation eq.(2.15) used at N 3 LO in other observables, we are able to report our predictions for fully differential distributions of the Higgs boson pair production. We have shown 6 differential kinematic distributions at √ s = 14 TeV in figure 17 as our illustrative examples, while the same differential cross sections at √ s = 13, 27, 100 TeV can be found in appendix D. These kinematics are the rapidity of the Higgs pair (up left panel of figure 17), the rapidity of a random Higgs boson (up right panel of figure 17), the transverse momenta p T of the harder (middle left panel of figure 17) and the softer Higgs (middle right panel of figure 17), the absolute rapidity difference |∆y| (low left panel of figure 17) and the azimuthal angle difference ∆φ (low right panel of figure 17) between the two Higgs particles. For the sake of clarity, we will only show the results of NLO mt (black), NNLO⊗NLO mt (dark-orange) and AN 3 LO⊗NLO mt (blue), where we have adopted the AN 3 LO calculations to approximate the N 3 LO differential cross sections.  1.2. The shape of the distribution is mainly driven by the partonic luminosity encoded in the PDF. The scale uncertainty band is reduced from NNLO⊗NLO mt to AN 3 LO⊗NLO mt by a factor of four.
Because the rapidity distributions of the leading-p T and subleading-p T Higgs bosons are sensitive to soft-gluon radiations, i.e. not IR safe at fixed orders, we instead show the  rapidity distribution of a random Higgs boson. The latter histogram is equivalent to the arithmetic mean of the former two histograms. Similar to the y hh distribution, the higherorder QCD corrections only change the shape slightly. The central region has a bit larger radiative corrections than the forward and backward regions. The difference is however quite insignificant, which is only at 1-2 percent level.  O(α 5 s ) corrections is evident from the obvious reduction of theoretical uncertainties. The differential cross sections in the transverse momenta of the leading-p T (harder) and the subleading-p T (softer) Higgs bosons can be found in the middle panels of figure 17. These two transverse momenta are identical at LO. Beyond LO, due to the presence of extra real radiations, the difference between the two emerges. It is quite often that the Higgs boson will pick up a larger p T if it recoils against the hardest real radiation. For this reason, the real emission topologies are dominant in the tail of the p T (h 1 ) distribution, which results in the growth of the scale uncertainties in the high p T (h 1 ) bins. The AN 3 LO⊗NLO mt scale uncertainty is +2% −5% at the bin p T (h 1 ) ∈ [800, 900] GeV, while those of NNLO⊗NLO mt and NLO mt are +7% −10% and +25% −19% respectively in the same bin. At low p T (h 1 ), the QCD radiative corrections become perturbatively unstable 8 due to the large logarithms of (p T (h 1 ) − p T (h 2 )) /µ 0 < p T (h 1 )/µ 0 → 0. The scale uncertainties of AN 3 LO⊗NLO mt are larger than NNLO⊗NLO mt in the first three bins. Such a pathological behaviour reflects the fact that more large logarithms due to the soft-gluon radiations appear in the higher order α s calculation. On the other hand, the subleading p T distribution receives quite uniform K factors at NNLO and AN 3 LO except the first bin, where the K factors are lower in the first bin than others. It has been shown in ref. [32] that the NLO QCD corrections are vanishing in the tail of the p T (h 2 ) distribution. This makes the NLO mt scale variation very small. However, we do not have an understanding in depth for such a behaviour at the moment. In the tail, we find that only AN 3 LO⊗NLO mt has the comparable size of the scale variation with NLO mt .
Finally, we are in the position to discuss the two kinematic correlation distributions between the Higgs boson pair. They are the rapidity difference ∆y and the azimuthal angle difference ∆φ in the low two panels of figure 17. The significance of the higher-order QCD corrections is slowly reduced from the two near Higgs boson (|∆y| ∼ 0) region to the region where the two Higgs particles are far away (i.e. |∆y| is large). This is because a large |∆y| usually corresponds to a large invariant mass of the Higgs pair m hh , where the latter is proportional to our hard scale. In particular, in Born kinematics, we have m hh = 2 m 2 h + p 2 T cosh ∆y 2 , where p T is the transverse momentum of an arbitrary Higgs boson.
The radiative corrections are dramatic in the ∆φ distribution. All the Born-like 2 → 2 events locate at ∆φ = π, as the two Higgs bosons are always in the back-to-back configuration in the transverse plane. All the contributions to the ∆φ < π regime must be from the events with at least one additional jet in the final states. In the bins of ∆φ < π, NLO mt , NNLO⊗NLO mt and AN 3 LO⊗NLO mt results correspond to the true LO, NLO and NNLO accuracy in α s . The K factor AN 3 LO⊗NLOm t NLOm t increases slightly from ∆φ = 0 to ∆φ = 0.8π and then drops quickly from ∆φ = 0.8π to ∆φ = π. The 9-point scale variations shift their central values by +45% −29% , +15% −17% and +10% −13% respectively in the first bin ∆φ ∈ [0, 0.05]π. The uncertainty reduction from NNLO to AN 3 LO is not as immense as in other cases. Since a small kick by soft gluon radiations will make the two Higgs bosons not being back-to-back anymore, a reliable prediction for the region ∆φ ∼ π can only be achieved after performing a resummation calculation.

Assessment of the top-quark mass approximations
Before we close the section, we will discuss how good are our top-quark mass approximations. Since the full NNLO and N 3 LO calculations with the full m t dependence are absent, the way of estimating the remaining m t uncertainties is not unique. One obvious way is to assess the missing m t uncertainties by trying different approximations. This has been taken at NNLO in ref. [54] even with the most advanced onethe FT approximation. In the inclusive cross sections, the FT approximation gives smaller predictions than other approximations, including the NNLO⊗NLO mt approximation, because of the additional m t contributions in the former. The difference is amplified a bit with the increasing of √ s. At NNLO, the difference between the FT approximation and the NNLO⊗NLO mt approximation is 5% at 13 TeV to 9% at 100 TeV. This is not surprising since the m t corrections become more important at larger energies. Given that the m t corrections are more-or-less orthogonal to the α s corrections, we expect the N 3 LO⊗NLO mt numbers in table 4 should be lowered by a similar amount after we applied the FT approximation at N 3 LO. Besides this normalisation, the shapes of NNLO⊗NLO mt and the FT approximation at NNLO are very close for y hh , p T (h 1 ), p T (h 2 ) and ∆φ distributions, while those for m hh are quite distinct. The deviation between the N 3 LO FT approximation and the N 3 LO⊗NLO mt scheme can be viewed as a way to assign the theoretical uncertainties from the missing top mass corrections. Such a difference is expected to be similar to what has been found at NNLO [54].
We can also follow the NLO discussions in ref. [32] to assess the goodness of our top-quark mass approximations in the differential distributions, where both the NLO mt (the results with the notation "NLO" in ref. [32]) and the NLO⊗LO mt (those with the notation "B-i, NLO HEFT" in ref. [32]) cross sections were computed. However, since we have already used the full NLO mt in our calculations, the remaining m t uncertainties are expected to be at least α s suppressed with respect to the estimations from NLO versus NLO⊗LO mt . The total cross sections are lowered by 14% (24%) at √ s = 14 (100) TeV from the NLO⊗LO mt approximation to the complete NLO mt calculations. Both the FT approximation at NLO and the NLO⊗LO mt results overestimate the true NLO QCD corrections at large m hh , p T (h 1 ), p T (h 2 ). On the other hand, the shapes of the rapidity distributions are quite similar between NLO⊗LO mt and NLO mt . In the former cases, the missing m t correction uncertainties will be underestimated by using the first approach described in the previous paragraph. Instead, a better way to assess the top mass corrections in N 3 LO⊗NLO mt is to multiply (dσ NLO mt − dσ NLO⊗LOm t ) with α s .

Summary
In the paper, we first carried out the N 3 LO QCD corrections to the Higgs boson pair production via ggF at high-energy hadron colliders in the infinite top-quark mass limit. We have shown that the corrections at this order are essential and quite remarkable due to the huge reduction of the scale uncertainties, which amount to a factor of four with respect to the known NNLO results. It paves the way for the precision theoretical studies of the Higgs potential at the percent level. Besides the total cross sections, we are also able to predict the various differential distributions at N 3 LO, where an approximation is used in the distributions other than the Higgs pair invariant mass distributions. In general, we have shown very good perturbative convergences in all distributions, and the scale uncertainties are in good control. Besides the SM case, we have also studied the N 3 LO impacts on the (differential) cross sections by varying the trilinear Higgs coupling λ hhh alone. The shapes are again found to be stable at N 3 LO with respect to those at NNLO. Based on these N 3 LO results, we include the important top-quark mass effects at O(α 4 s ) and O(α 5 s ) via three different approximations, where the full m t -dependent NLO calculations are taken from the public code [39,41]. The m t effects are indispensable for the realistic phenomenological applications. We take the (A)N 3 LO⊗NLO mt approximation as our best predictions in this paper. The most advanced FT approximation for the process, requiring the full differential knowledge, will be left for our future studies. The theoretical uncertainties are further improved by the inclusion of both the N 3 LO corrections and the finite m t corrections. The missing m t corrections are larger than the remaining scale uncertainties. Besides, there are several other additional uncertainty sources worthwhile being considered in order to improve the theoretical predictions further. They are the topquark mass scheme dependence [33], electroweak corrections, bottom quark effects and the parametric uncertainties (e.g. m t , α s and PDF).
As a follow-up paper of our previous short letter ref. [47], we have the opportunity to document all the technical details and validation materials here. In particular, we write down the analytic expressions of the one-loop amplitude and the new R 2 Feynman rules in the appendices. The NLO UFO model ready to be used in MG5_aMC is publically available and can be downloaded from http://feynrules.irmp.ucl.ac.be/wiki/HEFT_DH.

A Hard functions
The amplitudes for the Higgs boson pair production in the effective theory, g(p 1 ) + g(p 2 ) → h(p 3 ) + h(p 4 ), can be decomposed into two topologically distinct classes: Class-A with one effective vertex and Class-B with two effective vertices 9 , i.e.
where µ (p 1 ) and µ (p 2 ) are the polarisation vectors of the two initial gluons. The prefactor i v 2 in the above equation is chosen in order to recycle the same notations used in ref. [128]. The amplitudes for Class-A and Class-B can be decomposed into two Lorentz covariant and gauge invariant terms [29] M A/B,µν ab where the tensors are given by with p 2 T = (tû − m 4 h )/ŝ. The Mandelstam variables are defined aŝ For Class-A, after performing renormalisation in the MS scheme, we have where C g is the gluon structure function which has been calculated up to three loops [129,130]. For N 3 LO QCD corrections, we need the two-loop virtual correction to Class-B amplitudes. This was computed in [128], where the finite two-loop four-point amplitudes are obtained by subtracting the IR divergences following the method in ref. [131]. In our framework, a different subtraction method, namely the MS subtraction, is applied, and thus we have reconstructed the full amplitudes with IR poles in Class-B and then performed the renormalisation procedure according to the method in refs. [132,133]. As a result, we obtain the finite part including the real and imaginary contributions are needed in this paper. However, the explicit analytical results can not be found in the literature. In this work we calculated the one-loop amplitudes using FeynArts [134] and FIRE [135] packages, and the results read The dimensionless parameter y is defined as . These analytical expressions have been cross-checked against MadLoop [105,106] and the scale-dependent terms in M B, (2),f in i [128]. The analytic results of the two-loop amplitudes have been obtained in ref. [128] and are expressed in terms of the multiple polylogarithms, which can be evaluated numerically by the public Mathematica package PolyLogTools [136].
The hard functions of class-a, -b and -c are given by where we have averaged over the spins and colours of the two initial gluons and taken into account the symmetry factor 1 2 for the two identical Higgs bosons. Note that the renormalisation is performed at the amplitude level and there is no interference between the two Lorentz structures.

B The NLO model and Feynman rules for the rational R 2 terms
The NLO simulations in the MG5_aMC framework require the derivations of two necessary ingredients from the effective Lagrangian L eff in eq.(2.1) and the SM Lagrangian L SM on top of the information provided in a LO UFO model [112]. They are the UV counterterms to perform the one-loop renormalisation and the rational R 2 terms [137] originating from the integration of the (d − 4) parts of the loop integrands after decomposing their numerators into 4-dimensional and (d−4)-dimensional pieces, where d is the dimension of the spacetime in the dimensional regularisation.
The QCD UV renormalisation counterterms in the theory can be related to the renormalisations of the strong coupling α s and the wavefunctions of gluons and massless quarks. They are however identical to the QCD theory in the SM. Therefore, we will refrain from presenting them in the paper.
Similarly to the UV renormalisation, the computations of R 2 are also equivalent to those of tree-level amplitudes with a universal set of theory-dependent Feynman rules (see refs. [138][139][140][141][142] for QCD and electroweak corrections in the SM and refs. [143,144] for the beyond the SM cases 10 ). They can be derived once and for all (for each model) by just considering the one-particle irreducible one-loop Feynman diagrams. For di-Higgs production in the theory L eff + L SM , we have rederived the analytical expressions of R 2 Feynman rules for zero Higgs and one Higgs vertices by using an in-house Mathematica programme with the aid of FeynRules [145] and FeynArts [134] packages. They have been successfully validated against those in the literature [138,144] with the expressions: We have used the colour factors N c = 3, C F = N 2 c −1 2Nc = 4 3 , the Gell-Mann matrices t a in the fundamental representation of SU(N c ) group, the asymmetric structure constants f a 1 a 2 a 3 of SU(N c ), the colour charge g s = √ 4πα s , the parameter λ HV = 1(0) corresponding to dimensional regularisation (reduction) and the shorthand functions where T a is the colour matrix in the adjoint representation with its elements (T a ) bc = −if abc and the trace of T a into the trace of the Gell-Mann matrices are

C Renormalisation scale dependence
In the framework of SCET, the typical scales are hard, jet and soft scales in addition to a factorisation scale. In order to reproduce the fixed-order results from the resummation formula, all these scales are usually set to be equal (to the factorisation scale), i.e., there is only one scale in the expanded result. Since we want to investigate the scale uncertainties by varying factorisation and renormalisation scales independently, we must reconstruct the individual µ R and µ F dependence separately. In this appendix, we present details about the method we used to obtain the µ R dependence in the expanded results from transverse momentum resummation formula. As a by-product, we find a close relation between the contributions from class-a and class-b. Given that the N k LO cross section is scale (µ R = µ F = µ) invariant, we have The individual renormalisation and factorisation scale dependence is rebuilt from the evolution equation where the first term on the right hand is derived by expanding the transverse momentum resummation formula in the framework of SCET and the second term is given below. Since we use q T -subtraction to calculate the NNLO correction to the class-b diagrams, we focus on the scale dependence in this class. Firstly, we know that the total N 3 LO cross section is independent of the renormalisation scale at each fixed order, i.e., We have only calculated explicitly the results up to O(α 5 s ), so we omit higher-order terms. The first contribution on the right hand is known, The class-c cross section up to NLO QCD is scale invariant, As a consequence, the renormalisation group equation for class-b is derived, . (C.6) The ratio of C hh (µ R ) over C h (µ R ) can be expanded in terms of a s ≡ α s (µ R )/4π, with the coefficient δ 2 = 2 3 (16n f + 35) being scale independent and with β 0 = (11C A − 2n f )/3 and Then, eq. (C.6) turns out to be In the above equation, we have decomposed the class-b cross section as hh is the LO class-b cross section but has a close relation with the class-a cross section; see eq.(C.6). So eq. (C.11) indicates that the class-b cross section has a nonvanishing dependence on µ R only from two-loops. This is actually a consequence of the operator mixing as studied in ref. [146].

D Additional plots
We collect additional plots of 6 differential distributions from LO to AN 3 LO in the infinite top-quark mass limit and with top-quark mass effects in this appendix.   Figure 18. Same as in figure 14 but at √ s = 13 TeV.