A non-linear EFT description of $gg\to HH$ at NLO interfaced to POWHEG

We present the implementation of Higgs boson pair production in gluon fusion within a non-linear Effective Field Theory framework containing five anomalous couplings for this process. The code, available within the POWHEG-BOX-V2, includes full NLO QCD corrections with massive top quarks. All five couplings can be modified by the user. We show $m_{hh}$ distributions at seven benchmark points provided by an $m_{hh}$ shape analysis at NLO and showered $p_T^{hh}$ distributions resulting from an interface to Pythia-8 and Herwig-7.2.


Introduction
Precise measurements of the Higgs boson couplings to other particles and itself are among the main goals for the next phases of LHC and beyond. As the precision of the measurements increases, it is of great importance to have Standard Model predictions well under control, and to have reliable simulations of the effects of anomalous couplings. In fact, measurements of Higgs couplings to electroweak bosons and the top quark are already reaching a level where systematic uncertainties play an increasingly important role [1,2]. The trilinear Higgs-boson self-coupling c hhh still is rather weakly constrained, however the window of possible c hhh -values has been narrowed considerably in Run II [3,4].
Higgs boson pair production in gluon fusion in the SM has been calculated at leading order in Refs. [5][6][7]. The NLO QCD corrections with full top quark mass dependence became available more recently [8][9][10][11]. The NLO results of Refs. [8,9] have been combined with parton shower Monte Carlo programs in Refs. [12][13][14], where Ref. [14] allows the trilinear Higgs coupling to be varied. Before the full NLO QCD corrections became available, the m t → ∞ limit, sometimes also called "Higgs Effective Field Theory (HEFT)" approximation or "Heavy Top Limit (HTL)", has been used. In this limit, the NLO corrections were first calculated in Ref. [15] using the so-called "Born-improved HTL", which involves rescaling the NLO results in the m t → ∞ limit by a factor B FT /B HTL , where B FT denotes the LO matrix element squared in the full theory. In Ref. [16] 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 m t → ∞ approximation. In the m t → ∞ limit, the NNLO QCD corrections have been computed in Refs. [17][18][19][20]. The calculation of Ref. [20] has been combined with results including the top quark mass dependence as far as available in Ref. [21], and soft gluon resummation on top of these results has been presented in Ref. [22]. N 3 LO corrections have become available recently [23,24], where in Ref. [24] the N 3 LO results in the heavy top limit have been "NLO-improved" using the results of Refs. [12,14]. 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. [10,11]. For a more detailed description of the various developments and phenomenological studies concerning Higgs boson pair production we refer to recent review articles, e.g. Refs. [25][26][27][28].
The main purpose of this paper is to present an update of the public Monte Carlo event generator POWHEG-BOX-V2/ggHH, where the user can choose the values of five anomalous couplings relevant to Higgs boson pair production as input parameters. It is based on the implementation of the fixed-order NLO results [8,9], combined with a non-linear Effective Field Theory framework [29], in the POWHEG-BOX [30][31][32]. It builds on the code described in Ref. [14] which allows variations of the trilinear Higgs coupling (and the top Yukawa coupling) only. We also show results for seven benchmark points, which have been identified by an m hh shape analysis presented in Ref. [33], based on the full NLO calculation, and compare NLO effects to effects from anomalous couplings. Further, we show results matched to the PYTHIA-8 [34] and HERWIG-7.2 [35] parton showers to enable the assessment of parton-shower related uncertainties. This paper is organised as follows. In Section 2 we describe the theoretical framework and the definition of the anomalous couplings. In Section 3 we describe the code and the usage of the program within the POWHEG-BOX-V2. Section 4 contains the discussion of phenomenological results, before we conclude in Section 5. More detailed usage instructions are given in an appendix.

