NLO QCD Corrections to Higgs Pair Production including Dimension-6 Operators

New Physics that becomes relevant at some high scale $\Lambda$ beyond the experimental reach, can be described in the effective theory approach by adding higher-dimensional operators to the Standard Model (SM) Lagrangian. In Higgs pair production through gluon fusion, which gives access to the trilinear Higgs self-coupling, this leads not only to modifications of the SM couplings but also induces novel couplings not present in the SM. For a proper prediction of the cross section, higher order QCD corrections that are important for this process, have to be taken into account. The various higher-dimensional contributions are affected differently by the QCD corrections. In this paper, we provide the next-to-leading order (NLO) QCD corrections to Higgs pair production including dimension-6 operators in the limit of large top quark masses. Depending on the dimension-6 coefficients entering the Lagrangian, the new operators affect the relative NLO QCD corrections by several per cent, while modifying the cross section by up to an order of magnitude.


Introduction
With the discovery of the Higgs boson [1,2], its role has developed from the long-sought particle into a tool for exploring beyond the Standard Model (SM) physics [3], possibly paving the way into new physics (NP) territory. This is the more true, as to date we are lacking any direct evidence of physics beyond the SM (BSM). The Higgs boson itself behaves SM-like with its couplings to other SM particles being close to the predicted values, in particular the couplings to gauge bosons. In some NP models, however, the trilinear Higgs self-coupling can still deviate significantly from the SM expectations [4]. A means to describe BSM physics, that is realized at a scale well above the electroweak symmetry breaking scale, in a rather model-independent way is given by the Effective Field Theory (EFT) framework. Deviations from the SM are parametrized by higher-dimensional operators, which lead to modifications of the Higgs boson couplings to the other SM particles and to itself.
The trilinear Higgs self-coupling is accessible in double Higgs production [5][6][7][8], with the dominant production mechanism at the LHC provided by gluon fusion [9,10]. The leading order process is mediated by top and bottom quark triangle and box diagrams. As in single Higgs production [11,12], the next-to-leading order (NLO) QCD corrections to this process are important. They have first been obtained in the limit of large top quark masses [5]. While this approximation works quite well for single Higgs production, the uncertainties of the approximation are considerably larger for double Higgs production and even more in the case of differential distributions [13,14]. Top quark mass effects have been analysed in [15][16][17][18], and first results towards a fully differential NLO calculation have been presented in [17]. Recently, the next-to-next-to-leading order (NNLO) QCD corrections have been calculated in [19][20][21]. The authors of [22] have performed a soft gluon resummation at next-to-next-to-leading logarithmic order within the SCET approach.
Higher-dimensional operators relevant for Higgs pair production through gluon fusion have been discussed in [4,[23][24][25][26]. They lead to contributions that are different for the triangle and box diagrams mediating the pair production process. In this work we perform the computation of the NLO QCD corrections to gluon fusion into Higgs pairs including higher dimensional operators in the large top mass limit. Our result allows us to investigate the validity of an approximation applied in previous works. This approximation relies on the multiplication of the full leading order (LO) cross section by an overall K-factor, given by the ratio of the SM result for the NLO QCD cross section divided by the LO cross section, in the large top mass limit.
In the next section 2 we present the details of our calculation. This is followed by a numerical analysis in section 3. In section 4 we summarize and conclude.

Details of the calculation
Gluon fusion into Higgs pairs is mediated by top and bottom quark loops dominantly [9]. We compute the NLO QCD corrections in the heavy top quark limit and we neglect in the following in this framework the bottom quark loops, which only contribute with less than 1% [5,10].
If physics beyond the SM appears at some high-scale, NP effects can be parametrized in a rather model-independent way by introducing higher-dimensional operators. In case the Higgs boson is embedded in an SU (2)-doublet H the leading BSM effects are given by dimension-6 operators. 1 In the Strongly-Interacting-Light Higgs (SILH) basis the operators relevant for Higgs pair production are given by [27], where v is the vacuum expectation value v ≈ 246 GeV, M h = 125 GeV the Higgs boson mass, m W the W boson mass, y t the top Yukawa coupling constant, g s the strong coupling constant and G a µν the gluon field strength tensor. Note that we neglect CP-violating effects. An estimate of the size of the coefficientsc H ,c u ,c 6 andc g and the most important experimental bounds can be found in [28]. The first three operators in Eq. (2.1) modify the top Yukawa and the trilinear Higgs self-coupling with respect to the corresponding SM values, while the last operator parametrizes effective gluon couplings to one and two Higgs bosons not mediated by SM quark loops. The second operator furthermore introduces a novel two-Higgs two-fermion coupling [29].
In case the SU (2) L × U (1) Y is non-linearly realized and the physical Higgs boson h is a singlet of the custodial symmetry and not necessarily part of a weak doublet, the contributions relevant for our process are summarized by the non-linear Lagrangian [30] with α s = g 2 s /(4π). In contrast to the SILH parametrization, where the coupling deviations from the SM are required to be small, the non-linear Lagrangian is valid for arbitrary values of the couplings c i . From the SILH Lagrangian in the unitary gauge and after canonical normalization the relations between the SILH coefficients and the non-linear coefficients c i can be derived, leading to [4] with α 2 = √ 2G F m 2 W /π and G F denoting the Fermi constant. In the following we will give results for the non-linear parametrization and defer the SILH case to Appendix A.
In the low-energy limit of small Higgs four-momentum an effective Lagrangian valid for light Higgs bosons can be derived for the Higgs boson interactions. The effective Lagrangian can be used to compute the QCD corrections to Higgs pair production in the large top mass limit. From single-Higgs production it is known that the K-factor obtained in this limit approximates the exact value to better than 5%, when the full mass dependence is included in the LO cross section [11]. The low-energy approach has also been used to derive the QCD corrections to Higgs pair production [5]. Here the uncertainty of 10% induced in the K-factor is worse than in the single Higgs case [15][16][17][18]. 2 The effective Lagrangian for the Higgs couplings to gluons and quarks reads [5] The factor (1 + 11/4 α s /π) emerges from the matching of the effective to the full theory at NLO QCD. The Feynman rules for the effective couplings of two gluons to one and two Higgs bosons, based on the low-energy theorems [31,32], are given in Fig. 1. The generic diagrams contributing to gluon fusion into Higgs pairs at LO are depicted in Fig. 2. The LO partonic cross section can generically be cast into the form where µ R denotes the renormalization scale. The Mandelstam variables are given bŷ in terms of the scattering angle θ in the partonic center-of-mass (c.m.) system with the invariant Higgs pair mass Q and the relative velocity The integration limits at cos θ = ±1 arê (2.8) The form factors F ∆ and F 2 in F 1 and F 2 defined as contain the full quark mass dependence and can be found in [9]. In the heavy quark limit the form factors F ∆ , F 2 and G 2 approach and F 1 and F 2 simplify to We have furthermore introduced the abbreviations with (2.13) The terms proportional to c t , respectively c 2 t in F 1 and F 2 and in front of the form factor G 2 are the usual SM contributions including the modifications due to the rescaling of the top Yukawa coupling by c t . The contributions coming with c ∆ and c 2 originate from the effective two-gluon couplings to one and two Higgs bosons, while the term involving c tt is due to the novel two-Higgs two-top quark coupling.
We use the effective couplings to compute the NLO QCD corrections to Higgs pair production. They are composed of the virtual and the real corrections. Sample diagrams are shown (2.14) We obtain for the individual contributions of Eq. (2.14) where s denotes the hadronic c.m. energy and The Altarelli-Parisi splitting functions are given by [33], with N F = 5 in our case. We denote the factorization scale of the parton-parton luminosities dL ij /dτ by µ F . While the relative real corrections are not affected by the higher-dimensional operators, the virtual corrections are altered compared to the SM case because of the overall coupling modifications of the top Yukawa and the trilinear Higgs self-coupling and due to the additional contributions from the novel effective vertices. The coefficient C for the virtual corrections reads with the squared transverse momentum and the coefficients The third line in Eq. (2.22) and the terms proportional to c 1 in the second line originate from the third diagram with the two effective Higgs-two-gluon couplings in Fig. 3 (upper), and the remaining terms are due to the diagrams with gluon loops in Fig. 3 (upper). In the derivation of the coefficient C for the virtual corrections we have kept the full top quark mass dependence in the LO amplitude. The SM result is recovered for

Numerical Analysis
For the numerical analysis we have implemented the LO and NLO Higgs pair production cross sections including the contributions of dimension-6 operators, as presented in the previous section, in the Fortran program HPAIR [34]. We have chosen the c.m. energy √ s = 14 TeV and for comparison also the very high energy option √ s = 100 TeV. The Higgs boson mass has been set equal to M h = 125 GeV [35] and for the top and bottom quark masses we have chosen m t = 173.2 GeV and m b = 4.75 GeV, respectively. We have adopted the MSTW08 [36] parton densities for the LO 3 and NLO cross sections with α s (M Z ) = 0.13939 at LO and α s (M Z ) = 0.12018 at NLO.
In order to study the impact of the new couplings on the QCD corrections we show in Fig. 4 the K-factor, defined as the ratio of the NLO and LO cross sections, K = σ NLO /σ LO , where the parton densities and the strong couplings α s are taken at NLO and LO, respectively. Deviations with respect to the SM K-factor arise in the virtual corrections due to the second and third line in the formula Eq. (2.22) for the coefficient C. Additionally, the real corrections are affected because of the different weights in the τ integration due to the modified LO cross section. In Fig. 4 (upper) we have set all couplings but c g to their SM values, c 3 = c t = 1 and c tt = c gg = 0. We have varied c g away from its SM value c g = 0 in the range −0.15 ≤ c g ≤ 0.15. For illustrative purposes we have chosen a rather large range, that goes beyond current experimental limits obtained under certain assumptions [4,37]. In the lower plot we have instead set c g = 0 and varied c gg away from its SM value c gg = 0 in the range −0.15 ≤ c gg ≤ 0.15 [4]. As can be inferred from Fig. 4 the new contact interactions c g and c gg each vary the K-factor between 1.8 and 1.9 in the investigated range. Figure 4 (lower) shows that the maximal deviation from the SM K-factor is obtained for c gg = −0.15, where the K-factor practically becomes constant. It amounts to The impact on the total cross section, however, is much larger. Here we have at NLO max|σ cgg − σ SM |/σ SM = 5.8. While the effect of higher-dimensional operators on the K-factor is at the level of a few per cent, on the cross section itself it is enormous. In Fig. 4 (upper) the maximum is less pronounced, as c g modifies the coupling of a single Higgs boson to two gluons, which is attached to the Higgs propagator. At the c.m. energy of √ s = 100 TeV (not shown here) the K-factor is smaller and varies for non-zero c gg in the range 1.46-1.58. For c g the impact at very high energy is smaller, i.e. a per cent effect.
In Fig. 5 we have set all couplings to their SM values but allowed for the new contact interaction between two Higgs bosons and two top quarks parametrized by c tt , which we have varied between -1.5 and 1.5 [4]. The maximum K-factor is reached for c tt ≈ 0.7 where the LO cross section is minimized and we have The maximum deviation is smaller than for c gg . Note, furthermore, that the c tt value, for which the deviation is maximal, is much larger than for the c g and c gg variations discussed before. This is due to the normalization of these coupling factors as can be inferred from the Feynman rules in Fig. 1. At √ s = 100 TeV the K-factor varies between 1.49 and 1.59.
The value of the trilinear coupling is practically not constrained and we have allowed in Fig. 6 for a variation of c 3 between -10 and 10 [4], while setting all other couplings to their SM      values. The impact of c 3 is rather small. In the investigated range the maximum deviation of is reached for 5 < ∼ c 3 < ∼ 10. At √ s = 100 TeV the K-factor varies between about 1.42 and 1.51.
Finally, the variation of c t in the still allowed range 0.65 ≤ c t ≤ 1. 15 [37], while keeping all other couplings at their SM values, hardly changes the K-factor, and we do not show the corresponding results separately. In this parameter configuration c t enters in the virtual correction factor C, Eq. (2.22), through the terms proportional to c 1 and c 2 , and its effect in the nominator almost cancels against the one in the denominator.
The above discussion shows, that the new couplings affect the K-factor by only a few per cent. It has to be kept in mind though that we varied the couplings only one by one away from their SM values. The combination of various new couplings could have a larger impact on the NLO corrections.

Conclusions
We have computed the NLO QCD corrections to Higgs pair production including dimension-6 operators in the large top mass limit. The dimension-6 operators lead not only to coupling modifications of the SM Higgs couplings, but also induce effective gluon couplings to one and two Higgs bosons and a novel two-Higgs two-top quark coupling. The various contributions to Higgs pair production are affected differently by the QCD corrections. Depending on the relative size of the NP contributions, the K-factor is changed by several per cent in the parameter regions compatible with the LHC Higgs data. This small impact on the K-factor underlines the dominance of soft and collinear gluon effects in the QCD corrections. The inclusion of the QCD corrections in the gluon fusion process based on the effective theory approach to describe NP, is necessary for reliable predictions of the cross section.