Triple Higgs production at hadron colliders at NNLO in QCD

In this work we analyse the triple Higgs boson production cross section at hadron colliders via gluon fusion, and present for the first time the full set of QCD NNLO corrections in the heavy top limit. In order to account for finite top mass effects we perform two different reweighting procedures, and study the dependence of the result on the choice of the approximation. Combining the most accurate predictions available to date, we present the following result for the total NNLO cross section for triple Higgs boson production at a 100 TeV collider: σNNLO=5.56−6%+5%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\sigma}_{\mathrm{NNLO}}={5.56}_{-6\%}^{+5\%} $$\end{document}± 20% fb, where the first uncertainty is an estimate for higher order effects from scale variations, while the last one is an estimate for the missing finite top mass effects.


Introduction
The discovery of the Higgs boson in 2012 [1,2] marked a milestone for the particle physics community. Since then, great efforts have been made to measure its properties in order to determine if it corresponds to the Standard Model (SM) Higgs boson or if it opens a window to new physics beyond it. One of the crucial aspects of this quest is the determination of its self-couplings, which are directly related to the scalar potential that drives electroweak symmetry breaking, and can be modified in the presence of new-physics effects.
In the SM, the triple and quartic self couplings of the Higgs, λ 3 and λ 4 , are uniquely fixed by its mass m H through the enforcement of gauge symmetries and renormalisability of the theory, namely λ 3 = λ 4 = λ SM = m 2 H /(2v 2 ), where v ≈ 246 GeV is the Higgs vacuum expectation value. These couplings are directly accessible through double (for λ 3 ) and triple (for λ 4 ) Higgs boson production, although they might also be extracted indirectly from single Higgs production [3][4][5][6][7][8] and precision electroweak observables [9,10]. All of these measurements, however, are extremely challenging, and in particular already the determination of the triple self-coupling will prove very difficult at the high-luminosity phase of the LHC (see ref. [11] for a review). Unfortunately, the much smaller triple-Higgs production rate makes the determination of λ 4 prohibitive at the LHC, and even challenging at future colliders [12,13]. New physics effects, however, might change the prospects of measuring triple-Higgs production by inducing large enhancements of the cross section, and in this respect it is therefore desirable to provide precise predictions for the SM expectation.
As it occurs for single and double Higgs, the dominant triple Higgs production mechanism at hadron colliders is gluon fusion, mediated by a heavy-quark (mostly top quark) loop. This process, being gluon induced, is expected to present large QCD corrections. However, due to the massive loop appearing in the amplitude at Born level and the large number of external particles, the calculation of its higher order corrections is extremely JHEP03(2020)155 challenging. As a consequence, the exact cross-section for triple Higgs production is only known at leading order (LO) in QCD [14,15], while next-to-leading order (NLO) corrections are currently unknown. In order to provide precise predictions for triple Higgs production, approximate NLO corrections were calculated in ref. [16] within the heavy top limit (HTL) for the virtual amplitudes, while keeping the exact dependence on the top mass m t for the real emission diagrams (this approximation is denoted FT approx. ). In ref. [17] the NNLO virtual amplitudes where computed in the HTL, and the so-called soft-virtual approximation of the HTL NNLO corrections (NNLO SV ) was obtained [18]. Phenomenological results were presented by reweighting the NNLO SV cross-section by the exact LO result (i.e. with full m t dependence) [17], in what is known as the Born-improved (Bi) cross-section. The main goal of the present work is to extend the results of ref. [17] beyond the soft-virtual approximation, computing the complete set of NNLO QCD corrections for triple Higgs production within the HTL. In addition to that, we include partial finite-m t effects by taking into account the NLO FT approx. results of ref. [16], thus providing a final estimate of the cross section that combines the most advanced results available to date.
The paper is organised as follows. In section 2 we examine the structure of the LO triple Higgs production cross-section and its dependence on λ 3 and λ 4 . Then, in section 3 we complete the calculation of the NNLO corrections in the HTL by adding the real emission corrections to the results presented in ref. [17]. Then, in section 4 we present the phenomenological results for triple Higgs production at the LHC and future hadron colliders. We estimate the dependence on the reweighting method by comparing the Born-improved result with a modified prescription (introduced in ref. [19]), that we call dynamically Bornimproved. Finally, we combine our result with the one presented in ref. [16] to report our best prediction, and in section 5 we present our conclusions.