Anomalous couplings in Higgs boson pair production
The calculation builds on the ones presented in Refs. [14,29] and therefore will be described only briefly here. We work in a non-linear EFT framework, sometimes also called Electroweak Chiral Lagrangian (EWChL) including a light Higgs boson [36,37]. It relies on counting the chiral dimension of the terms contributing to the Lagrangian [38], rather than counting the canonical dimension as in the Standard Model Effective Field Theory (SMEFT). In this way, the EWChL is also suitable for describing strong dynamics in the Higgs sector. Applying this framework to Higgs boson pair production in gluon fusion, keeping terms up to chiral dimension four, 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 leading-order diagrams are shown in Fig. 1. In Ref. [29] the NLO QCD corrections were calculated within this framework, and NLO results were presented for the twelve benchmark points defined in Ref. [39].  131.06 ± 0.08 2.28 3.98 Table 1: NLO benchmark points derived in Ref. [33]. The values for the cross section are given at √ s = 14 TeV.
In Ref. [33], shapes of the Higgs boson pair invariant mass distribution m hh were analysed in the 5-dimensional space of anomalous couplings using machine learning techniques to classify m hh -shapes, starting from NLO predictions. This method of clustering leads us to define seven new benchmark points, which we use here to discuss our phenomenological results. For convenience we repeat the benchmark points here in Table 1, together with the corresponding values for the cross section at √ s = 14 TeV. There are different normalisation conventions for the anomalous couplings in the literature. In Table 2 we summarise the conventions commonly used.
Eq. (2.1), i.e. L of Ref. [29] Ref. [39] Ref. [40] c We also give the relation to the corresponding parameters in the SMEFT, using the following Lagrangian based on the counting of canonical dimensions: where we follow the conventions used in [40,41], except forc g which differs by inclusion of the weak coupling g 2 :c g Ref. [40] = g 2c . The term proportional toc ug denotes the chromomagnetic operator, which does not contribute at the order in the chiral counting we are considering here (dχ ≤ 4), because it gets an additional loop suppression factor 1/16π 2 due to the fact that dimension-6 operators involving field strength tensors (such as σ µν G µν ) can only be generated through loop diagrams [29,42]. The remaining coefficientsc i in Eq. (2.2) can be related to the couplings of the physical Higgs field h and compared with the corresponding parameters of the chiral Lagrangian (2.1). After a field redefinition of h to eliminatec H from the kinetic term one finds [40,43] 3 Description of the code

Structure of the code
The code is an extension of the one presented in Ref. [14] to include the possibility of varying all five anomalous couplings rather than only the trilinear Higgs coupling.
For the virtual two-loop corrections, we have built on the results of the calculations presented in Refs. [8,9]. The real radiation matrix elements were implemented using the interface between GoSam [44,45] and the POWHEG-BOX [32,46]. The extra matrix elements occurring in the EFT framework have been generated by GoSam via a model file in UFO format [47] which has been developed in Ref. [29], derived from the effective Lagrangian in Eq. (2.1) using FeynRules [48]. The framework presented in Ref. [14] to interface the two-loop virtual contribution in POWHEG is generalised in the following way: instead of a second-order polynomial (as for variations of c hhh only), at NLO we can write the squared matrix element for variations of all five anomalous couplings as in Eq. (3.1), following Refs. [29,39,43].
For the Born-virtual interference term, we produce grids using 6715 points (5194 points at √ s = 14 TeV and 1521 points at 100 TeV) for 23 linearly independent sets of couplings. This enables us to derive, for each phase-space point, the coefficients a 1 , . . . , a 23 by interpolation. Once the user has chosen a set of anomalous couplings, the 23 grids are combined into one using Eq. (3.1). This step is performed only once, in the first POWHEG parallel stage. The Born and real contributions are evaluated exactly for the chosen anomalous couplings without relying on a grid or interpolation. Note that the a i coefficients of Eq. (3.1) are not equal to the A i coefficients of Ref. [29], which are derived for the (normalised) cross-section and not for the Born-virtual interference term. The original phase space points were produced with SM couplings, therefore some regions, for example the low m hh -region, are less populated than the peak of the m hhdistribution in the SM. The statistical uncertainty on the input data induces a systematic uncertainty on the 23 grids. We have checked the relative size of the uncertainties of the virtual corrections in each bin of the m hh distributions for all seven benchmark points. This uncertainty is below 2% throughout the whole m hh range, except for the first bin. This bin is poorly populated in the SM and therefore the uncertainties in this bin are larger for coupling configurations where the low-m hh region is very different from the SM case or where the relative size of the virtual contribution is large. Hence we find uncertainties in the first m hh bin of about 6% for all benchmarks except for benchmark 5, which has a 12% uncertainty in this bin.

