Two-loop corrections to the triple Higgs boson production cross section

In this paper we compute the QCD corrections for the triple Higgs boson production cross section via gluon fusion, within the heavy-top approximation. We present, for the first time, analytical results for the next-to-leading order corrections, and also compute the soft and virtual contributions of the next-to-next-to-leading order cross section. We provide predictions for the total cross section and the triple Higgs invariant mass distribution. We find that the QCD corrections are large at both perturbative orders, and that the scale uncertainty is substantially reduced when the second order perturbative corrections are included.


Introduction
After the discovery of a Higgs boson [1,2] at the Large Hadron Collider (LHC), it is crucial to study its properties in order to determine whether it is indeed the particle predicted by the Standard Model (SM) or not. Besides its couplings to fermions and gauge bosons, which are so far compatible with the SM expectations, it is of great interest to determine the Higgs boson self-couplings, which will allow to shed light on the scalar potential, and therefore the electroweak symmetry breaking mechanism.
The Higgs boson trilinear and quartic self-couplings λ 3 and λ 4 can be studied in hadron colliders via double and triple Higgs production, respectively (see Ref. [3] for an alternative method based on single Higgs production). The SM expectations for these processes, corresponding to λ 3 = λ 4 = m 2 H /(2v 2 ), being v 246 GeV the Higgs vacuum expectation value and m H its mass, are very low. For a collider energy of 14 TeV, the leading order (LO) predictions for the double and triple Higgs production cross sections are of O(20 fb) and O(0.05 fb). As a consequence, in a SM-like scenario, a determination of the trilinear coupling will be very challenging at the LHC, while the measurement of the quartic coupling via triple Higgs boson production will be at best relegated to a future collider [4,5]. Of course, the situation can be largely modified in the presence of new physics scenarios for the Higgs sector.
As it also happens for single and double Higgs production, the triple-Higgs final state is mainly produced in the SM via gluon fusion, mediated by a heavy quark (mainly top) loop. For this production mechanism, the corresponding cross section is only known at LO in perturbation theory. However, the QCD corrections are expected to be large, as it also happens for the other gluon initiated loop-induced processes mentioned above. Unfortunately, their computation is very difficult. For instance, the next-to-leading order (NLO) corrections for double Higgs production (a simpler process with one particle less in the final state) were not known until very recently [6]. In the absence of a full NLO calculation, and in order to provide an estimate of the size of the perturbative corrections, approximate NLO predictions were obtained in Ref. [7], where only the exact real corrections were included.
In this paper, we present the first calculation of the QCD perturbative corrections for triple Higgs production within the Higgs effective field theory (HEFT). Within this framework, which formally corresponds to the large top quark mass limit of the SM, the Higgs bosons couple directly to gluons via an effective Lagrangian. This approach has been successfully used to compute the QCD corrections for single and double Higgs production. Motivated by this, we apply it to compute the NLO corrections and the next-to-next-to-leading order (NNLO) soft and virtual contributions for the total cross section and the triple Higgs system invariant mass distribution.
This work is organized as follows: in Section 2 we present the virtual corrections up to NNLO, and later in Section 3, after combining with the corresponding real corrections, we present the NLO and NNLO partonic cross sections, the latter within the soft-virtual approximation. In Section 4 we present the numerical results for the LHC and future colliders, and compare our predictions with the other NLO approximation available. Finally, in Section 5 we present our conclusions.

