Higgs boson pair production via gluon fusion at N 3 LO in QCD

Long-Bin Chen, ∗ Hai Tao Li, † Hua-Sheng Shao, ‡ and Jian Wang § School of Physics and Electronic Engineering, Guangzhou University, Guangzhou 510006, China Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, 87545, USA Laboratoire de Physique Théorique et Hautes Energies (LPTHE), UMR 7589, Sorbonne Université et CNRS, 4 place Jussieu, 75252 Paris Cedex 05, France School of Physics, Shandong University, Jinan, Shandong 250100, China We present next-to-next-to-next-to-leading order (NLO) QCD predictions for the Higgs boson pair production via gluon fusion at hadron colliders in the infinite top-quark mass limit. Besides the inclusive total cross sections at various collision energies, we also provide the invariant mass distribution of the Higgs boson pair. Our results show that the NLO QCD corrections enhance the next-to-next-to-leading order cross section by 3.0% (2.7%) at √ s = 13 (100) TeV, while the scale uncertainty is reduced substantially below 3% (2%). We also find that a judicious scale choice can significantly improve the perturbative convergence. For the invariant mass distribution, our calculation demonstrates that the NLO corrections obviously improve the scale dependence but almost do not change the shape.

Introduction -The discovery of the Higgs boson [1,2] marks the completion of the standard model (SM) of particle physics and the start of a new era for the physics searches at the LHC. The next primary goal of the LHC is to precisely pin down its interactions with other SM particles or itself. In particular, the precision study of the Higgs potential is ultimately crucial for understanding the electroweak symmetry breaking mechanism. At the LHC, it has been found that its interaction couplings with massive gauge bosons and fermions agree with the SM expectations [3][4][5][6], while there are only quite weak constraints on the trilinear Higgs self-coupling [7,8]. However, in the future, the experimental probe will be significantly improved as the increase of the integrated luminosity and the collision energy, and/or by employing novel analyzing methods [9]. Theoretically, there indeed exist beyond-the-SM (BSM) models in which the trilinear Higgs self-coupling deviates from the SM value by about 100% while the Higgs couplings to gauge bosons and fermions are almost SM-like [10]. Therefore, the precise measurement of the trilinear Higgs self-coupling at the LHC would be of paramount importance to fully explore the elusive BSM signals.
The direct manner to probe the trilinear Higgs selfcoupling at the LHC is via the Higgs boson pair production. Like the single Higgs case, the gluon-gluon fusion (ggF) channel is dominant, while other channels like vector-boson fusion (VBF) are at least one order of magnitude lower in their yields [11]. Similarly to the cross section of single Higgs boson production, the ggF di-Higgs cross section is plagued with large theoretical uncertainties, dominated by the QCD scale uncertainty [12] and the top-quark mass scheme dependence [13,14]. The computations of the cross section have been carried out both in the infinite top-quark mass limit and with full top-quark mass dependence.
On the other hand, there are also many attempts to go beyond the infinite top-quark mass approximation. The full top-quark mass dependence was first included in the real-emission part at NLO [11,22]. The NLO virtual corrections, involving multi-scale two-loop integrals since the LO is already a loop-induced process, have been evaluated by expansion in the heavy top-quark mass limit up to O(1/m 12 t ) [23][24][25], in the small top-quark mass limit [26,27], and in terms of a small Higgs transverse momentum [28] or a small Higgs mass [29]. Recently, the expansion of the three-loop virtual corrections in the heavy top-quark mass limit has been presented [30]. Finally, the full NLO QCD corrections including exact dependence on the top-quark mass were computed numerically by two groups [14,[31][32][33], either by using a quasi-Monte Carlo method [34,35] or via a direct Monte Carlo integration by Vegas. The matching to parton showers has also been carried out [36][37][38].
Although the infinite top-quark mass approximation is usually insufficient for the corresponding phenomenology studies, a standard way of improving the theoretical prediction on the ggF di-Higgs cross section is to use the lower-order result with full top-quark mass dependence, and to augment it with higher-order corrections in the infinite top-quark mass limit [39].
In this Letter, we provide the first next-to-next-tonext-to-leading order (N 3 LO) perturbative QCD predictions for the Higgs boson pair production via gluon fusion at hadron colliders in the infinite top-quark mass limit. This result becomes one of a few highest-precision computations for scattering processes relevant at the LHC. The existing calculations performed at N 3 LO include the inclusive cross sections of the ggF [40,41], VBF [42] and bottom-quark fusion [43] of single Higgs boson produc-tion, as well as the VBF of di-Higgs production [44]. Some differential distributions approximated at the same order are also known for the ggF of single Higgs boson production [45][46][47]. In our calculation, we will provide both the inclusive cross sections and the invariant-mass distribution of the Higgs boson pair at N 3 LO, where the latter is the first exact N 3 LO differential distribution for the ggF channel.
Theoretical framework -In the infinite top-quark mass limit, the effective Lagrangian describing the coupling of two gluon field strength tensors with one or two Higgs bosons reads where the vacuum expectation value of the Higgs field v can be related to the Fermi constant by v = √ 2G F −1/2 .
C h and C hh are the Wilson coefficients by matching the full theory to the effective theory, which start from O(α s ) and have been calculated up to O(α 4 s ) [48][49][50][51][52][53][54][55][56]. The ggF di-Higgs cross section can be organized according to the number of the effective vertex insertions in the amplitude squared, where three representative Born cut-diagrams are shown in Fig. 1. They are the ones with two, three and four effective vertices, and are denoted as class-a, -b and -c respectively in the following context. Accordingly, the (differential) cross section can be splitted into three parts, Their contributions to various α s orders are tabulated in Table 1. At N 3 LO in α s , we need to calculate the class-a (class-b and class-c) contribution at N 3 LO (NNLO and NLO). Because of the similar topologies, the class-a part can be obtained from the calculation of a single Higgs boson production. They can be related by where λ is the Higgs self-coupling and the function f h→hh accounts for the phase space difference between the single and double Higgs boson production, TABLE I: 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 of order counting applies to class-c. and σ h (m h → m hh ) denotes the cross section calculated using iHixs2 [57] after replacing the Higgs boson mass with the invariant mass of the Higgs boson pair in the code. Such a method has also been used in the earlier NNLO calculation of the ggF di-Higgs production in Ref. [16].
The class-b part can be obtained through the q Tsubtraction method [58], in which we divide the cross section into two parts, where p hh T represents the transverse momentum of the Higgs pair system. An artificial cutoff parameter p veto T is introduced in order to deal with the infrared divergences. In the first part, the transverse momentum of the Higgs pair system is required to be less than the cutoff parameter. With a sufficiently small p veto T , we can safely ignore all the power-suppressed terms in p veto T . In such a case, the cross section admits a factorization form which can resum the large logarithms ln n (m hh /p veto T ) to all orders in α s . In the following context, we will use the softcollinear effective theory [59][60][61][62][63], in which the factorized cross section can be written as a convolution of the transverse momentum dependent (TMD) beam function, soft function and hard function [64]. The rapidity divergences appearing in the calculation of the TMD beam function and the soft function need an additional regulator besides the dimensional regularization [65][66][67][68]. However, the final physical cross section is independent of such a regulator. The two-loop analytical results for these ingredients can be found in [69][70][71][72][73]. The NNLO hard function can be obtained by combining the two-loop amplitudes calculated in [74] and one-loop amplitudes we calculate analytically. The final results are expressed in terms of multiple polylogarithms, evaluated by the public Mathematica package PolyLogTools [75]. We have set up a streamline to combine the various components together in the computations of the NNLO differential cross sections of W hh [76] and Zhh [77] associated production processes. As opposed to the quark anti-quark initial states in the previous calculations, we extend our program to the gluon-gluon initial states in this work.
In the second part of class-b in Eq. (5), the transverse momentum of the Higgs pair system is imposed to be larger than the cutoff parameter p veto T . In such a case, there must be an additional jet in accompany with the Higgs pair. Therefore, in order to have NNLO cross section of class-b, we only need to calculate the NLO corrections to hh plus a jet, of which the underlying Born is represented for example by Fig.1(b) but with an additional gluon emission. In this work, we use the Mad-Graph5 aMC@NLO [78] framework to perform such calculations. The two Wilson coefficients are also expanded in a series of α s . Since the contribution of this class is from the interference between the amplitudes with only one effective vertex insertion and with two effective vertices, one has to organize these coefficients and amplitudes in an appropriate way. Thanks to the recent development [79] to handle mixed-order scenarios, we are able to obtain the results order-by-order in α s . To calculate the one-loop amplitudes automatically, we prepare the model files by using FeynRules [80], FeynArts [81] and an in-house Mathematica program, which has been validated in [82,83]. The counter-terms, especially the rational R 2 terms, have been extensively checked with the results in the literature [84,85]. The tensor integrals appearing in the one-loop amplitudes are evaluated by MadLoop [78,86] equipped with Collier [87], while the real emission contribution is computed with the module MadFKS [88,89] with the FKS subtraction method [90,91]. We want to stress that the inclusion of the contribution from class-b is indispensable in the sense that it not only contributes to the same order in α s but also cancels the remaining scale dependence in classa at N 3 LO (details shown in the supplemental material). Finally, since the NLO cross sections of class-c can be obtained with full-fledged methods, we refrain ourselves from presenting details about them, but they have been routinely included in our final results.
We have performed many cross checks and validations in our calculations. All the terms except for O(α 5 s ) terms of class-a and class-b listed in Table I have been cross checked at least by two independent calculations at the inclusive total cross section level. Specifically, we have reproduced the cross section of a single Higgs boson production up to NNLO in iHixs2 by using our program. This agreement can check our implementations of the two-loop beam and soft functions, as well as the calculation of one-loop amplitudes with one effective vertex. In addition, we have calculated the NLO and NNLO corrections to Higgs pair production in the infinite top-quark mass limit, and found agreement with HPair2 [12,13] and Ref. [18], respectively. This helps to check Eq.(3) and the calculation of one-loop amplitudes with two effective vertices. These nontrivial checks already ensure the correctness of many components of our calculations. For the O(α 5 s ) term of class-a, we simply used iHixs2 by employing Eq.(3). Such a program has been validated with the Higgs pair cross sections from LO to NNLO, which makes us convinced that the O(α 5 s ) piece of classa is correct. For the remaining O(α 5 s ) part of class-b, we carefully checked the various pieces that are used in our calculation. In particular, we have checked the scale dependence of the finite part in the two-loop amplitudes with two effective vertices [74] by the renormalization group equation that the hard function should satisfy. The one-loop amplitude can also been extracted from the scale-dependent part of the two-loop amplitudes, and it has been compared against the analytical result we calculated with fire [92] and to the numerical result from MadLoop. Again, we find perfect agreements. Moreover, we have checked the independence of the final NNLO results for class-b on the values of p veto T over the range from 4 GeV to 20 GeV (see the supplemental material).
Results -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.2 GeV. We use the PDF4LHC15 nnlo 30 PDF [93][94][95][96] provided by LHAPDF6 [97], and the associated strong coupling α s . The default central scale is chosen to be the invariant mass of the Higgs pair divided by 2, i.e. µ 0 = m hh /2, and the scale uncertainty is evaluated through the 9-point variation of the factorization scale µ F and the renor-  We present the inclusive total cross sections (from LO to N 3 LO) of the Higgs boson pair production at different center-of-mass energies in Table II and Fig. 2. Similarly to the single Higgs case, the QCD higher-order corrections are prominent. The NLO corrections increase the LO cross section by 87% (85%) at √ s = 13 (100) TeV. The NNLO corrections improve the NLO cross section further by 18% (16%), reducing the scale uncertainty by a factor of 2 to 3 to be below 8%. Finally, the N 3 LO corrections turn out to be 3.0% (2.7%), which lies well within the scale uncertainty band of the NNLO result. Now, the scale uncertainty at N 3 LO is less than 3% (2%), with another significant reduction of 2-3 times. For the purpose of the comparison, the PDF parameterization uncertainty at 13 TeV amounts to ±3.3%, which is larger than the current scale uncertainty. Such an improvement can be more clearly seen in Fig. 3, where we have varied the scale by a factor of four around the default choice with imposing µ R = µ F . The plot illustrates the importance of the choice of scales in a lower order perturbative calculation. If one chooses a scale to be larger than m hh , the higher-order QCD corrections are very sizable. Instead, if one chooses a judicious scale between m hh /4 and m hh /3, the perturbative convergence of the α s series is very good from NLO to N 3 LO. It also indicates that the commonly used central scale m hh /2 is not optimal.
Besides the inclusive total cross section, we are also able to obtain the exact N 3 LO results for a differential distribution, i.e., the invariant mass m hh distribution shown in Fig.4. As in the total cross section case, the inclusion of the N 3 LO corrections dramatically stabilizes the perturbative calculation of the invariant mass differential distribution. It can also be seen that the higherorder QCD corrections do not change the peak position, and the K factor of N 3 LO over NNLO is almost flat over a large region of m hh . The N 3 LO result with small scale uncertainty is completely enclosed within the NNLO uncertainty band. Such a feature consolidates that the perturbative expansion of this differential cross section in a series of α s is convergent up to this order.
Conclusion -We have calculated the N 3 LO QCD corrections to the Higgs boson pair production via gluon fusion at hadron colliders in the infinite top-quark mass limit. We find that the total cross section at N 3 LO increases by 3.0% (2.7%) at √ s = 13 (100) TeV with respect to the NNLO result under the central scale choice µ 0 = m hh /2. The scale uncertainty has been significantly improved at N 3 LO compared to the previous result at NNLO, which is now less than 3% (2%). In contrast, the PDF uncertainty is ±3.3% at the 13 TeV LHC. Moreover, we have computed the invariant mass distribution at N 3 LO for the first time, and almost constant improvement has been found. The perturbative series of both the total inclusive cross section and the invariant mass distribution are found to be convergent up to this order. In the future, for the phenomenological applications, we plan to combine our N 3 LO calculations in the infinite top-quark mass limit with the NLO results including exact top-quark mass dependence.

Renormalization scale dependence
In this Supplemental material, we present the method to obtain the µ R dependence in the framework of soft-collinear effective theory (SCET), where the renormalization scale is usually set to be the same as the factorization scale, while they can be distinguished clearly in fixed-order calculations.
The cross section is scale ( The individual renormalization and factorization scale dependence can be obtained through where the first part on the right hand can be predicted with q T -subtraction method in the framework of SCET and the second part is obtained by requiring the renormalization scale independence of the total cross section. The N 3 LO cross section for Higgs pair production is renormalization scale invariant up to O(α 6 s ) corrections, i.e.
For class-a, the differential equation is where σ h has the expansion σ h = σ h ∝ α 2+i s . Up to N 3 LO, we need the NLO QCD corrections to class-c cross section which is standalone and scale invariant. Therefore, for class-b, the renormalization group equation is The ratio of C hh over C h can be written in a series of a s ≡ α s (µ R )/4π, where the coefficient δ 2 is scale independent. Therefore, we obtain with β 0 = (11C A − 2n f )/3 . Then Eq. (5) becomes The scale-invariance violation terms in Eq. (8) start from NNLO corrections to class-b Higgs pair production.