Usage of the code
The code can be found at the web page http://powhegbox.mib.infn.it under User-Processes-V2 in the ggHH process directory. An example input card (powheg.input-save) and a run script (run.sh) are provided in the testrun folder accompanying the code. In the following we only describe the input parameters that are specific to the process gg → HH including five anomalous couplings. The parameters that are common to all POWHEG-BOX processes can be found in the POWHEG manual V2-paper.pdf in the POWHEG-BOX-V2/Docs directory.

Running modes
The code contains the SM NLO QCD amplitudes with full top quark mass dependence. A detailed description of the different approximations can be found in Ref. [9]. For the Standard Model case as well as for BSM-values of the trilinear Higgs coupling c hhh , the code can be run in four different modes, either by changing the flag mtdep in the POWHEG-BOX run card powheg.input-save, or by using the script run.sh [mtdep mode]. If all five anomalous couplings are varied, there is only the possibility of either calculating at • NLO with full top quark mass dependence, or • LO (setting bornonly=1) in either the full theory or in the m t → ∞ limit.
In more detail, the following choices are available: mtdep=0: all amplitudes are computed in the m t → ∞ limit (HTL). This option is only available at NLO in the SM case or if only c hhh is varied, or at LO.

Input parameters
The bottom quark is considered massless in all four mtdep modes. The Higgs bosons are generated on-shell with zero width. A decay can be attached through the parton shower in the narrow-width approximation. However, the decay is by default switched off (see the hdecaymode flag in the example powheg.input-save input card in testrun).

Phenomenological results
We present results calculated at a centre-of-mass energy of √ s = 14 TeV using the PDF4LHC15 nlo 30 pdfas [49][50][51][52] parton distribution functions interfaced to our code via LHAPDF [53], along with the corresponding value for α s . 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. Jets are clustered with the anti-k T algorithm [54] as implemented in the FastJet package [55,56], with jet radius R = 0.4 and a minimum transverse momentum p jet T,min = 20 GeV. The scale uncertainties are estimated by varying the factorisation/renormalisation scales µ F , µ R , where the bands represent 3-point scale variations around the central scale µ 0 = m hh /2, with µ R = µ F = c µ 0 , where c ∈ {0.5, 1, 2}. For the case c hhh = c SM hhh = 1 we checked that the bands obtained from these variations coincide with the bands resulting from 7-point scale variations. In Fig. 2a we show the Higgs boson pair invariant mass distributions for benchmark points 1, 2 and 3 compared to the SM case. The magnitudes of the cross sections are similar, due to the fact that the benchmark points were defined with the constraint that the total cross section should not exceed 6.9 × σ SM at 13 TeV [33]. Nonetheless, the curves are normalised to the SM cross section, such that only shape differences appear in the figure. Benchmark point 3 has a value of c hhh where the destructive interference between box-and triangle-type contributions is large, which leads to the dip in the m hh spectrum, while the tail is enhanced due to non-zero c ggh and c gghh values. Benchmark point 1 shows the largest enhancement of the very low m hh region, even though its value for c hhh is smaller than the one for benchmark point 2. This behaviour can be attributed to the interplay with the nonzero value of c tt , as can be concluded from the analysis in Ref. [33]. In Fig. 2b the m hh distribution for benchmark points 5, 6 and 7 is shown, normalised to the SM cross section. Benchmark point 5 shows a narrow dip below m hh = 2m t , which would not be present for c hhh = 3.95 if all other couplings were SM-like. In fact, from the analysis in Ref. [33] it can be inferred that the negative c gghh value in combination with c hhh = 3.95 is causing this dip in the shape.  In Fig. 3a we consider benchmark point 4 compared to the full NLO SM as well as in the Born-improved m t → ∞ limit, matched to PYTHIA-8 in all cases. The curves for the Born-improved HTL SM case and for benchmark point 4 are normalised to the SM cross section. Even though the m t → ∞ approximation shows an enhanced tail compared to the full SM, the enhancement of the tail in the case of benchmark 4 is much more pronounced. The situation is different for the p hh T distribution, shown in Fig. 3b. For this observable, the results for benchmark 4 and the Born-improved m t → ∞ approximation are very close. This fact again shows the importance of the full NLO corrections in order to clearly identify new physics effects. In both Fig. 3a and Fig. 3b, in order to obtain the scale uncertainty bands, the variation curves were normalised by the ratio of the central-scale prediction to the SM cross section. Thus the bands have the same relative size as in an unnormalised plot. We also investigated a different option to produce the scale bands for the normalised cross section, where the scale uncertainties are not normalised by the ratio σ SM /σ(µ 0 ), but rather by σ SM /σ(c µ 0 ), c ∈ {0.5, 1, 2}, i.e. by their own cross section at the considered scale. For the m hh distribution this type of normalisation makes the scale bands disappear within the statistical uncertainties. This is because for central scale choice µ 0 = m hh /2 the scale variations do not introduce significant shape changes for this observable. For the p hh T distribution the situation is different, firstly because the central scale choice is not aligned with the observable, and secondly because the tail of the p hh T distribution is dominated by HH+jet events, which are the leading order in this channel and therefore show larger scale uncertainties. Indeed we observe that when normalised to their own cross-section prediction, the scale uncertainty bands differ from the central prediction by ±(3-4)% for p hh T 200 GeV, and up to ∓18% at p hh T = 600 GeV. This is to be compared to an overall scale uncertainty of 30-40% at p hh T = 600 GeV with the default normalisation method. To confirm this interpretation, we also investigated scale-varied differential cross sections for the distribution of the transverse momentum of one (any) of the Higgs bosons p h T , which is an observable that is not aligned with the choice of scale µ = m hh /2, but which does get contributions from genuine radiative corrections at NLO. When normalised to their own cross section, the scale-varied predictions differ by ±1% from the central prediction at low-to-moderate p h T , and grow up to ∓(8-10)% in the tail at p h T = 600 GeV. As expected, the scale uncertainties are non-vanishing but still smaller than for the p hh T distribution. In comparison, the full (unnormalised) scale uncertainties are of the order of ±(20-25)% across the p h T range. In Fig. 4 we compare NLO predictions matched to different parton showers, namely PYTHIA-8 and two HERWIG-7.2 parton showers (the angular-orderedq and the dipole shower), to the fixed-order case, (a) for benchmark point 4, and (b) for the SM case. We observe that the enhancement of the tail with POWHEG+PYTHIA-8 present in the SM case is much less pronounced for benchmark 4, where the POWHEG+PYTHIA-8 result also touches onto the fixed-order result at large p hh T . This behaviour is most likely due to the non-zero values for c ggh and c gghh for benchmark 4, which cause the tail of the distribution already to be harder than in the SM, such that additional hard radiation created by the shower has a lower relative impact.