Virtual corrections up to NNLO
In this section we present the one and two-loop corrections to the triple-Higgs boson production cross section in hadronic collisions via gluon fusion. As was stated before, we work within the HEFT were the Higgs bosons couple directly to gluons via the effective Lagrangian and where the matching coefficients can be expanded in powers of the strong coupling α S as The three coefficients needed for our calculation are known up to fourth order in their perturbative expansion [8][9][10][11][12][13][14].
For the generation of the relevant Feynman diagrams we employed qgraf [15], while the manipulation of the resulting amplitudes was performed with in-house routines written for Mathematica. Finally, we reduced the result into master integrals using the algorithm Fire [16]. The infrared divergent results were handled using dimensional regularization with D = 4 − 2 dimensions.
The virtual corrections to the partonic cross section can be written in terms of the squared matrix element asσ v = 1 2s where we include the flux factor, the average over initial state colors and helicities, and the factor 1/3! arising from the identical particles in the final state. Here dPS 3 represents the three particle phase space. Expanding in powers of the strong coupling, we have Exploiting the well known one and two-loop infrared behaviour of QCD amplitudes [17][18][19], we can write the renormalized NLO and NNLO virtual corrections as where I (1) and I (2) represent the one and two-loop insertion operators defined, for instance, in Ref. [17].
The D dimensional LO cross section can be written as and where the coefficient C 3H LO is defined as Here s ij...k = (p i + p j + · · · + p k ) 2 . For simplicity, we have set Γ H = 0 (the numerical effect due to the Higgs width is negligible), in which case the C 3H LO coefficient is a real number. The one and two-loop infrared-regulated parts can be organized in the following way: The contributions labelled F arise from diagrams with only one HEFT operator insertion. The ones in R originate from the interference between amplitudes with two HEFT operator insertions and amplitudes with only one insertion. On the other hand, contributions in T arise from the interference between diagrams with three operator insertions and the LO, and the ones in V come from the square of amplitudes with two insertions. Finally, contributions to S have their origin on the difference between the NNLO QCD corrections to the effective vertices Hgg, HHgg and HHHgg. In Figure 1 we show illustrative examples of the different Feynman diagrams involved in the calculation of the virtual corrections. As already mentioned, since we adopted the Higgs zero-width approximation, both C 3H LO and C 2H LO are real numbers. Beyond that limit, there is also a numerically negligible contribution proportional to Im(C 2H LO )Im(C 3H LO ), which we will ignore. We start by presenting the one-loop corrections. For simplicity, we set µ R = µ F = Q through-out the rest of the work, being Q the invariant mass of the triple-Higgs system. We find that where N f represents the number of light partons, ζ n stands for the Riemann zeta function, and with t ij = (p i − p j ) 2 , and where we have defined Notice that, given that these contributions enter in the two-loop result multiplied by a double pole, we need their expansion up to O( 2 ).
We start now with the NNLO results. For the coefficient F we find where m t stands for the top quark mass. The function R 3H can be written at this order as where For S where the NNLO corrections to the effective vertices between Higgs bosons and gluons imply [8][9][10][11][12][13][14] For the function T 3H we write where we have defined where the last term indicates that the momenta p 1 and p 2 have to be exchanged. Finally, the function V 3H can be expressed as with the following definitions, and Here the Lorentz invariants are defined as with p i+j = p i + p j .
With the above results we complete the presentation of the NNLO virtual corrections for the triple Higgs boson production cross section. It is worth to mention that some of the expressions obtained can be directly related to their analogous in double Higgs production. This is of course the case of the F contributions, which arise directly from the gluon form factor, and take exactly the same value for both processes. Less trivially, the functions r (1) and r (2) defined in Eqs. (10) and (14) are equal to the coefficients R (1) and R (2) defined in Eqs. (11) and (15) of Ref. [20], provided that the limit m 1 = m 2 = m H is taken in the former.

NLO and NNLO SV partonic cross sections
We present here the partonic cross sectionσ at NLO and NNLO, obtained by combining the results from the previous section with the corresponding real corrections (computed within the soft approximation for the NNLO case). We recall that the inclusive hadronic cross section can be obtained from the partonic result as where √ s H represents the collider center-of-mass energy. The parton densities are denoted by f i/h (x) and the subscripts i, j label the type of massless partons. For simplicity, the dependence on the factorization and renormalization scales is always understood.
The partonic cross sectionσ can be expanded in powers of the strong coupling α S . Up to NNLO, we have where the LO cross section in the HEFT takes the form and where C 3H LO is defined in Eq. (7). At LO only the gluon initiated subprocess contributes, and therefore we have where x = Q 2 /s, being √ s the partonic center-of-mass energy.
For the calculation of the NLO results, we exploited the relation that can be established between some contributions to the triple-Higgs boson cross section and the single-Higgs boson one, as was already discussed for the double-Higgs case in Ref. [21]. We find the following results, where the plus distributions D i (x) are defined as usual, and the coefficient R 3H is defined in Eq. (9). Since we are dealing in this section with finite quantities, the = 0 limit can be taken for this coefficient.
The results above complete the NLO calculation within the HEFT which, to the best of our knowledge, has not been presented before. For the NNLO SV cross section, we make use of the universal formula derived in Ref. [22]. The soft and virtual contributions are those proportional to the delta function δ(1 − x) and the plus distributions D i (x), which in Mellin space correspond to constants and threshold enhanced logarithms. These contributions only appear in the gluon initiated partonic channel, and they can be expressed as η (2) gg(SV) = δ(1 − x) 11 18 where the finite reminders of the one and two-loop virtual corrections dσ fin are defined in Eq. (8).
As it was already observed in Refs. [22,23], the SV approximation yields better results when defined in Mellin space. To this end, we need to compute the following N -moments, η (2) gg(SV),N = η , where we additionally take the large-N limit on the resulting expression for the Mellin transform η (2) gg(SV),N , retaining only the logarithmically enhanced and constant terms. Its explicit expression can be easily derived from the results in Ref. [22]. Finally, the NNLO contribution to the physical cross section in the SV approximation can be obtained by Mellin inversion, where the constant C M P defining the contour of integration is on the right of all the singularities of the integrand, as defined in the Minimal Prescription in Ref. [24].