The amplitude at LO
In this section we will examine the structure of the LO amplitude and cross-section. The Born amplitudes needed for the numerical calculation were obtained using Recola2 [20]. For the parton distribution functions we adopted the MMHT2014 [21] set interfaced via LHAPDF [22], while the CUBA [23] library was used to perform the numerical integration. The values implemented for the physical input parameters are G F = 1.16656 × 10 −5 GeV −2 for the Fermi constant, m H = 125 GeV, m t = 173.2 GeV and Γ H = Γ t = 0 for the masses and widths of the Higgs boson and the top quark, respectively, and α S (m Z ) = 0.135 for the strong coupling constant at LO, as provided by the MMHT2014 set. Throughout this work, the on-shell top quark mass scheme is used. All the plots in this section correspond to a collider centre of mass (CM) energy of 100 TeV, although we explicitly checked that all our conclusions also hold at 14 and 27 TeV.
For triple Higgs production, the relevant diagrams (modulo permutations of the final state particules) are shown in figure 1. We can split them in four different categories: pentagons (P), boxes (B) and two triangle contributions (T 1 and T 2 ), each one of these with a specific dependence on the parameters κ 3 = λ 3 /λ SM and κ 4 = λ 4 /λ SM that parametrise JHEP03(2020)155 departures of the self couplings from the SM expectations, As in the case of double Higgs production [24], there are only two independent helicity configurations of the initial gluons, that we call according to the value of total spin along the collision axis. The Spin 2 configuration vanishes in the limit m t → ∞, while the Spin 0 configuration remains. We observe in figure 2 that the contribution from the Spin 2 piece is rather small (below 5% of the total cross-section). Also, we notice that triangle contributions T only contribute to the Spin 0 helicity configuration, and therefore the Spin 2 configuration is not sensitive to the quartic Higgs self coupling κ 4 .
It is interesting to observe the share of the cross-section from each topological contribution and their corresponding interferences. In figure 3 we show different contributions of the Spin 0 component to the invariant mass distribution of the triple Higgs system. In the left panel we plot the interference structure between the P and the B diagrams, similar to the one usually presented for double Higgs production between the box and triangle contributions. This pattern can be better understood if we parametrise the amplitudes in terms of quark loop form factors (factoring out the Higgs boson couplings and propagators)  and look at their behaviour in the HTL: where p 1 and p 2 is the four-momenta of the initial gluons, and p i (i = 3, 4, 5) the ones of the outgoing Higgs bosons, Q 2 = (p 1 +p 2 ) 2 and s ij = (p i +p j ) 2 are the squared invariant masses of the triple Higgs boson system and of the different pairs (ij), respectively, and (ij) is a sum over the three different pairs {(34), (45), (53)}. The dominant contributions come from the pentagon and the box diagrams, as they are less suppressed by Higgs propagators, although the difference in sign between the F P with F B is responsible for a large negative interference between them. This effect is particularly strong at threshold, and therefore this region presents an enhanced sensitivity to any BSM effect that can alter the delicate cancellation between diagrams present in the SM. In the right panel of figure 3 we present, in solid lines, the contribution originated from different powers of κ 4 . The quadratic dependence arises from the T 4 triangle diagram, which is suppressed by a Higgs propagator (∼ Q −2 ) making it negligible, with a contribution at the peak of ∼ 1.5%. The linear term has two contributions, plotted as dashed lines, that cancel each other to a large extent. This cancellation is due to the sign difference  of the form factors in eq. (2.2), making the κ 4 linear contribution about a 15% of the κ 4 -independent one. Because of the different dependence on the self-couplings of the contributions presented in eq. (2.1), a departure from the SM value κ 3,4 = 1 might spoil the destructive interference patterns depicted in figure 3. This can be seen clearly in figure 4, where a large sensitivity to the κ 3 coupling is present particularly around the threshold, where the cross-section can be more than 30 times larger than the SM expectation. Indeed, deviations from κ 3 = 1 spoil the cancellation between the diagrams with most significant contributions, P and B. In the case of κ 4 , because it only affects the 2Re((P + B)T * 4 ) contribution, which is suppressed with respect to |P + B| 2 , the regions of larger sensitivity are those in which the latter is small, namely the production threshold (for κ 3 = 1) and the tail of the invariant mass distribution.
We can also analyse the dependence of the inclusive cross-section on the couplings κ 3,4 , by performing variations with respect to the SM value, as shown in figure 5. For illustrative purposes, both coupling modifiers are varied in the range κ i ∈ [−1; 3]. While the dependence on κ 3 is large (e.g. σ > 8σ SM for κ 3 ∼ −1 and κ 4 = 1) and quadratic JHEP03(2020)155 contributions become noticeable, the dependence on κ 4 is rather small (with departures of |σ/σ SM − 1| < 28% in the range under study), which is compatible with the right panel of figure 3 that shows the |T 4 | 2 contribution to the invariant mass distribution.