Conclusions
We have presented a publicly available implementation of Higgs boson pair production in gluon fusion within an Effective Field Theory framework, calculated at full NLO QCD, in the POWHEG-BOX-V2. The code allows five anomalous couplings relevant for di-Higgs production to be varied and offers the possibility to produce fully differential final states. We have also investigated the behaviour of the shape of the invariant mass distribution of the Higgs boson pair, m hh , for seven benchmark points representing characteristic m hh -shapes based on an NLO analysis [33]. In addition, for one of the benchmark points (benchmark 4), characterised by an enhanced tail in the m hh distribution, we carried out a comparison to the full NLO SM result as well as to the m t → ∞ approximation, including scale uncertainties, for both the m hh and the p hh T distributions. We found that the two distributions show different characteristics concerning the distinction of the BSM curve from the SM (and approximate SM) curves: while in the m hh distribution the enhanced tail of benchmark 4 is clearly outside the uncertainty bands of the SM predictions, in the p hh T distribution the bands for benchmark 4 and the SM in the m t → ∞ approximation overlap. This again demonstrates the importance of using full NLO predictions, particularly in trying to resolve partly degenerate directions in the space of anomalous couplings. We further produced results for the p hh T distribution matched to three different parton showers: PYTHIA-8 and two different HERWIG-7.2 parton showers. The two HERWIG-7.2 parton showers show very similar results. From previous SM results, it is known that loop-induced processes like gg → HH in POWHEG matched to PYTHIA-8 can show a harder tail than with other parton showers. However, in the case of benchmark point 4, PYTHIA-8 produces less additional hard radiation than in the SM case, such that the PYTHIA-8 and HERWIG-7.2 results are much more similar. Our studies show that the behaviour of higher-order effects known from the SM does not necessarily carry over to the BSM case, such that precise predictions for both cases are necessary to clearly identify new physics effects. With our code, fully exclusive studies of anomalous couplings in Higgs boson pair production at NLO QCD are possible.

