Differential Higgs Boson Pair Production at Next-to-Next-to-Leading Order in QCD

We report on the first fully differential calculation for double Higgs boson production through gluon fusion in hadron collisions up to next-to-next-to-leading order (NNLO) in QCD perturbation theory. The calculation is performed in the heavy-top limit of the Standard Model, and in the phenomenological results we focus on pp collisions at 14 TeV. We present differential distributions through NNLO for various observables including the transverse-momentum and rapidity distributions of the two Higgs bosons. NNLO corrections are at the level of 10%-25% with respect to the next-to-leading order (NLO) prediction with a residual scale uncertainty of 5%-15% and an overall mild phase-space dependence. Only at NNLO the perturbative expansion starts to converge yielding overlapping scale uncertainty bands between NNLO and NLO in most of the phase-space. The calculation includes NLO predictions for pp ->HH+jet+X. Corrections to the corresponding distributions exceed 50% with a residual scale dependence of 20%-30%.


Introduction
The discovery of the Higgs boson [1,2] during Run I of the Large Hadron Collider (LHC) opened the door towards direct tests of electroweak symmetry breaking. To this end the search for the production of Higgs boson pairs is one of the main goals of ongoing and future runs of the LHC. Only this production mode allows for direct tests of Higgs trilinear self-couplings, whose knowledge in turn is necessary to reconstruct the scalar potential responsible for electroweak symmetry breaking.
As it is the case for single Higgs boson production, the dominant production mode for Higgs boson pairs in the Standard Model (SM) proceeds at hadron colliders via gluon fusion, mediated by heavyquark loops. At the leading order (LO) in QCD [3][4][5] there are two interfering production mechanisms: either the two Higgs bosons couple directly to a heavy-quark loop via a box diagram, gg → HH, or they couple via the trilinear Higgs coupling λ to an off-shell Higgs boson, which in turn is produced via a triangular loop, similarly to single Higgs boson production, gg → H * → HH.
Due to the loop suppression of the LO process and additional large accidental cancellations between "triangle" and "box" contributions in the scattering amplitude [3], signal rates for Higgs boson pair production are only at the level of few fb at 8 TeV and of tens of fb at 13 − 14 TeV. These small rates, together with large irreducible backgrounds in the relevant bbγγ [6][7][8], bbττ [7][8][9], bbW + W − [10] and bbbb [11,12] final states, pose a true challenge to experimental searches for Higgs boson pair production. Consequently, currently only upper limits exist, constraining the cross section for Higgs boson pair production to the level of about 70 times the SM prediction [13][14][15][16][17][18]. However, the production cross section can significantly be altered by new-physics effects, for example due to new loop contributions [19], altered tth or novel tthh couplings [20], or due to new resonances [21]. Thus, future measurements (together with precision predictions) on Higgs boson pair production, and in particular ratios of cross-section measurements [22], do not just serve as stringent tests of electroweak symmetry breaking in the SM, but might also open the door towards physics beyond the SM.
Given the loop-induced nature of the scattering process for Higgs boson pair production, higher orders in perturbation theory are extremely difficult to calculate. Only very recently a complete nextto-leading order (NLO) calculation became available [23], where the required multi-scale two-loop scattering amplitudes have been evaluated via numerical integration. Prior to the pioneering work of Ref. [23] the possible impact of the NLO corrections has been studied in Refs. [24][25][26][27][28] via asymptotic expansions of the two-loop virtual diagrams in the inverse top-quark mass.
Assuming the top quark to be heavy and all other quarks to be massless, an effective theory can be formulated, introducing a tree-level coupling of gluons and Higgs bosons. In this heavy-top limit NLO corrections to Higgs boson pair production have been presented in Ref. [29], where a rescaling with the exact Born cross section was performed. The obtained NLO corrections in the so-called Bornimproved heavy-top limit increase the total cross section by about 100% with sizable remaining scale uncertainties. In contrast, the exact NLO computation presented in Ref. [23] yields a result that is smaller by 14%.
In order to further improve on these predictions, and in particular to reduce the remaining scale uncertainties, the effective theory allows for the computation of perturbative corrections beyond NLO. To this end, employing the amplitudes derived in Ref. [30] (supplemented by Ref. [24]) a calculation of Higgs boson pair production at next-to-next-to-leading order (NNLO) accuracy in the heavy-top approximation was presented in Ref. [31]. At the inclusive level the NNLO corrections increase the cross section by about 20% with respect to the NLO prediction, leaving scale uncertainties at the level of 8% − 10%. Besides inclusive NNLO cross sections the calculation of Ref. [31] offers differential predictions in the invariant mass of the produced Higgs boson pair indicating a rather mild phase-space dependence of the NNLO corrections. Still working in the heavy-top limit, soft-gluon resummation up to next-to-next-to-leading logarithmic (NNLL) accuracy has been carried out, and the results were matched to the NLO [32] and the NNLO [33] fixed-order computations. At NNLL+NNLO accuracy the theoretical uncertainties on the inclusive cross section due to QCD effects are reduced to about 5% [33]. Furthermore, extending the SM with additional dimension-6 operators [34], NLO corrections in the heavy-top limit have been presented in Ref. [35]. Notably, in Refs. [36,37] a reweighting technique has been presented, allowing to combine exact one-loop real corrections of Higgs boson pair production with the corresponding virtual contributions, the latter computed in the effective theory. On the other hand, the real corrections, which imply the evaluation of one-loop amplitudes with an extra parton in the final state, have been computed in an exact way and used to merge LO samples for HH + 0, 1 jets [38,39] in order to obtain more reliable exclusive distributions.
In this paper we extend the calculation of Ref. [31] providing fully differential NNLO predictions for Higgs boson pair production in the heavy-top approximation of the SM via a flexible Monte Carlo implementation. The calculation is based on the combination of the q T subtraction formalism [40] with the Monte Carlo framework Munich † , supplemented by tree and one-loop amplitudes from Open-Loops [41]. Employing these tools we provide NNLO predictions for various kinematic distributions that are relevant for searches and precision measurements of Higgs boson pair production at the LHC. The calculation includes NLO predictions for pp → HH + jet + X. Corresponding differential distributions are studied in detail. In our study we refrain from a reweighting using exact LO or NLO matrix elements or cross sections. Such a reweighting should eventually be performed employing the results of Ref. [23]. For the time being we focus on the differential NNLO/NLO correction factors obtained in the effective theory, which are the main result of our paper. Such correction factors provide valuable information that can directly be applied to any pp → HH + X NLO prediction in different approximations.
The paper is organized as follows. In Sect. 2 we introduce the heavy-top limit for multi-Higgs production at higher orders in perturbation theory together with the technical ingredients of our calculation. Numerical results are presented in Sect. 3, and in Sect. 4 we summarize our results.