NNLO corrections
After discussing the different contributions to the LO cross-section in the previous section, we will present the results for the full NNLO corrections in the HTL.
In the HTL the Higgs bosons couple directly to gluons via the effective Lagrangian where the matching coefficients can be expanded in powers of the strong coupling constant α s as where the expansion is known up to fourth order for X = H, HH, HHH [25][26][27][28][29][30][31]. The computation of the complete NLO corrections and of the virtual amplitudes up to NNLO, both in the HTL, was presented in ref. [17]. The latter were used to construct the soft-virtual approximation (based on the results from ref. [18]) to provide a phenomenological NNLO SV cross-section. In the present work we computed the real emission amplitudes to obtain a full NNLO cross-section in the HTL, including all partonic channels in addition to gluon-fusion. To this end, we exploited the known relation between the single Higgs boson cross-section and some contributions to the multiple Higgs one. This is done in the same fashion as the calculations for double Higgs production [32] and its extension to the dimension 6 Standard Model Effective Theory and Higgs Effective Theory [19].

JHEP03(2020)155
Contributions involving only one HTL operator at the amplitude level (i.e. a single effective vertex between Higgs bosons and gluons) are totally equivalent, apart for an overall normalization and the corresponding matching coefficient, to those from single Higgs production. In terms of the degree of difficulty, these contributions are truly at the NNLO level, but the results can be borrowed from single Higgs production, exploiting the above mentioned similarities. On the other hand, when considering diagrams with more than one HTL operator insertion (actually, their interference with the ones with just one insertion), as each of them carries an extra coupling C X = O(α s ), only tree level configurations can appear to NLO, while at NNLO accuracy the only possible contributions arise from one-loop and single real emission diagrams. Because of this simplification, their infrared divergences can actually be handled with standard NLO procedures.
In our calculation, the FeynArts [33] and FeynCalc [34,35] packages were used in Mathematica to compute the amplitudes. The cancellation of their infrared singularities with the ones present in the virtual amplitudes presented in ref. [17] and the absorption of the remaining ones into the evolution of the parton distribution functions was performed following the FKS [36,37] approach. In this work we neglect effects coming from the Higgs width, and set it equal to zero throughout the calculation. We present below the final result with a notation suitable for the discussion on the reweighting following in section 3.1. More details on the derivation can be found in the appendix.
As usual, the cross-section can be written as where Q is the invariant mass of the triple Higgs system, √ s H the collider centre of mass energy, and i, j are the labels for the massless partons inside the hadrons h 1 and h 2 with respective parton density f (i,j)/(h 1 ,h 2 ) . Here the dependence on the factorisation and the renormalisation scales is implicitly understood. The partonic cross-sectionσ is computed order by order as an expansion in the strong coupling α S , such that up to NNLO we write it as where the LO amplitude M 3H can generically be written as and The perturbative coefficients η (0,1,2) ij , within the HTL approximation and its extensions (reweightings) described in section 3.1, are up to NNLO where η H,(n) ij , n = 0, 1, 2, are the corresponding QCD corrections for single Higgs production that can be found in ref. [38] (and coincide with the results in refs. [39,40]), the renormalisation and factorisation scales were set to µ F = µ R = Q, k and l label the final state Higgs bosons, and (kl) denotes the sum over distinct pairs of them. C H LO and C 2H LO are (up to a normalisation factor) the LO amplitudes for single and double Higgs production 11) and the coefficients R 3H and V 3H are the finite remainders of the virtual corrections at NNLO presented in ref. [17], with the only difference being that we use (as we will discuss later) the general definition of eq. (2.2) for C 2H LO . If we express this coefficient in the HTL, the expressions are identical to those in ref. [17].
The only missing ingredient, ρ ij , corresponds to the finite remainder of the real emission corrections to the diagrams with more than one HTL operator insertion (which are already included in η H(2) ij ). These only contain diagrams with a single parton emission whose divergences were regulated using dimensional regularisation in D = 4 − 2 dimensions in the FKS [36,37] framework. After subtracting the singularities, we can write the remainder as that contains a soft-collinear term ρ (sc) ij and a regular term ρ (r) ij whose explicit expressions, together with the definition of the standard plus-prescription, are given in the appendix.