A Appendix
A.1 Running with full top quark mass dependence (mtdep=3) In this appendix we give some further details about the running mode with full top quark mass dependence. The two-loop virtual amplitudes in the NLO calculation with full top quark mass dependence are computed via a grid which encodes the dependence of the virtual two-loop amplitude on the kinematic invariantsŝ andt [12]. We emphasize that the numerical values m H = 125 GeV and m t = 173 GeV are hardcoded in this grid and therefore should not be changed in powheg.input-save when running in the mtdep=3 mode. The grid is generated using python code and is directly interfaced to the POWHEG-BOX fortran code via a python/C API. In order for the grid to be found by the code, the files (events.cdf, createdgrid.py, Virt full *E*.grid) from the folder Virtual need to be copied into the local folder where the code is run. Instead of copying the files, we suggest to create a symbolic link to the needed files. All this is done automatically if you use the script run.sh. To do this manually: assuming the code is run from a subfolder (e.g. testrun) of the process folder, the link can be created in this subfolder as follows: ln -s ../Virtual/events.cdf events.cdf ln -s ../Virtual/creategrid.py creategrid.py for grid in ../Virtual/Virt full *E*.grid; do ln -s $grid; done Once the links are in place, the code can be run with mtdep=3 as usual. The python code creategrid.py will then combine the virtual grids generated with the 23 combinations of coupling values to produce a new file Virt full *E*.grid corresponding to the values of c hhh , c t , c tt , c ggh , c gghh defined by the user in the powheg.input-save file. The python code for the grid relies on the numpy and sympy packages, which the user should install separately. When building the ggHH process the Makefile will find the embedded python 3 library via a call to python3-config, which the user should ensure is configured correctly and points to the correct library. Note that on some systems the python/C API does not search for packages (such as numpy and sympy) in the same paths as the python executable would, the user should ensure that these packages can be found also by an embedded python program. To ensure that the linked files are found, we recommend to add the run subfolder to PYTHONPATH.

A.2 Powheg input and run scripts
The run.sh script in the testrun folder allows the different stages of POWHEG to be run easily. By typing ./run.sh without any argument a menu with the 4 mtdep running modes described above is shown. For all mtdep running modes, run.sh will make the code go through the various steps (parallel stages) of the calculation: parallelstage=1: generation of the importance sampling grid for the Monte Carlo integration; parallelstage=2: calculation of the integral for the inclusive cross section and an upper bounding function of the integrand; parallelstage=3: upper bounding factors for the generation of radiation are computed; parallelstage=4: event generation, i.e. production of pwgevents-*.lhe files.
Please note: if you use the script run.sh [mtdep], the value for mtdep given as an argument to run.sh will be used, even if you specified a different value for mtdep in powheg.input-save. After running parallelstage=4, the LHE files produced by POWHEG can be directly showered by either PYTHIA-8 or HERWIG-7.2. We provide a minimal setup for producing parton-shower matched distributions in test-pythia8, respectively test-herwig7. Both the angular-ordered and the dipole-shower implemented in HERWIG-7.2 can be used by changing the showeralg flag to either default or dipole in HerwigRun.sh. Further, we should point out that POWHEG offers the possibility to use a damping factor h = hdamp of the form [57,58] F = h 2 (p hh T ) 2 + h 2 , where p hh T is the transverse momentum of the Higgs boson pair, to limit the amount of hard radiation which is exponentiated in the Sudakov form factor. The setting F ≡ 1, corresponding to hdamp= ∞, results in quite hard tails for observables like p hh T [12,14]. Changing the damping factor F by setting the flag hdamp to some finite value in the input card softens the high transverse momentum tails. Varying hdamp allows shower uncertainties to be assessed within the POWHEG matching scheme. However, hdamp should not be so low that it starts to cut into the Sudakov regime. In fact, a too low value for hdamp could spoil the logarithmic accuracy of the prediction. For this reason we suggest not to choose values for hdamp below ∼ 200. Our default value is hdamp=250.