Higgs boson pair production through NNLO in the heavy-top limit
In the heavy-top approximation effective tree-level couplings between gluons and Higgs bosons are introduced via the effective Lagrangian [29,42,43] where v 246 GeV is the vacuum expectation value of the Higgs field. In this effective Lagrangian only couplings relevant for our calculation are shown, while in general within this effective theory there are also further couplings for any number of Higgs bosons to gluons. The matching coefficients C H and C HH can be expanded in powers of α S via the following parametrization, The perturbative expansion for both coefficients is known up to O(α 3 S ) and reads [24,42,44,45] where n F is the number of light quarks, µ R the renormalisation scale and m t the pole mass of the (heavy) top quark. As can be seen from Eq.
As already discussed, in the full theory Higgs boson pair production at LO is governed by "box" and "triangle" contributions. In the heavy-top limit the corresponding scattering amplitudes manifest as tree-level diagrams with one double-Higgs and one single-Higgs operator insertion, respectively. At NLO, in the perturbative expansion we have the usual real and virtual contributions, where the former includes gluon and quark bremsstrahlung, and the latter are given by one-loop corrections to the diagrams mentioned before. However, at the same order of perturbation theory there is an additional contribution with Born-level kinematics, originating from amplitudes with two single-Higgs operator insertions in interference with the LO amplitude [29]. In the full theory such contributions correspond to reducible double-triangle two-loop diagrams.
This pattern also appears at higher orders, and in particular the NNLO virtual contributions have to include both two-loop corrections to amplitudes with one single-or double-Higgs operator insertion, and one-loop corrections to amplitudes with two single-Higgs operator insertions [30]. These NNLO virtual contributions have to be combined via an appropriate subtraction scheme with double-real and real-virtual contributions of the same perturbative order. Similarly to what was discussed before, the real-virtual contributions, i.e. the virtual amplitudes for HH + jet production, have to be extended to include two single-Higgs operator insertions in interference with the corresponding tree-level amplitude. More details on the technical implementation of such double-operator insertions in our calculation are given in Sect. 2.3.