Reweighting
The results presented in this section complete the full NNLO corrections to triple Higgs production in the HTL when the quantities F P , F B and F T are expressed in this limit (see eq. (2.2)). In order to retain part of the m t dependence, this expression is usually reweighted by using the exact result for M 3H in eq. (3.4), while keeping the coefficients η (0,1,2) ij in the HTL. This procedure is usually referred to as the Born-improved (Bi) NNLO cross-section, and its accuracy has been studied for double Higgs production [41][42][43] at NLO finding that it overestimates the exact inclusive cross-section by a 32% at collider energies of 100 TeV (and around 16% at 14 TeV). This overestimation is enhanced in the tail of the Higgs pair invariant mass distribution.
In order to parametrise the dependence of the result on the reweighting procedure, we also considered what we call a dynamically-Born-improved (dBi) NNLO cross-section, which we obtain by using the full dependence on the kinematics of the outgoing particles of F T , F B and F P in the definitions of P, B, C H LO and C 2H LO , and therefore also in the definitions of R 3H and ρ (2) ij . In this way, we reweight the HTL insertion operators with the respective LO amplitude diagram by diagram (e.g. the (H)Hgg vertex with the (double) single Higgs production amplitude ∼ C (2)H LO ). This reweight cannot be applied in a straightforward way, as the form factor F B is not always defined for the kinematics characterising the HTL vertex. This happens in diagrams with two HTL vertices connected via an off-shell gluon (e.g. a box and a triangle loop). To fix this problem, we modify the kinematics in a way such that we preserve the momenta of the outgoing Higgs bosons, while redefining the momenta of the gluons entering the vertex to be on-shell.
For the triangle form factor F T , we evaluate it at the invariant mass of the outgoing Higgs boson, just as expressed in the eqs. (3.10) and (3.11). For the box form factor F B appearing in eq. (3.11), we need to define a prescription that corrects the momenta of the initial gluons to compensate for the recoil of all other particles not involved in the HTL vertex. Lets recall we are labelling p 1 and p 2 as the momenta of the incoming partons, and p (3,4,5) the momenta of the outgoing Higgs Bosons. For the vertex with outgoing Higgses {i, j}, we define q := p i + p j , M 2 := q 2 , and q µ T as the transverse component of q with respect to p 1 and p 2 , and then define the momenta of the gluons entering the form factor as k 1 and k 2 := q − k 1 in the following way 14) where the different choices of ξ ∈ [0, 1/2] define different consistent prescriptions to account for the Higgs pair recoil. This prescription corresponds to the one presented in ref. [44] in the context of transverse-momentum resummation, if one chooses k 1T k 1T k 1T = ξq T q T q T for the eqs. (25) and (26) therein. In particular, if ξ = 0 the transverse recoil is compensated by JHEP03(2020)155 ). Let us emphasise that this prescription leaves unchanged the momenta of all the Higgs bosons involved, including their virtuality, and its main purpose is to redefine the momenta of the virtual gluons entering Box diagrams, so that these are on-shell and F B is well defined.
With this prescription, the dynamically Born improved approximation is well defined and we will denote it as dBi ξ . In particular, we will present results for dBi 0 and dBi 1/2 as benchmarks to show the dependence on the choice of ξ, and we will show that this dependence is numerically negligible at NLO and NNLO.

