Anomalous couplings in Higgs-boson pair production at approximate NNLO QCD

We combine NLO predictions with full top-quark mass dependence with approximate NNLO predictions for Higgs-boson pair production in gluon fusion, including the possibility to vary coupling parameters within a non-linear Effective Field Theory framework containing five anomalous couplings for this process. We study the impact of the anomalous couplings on various observables, and present Higgs-pair invariant-mass distributions at seven benchmark points characterising different $m_{hh}$ shape types. We also provide numerical coefficients for the approximate NNLO cross section as a function of the anomalous couplings at $\sqrt{s}=14$ TeV.


Introduction
The Higgs sector is being explored at the LHC with impressive success, having established the coupling of the Higgs boson to all the electroweak gauge bosons [1], the top [2,3] and bottom quarks [4,5], the tau lepton [6,7] and recently moving towards indication for Higgs couplings to muons [8,9]. In addition to these important results, it is also crucial to study Higgs boson self-couplings, which provide a way to explore the potential that drives the electroweak symmetry breaking mechanism. The trilinear Higgs boson self-coupling is still rather weakly constrained, even though the limits have improved considerably in Run II [10][11][12][13]. The precision of the measurements of Higgs boson couplings to other particles and itself will increase substantially in the high-luminosity phase of the LHC [14,15]. Therefore, simulations of the effects of anomalous couplings in the Higgs sector also need to achieve rather small uncertainties.
Higgs-boson pair production in gluon fusion is a process with direct access to the trilinear Higgs boson self-coupling c hhh . An established deviation from its Standard Model (SM) value would be a clear sign of new physics. However, if c hhh is different from the value predicted by the SM, it is very likely that other Higgs couplings are also modified, such that several anomalous couplings, entering the process gg → HH dominantly within an Effective Field Theory (EFT) framework, should be studied simultaneously. Furthermore, in order to get reliable predictions, higher-order QCD corrections need to be taken into account, not only for the SM cross section but also in the case anomalous couplings are included.
The loop-induced process gg → HH at leading order has been calculated in Refs. [16][17][18]. Before the full next-to-leading order (NLO) QCD corrections became available, the m t → ∞ limit ("Heavy Top Limit, HTL"), sometimes also called "Higgs Effective Field Theory (heft)" approximation, has been used. In this limit, the NLO corrections were first calculated in Ref. [19] using the so-called "Born-improved HTL", which involves rescaling the NLO results in the m t → ∞ limit by the LO result in the full theory. In Ref. [20] an approximation called "FT approx " was introduced, which contains the real radiation matrix elements with full top-quark mass dependence, while the virtual part is calculated in the Born-improved HTL approximation. The NLO QCD corrections with full top-quark mass dependence became available more recently [21][22][23][24]. The NLO results of Refs. [21,22] have been combined with parton shower Monte Carlo programs in Refs. [25][26][27][28], where Ref. [28] introduced the possibility of varying five Higgs couplings.
In the m t → ∞ limit, the next-to-next-to-leading order (NNLO) QCD corrections have been computed in Refs. [29][30][31][32][33]. The calculation of Ref. [33] has been combined with results including the top-quark mass dependence as far as available in Ref. [34], defining an NNLO FTapprox result which contains the full top-quark mass dependence at NLO as well as in the double real radiation part. Soft gluon resummation combined with these results has been presented in Ref. [35]. N 3 LO corrections are also available [36,37], where in Ref. [37] the N 3 LO results in the HTL have been "NLO-improved" using the results of Refs. [25,27]. In Ref. [38], NNLO results in the Born-improved heavy top limit including the effect of anomalous couplings were presented, see also Ref. [39] for the NLO case.
The scale uncertainties at NLO are still at the 10% level, while they are decreased to about 5% when including the NNLO corrections and to about 3% at N 3 LO in the "NLO-improved" variant. The uncertainties due to the chosen top mass scheme have been assessed in Refs. [23,24,40].
In this work we present a study of the anomalous couplings relevant to the process gg → HH within a non-linear EFT operator expansion, at approximate NNLO in QCD (which we will denote by NNLO ). It builds on the results presented in Ref. [38], and extends them to include the full NLO QCD corrections from Ref. [41]. Thus NNLO contains NNLO results in the heavy top limit, where exact Born expressions have been used whenever the higher-order corrections in the HTL factorise, as well as NLO corrections with full top-quark mass dependence.
In particular, we provide coefficients for all the possible combinations of anomalous couplings that can occur at NNLO for the total cross section, analogous to what has been provided at NLO in Ref. [41]. These coefficients can then be used to reconstruct the cross section for any combination of coupling values. Furthermore, we show results for seven benchmark points characteristic of certain shape types of the Higgs-boson pair invariant-mass distribution m hh , which have been identified by an NLO shape analysis presented in Ref. [42]. This paper is organised as follows. In Section 2 we describe the theoretical framework and the definition of the anomalous couplings. Section 3 contains the phenomenological results, including total cross sections and differential distributions for the benchmark points, and the description of the fitting procedure for the coefficients of the anomalous couplings, together with the corresponding results. Finally, we conclude in Section 4.