q T subtraction
In order to handle infrared singularities in the NNLO calculation, we apply the q T subtraction formalism [40]. In this approach the separation between genuine NNLO singularities, located where the transverse momentum of the Higgs pair, q T,HH , is zero, from NLO-like singularities in the HH + jet contribution is explicit. As a consequence, the contribution dσ HH+jet NLO in the q T subtraction formula, can be evaluated using any well-established subtraction procedure at NLO. The remaining divergence in the limit q T,HH → 0 is cancelled by the process-independent counterterm dσ CT NNLO . The implementation is fully general, and it is based on the universality [46] of the hard-collinear coefficients H HH pp → HH + X @ 14 TeV pp → HH + X @ 14 TeV pp → HH + X @ 14 TeV pp → HH + X @ 14 TeV Figure 1: Dependence of the pp → HH + X cross sections at 14 TeV on the q T -subtraction cut, r cut , for both NLO (left plot) and NNLO (right plot) results. NLO results are normalized to the r cutindependent NLO cross section computed with Catani-Seymour subtraction, and the NNLO results are normalized to the r cut -independent inclusive NNLO cross section calculated in the framework of Ref. [33]. The blue band indicates the NNLO result from q T subtraction in the limit r cut → 0, with an approriate extrapolation-error estimate.
appearing in the first term on the right hand side of Eq. (4). The general structure of these coefficients for gluon-initiated processes has been presented in Refs. [47,48]. Their process dependence is embodied in a single perturbative hard factor which is obtained from the two-loop virtual correction, derived for this process in Ref. [30], through an appropriate subtraction procedure [46].
The difference in the square bracket in Eq. (4) is formally finite as q T → 0, but each term separately exhibits logarithmic divergences in this limit. In practice a small technical cut, r cut , needs to be applied on r ≡ q T /Q, where Q is typically chosen as the invariant mass of the final-state system (so here Q = m HH ). After cancellation of these logarithms between the real contribution dσ HH+jet NLO and the counterterm, the remainder shows a very slight r cut dependence below about r cut = 1%; we thus use the finite-r cut results to extrapolate to r cut = 0, and conservatively assign an additional extrapolation error to our results. We verified in detail that for Higgs boson pair production the NNLO result is indeed very stable when varying the cut parameter. More precisely, in the range below r cut = 1%, the variation in the NNLO cross section is of O(0.1%), and our extrapolated result is in good numerical agreement with the analytic result of Ref. [33] (see Fig. 1).

Tree and one-loop amplitudes from OpenLoops
All tree and one-loop amplitudes, i.e. in particular the one-loop amplitude for the dσ HH+jet NLO contribution, are provided by the publicly available ‡ OpenLoops amplitude generator [41], which is based on a fast numerical recursion for the generation of tree and one-loop scattering amplitudes [49].
In order to extend OpenLoops to one-loop corrections in the heavy-top limit of the SM, all relevant Feynman rules for single-Higgs [50] and double-Higgs [35,51] production have been implemented ‡ The publicly available OpenLoops process library includes all relevant matrix elements to compute NLO QCD corrections, including colour-and helicity-correlations and real radiation as well as loop-squared amplitudes, for more than a hundred LHC processes. Amplitudes for Higgs boson pair production (+1 jet) at NLO in the heavy-top limit have been made available together with this publication, while amplitudes for pp → H, pp → Hj, pp → Hjj and pp → Hjjj have been publicly available already for some time.
in the framework of the numerical open-loops recursion including UV renormalisation and the rational contributions of type R 2 [52]. Combined with the OPP reduction method [53] implemented in Cut-Tools [54] and the scalar one-loop library OneLOop [55] the employed recursion permits to achieve very high CPU performance and a high degree of numerical stability. The small fraction of numerically unstable matrix elements is automatically detected and rescued through re-evaluation in quadruple precision.
Technically, the effective field theory of Eq. (1) introduces various features which do not appear in the Standard Model. Most notably, the Feynman rules for the dimension-5 and -6 operators G µν G µν H(H) introduce, apart from 5-and 6-point vertices, the Lorentz structure p 1 p 2 g µν − p ν 1 p µ 2 (where p µ 1 and p ν 2 are the gluon momenta) which, if present in a loop, raises the tensor rank of the amplitude by 2. In the calculation of HH(+jet) production at one-loop level such operators enter only once, leading to a tensor rank up to one higher than the number of loop propagators. The reduction of such amplitudes with one "additional" tensor rank is supported by CutTools. Furthermore, the O(α S ) contributions to the matching coefficients C H and C HH must be included. Considering the order of coupling powers it is natural to treat these contributions as counterterms. As discussed above, at the same order of perturbation theory as the one-loop scattering amplitudes for HH(+jet) production, contributions from two single-Higgs operator insertions at tree-level in interference with the LO tree-level amplitude with one double-Higgs operator insertion have to be considered. Similarly to the O(α S ) contributions to C H and C HH , these contributions are included via dedicated O(α S ) pseudo-counterterms.