Phenomenological results
We present in this section our numerical predictions for the LHC and future hadron colliders. The NNLO corrections are computed within the soft-virtual (SV) approximation, which has proven to be an excellent estimation to the full HEFT result for other gluon-initiated processes, like single and double Higgs boson production [20,25]. In particular, for the triple Higgs production cross section, we find that at NLO the SV approximation differs from the full HEFT result by less than 2.5%.
In order to partially retain the dependence on the top quark mass, we normalize the QCD corrections computed in the HEFT with the exact LO result, differentially in the triple Higgs system invariant mass. The full (loop-induced) LO cross section was computed using Mad-Graph5 aMC@NLO [26,27].
For the numerical implementation we set the values m H = 125 GeV and m t = 173.2 GeV for the Higgs boson and top quark masses. We use the MMHT2014 sets [28] for the parton flux and strong coupling, at the corresponding order for the LO, NLO and NNLO predictions. For the renormalization and factorization scales we use two different central scale choices, µ R = µ F = µ 0 with µ 0 = Q and µ 0 = Q/2. As usual, the theoretical uncertainty arising from the missing higher orders is estimated by varying these scales independently in the range [µ 0 /2; 2µ 0 ], with the constraint 0.5 < µ R /µ F < 2.
In Figure 2 we show the triple Higgs system invariant mass distribution for a collider energy of 14 TeV, for the two central scale choices. As can be seen in the lower panel, for both choices the NLO corrections turn out to be large, with almost flat K-factors which approximately take the values 1.8 and 1.6 for µ 0 = Q and µ 0 = Q/2, respectively. The relative scale uncertainty is reduced at NLO, but still remains rather large, ∼ 33% for both scales. It is worth to notice that the scale variation at LO fails to anticipate the size of the higher order corrections, as it occurs for single and double Higgs production as well.
The NNLO corrections are still very sizeable, specially for µ 0 = Q, where they represent an increase of about 28% with respect to the NLO. Corrections are more moderate for µ 0 = Q/2, increasing the NLO result by about 13%, and with a large overlap between the corresponding uncertainty bands. The K-factors have a mild dependence on the invariant mass of the system, showing a small increase towards larger values of Q. For both scale choices, the total scale uncertainty is substantially reduced at NNLO, being of about 15% and 12% for µ 0 = Q and µ 0 = Q/2, respectively. It is worth pointing out that at NNLO both scale choices give fully compatible results, with a difference between the central values smaller than 4%. In Figure 3 we show the triple Higgs production cross section at the different accuracy levels as a function of the collider energy. Also, we present in Table 1 the total cross section for E cm = 13, 14, and 100 TeV. We can observe that both NLO and NNLO corrections are very sizeable in the whole range under study, taking lower values for higher energies. Once again, the overlap between the NLO and NNLO predictions is larger for µ 0 = Q/2. Figure 3, we can also observe that the scale uncertainty is substantially reduced once the NNLO corrections are included, independently from the collider energy under consideration. Specifically, at 13 TeV and for µ 0 = Q/2 the total uncertainty goes from 59% to 33% and 12% when going from LO to NLO and NNLO. The analogue uncertainties at 100 TeV are 37%, 26% and 14%, where we can also observe this reduction. Similar results are obtained with µ 0 = Q.

From the results in
For completeness, we show in Figure 4 the dependence of the NNLO total cross section on the value of the Higgs self-couplings. We do not intend to perform a full analysis for non-SM scenarios, but just to illustrate the sensitivity of this observable to λ 3 and λ 4 arround their SM expectations (in particular, in the range λ/λ SM ∈ [0; 2]), which is of course relevant for a future measurement. The results correspond to a center of mass energy of 100 TeV, and the scale choice µ R = µ F = Q/2. We can observe that, unfortunately, the dependence on λ 4 is rather small, and that the situation slightly improves for λ 3 > λ SM , while the sensitivity is even lower for λ 3 < λ SM . It is worth to mention that the dependence of the ratio between the NNLO and LO cross sections on the value of the self-couplings is quite small, finding almost flat K-factors in the whole range under analysis.   Finally, we want to evaluate the applicability of the HEFT for the computation of the QCD corrections for triple Higgs production. Formally corresponding to the large top quark mass limit, this approximation completely fails to reproduce the known LO result. However, it should be more reliable for the computation of the radiative corrections, once the exact LO is used to normalize the latter. In fact, for the double Higgs production cross section (where the main contribution comes from the region in which the invariant mass of the di-Higgs final state is larger than the threshold 2m t ) it has been shown that the NLO HEFT prediction overestimates the full NLO result by a 14% and 24% for E cm = 14 and 100 TeV [6,29] (and with deviations of the same order for the shape of the invariant mass distribution).
Of course, for triple Higgs production the exact NLO result is not available. However, in order to estimate the level of accuracy of the HEFT we can rely on the approximate NLO results presented in Ref. [7], where the exact real corrections were included via a reweighting technique. This approximation, for instance, improves the HEFT result for double Higgs production, overestimating the 14 and 100 TeV total cross sections by only 4% and 6% respectively.
Using the same setup of Ref. [7] (in particular the same PDF sets and strong coupling), we find that the NLO HEFT result overestimates their approximate NLO prediction by about 7% and 9% for E cm = 14 and 100 TeV respectively. This relatively small deviation, combined with the good level of accuracy shown by the approximation of Ref. [7] for the di-Higgs production cross section, indicates that, within its expected limitations, the (Born-normalized) HEFT can be used to gauge the size of the QCD higher order corrections for triple Higgs boson production. At NLO, we estimate the uncertainty of the HEFT approach to be of O(20%).