Phenomenological results
We present results for the NNLO cross-section for triple Higgs production, reweighted using the Bi, dBi 0 and dBi 1/2 prescriptions. The set-up is the same as the one used in section 2. We use the MMHT2014 [21] set of parton distributions, at the corresponding order in the strong coupling constant for each contribution (LO, NLO, NNLO). For phenomenological purposes, we compute this for collider energies of 100 TeV and 27 TeV that are relevant for physics at the Future Circular Collider and High-Energy LHC, respectively. To estimate the theoretical uncertainty arising from the missing higher orders in the perturbative series, we perform an independent variation of the factorisation and renormalisation scales in the range [µ 0 /2, 2µ 0 ], with the constrain 0.5 < µ R /µ F < 2. The choice of the central values µ 0 used in this work were µ 0 = Q/2 and µ 0 = Q, where Q is the invariant mass of the triple Higgs system.
In figure 6 we see the cross-section computed in the dinamically Born improved approximation up to different orders, as well as the K factor defined as usual, K = dσ / dσ LO . As seen also within the soft-virtual approximation [17], the cross-section begins to stabilise only from NNLO. The K factors are rather flat at the peak of the invariant mass distribution, with values around 1.7 and 1.8 for collider CM energies of 100 and 27 TeV respectively, while the NNLO K factors present a suppression in the tail. Due to this suppression, the entire NNLO band falls inside the scale variation of the NLO, suggesting that the perturbative series is more stable in this region. The total scale uncertainty is reduced from 37% to 27% and to 11% when going from LO to NLO to NNLO, at 100 TeV with a central scale choice of µ 0 = Q/2. For 27 TeV the reduction is similar, from 48% to 30% and then to 12%.
To measure the effect of the reweighting, we can compare the different approximations Bi, dBi 0 and dBi 1/2 and observe the corresponding effect on the invariant mass distribution. In figure 7 we see that, although at NLO the three approximations are almost completely compatible, at NNLO there is a significant decrease in the tail of the distribution when using a dBi ξ instead of the Bi approximation, a discrepancy that is even bigger than the scale variation in this region. We also notice that the dependence on ξ of the dBi ξ approximation is phenomenologically negligible. We know that for double Higgs production the Bi approximation overestimates the tail of the distribution at NLO respect to the exact calculation [41][42][43]. If this effect holds also for triple Higgs production, we can expect the JHEP03(2020)155  dBi approximation to provide more reliable predictions, as it predicts a smaller tail of the distribution just by reweighting each contribution by the associated form factor instead of using the full amplitude.
In order to understand why the discrepancies in the tail of the invariant mass distribution between the Bi and dBi prescriptions arise only at NNLO, lets recall what are the main differences between the two reweighting procedures. The Bi reweights all amplitudes by the Born amplitude, including those that have more than one HTL vertex. The dBi only does this to amplitudes containing a single HTL vertex, while applying a different prescription for those amplitudes with more than one HTL vertex. In this way, the Bi JHEP03(2020)155 reweighting procedure increases the relative significance of the amplitudes with many HTL vertices, respect to the dBi. At NLO such amplitudes appear only at tree level, due to the power counting of the HTL vertices, making the discrepancy phenomenologically negligible. At NNLO, such diagrams are enhanced by real emission corrections, making the discrepancy at the tail of the invariant mass distribution more noticeable.
In table 1 we present the results for the inclusive cross-section obtained for different collider energies, choices of central scales µ 0 , reweighting procedures and orders in the perturbative expansion. The corresponding K factors for the dBi results are presented in table 2. When comparing our results with the ones obtained in ref. [17] in the soft-virtual approximation, we find that although the SV result differs only in about 1% from the full NNLO in the HTL, when using the Bi reweight the difference grows to a 2.5% and 4.4% increase in the inclusive cross-section at 14 TeV and 100 TeV, respectively. This larger difference is due to the fact that in the HTL the SV approximation is slightly smaller than the complete NNLO for small invariant masses, but compensates at large invariant masses, resulting in a small difference in the inclusive cross section. After the reweighting procedure, the region of large invariant masses is suppressed, and therefore this accidental JHEP03(2020)155  Table 2. NLO and NNLO K factors for the inclusive triple Higgs boson production at different collider energies. These defined as the quotient between the dBi and the LO cross section. The results are shown for central scale values of Q (top) and Q/2 (bottom). The dependence on ξ in the dBi ξ reweight procedure is below the per-mill level and therefore omitted. The last row shows our best available prediction for the different collider energies.
compensation is reduced, increasing the difference in the inclusive cross sections between SV and complete NNLO. The different reweighting procedures produce a small difference ranging from ∼ 0.7% at 14 TeV up to ∼ 1.3% at 100 TeV. Although in figure 7 we see that the different reweights lead to discrepancies larger than the theoretical uncertainties in the tail of the distribution, since this region has a small impact on the inclusive cross-section the results for the total cross-section are consistent within the theoretical uncertainties. Of course, this small difference at the total cross section level can only be a lower bound on the expected finite top mass effects, and from the results obtained at NLO within the FT approx. (which are ∼ 10% smaller than the dBi prediction) it is clear that they are expected to be much larger. What we can conclude from this exercise therefore, is that the systematic uncertainties related to the choice of the reweighting procedure (among the choices presented here) is expected to be marginal compared to the full size of the finite-m t effects.
In order to provide the best possible estimate of the triple-Higgs production cross section, it becomes necessary to include the partial finite-m t effects obtained in ref. [16] within the FT approx . To this end, we use the predictions presented therein and in ref. [45] for the total cross section, and encode the finite mass effects in the parameter δ t defined by and we define our best prediction as This procedure is similar to the prescription that was implemented in ref. [45] for double-Higgs production. The values that we obtain for δ t at the different collider energies are δ t = −0.107, −0.073 and −0.146 for 14, 27 and 100 TeV, respectively. 1 The corresponding cross sections and K factors are presented in tables 1 and 2, respectively.