Higgs-boson pair production within an EFT framework
The calculation builds on the ones presented in Refs. [28,38,41] and therefore the methods will be described only briefly here, focusing on the new aspects.
We work in a non-linear EFT framework, sometimes also called Electroweak Chiral Lagrangian (EWChL) including a light Higgs boson [43,44] or HEFT (Higgs Effective Field Theory), not to be confused with the heavy top limit, which is sometimes also called heft. It relies on counting the chiral dimension of the terms contributing to the Lagrangian [45], rather than counting the canonical dimension as in the Standard Model Effective Field Theory (SMEFT). As a consequence, the EWChL is also suitable to describe strong dynamics in the Higgs sector. Applying this framework to Higgs-boson pair production in gluon fusion, keeping terms up to, and including, chiral dimension d χ = 4, we obtain the effective Lagrangian relevant to this process as In the EWChL framework there are a priori no relations between the couplings. In general, all couplings may have arbitrary values of O(1). The conventions are such that in the SM c t = c hhh = 1 and c tt = c ggh = c gghh = 0. The EWChL coefficients can be related [39] to those in the SMEFT at Lagrangian level, however how to treat double insertions of operators and squared dimension-6 terms at cross section level is less straightforward when attempting to relate the two EFT frameworks. The diagrams at leading order in QCD and chiral dimension four are shown in Fig. 1. There are different normalisation conventions for the anomalous couplings in the literature. In Table 1 we summarise some conventions commonly used. For the relation Eq. (2.1), i.e. L of Ref. [41] Ref. [46] Ref. [39] c to the corresponding parameters in the SMEFT we refer to Ref. [28].
In Ref. [41] the NLO QCD corrections were calculated within this framework. These corrections have been implemented into the code ggHH [27,28], which is a public code available at [47], where the real radiation matrix elements were implemented using the interface between GoSam [48,49] and the POWHEG-BOX [50][51][52] and the virtual twoloop corrections use the results of the calculations presented in Refs. [21,22]. In turn, approximate NNLO QCD corrections in a similar EFT framework have been computed in Ref. [38], working in the heavy top limit improved by inserting LO form factors with full top mass dependence. The results described in this paper are obtained by performing an additive combination of the ones presented in Refs. [38,41], i.e. keeping the full top mass dependence up to NLO and adding the approximate results for the genuine NNLO piece. We will denote our results by NNLO .
We can describe the dependence of the NNLO cross section on the five anomalous couplings in terms of 25 coefficients a i , following Refs. [41,46,53]: While at LO only the first 15 coefficients contribute, at NLO 23 coupling combinations occur and at NNLO 25 combinations are possible.
It is worth stressing that the functional dependence described in Eq. (2.2), which holds for the exact NNLO cross section, is also valid for the approximate results presented in this work. To this end, the fact that the combination of the full NLO and approximate NNLO is done through an additive approach becomes crucial, since e.g. a simple multiplicative rescaling would not respect the functional form of Eq. (2.2) (and, in particular, due to the new combinations appearing at NNLO it would be ill-defined for certain values of the anomalous couplings).

Phenomenological results and coupling coefficients
Our results are calculated at a centre-of-mass energy of √ s = 14 TeV. The parton distribution functions PDF4LHC15 [54][55][56][57] have been used, interfaced via LHAPDF [58], along with the corresponding value for α s (µ). We use NLO parton densities and strong coupling evolution for the LO and NLO predictions, and NNLO PDFs and strong coupling evolution for the NNLO calculation. The masses of the Higgs boson and the top quark have been fixed to m h = 125 GeV, m t = 173 GeV and their widths have been set to zero. The top-quark mass is renormalised in the on-shell scheme.
The scale uncertainties are estimated by varying the factorisation/renormalisation scales µ F , µ R , where the bands represent 3-point scale variations around the central

Total cross section
In Ref. [42], shapes of the Higgs-boson pair invariant-mass distribution m hh were analysed in the 5-dimensional coupling parameter space, using machine learning techniques to classify m hh -shapes from NLO predictions. This procedure led to seven shape characteristic benchmark points, for which we also show NNLO results. For convenience we repeat the benchmark points in Table 2, together with the corresponding values for the NLO and NNLO cross sections at √ s = 14 TeV.
The NNLO result is rescaled by a constant factor of 0.958 in order to match the NNLO FTapprox total cross section of Ref. [34] when the SM limit is taken. 1 We note that all other results in this paper (apart from the invariant-mass distributions in Fig. 2), including the parametrisation presented in the following section, are insensitive to this additional factor, which cancels out in the ratio σ/σ SM .
From the results in Table 2 we can observe that the size of the QCD corrections, described by the NNLO K-factor, have a sizeable dependence on the EFT parameters. In the particular case of the seven benchmark points under consideration, the K-factor can differ from the SM one by almost 40%. For all the EFT points in Table 2   observe a sizeable reduction in the scale uncertainties once the NNLO corrections are included.