Validation
One-loop amplitudes for single Higgs boson production plus up to two jets have been extensively validated against the results of Refs. [56][57][58][59][60][61][62] implemented in Sherpa [63] and MCFM [64] (via the corresponding implementation in the POWHEG-BOX [64,65]). Due to the lack of publicly available alternatives, the validation of the one-loop amplitudes for Higgs boson pair production plus jets had to rely on various internal cross checks.
We performed a calculation of pp → H + X up to NNLO in the heavy-top limit in the same framework as employed for pp → HH +X, and compared against the results obtained with the inclusive analytical codes of Refs. [66][67][68], where agreement well beyond the per mill level was found. Due to the similarity of the two processes this serves as a strong cross check for many technical ingredients of the calculation presented here.
In order to validate all ingredients of the computation of Higgs boson pair production in the heavy-top limit presented in this paper, the LO, NLO and NNLO inclusive cross sections computed in Ref. [31] have been reproduced at the per mill level § (see Fig. 1). Additionally, mutual agreement has been found for the invariant mass distribution m HH up to NNLO comparing against the results of Ref. [31]. Furthermore, the NLO results have been computed in the q T subtraction formalism and also employing the dipole subtraction framework [69,70] within Munich, where again we found mutual agreement far beyond the per mill level (see again Fig. 1). § In [31] the relation C