JHEP03(2020)155
Before providing our final results, we address the issue of the remaining uncertainties associated to finite-m t effects. There is in principle no reason to expect these effects to be smaller than the corresponding ones present in double Higgs production; in fact the typically larger invariant masses involved might point to the opposite direction. Of course, the more complicated structure of the triple Higgs production amplitudes, involving additional topologies, might lead to accidental cancellations of these effects that cannot be predicted at this point. In the absence of any prediction with full top mass dependence beyond the LO, one can only rely on approximated results in order to estimate the associated uncertainties at NNLO. The best available approximation at NLO, the FT approx. , differs from the Born-improved NLO at O(10%) at 14 and 27 TeV, and O(15%) at 100 TeV. The size of this difference can be a good estimation of the missing finite top mass effects, and indeed this is the case for the double-Higgs production cross section. In order to provide a conservative estimation, and having in mind the possibly worse situation in triple Higgs as compared to Higgs pair production, we estimate the uncertainty of our prediction to be of ±15% at 14 and 27 TeV, and ±20% at 100 TeV.
Given that the full dependence on the top mass is only retained at LO, it is not possible to perform a complete analysis on the top mass scheme uncertainty (which was found to be large at NLO in the case of double Higgs [43]). Nevertheless, a parametric variation of the default on-shell value m t = 173.2 GeV used along this work to the correspondent one in the MS scheme, m t (m t ) = 163.6 GeV (using a three-loop conversion between schemes), shows a decrease of about 25% in the cross section, which indicates an uncertainty in line with the one estimated for the finite top mass effects.
Compiling all the ingredients described in the last paragraphs, we arrive therefore to the following final prediction for the triple Higgs production cross section:

