Higgs boson pair production in non-linear Effective Field Theory with full $m_t$-dependence at NLO QCD

We present a calculation of the NLO QCD corrections to Higgs boson pair production within the framework of a non-linearly realised Effective Field Theory in the Higgs sector, described by the electroweak chiral Lagrangian. We analyse how the NLO corrections affect distributions in the Higgs boson pair invariant mass and the transverse momentum of one of the Higgs bosons. We find that these corrections lead to significant and non-homogeneous K-factors in certain regions of the parameter space. We also provide an analytical parametrisation for the total cross-section and the $m_{hh}$ distribution as a function of the anomalous Higgs couplings that includes NLO corrections. Such a parametrisation can be useful for phenomenological studies.


Introduction
Exploring the Higgs sector and the mechanism of electroweak symmetry breaking is one of the primary goals for the current and future LHC program as well as other planned experiments. While some of the properties of the Higgs boson, like its mass, spin and couplings to electroweak bosons, have been measured meanwhile impressively well [1], other parameters, like the couplings to (light) fermions, and in particular the self-coupling, are still largely unconstrained and leave room for physics beyond the Standard Model, see e.g. ref. [2] for a recent review.
In the Standard Model (SM) the strength of all Higgs boson couplings is predicted; however, effects of physics beyond the Standard Model (BSM) may lead to deviations which, once firmly established, are a clear sign of New Physics. Since higher-order QCD corrections are known to be important in Higgs boson production processes, they need to be taken into account to improve the sensitivity to New Physics effects.
Given the energy gap between the electroweak scale at v 250 GeV and a New Physics scale Λ which is supposed to be in the TeV range, it is natural to parametrise the BSM effects in a model-independent way in an Effective Field Theory (EFT) framework. Such a -1 -
In the SM, Higgs boson pair production has been calculated at leading order in refs. [45][46][47]. As it is a loop-induced process, higher order calculations with full top quark mass dependence involve multi-scale two-loop integrals. Therefore, the NLO calculations until recently have been performed in the m t → ∞ limit [48] also called HEFT ("Higgs Effective Field Theory"), 1 and then rescaled by a factor B FT /B HEFT , B FT denoting the leading order matrix element squared in the full theory. This procedure is called "Born-improved HEFT" in the following. In refs. [49,50], an approximation called "FT approx " was introduced, which contains the full top quark mass dependence in the real radiation, while the virtual part is calculated in the HEFT approximation and rescaled at the event level by the re-weighting factor B FT /B HEFT .
In addition, the HEFT results at NLO and NNLO have been improved by an expansion in 1/m 2 t in refs. [51][52][53][54]. The NNLO QCD corrections in the heavy-top limit have been computed in refs. [52,[55][56][57], and they have been supplemented by an expansion in 1/m 2 t in ref. [53] and by threshold resummation, at NLO+NNLL in ref. [58] and at NNLO+NNLL in ref. [59], leading to K-factors of about 1.2 relative to the Born-improved HEFT result.
The full NLO corrections, including the top quark mass dependence also in the virtual two-loop amplitudes, have been calculated in ref. [60]. Phenomenological studies at 14 TeV and 100 TeV, including variations of the Higgs boson self-coupling, have been presented in ref. [61]. The full NLO calculation was supplemented by NLL transverse momentum resummation in ref. [62]. It also has been matched to parton shower Monte Carlo programs [63,64], where the matched result of ref. [63] is publicly available within the POWHEG-BOX-V2 framework.
Recent work also includes a combination of an analytic threshold expansion and a largem t expansion together with a Padé approximation framework [65], and analytic results based on a high energy expansion for the planar part of the two-loop amplitude [66]. 1 Sometimes the electroweak chiral Lagrangian with a light Higgs boson is also referred to as Higgs Effective Field Theory (HEFT) in the literature. The two EFTs are unrelated and should be carefully distinguished. Here we employ the term electroweak chiral Lagrangian for the non-linear EFT of physics beyond the SM, and reserve the expression HEFT for the heavy-top limit in Higgs interactions.  , top quark mass effects have been incorporated in the NNLO HEFT  calculation, including the full NLO result and combining one-loop double-real corrections  with full top mass dependence with suitably reweighted real-virtual and double-virtual contributions evaluated in the large-m t approximation [67].
Within a non-linear EFT framework, higher order QCD corrections have been performed in the m t → ∞ limit. The NLO QCD corrections have been calculated in ref. [68], recently also supplemented with the case of CP-violating Higgs sectors [69]. The NNLO QCD corrections in the m t → ∞ limit including dimension 6 operators have been presented in ref. [70]. These calculations found rather flat K-factors, which however could be an artefact of the m t → ∞ limit. One of the main goals of the present paper is to investigate whether this feature is preserved once the full top quark mass dependence is taken into account. We calculate the NLO QCD corrections to Higgs boson pair production in gluon fusion within the non-linear EFT framework, retaining the full top quark mass dependence, based on the numerical approach developed in ref. [60]. In order to quantify the different effects of the five operators and corresponding couplings that can lead to deviations from the SM in the Higgs sector, we give results for the total NLO cross section parametrised in terms of 23 coefficients of all possible combinations of these couplings, as introduced at LO in refs. [31,71]. We also show differential distributions for 12 benchmarks points which should be characteristic for clusters of BSM scenarios. Such clusters were identified in refs. [71][72][73] at leading order and represent partitions of the BSM parameter space according to the shape of the differential distributions. We demonstrate that there are regions where the NLO corrections lead to substantial and non-homogenous K-factors and provide numbers for the parametrisation of the NLO cross section, which can be used in subsequent phenomenological studies. This paper is organised as follows. In section 2, we explain the framework of the calculation. In particular, we introduce the Higgs-electroweak chiral Lagrangian and describe how it is applied to Higgs boson pair production, including the NLO QCD corrections. Section 3 is dedicated to the phenomenological results. We provide a parametrisation of the NLO cross section in terms of coefficients of all combinations of couplings occurring in the NLO cross section. Based on this parametrisation we show heat maps both at LO and at NLO, where we vary two couplings while keeping the others fixed to the SM values. Then we give results for total cross sections and differential distributions at twelve benchmark points and discuss their implications before we conclude. An appendix explains the conventions used for the tables containing the differential coefficients of the couplings in the Higgs boson pair invariant mass distribution. The values are available in csv format as ancillary files to the arXiv submission and the JHEP publication. A further appendix compares the treatment of Higgs-pair production in the Higgs-electroweak chiral Lagrangian and in SMEFT.

The Higgs-electroweak chiral Lagrangian
In the present analysis, we will describe the potential impact of physics beyond the Standard Model through the electroweak chiral Lagrangian including a light Higgs boson [19,21,74].

JHEP09(2018)057
This framework provides us with a consistent effective field theory (EFT) for New Physics in the Higgs sector, as we will summarise in the following.
To leading order the Lagrangian is given by The first line is the unbroken SM, the remainder represents the Higgs sector. Here h is the Higgs field and U = exp(2iϕ a T a /v) encodes the electroweak Goldstone fields ϕ a , with T a the generators of SU (2). v is the electroweak vacuum expectation value, P ± = 1/2±T 3 , and The trace of a matrix A is denoted by A . The left-handed doublets of quarks and leptons are written as q L and l L , the right-handed singlets as u R , d R , e R . Generation indices are omitted. In the Yukawa terms the right-handed quark and lepton fields are collected into can arise at every order in the Higgs field h n , in addition to the usual Yukawa matrices Y u,d,e . The h-dependent functions are In the limit where and all other couplings f U,n , f V,n , Y (n) f equal to zero, the Lagrangian in (2.1) reduces to the usual SM. For generic values of those parameters, the Lagrangian describes the SM with arbitrary modifications in the Higgs couplings. While the deviations of these couplings from their SM values could, in principle, be of order unity, the parametrisation in (2.1) remains relevant as long as the anomalous Higgs couplings are the dominant New Physics effects at electroweak energies. Employing (2.1), we assume that this is the case. Such a hypothesis remains to be tested experimentally. We emphasise, however, that the assumption is well motivated by the current status of Higgs coupling measurements. These still allow deviations from the SM of 10 -20% or more, considerably larger than the typical precision of 1% reached in the electroweak gauge sector [75]. A useful property of the Lagrangian (2.1) is therefore that it allows us to concentrate on anomalous Higgs couplings in a systematic way [22,76].

JHEP09(2018)057
In fact, the intuitive picture of introducing (2.1) as the SM with modified Higgs couplings can be formulated as a consistent EFT. Because of the need to write the modified Higgs couplings in a gauge-invariant way, the Higgs field has to be represented as an electroweak singlet h, independent of the Goldstone matrix U = exp(2iϕ a T a /v). The latter transforms as U → g L U g † Y under the SM gauge group. The symmetry is non-linearly realised on the Goldstone fields ϕ a . The Lagrangian (2.1) is then nonrenormalisable (in the traditional sense) as it contains interaction terms of arbitrary canonical dimension. The EFT is therefore not organised by the canonical dimension of operators, but rather by chiral counting in analogy to the chiral perturbation of pions in QCD. Chiral counting is equivalent to an expansion in loop orders L, which can be conveniently counted by assigning chiral dimensions d χ ≡ 2L + 2 to fields and weak couplings. This assignment is simply 0 for bosons, and 1 for each derivative, fermion bilinear and weak coupling: Here A µ represents a generic gauge field, ϕ the Goldstone bosons, and h the Higgs scalar. g denotes any of the SM gauge couplings g, g , g s , and y any other weak coupling, such as the Yukawa couplings or the square-roots of the parameters f V,n in the Higgs potential. Based on this counting, the leading-order expression (2.1) can be constructed from the SM field content and symmetries as the most general Lagrangian of chiral dimension 2. Leading processes are described by tree-level amplitudes from (2.1). Next-to-leading order effects come from one-loop contributions of (2.1) and from tree-level terms of the NLO Lagrangian L 4 . Both are considered to be of 'one-loop order', or chiral dimension d χ = 4.
We next apply this framework to Higgs-pair production gg → hh. Since this process is loop-induced, at leading order both one-loop diagrams built from the LO interactions, as well as tree contributions from the NLO Lagrangian have to be taken into account. The relevant terms from the effective Lagrangian L 2 + L 4 are given by [77] The first three couplings, c t , c tt , c hhh , are from L 2 , the Higgs-gluon couplings c ggh and c gghh from L 4 [19,76]. To lowest order in the SM c t = c hhh = 1 and c tt = c ggh = c gghh = 0.
In general, all couplings may have arbitrary values of O(1). Note that we have extracted a loop factor from the definition of the Higgs-gluon couplings. The leading-order diagrams are shown in figure 1. All diagrams are at the same order in the chiral counting (chiral dimension 4, equivalent to one-loop order). They illustrate the interplay between leading order anomalous couplings (black dots) within loops, and next-to-leading order terms (black squares) at tree level. All the five couplings defined in (2.6) appear in figure 1. In the following section we discuss the extension of this analysis to the next order in QCD.

Calculation of the NLO QCD corrections
Within the framework of the electroweak chiral Lagrangian, the calculation of the gg → hh amplitude can be extended to the next order in the loop expansion, that is to two-loop  order, or chiral dimension 6. In full generality, this would require to also include two-loop electroweak corrections and local terms from the Lagrangian at chiral dimension 6. The latter introduce additional couplings, parametrising subleading new-physics effects. Such effects are beyond the experimental sensitivity in the foreseeable future, given that even the determination of the LO couplings in (2.6) remains a substantial challenge. On the other hand, radiative corrections from QCD are known to be very important for gg → hh and similar processes.
For this reason, we extend the calculation of gg → hh to the next order in the non-linear EFT, but restrict the NLO corrections to the effects from QCD. Within the systematics of the EWChL this approximation corresponds to including those corrections at chiral dimension 6 that come with a relative factor of the QCD coupling g 2 s . This procedure is consistent without introducing further anomalous couplings, beyond the ones in (2.6), because this effective Lagrangian is renormalisable with respect to QCD [22]. Since the LO amplitude for gg → hh scales as ∼ g 2 s , the NLO virtual corrections of interest to us comprise all the diagrams at two-loop order carrying a factor of g 4 s . They exist as two-loop, one-loop and tree topologies, as illustrated in figures 2, 3 and 4, respectively.
In addition, real emission diagrams at O(g 3 s ) have to be included as shown in figure 5. To further clarify our approximation with respect to the full chiral expansion at NLO, we give in figure 6 a few examples of higher-order effects that are consistently neglected in our scheme.
Example (a) shows a correction from electroweak-boson exchange. It is of two-loop order, but scales as g 2 s g 2 , rather than g 4 s . It is not a NLO QCD effect and we neglect it here. Similarly, the one-loop topology in (b) counts as two-loop order, but scales only as g 2 s c hhhh , with c hhhh the (anomalous) quartic Higgs coupling.
In example (c) we consider an anomalous top-gluon coupling of the form Q ttG = y t g stL σ µν G µν t R , where the top Yukawa coupling reflects the change in chirality. This  operator is therefore (at least) of chiral dimension 4 (one-loop order) and the diagram in figure 6 (c) of two-loop order, but again not of order g 4 s . Since (2.1) assumes that the top quark is weakly coupled to the (possibly strongly interacting) new-physics sector, it is more likely that the operator comes with further weak couplings from t L and t R and thus carries chiral dimension 6. In this case, diagram (c) is of three-loop order and clearly negligible. The effect of the chromomagnetic operator on single Higgs boson production has been calculated recently in the context of SMEFT in ref. [78].
Example (d) illustrates the effect of a local Higgs-gluon interaction of chiral dimension 6, which enters at two-loop order as a tree-level topology. A possible operator would be g 2 s G a µν G a,µν ∂ λ h∂ λ h. However, this effect, although of two-loop order, does not scale as g 4 s .
Finally, we may have an operator g 3 s f abc G a µν G b,ν λ G c,λµ h, also of chiral dimension 6. Diagram (e) then amounts to a two-loop order interaction with real emission, which is beyond our approximation.
At the technical level, the NLO QCD corrections have been calculated building on the setup described in refs. [60,61], summarised briefly below.    -8 - (c) Figure 6. Higgs-pair production in gluon fusion at NLO: examples for contributions that are consistently neglected within our approximation. The dotted square indicates a local term at chiral dimension 6 (two-loop order). See text for further explanation.

Virtual corrections
The virtual part of order α 3 s consists of genuine two-loop diagrams as well as one-loop and tree-level diagrams, see figures 2, 3 and 4.
For the two-loop part, we made use of the numerical results for the two-loop virtual diagrams in the Standard Model (SM) by dividing them into two classes: diagrams containing the Higgs-boson self-coupling ("triangle-type"), and diagrams without ("box-type"). The tthh coupling generates new two-loop topologies, see e.g. the second line of figure 2. The results for these diagrams however can be obtained from the SM triangle-type diagrams by multiplying them with the inverse Higgs boson propagator and rescaling the couplings, i.e. multiplying with c tt /c hhh . The other two-loop diagrams occurring in our calculation have the same topologies as in the SM and therefore can be derived from the SM results by rescaling of the couplings c t and c hhh .
The one-loop part containing the Higgs-gluon contact interactions has been calculated in two ways: first, using GoSam [79,80] in combination with a model file in ufo format [81], derived from an effective Lagrangian using FeynRules [82], and second analytically as a cross-check.
As we are only considering QCD corrections, the renormalisation procedure is the same as in the SM and is described in ref. [61].

Real radiation
The real corrections consist of 5-point one-loop topologies with closed top quark loops as well as tree-level diagrams, see figure 5. Both classes of diagrams have been generated with GoSam and arranged such that interferences between the two classes are properly taken into account.
In order to isolate the singularities due to unresolved radiation, we use the same framework as in ref. [61], i.e. we use the Catani-Seymour dipole formalism [83], combined with a phase space restriction parameter α as suggested in ref. [84].
The various building blocks are assembled in a C++ program and integrated over the phase space using the Vegas algorithm [85] as implemented in the Cuba library [86].

Parametrisation of the total cross section
To parametrise the deviations of the total cross section from the one in the SM, we write the LO cross section in terms of the 15 coefficients A 1 , . . . , A 15 , following refs. [31,71].
At NLO the coefficients A 1 , . . . , A 15 are modified and new terms appear. We find:

Validation of the calculation
To validate our results, we have compared the Born-improved NLO HEFT results calculated with our setup with the ones from ref. [68], where we find agreement if we use µ r = µ f = m hh and MSTW2008 [87] PDFs at LO/NLO for the LO/NLO calculation, along with the corresponding α s value. 2 We also have cross-checked the results by using two independent codes, where the only common parts are the ufo model files and the SM virtual two-loop corrections.
In addition, we have compared the leading order distributions, benchmark points and fits of the coupling coefficients in the total cross section (see eq. (2.7)) with the ones given in refs. [31,71,72]. We find agreement with ref. [31] for all A i coefficients at the 1% level. Comparing to refs. [71,72], we systematically find values that differ by 15-20% for coefficients linear in c ggh and by ∼ 40% for the coefficient quadratic in c ggh . We also compared our results with the distributions shown in refs. [71,72], finding agreement for all benchmark points except for benchmark point 8. While in refs. [71,72] a dip in the leading order distribution is found for benchmark point 8, we find no such dip. This is why we chose a different point of cluster 8 which does show a dip, and which we call 8a.

Phenomenological results
In this section we present numerical results for benchmark points which were identified in ref. [71] to represent partitions of the BSM parameter space according to characteristic shapes of differential distributions, in particular the Higgs boson pair invariant mass distributions. All our results are for a centre-of-mass energy of √ s = 14 TeV. The results were computed using the PDF4LHC15 nlo 100 pdfas [88][89][90][91] parton distribution functions interfaced via LHAPDF [92], along with the corresponding value for α s (µ), with α s (M Z ) = 0.118. The masses of the Higgs boson and the top quark have been JHEP09(2018)057 set to m h = 125 GeV and m t = 173 GeV (pole mass), respectively. The widths of the top quark (and the Higgs boson) have been set to zero. Bottom quarks are treated as massless and therefore are not included in the fermion loops. The scale uncertainties are estimated by varying the factorisation scale µ F and the renormalisation scale µ R around the central scale µ 0 = m hh /2, using the envelope of a 7-point scale variation. The latter means that we use µ R,F = c R,F µ 0 , where c R , c F ∈ {2, 1, 0.5}, and consider each combination except the two extreme ones c R = 0.5, c F = 2 and c R = 2, c F = 0.5. In the SM case, the combinations c R = c F = 0.5 and c R = c F = 2 always coincided with the envelope of the 7 combinations to vary c R , c F .

NLO cross sections and heat maps
In this section we will provide results for the coefficients defined in eqs. (2.7) and (2.8), i.e. for the expression We evaluated the coefficients in two different ways: determination via projections and performing a fit, finding agreement of the results within their uncertainties. The results of the projection method, including uncertainties, are summarised in table 1.
In the following we show heat maps for the ratio σ/σ SM , based on the results for A 1 , . . . , A 23 . For the fixed parameters the SM values are used. Further we use σ LO SM = 19.85 fb, σ NLO SM = 32.95 fb. The couplings are varied in a range which seems reasonable when taking into account the current constraints on the Higgs coupling measurements [1,93,94], as well as recent limits on the di-Higgs production cross section [95][96][97][98].
In figure 7 we display heat maps where the anomalous coupling c tt is varied in combination with the Higgs-gluon contact interactions c gghh and c ggh . We show the ratio to the SM total cross section both at LO and at NLO. We can see that the NLO corrections can lead to a significant shift in the iso-contours. It also becomes apparent that the cross sections are more sensitive to variations of c tt than to variations of the contact interaction c ggh . Figure 8 shows variations of the triple Higgs coupling c hhh in combination with c ggh and c tt . We observe that the deviations from the SM cross section can be substantial, and again we see a rapid variation of the cross section when changing c tt .
In figure 9 we display variations of c t versus c hhh , and variations of c t versus c tt . We see that values of c t around 2.0 in combination with large negative values of c hhh can enhance the cross section by two orders of magnitude. Current experimental limits suggest that the total cross section for Higgs boson pair production does not exceed about 13-24 times the SM value, assuming a SM-like shape in the distributions [96,97]. Together with the prospects that c t will be increasingly well constrained in the future, e.g. from measurements -  Table 1. Results for the coefficients defined in eq. (3.1). The uncertainties are obtained from the uncertainties on the total cross sections entering the projections, using error propagation which neglects correlations between these cross sections. of ttH production [99,100], this should allow to constrain some of the parameter space for c hhh . 3 Figure 10 shows variations of c gghh versus c ggh and c t . We observe that the impact on the NLO corrections is milder in this case.
In figure 11 we show the K-factors as a function of the coupling parameters, with the others fixed to their SM values. It shows that the rather flat K-factors which have been found [68,70] in the m t → ∞ limit (flat with respect to variations of one of the coupling parameters) show a much stronger dependence on the coupling parameters once the full top quark mass dependence is taken into account.

Cross sections and distributions at several benchmark points
In the following we will show results for the benchmark points defined in ref. [71], except for benchmark point 8, where we choose a different one (denoted as "outlier" number 5 for cluster 8 in ref. [72]) which has a more characteristic shape, and which we call 8a. The conventions for the definition of the couplings between our Lagrangian, given in eq. (2.6), and the one of ref. [71] are slightly different. In table 2 we list the conversion factors to translate between the conventions.
The benchmark points translated to our conventions are given in table 3.

Total cross sections
We first show the values for the total cross sections, together with their statistical uncertainties and the uncertainties from scale variations. We should point out that the cross sections for benchmark points B 3 , B 4 and B 12 are larger than the limits measured in the bbγγ decay channel [97,98]. However, within the same cluster [72], i.e. the set of couplings -14 -

JHEP09(2018)057
EWChL eq. (2.6) Ref. [71] c hhh κ λ c t κ t c tt c 2 c ggh  which lead to a similar shape of the m hh distribution, one can easily find combinations of couplings where the value of the total cross section is below the experimental exclusion bound. For example, taking the point c hhh = 1, c t = 1, c tt = 0, c ggh = 4/15, c gghh = −0.2 in cluster 4 leads to a cross section of about 1.8 times the SM cross section, still far from being excluded, see figure 16. The large differences in the statistical uncertainties for the different benchmark points are due to the fact that the results for the virtual two-loop part are based on rescaling of the SM numerical results, which are distributed differently in the phase space. Therefore the statistical uncertainties are largest for benchmark points where the distribution in phase space is very different from the SM case. For example, benchmark 10 has a large differential cross section at low m hh values, where the SM statistics is very low. This translates into the large statistical uncertainty for benchmark 10.

m hh and p T ,h distributions
Now we consider differential cross sections for the 12    in red, and compare it to the two approximations "Born-improved NLO HEFT" (purple) and FT approx (green). The leading order (yellow) as well as the SM results are also shown (blue NLO, black LO). The lower ratio plot shows the ratio of the two approximate results to the full NLO result. The upper ratio plot shows the differential BSM K-factor, i.e. NLO BSM /LO BSM , both evaluated with the same PDFs. cross section is about 6 times the SM cross section, and the shape of the m hh distribution is completely different from the SM. In fact, one can show analytically that the LO cross section in the m t → ∞ limit exactly vanishes near m hh = 364 GeV, which relates to the dip in the distribution. The huge enhancement at low m hh values is due to the large value of c hhh . Figure 13, corresponding to benchmark 2, shows a very different behaviour. The result is very much suppressed in the region where the SM shows a peak, while there is a large enhancement in the tail of both the m hh and the p T,h distributions. The enhancement in the tail is mainly due to the nonzero c gghh value, as the amplitude proportional to c gghh grows likeŝ [31]. We also notice that the approximations "Born-improved NLO HEFT" and FT approx cannot describe the pattern around the 2m t threshold, where the nonzero value of c tt seems to play a significant role. The K-factor for benchmark 2 is very nonhomogeneous around the dip in the m hh distribution, and can reach up to a factor of three. This is a clear example where rescaling the LO result with a K-factor obtained from higher order calculations in the HEFT approximation would lead to very different results.
Benchmark point 3, shown in figure 14, has the same values for c hhh and c t as benchmark point 2 (the SM values), but the distributions show a very different behaviour. As in the SM, there is a peak around the 2m t threshold, but the cross section is largely enhanced, not only in the peak region. As mentioned above, with a total cross section of about 32 times the SM NLO cross section, this parameter point is above the current limit deduced at 95% CL from the measured pp → HH → γγbb cross section [97,98].
Benchmark point 4, shown in figure 15, has negative values for c hhh and c tt , a slightly encreased Yukawa coupling c t , and no Higgs-gluon contact interactions. This combination removes the destructive interference between different types of diagrams present in the SM, and therefore leads to a very large cross section. The differential K-factor is about 2, as for the other benchmarks, and rather constant over the whole m hh range (whereas for benchmark 2, the differential K-factor is far from being homogeneous). Benchmark 4 is  the one with the largest cross section of all the considered benchmark points, with a total cross section of about 270 times the SM one. This point in parameter space is excluded already. Therefore, in figure 16, we also show results for another point from cluster 4, defined by c hhh = 1, c t = 1, c tt = 0, c ggh = 4/15, c gghh = −0.2, which leads to a similar shape as benchmark point 4, but to σ/σ SM = 1.8, and hence is not yet excluded. This parameter point also has the interesting feature that the distributions for NLO SM and LO BSM almost coincide. However, there is no degeneracy with the SM distribution once the NLO corrections are taken into account. Figure 17 shows distributions for benchmark point 5, where c tt is zero and c hhh and c t are as in the SM, while the Higgs-gluon interactions are nonzero. Similar to benchmark point 2, we observe a dip near m hh = 350 GeV, but the LO HEFT amplitude does not vanish there. The total cross section for benchmark point 5 is very similar to the SM one. This is an example where differential measurements are crucial to establish a clear BSM disappear completely, leading to a total cross section which is about 6.7 times larger than the one for benchmark 6, and a large enhancement of the low m hh and low p T,h regions. The distributions also show that the full top quark mass dependence in the "triangle-type" diagrams containing c hhh , which dominate the low m hh region, seems to play a significant role, as the full NLO result is quite different from the approximate results.
Benchmark point 8a, displayed in figure 20, again shows a characteristic dip just before the 2m t threshold. It is also an example where the total cross section is very similar to the SM one, but the shape of both the m hh and the p T,h distributions clearly discriminates the SM from the BSM case.
Benchmark point 9, displayed in figure 21, shows a large enhancement in the tails of the distributions, similar to benchmarks 2 and 3, which can be attributed mainly to the rather large value of c gghh , in combination with a non-zero value of c tt .
For benchmark point 10, shown in figure 22, the large value of c hhh = 10 completely dominates the shape, leading to a large enhancement in the low m hh and p T,h regions. With a value for the total cross section which is about 17 times larger than the SM cross section, benchmark point 10 is still allowed by the limits given by CMS [97], where separate limits for the various benchmark points are given.
Benchmark point 11, displayed in figure 23, has the same value for c hhh as benchmark 6, which is the one where the destructive interference would be maximal if all other couplings are kept SM-like. However, the destructive interference is compensated by the large and non-zero values of c ggh and c gghh , such that the total cross section for benchmark 11 is about 5 times larger than the SM cross section. In view of the fact that this benchmark point is dominated by the Higgs-gluon contact interactions parametrised by c ggh and c gghh , it is not a surprise that the approximations FT approx and Born-improved HEFT agree quite well with the full calculation, as all three curves have these contributions in common, while the part which differs is damped by the destructive interference.
Benchmark point 12, shown in figure 24, has all couplings SM-like except c tt = 1 and c hhh , where for the latter an extreme value of c hhh = 15 is chosen, leading to a cross section about 100 times larger than the SM cross section. This scenario is already ruled out by current LHC measurements.
All the distributions show that the NLO K-factors are large, being about a factor of two or larger. Therefore it is essential to take NLO corrections into account. The approximations where the top quark mass dependence is only partly taken into account also differ substantially in the shape from the full result for some of the benchmark points, which emphasises the importance of including the full top quark mass dependence.
In figure 25, we show the LO and NLO scale variation bands for benchmark point 5. This benchmark point is an example where the scale variation band of the 7-point scale variation mainly decreases the differential cross section over almost the whole m hh range, where the upper limit of the scale variation band is mostly given by the combination µ F = µ 0 /2, µ R = µ 0 , for some of the bins also by µ F = µ 0 , µ R = 2µ 0 . In the SM, the upper limit of the 7-point scale variation band is given by µ F = µ R = µ 0 for all bins of the m hh  distribution. We further notice that LO and NLO scale variation bands do not overlap for the m hh distribution. However, this feature is also present in the SM.

Discussion of the benchmark points
Attempting a more global view of the behaviour of the m hh distribution as a function of the five BSM parameters, we can identify the following patterns: a dip in the m hh distributions is present for benchmark points 1, 2, 5, 6 and 8a. The presence of a non-zero value for c tt or c ggh is a characteristic feature of many parameter space points that show a dip in the m hh distribution, but this is not a necessary condition for the presence of the dip. For instance, points with c hhh 2.5c t and the other couplings vanishing also show such a dip. For the subset (1, 2, 6) of the above points there is a m hh value where the LO amplitude in the m t → ∞ limit exactly vanishes, which is a feature that can cause the dip. The low m hh region is enhanced for benchmark points 1, 6, 7, 10, 11, 12, which is mainly due to the large value of c hhh , as the matrix element squared proportional to c 2 hhh for largeŝ behaves like m 2 h /ŝ log 2 m 2 t /ŝ [31] and therefore dominates at low values ofŝ. The term proportional to c 2 tt for largeŝ behaves like log 2 m 2 t /ŝ and seems to partially cancel the logarithmic terms from c hhh , such that benchmark 4 has a SM-like shape even though the absolute value for c hhh is large. The matrix element squared proportional to c gghh grows likeŝ, this is why for benchmark points which have large values of c gghh , the tail of the m hh distribution is enhanced.
In order to assess the effect of a variation of c tt while the other couplings are fixed to their SM values, we show m hh distributions for the c tt values c tt = −2, −1, 0, 1, 2 in figure 26. The minimum of the cross section is at c tt ∼ 0.25. We observe that the enhancement of the cross section as |c tt | increases is growing more rapidly for negative values of c tt , see also figure 8b. The shape changes compared to the SM are most pronounced in the low m hh region.

Conclusions
We have calculated the NLO QCD corrections with full m t dependence to Higgs boson pair production within the framework of the electroweak chiral Lagrangian, a non-linearly realised Effective Field Theory in the Higgs sector, which allows to focus on anomalous Higgs boson properties. This restricts the BSM parameter space to five possibly anomalous couplings, c hhh , c t , c tt , c ggh and c gghh .
We gave a parametrisation of the total NLO cross section and of the m hh distribution in terms of 23 coefficients of all combinations of these couplings, and also showed iso-contours of LO and NLO cross section ratios σ/σ SM for two-dimensional projections of the parameter space. These studies showed that the cross sections are very sensitive to variations of c tt , the effective tthh coupling, and that the K-factors can be large and non-uniform as the anomalous couplings are varied.

JHEP09(2018)057
We have also shown differential cross sections for m hh and p T,h at several benchmark points which exhibit characteristic shapes of the distributions. The differential K-factors for the m hh distributions are of the order of two, but can reach up to three and can be very non-uniform over the m hh range. This means that a rescaling of the LO distribution with a global K-factor can be rather misleading.
Some combinations of couplings lead to a huge enhancement of the cross section, others lead to a total cross section which is nearly degenerate to the SM one, but the corresponding m hh and p T,h distributions have a shape which is very different from the SM one, and therefore should have discriminating power even with low statistics, which emphasises the importance of measuring distributions.
Our analytical parametrisation of the total NLO cross section and of the m hh distribution in terms of all possible combinations of anomalous couplings should open the door to further studies of the considered BSM parameter space and lead to refined limits on anomalous Higgs boson couplings in the not too distant future.
to the relation in (B.4). It should, however, be kept in mind that the extraction of Higgs couplings ultimately requires a global analysis, where other Higgs-related processes are also taken into account. In this case, in particular for the important single-Higgs observables, the EWChL has overall fewer parameters to describe leading NP effects [76]. It has the advantage to focus on anomalous Higgs properties and to naturally allow for large deviations from the SM in the Higgs sector.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.