Differential results and heat maps
In Fig. 2 we present the di-Higgs invariant-mass distribution for the SM and for the seven shape benchmarks listed in Table 2, both at NLO and NNLO .
In Fig. 2 we can observe a clear reduction of the scale uncertainties between NLO and NNLO . Furthermore, some of the benchmark points show a rather non-uniform differential NNLO /NLO K-factor, where in particular the region close to the Higgsboson pair production threshold is enhanced compared to NLO, while in the SM case the enhancement is stronger in the tail of the distribution. Note that the c hhh value for both benchmark 3 and 4 is close to the value c hhh 2.5 where the cross section goes through a minimum as a function of c hhh only, and shows a dip in the m hh distribution if the other couplings are SM-like. However, for benchmark 4, this dip is not present because the values of the other couplings destroy the strong cancellations that lead to the dip.   Table 2, at NLO (blue) and NNLO (red).
The low value of c t for benchmark 4 has the effect of reducing the contribution of the box-type diagrams to the cross section, such that the triangle-type contributions, dominating at low m hh , play a larger role, even more so at NNLO due to the dominance of the soft-virtual contributions in this kinematic region. Note also that the coefficients a 24 and a 25 in Eq. (2.2), which are only present at NNLO, contain c ggh to the powers 4 and 3, respectively, therefore the influence of this coupling has a different functional dependence at NNLO, which can lead to shape distortions.
In Fig. 3 (Fig. 4) we present the NNLO total cross section and K-factor, both normalised to the corresponding SM value, in the c hhh -c ggh (c hhh -c tt ) plane. 2 Here the K-factor is defined as the ratio between the NNLO and LO predictions. (3.1) The chosen range for each coefficient is motivated by current experimental constraints [12,13,[59][60][61][62]. However, since these limits are not based on a simultaneous fit of all the couplings under consideration here, we have conservatively increased the allowed range   Fig. 3 for the c hhh − c tt plane. of variation. 3 We note that in the case of c gghh the interval in Eq. (3.1) is chosen arbitrarily, since there are no experimental constraints on its value.
We can compare the curves in Fig. 5 (left) with the best experimental limit on the total di-Higgs production cross section of 4.1 times its SM value at 95% confidence  Figure 5: Total cross section, normalised to its SM value, as a function of each of the anomalous couplings. The horizontal line corresponds to the best current experimental limit [13]. The panel on the right shows the ratio to the NLO curves, that is (σ/σ SM ) NNLO /(σ/σ SM ) NLO .
level [13], indicated by the yellow horizontal line. From this comparison we conclude that, with the only exception of c ggh , all of the anomalous couplings can, within the limits given in Eq. (3.1), generate variations in the di-Higgs cross section which are larger than the current experimental limit. While Eq. (3.1) describes the allowed region in the EFT parameter space only in a qualitative way, the results in Fig. 5 (left) clearly indicate that, with the present level of precision achieved by the LHC, a simultaneous variation of all couplings is needed for a meaningful EFT analysis.
The right panel of Fig. 5 shows the ratio of the NNLO and NLO results for the same coupling modifications that are present in the left panel of the figure. As we are comparing perturbative predictions for the ratio σ/σ SM , we can expect a reduced impact from higher-order corrections, since their effect will partially cancel between numerator and denominator (in particular, the quantity presented in the right panel of Fig. 5 is, by definition, equal to 1 in the SM limit). In fact, we can observe that for some of the curves displayed in Fig. 5 the difference between NNLO and NLO is below 1%. Other coupling variations (c ggh and c gghh ) feature a larger deviation from the NLO prediction, the difference going up to 8% in the case of c gghh for the considered range. We note that larger deviations between the NNLO and NLO predictions for σ/σ SM can be obtained when all couplings are allowed to vary simultaneously.

Fitting procedure and results for the coupling coefficients
In order to obtain a parametrisation of the NNLO total cross section as a function of the anomalous couplings, we perform a non-linear fit of the a i coefficients in Eq. (2.2). The fit is carried out with the software Mathematica, and is based on the values of the NNLO cross section computed for 43 different points in the EFT parameter space. The list of the points is shown in Appendix A. We repeat the fit for the three scale choices of µ R and µ F described in the previous section. The result of the fits is shown in Table 3.
We perform an inverse-variance weighted fit, using the uncertainty of the NNLO cross sections. This uncertainty is dominated by the statistical uncertainty associated with the determination of the top-mass-dependent two-loop virtual corrections, which are calculated numerically from a finite number of 6715 phase-space points, while the numerical uncertainties from the Monte Carlo integration are considerably smaller. Because the phase-space points were originally generated to populate the m hh distribution for SM values of the couplings, the phase-space cover is suboptimal in regions where the differential cross section is much larger than in the SM case. In order to estimate this uncertainty, we propagate the statistical uncertainty stemming from the virtual contribution to the total differential cross section in m hh , for all of the 43 EFT points used in the fit, see Tables 4 and 5. Further details about the choice of the EFT points are given in Appendix A.
This uncertainty of the NLO contribution is ≤ 2% in all bins, except in the first m hh bin. There it can reach ∼ 20% for a few points, and is otherwise around 4 − 12% for most points. Thus the uncertainty on the total cross section is typically at the percent level for the majority of the points used in the fit. If larger uncertainties are found for some points in the EFT parameter space they feature either: • large numerical cancellations between different contributions at NLO, • an enhancement of the low-m hh region (and to a lesser extent of the tail for benchmark points with non-zero values of c ggh and c gghh ). These regions, being suppressed in the SM case, are less densely populated by the numerical grid used to compute the NLO virtual corrections, • or a proportionally large contribution from the virtual corrections to the total cross section.
In order to check the robustness of our fit and to estimate the uncertainties associated to it, we produce 1000 replicas of the 43 cross section values entering the fit,  by randomly generating Gaussian variations around their central value, and using the corresponding uncertainty as the standard deviation. Based on these replicas, we in turn obtain 1000 replicas of the fitted cross section. Finally, we perform a scan in the EFT parameter space, randomly generating 10000 points in the range indicated in Eq. (3.1).
We study the variation in the cross section by computing its standard deviation over the 1000 fit replicas, for each of the 10000 points in the EFT parameter scan. We note that this standard deviation can be regarded as an estimate of the final uncertainty in the cross section, encompassing both the uncertainties in the 43 original values of the cross section as well as the uncertainties coming from the fitting procedure. We have found these standard deviations to be at the percent level. Specifically, for 68% (95%) of the points in the EFT scan this deviation is below 2.0% (5.8%).
Based on this analysis, we conclude that the uncertainties in the cross sections that are obtained by using Eq. (2.2) together with the results in Table 3 are at the O(5%) level. We note that this uncertainty is of the same order of magnitude as the one in the cross section values used to perform the fit. In other words, the fitting procedure does not produce a substantial increase in the uncertainty.
A more detailed estimation of the numerical uncertainties is beyond the scope of this paper. We note, however, that our O(5%) estimation is mostly below the scale uncertainties at NNLO , as can be seen from Table 2, and also below uncertainties related to the renormalisation scheme dependence of the top-quark mass, which are estimated to be between +4% and −18% at total cross section level in the SM case [40].

Conclusions
We have presented a combination of NLO predictions with full top-quark mass dependence with approximate NNLO corrections for Higgs-boson pair production in gluon fusion. Anomalous couplings are included in the framework of an Effective Field Theory where the dominant operators contributing to this process are parametrised by five anomalous couplings, c hhh , c t , c tt , c ggh and c gghh . Our combination, which we call NNLO , is based on an additive approach of the full NLO corrections and the NNLO corrections in the heavy top limit improved by inserting LO form factors with full top-quark mass dependence.
We have given NNLO results including scale uncertainties at √ s = 14 TeV at seven benchmark points, for the total cross section as well as for the m hh distribution. We observe a strong reduction of the scale uncertainties compared to the NLO case, by a factor of 2-3, depending on the benchmark point. The differential results also show shape distortions at NNLO for some of the benchmark points, in particular close to the Higgs-boson pair production threshold.
Parametrising the cross section in terms of all possible coupling combinations leads to 25 coupling coefficients a i , of which two combinations occur at NNLO which are not present at NLO (while 10 new combinations occur comparing LO and NNLO). Based on the above parametrisation of the cross section in terms of coupling combinations, we have produced heat maps for slices of the coupling-parameter space, for both the cross section and the K-factor (normalised to their SM values), where we varied the anomalous couplings in a range motivated by current constraints. The feature that the cross section is much more sensitive to variations of c hhh and c tt than of c ggh , found already at NLO, also can be seen at NNLO . Furthermore, we presented individual variations of couplings and compared them to the current best limit on the total cross section, which leads to the conclusion that simultaneous variations of all couplings are needed for a meaningful EFT analysis. We also provided the values for the a i coefficients at three different scale choices, thereby allowing fast and flexible studies of anomalous coupling variations at NNLO level.