Conclusions
In this work we presented for the first time the complete set of NNLO corrections to triple Higgs boson production at hadron colliders in the heavy top limit. To partially retain finite top mass effects, two different reweighting procedures have been implemented: the usual Born-improved approximation (Bi), and a new procedure that we call dynamically Born improved approximation (dBi). Both procedures coincide at LO, and from their difference at higher orders we infer the dependence of the result upon the reweighting procedure. Overall, we found that the invariant mass distribution is sensitive to the reweighting procedure only in the tail, where the cross-section is already small, while for the inclusive cross-section the dependence on this procedure is O(1%), which falls inside the scale variation uncertainties of O(7%).
In order to provide a prediction based on the most advanced results available in the literature, we combined our dBi-reweighted NNLO results with the NLO predictions obtained JHEP03(2020)155 within the FT approx . From the differences between the available NLO approximations we estimated the size of the missing finite top mass effects. Based on this, our final prediction for the triple Higgs production cross section at a 100 TeV collider is σ NNLO = 5.56 +5% −6% ±20% fb.

Acknowledgments
This work is supported in part by the Swiss National Science Foundation (SNF) under contracts 200020 169041 and IZSAZ2 173357, by MINCyT under contract SUIZ/17/05, by Conicet and by ANPCyT. We thank Jean-Nicolas Lang for his support with the Recola2 library.

A Real corrections
In this appendix we present the results for the corrections arising from real emission in diagrams with more than one HTL insertion operator. Because each HTL operator insertion carries a α s factor, tree level diagrams with two operator insertions are already NLO, so their real radiation corrections correspond to a single parton emission. In the same way, diagrams with three operator insertions appear only at NNLO and their real emission corrections are of higher order, and therefore not considered. The single real emission amplitudes present divergences when the emitted parton becomes unresolved, some of which will cancel against the divergences present in the corresponding loop corrections calculated in ref. [17] and the rest have to be absorbed in the NLO evolution of the parton distribution functions. In order to perform such cancellations, we used an FKS approach in D = 4 − 2 dimensions.
The key idea of the FKS method is to divide the phase space of the real corrections dPS 4 into soft, collinear and regular regions. To do so, we express it in terms of the LO phase space dPS 3 plus the dependence on the emitted particle where M r is the amplitude for the real emission process we are considering, s is the invariant mass of the incoming partons, dPS is the Born phase space evaluated at s → xs and Q 2 is the invariant mass of the triple Higgs system such as x = Q 2 /s. The emitted parton therefore has a momentum fraction of (1 − x)s, and we define its orientation with respect to one of the initial partons with the azimuth φ and the cosine of the polar angle y.