Results
In the following we present predictions for Higgs boson pair production at the LHC including perturbative fixed-order corrections up to NNLO in the heavy-top limit. Inclusive results will be presented for centre-of-mass energies of √ s = 13 TeV and √ s = 14 TeV, while at the differential level we restrict ourselves to √ s = 14 TeV. SM input parameters are chosen according to the recommendations of [71], which in particular implies v = 246.2 GeV, m t = 173.2 GeV and Here, the top-quark mass does only enter via the NNLO contributions to the matching coefficients, as given in Eq. (3). For the calculation of hadron-level cross sections we employ the PDF4LHC15 [72] parton distribution functions (PDFs), and use the corresponding NLO PDFs for our LO and NLO predictions and NNLO PDFs for the NNLO predictions. ¶ Couplings are evaluated using the running strong coupling provided by the respective PDFs. All light quarks, including bottom quarks, are treated as massless particles, i.e. n F = 5, while the top quark does not contribute explicitly in the employed heavy-top limit. To define jets, we employ the anti-k T jet clustering algorithm [74] with R = 0.4 and require p T j > 30 GeV and |η j | < 4.4. In all results the renormalisation scale µ R and factorisation scale µ F are set to µ R,F = ξ R,F µ 0 , with µ 0 = m HH /2 and  and NNLO corrections are sizable, and only at NNLO the perturbative convergence becomes manifest with overlapping scale uncertainty bands between the NLO and NNLO predictions in most of the considered phase-space regions. At the same time theoretical uncertainties estimated by the scale variations described above are approximately halved when going from NLO to NNLO. The NNLO distributions shown in Figs. 5-7 are effectively only of next-to-leading order as they are either trivial or not defined at LO. They can be considered as a computation of HH + jet at NLO. Nevertheless, in the following discussion we always denote the highest considered perturbative order as NNLO (with respect to pp → HH + X).
Differential distributions in the transverse momentum and the rapidity of the two Higgs bosons, ordered by their hardness in p T , are shown in Fig. 2 and Fig. 3, respectively. NNLO corrections are overall at the level of 10% − 25% with a rather mild phase-space dependence. In particular, in the transverse-momentum distribution of the harder Higgs boson, p T,H 1 , the NNLO corrections slowly increase as p T,H 1 increases, while for the softer Higgs boson the corrections are to a large extent independent of p T,H 2 , except for the very low p T,H 2 region. As for the rapidity distributions, the NNLO effect is largely constant and at O(20%) for both y H 1 and y H 2 .
In the left plot of Fig. 4 we show predictions for the invariant-mass distribution of the produced Higgs boson pair, m HH . NNLO corrections are at the level of 18% with respect to NLO and hardly show any phase-space dependence. NNLO predictions in m HH have already been presented in Ref. [31], and the corresponding results are overlaid in Fig. 4 (left). In the computation of Ref. [31] IR singularities are analytically cancelled, thereby leading to negligible numerical fluctuations in the shown distribution. Within statistical uncertainties the results obtained from the two completely independent implementations agree perfectly.  In the right plot of Fig. 4 predictions for the rapidity of the Higgs boson pair, y HH , are presented. Again, we observe a mild phase-space dependence, with increasing NNLO corrections only for large rapidities. In all distributions in Figs. 2-4, NNLO scale uncertainties are reduced to the level of ±(5% − 12%), compared to ±(15% − 20%) at NLO. In Fig. 5 we show distributions in the transverse momentum of the Higgs boson pair, p T,HH , and of the hardest jet, p T,j 1 . At NLO (which is effectively LO for non-vanishing transverse momenta) these two distributions are directly related, and in both distributions scale uncertainties reach almost 50%. The NNLO effect is larger on the p T,HH distribution than on the p T,j 1 distribution, reaching 80% for p T,HH ≈ 200 GeV compared to 60% for p T,j 1 ≈ 200 GeV. The NLO nature of these NNLO corrections is furthermore reflected by sizable scale uncertainties at the level of 30% − 40%. In the limit p T,HH → 0 the perturbative expansion fails due to the appearance of large logarithmic terms of the form log n (p T,HH /m HH ). Here, a proper resummation of such terms is required in order to achieve a reliable theoretical prediction.
At LO the two Higgs bosons are always produced back-to-back. However, at higher orders additional QCD radiation allows for a non-trivial angular separation between the two Higgs bosons. In Fig. 6 we show the corresponding distribution in the azimuthal angle between the two Higss bosons, ∆φ HH . In our fixed-order approach, NNLO corrections are large and positive in the back-to-back configuration, where they are driven by soft-gluon emission, and jump to negative values for ∆φ HH π, due to the mis-cancellation between real and virtual contributions. In this region of phase-space, large logarithmic terms should again be resummed for achieving a reliable theoretical prediction. Configurations at small angles, i.e. ∆φ HH → 0, are driven by hard gluon emission, and NNLO corrections are at the level of 60% with respect to NLO.
Finally, in Fig. 7 we investigate corrections to the ∆R separation between the two Higgs bosons and the hardest jet. Overall corrections to these observables are moderate at the level of 20% − 40% with largely overlapping uncertainty bands between NNLO and NLO. However, for small ∆R H 1 j 1 separations, due to the ordering of the Higgs bosons according to their transverse momenta the entire phase-space opens up only at the NNLO level, inducing sizable correction factors at the NLO boundary ∆R H 1 j 1 π/2.

Summary and Outlook
In this paper we have presented the first fully differential calculation for double Higgs boson production through gluon fusion in hadron collisions. We worked in the heavy-top limit, and presented results for differential distributions through NNLO QCD for various observables including the transverse-momentum and rapidity distributions of the two Higgs bosons. NNLO corrections amount to about 10% − 25% with respect to the NLO prediction and are mildly dependent on the kinematics. The residual scale uncertainty at NNLO is about 5% − 15%. Only at NNLO the perturbative expansion starts to converge, and the uncertainty bands obtained through scale variations at NLO and NNLO overlap in most of the phase-space. The calculation includes NLO QCD predictions for pp → HH + jet + X. Corrections to the corresponding distributions exceed 50% with a residual scale dependence of about 20% − 30%.
The calculation presented here is based on the combination of the q T subtraction formalism [40] with the Monte Carlo framework Munich, supplemented by tree and one-loop amplitudes from Open-Loops [41]. This framework, to be integrated in the new numerical program Matrix , which is currently under development, allows for an extremely flexible implementation.
In the present paper we have limited ourselves to strictly work in the heavy-top limit of the SM. With the exact NLO virtual contributions available since recently, a combination of the exact results at NLO accuracy with the NNLO calculation in the heavy-top limit should be performed in the future. This combination, together with the inclusion of the Higgs boson decays, will facilitate realistic phenomenological studies at an unprecedented level of precision, as required for future measurements of the Higgs trilinear coupling.