Combining chromomagnetic and four-fermion operators with leading SMEFT operators for gg → hh at NLO QCD

: We present the calculation of the contribtuions from the chromomagnetic and four-top-quark-operators within Standard Model Effective Field Theory (SMEFT) to Higgs boson pair production in gluon fusion, combined with QCD corrections that are at NLO with full m t -dependence for the leading operators. We study the effects of these operators on the total cross section and the invariant mass distribution of the Higgs-boson pair, at √ s = 13 . 6 TeV. These subleading operators are implemented in the generator ggHH SMEFT , in the same Powheg-Box-V2 framework as the leading operators, such that their effects can be easily studied in a unified setup.


Introduction
Where is New Physics?If it resides at energy scales well separated from the electroweak scale, our ignorance about its exact nature can be parametrised within an Effective Field Theory (EFT) framework [1][2][3][4].
Here we will focus on Higgs boson pair production in gluon fusion, combining the NLO QCD corrections with full top quark mass dependence with anomalous couplings within SMEFT.The full NLO QCD corrections have been calculated in Refs.[38][39][40][41], based on numerical evaluations of the two-loop integrals entering the virtual corrections.The results of [38] have been implemented into the Powheg-Box-V2 event generator [42][43][44], first for the SM only [45], then also for κ λ variations [46] as well as for the leading operators contributing to this process in non-linear EFT (HEFT) [47,48] and SMEFT [13].Recently, the NLO QCD corrections obtained from the combination of a p T -expansion and an expansion in the high-energy regime have been calculated analytically and implemented in the Powheg-Box-V2 [49], allowing to study top mass scheme uncertainties in an event generator framework.
In Ref. [50] the combination of NNLO corrections in an m t -improved heavy top limit (HTL) has been performed including anomalous couplings, extending earlier work at NLO in the m t -improved HTL [51,52].The work of [50] has been combined with the full NLO corrections within non-linear EFT of Ref. [47] to provide approximate NNLO predictions in Ref. [53], dubbed NNLO ′ , which include the full top-quark mass dependence up to NLO and higher order corrections up to NNLO in the m t -improved HTL, combined with operators related to the five most relevant anomalous couplings for the process gg → hh.Recently, full NLO electroweak corrections have been computed in Ref. [54], following the emergence of previous partial results, i.e. the full NLO electroweak corrections in the large-m t limit [55], the NLO Yukawa corrections in the high-energy limit [56] and Yukawa corrections in the (partial) large-m t limit [57].
In this paper, we investigate the effect of two classes of operators contributing at dimension-6 level to the process gg → hh, which however are suppressed by loop factors compared to the leading operators considered in Ref. [13] when the potential UV completion is assumed to be a weakly coupling and renormalisable quantum field theory.These are the chromomagnetic operator and 4-top-operators.As has been shown in Ref. [58] for the case of single Higgs production, the latter are intricately related since they are individually γ 5 -scheme dependent, the scheme dependence only dropping out when they are consistently combined in a renormalised amplitude.Apart from the γ 5 continuation scheme, other sources of scheme differences in bottom-up SMEFT calculations also have been studied recently [29,59,60].
The subsequent sections are organised as follows: in Section 2, we describe these contributions and their scheme dependence in detail.Their implementation into the POWHEG ggHH SMEFT generator is described in Section 3, together with instructions for the user how to turn them on or off.Section 4 contains our phenomenological results, focusing on the effects of these newly included operators on the total cross section and on the Higgs boson pair invariant mass distribution, before we summarise and conclude.

Contributions of the chromomagnetic and four-top operators
In this section we describe our selection of contributing operators.Subsequently we recapitulate the power counting scheme for SMEFT and discuss the new contributions in detail, which will be identified as subleading based on the minimal assumption of a weakly coupling and renormalisable UV theory.Any bottom-up EFT is defined by its degrees of freedom, the imposed symmetries and a power counting scheme.Since SMEFT builds upon the SM, the above specifications are given by the field content and gauge symmetries of the SM and the main power counting, which relies on the counting of the canonical (mass) dimension.Due to strong experimental constraints it is common to exclude baryon and lepton number violating operators, hence only operators of even dimension are considered.Therefore, the dominant contributions are expected to be described by dimension-6 operators, on which we focus our attention in this paper.To further cut down the number of operators, 1we impose an exact flavour symmetry U (2) q × U (2) u × U (3) d in the quark sector for a first investigation, which forbids chirality flipping bilinears involving light quarks (bquarks included) and right-handed charged currents [8,61,62].This effectively makes the CKM matrix diagonal and sets all fermion masses and Yukawa couplings to zero, with the top quark as the only exception, thus being well compatible with a 5-flavour scheme in QCD which we employ.In addition, this flavour choice reflects the expected prominent role of the top quark in many BSM scenarios and could be a starting point for a spurion expansion as in minimal flavour violation [61,63].
We also neglect operators whose contributions involve only diagrams with electroweak particles propagating in the loop.In principle, electroweak corrections and such electroweak-like operator contributions can be of the same order in the power counting as the subleading contributions studied in this paper.In addition, the close connection between operators of class ψ 2 ϕ 2 D of Ref. [5] and C tG , observed by the structure of the γ 5 -scheme dependence in Ref. [58], demonstrates that our subset does not fully comprise a consistent subleading order in a systematic power counting.Nevertheless, we expect it to be useful to investigate the sensitivity of the process gg → hh to the chromomagnetic operator and 4-top operators in the presented form, especially since even in the simpler case of the SM, full electroweak effects to gg → hh have only become available very recently [54].
With these restrictions, all dimension-6 CP even operators that contribute to gg → hh are given by where σ µν = i 2 [γ µ , γ ν ] and φ = iσ 2 ϕ is the charge conjugate of the Higgs doublet.For the covariant derivative, we use the sign convention in order to be compatible with FeynRules [64,65] conventions and tools relying on UFO [66,67] models.The first two lines in Eq. (2.1) comprise the leading EFT contribution which has been studied in Ref. [13].For convenience of the reader and later reference, we show the Born-level diagrams related to those operators in Fig. 1 The chromomagnetic operator and the 4-top operators of Eq. (2.1) together form the subleading contribution that will be the focus of this work.Below the scale of electroweak symmetry breaking, and after performing a field redefinition for the physical Higgs field in unitary gauge [13], the relevant interaction terms of the Lagrangian have the form ) which is valid up to O(Λ −4 ) differences.Here v denotes the full vacuum expectation value including a higher dimensional contribution of C H Λ 2 and 3 where y t is the top-Yukawa parameter of the dimension-4 Lagrangian.
In the following, we will briefly comment on the notions of 'leading' and 'subleading' we have used above.In SMEFT, the operators are ordered by their canonical dimension, i.e. the expansion is based on powers in E/Λ.However, in a perturbative expansion, in particular in the combination of EFT expansions with expansions in a SM coupling, loop suppression factors also play a role.Therefore, a classification of operators into potentially tree-level induced and necessarily loop-generated operators [69], the latter thus carrying an implicit loop factor L = (16π 2 ) −1 , leads to a more refined counting scheme, corroborated by observations from renormalisation and the cancellation of scheme-dependent terms [58].The same loop factors can be derived by supplementing the SMEFT expansion by a chiral counting of operators [70], see also [71].Such a classification can only be made when making some minimal UV assumptions, which are however quite generic, assuming renormalisability and weak coupling of the underlying UV complete theory 4 .Therefore, under these assumptions and if the Wilson coefficients C i in the SMEFT expansion are considered to be of similar magnitude, it makes sense to expand in C i × 1/Λ a × 1/(16π 2 ) b .Fixing a = 2 (dimension-6 operators) we call the operator contributions with b = 0 'leading' and those with b > 0 'subleading'.The above factors are to be combined with explicit loop factors 1/(16π 2 ) c from the SM perturbative expansion.Nonetheless, one has to keep in mind that this approach does not cover UV effects in full generality and that this classification is not invariant under field redefinitions and thus necessarily basis dependent [3].
Applying those rules to the Born contributions of Fig. 1 and collecting loop factors of QCD origin together with associated powers of g s leads to M Born ∼ O ((g 2 s L)Λ −2 ).Here we identify both types of contributions: explicit diagrammatic loop factors combined with tree-generated operator insertions (first line, grey dots, b = 0, c = 1 in the above classification), and tree diagrams combined with implicitly loop-generated operators (second line, grey squares, b = 1, c = 0 in the above classification).The power counting of the subleading contributions is addressed in Sections 2.1 and 2.2.
At cross section level, we therefore have where and contain the contributions with a single insertion of C tG and/or 4-top operators.Values inside the square brackets in Eqs.(2.7) and (2.8) denote the order in power counting of the respective contribution at cross section level.
In the subsequent parts of this section, we discuss the structure of the contributions to the amplitude which involve single insertions of the chromomagnetic operator and the 4-top operators of eq.(2.1).All relevant diagrams were generated with QGraf [72] and the calculation was performed analytically using FeynCalc [73][74][75].UV divergences are absorbed in a mixed on-shell-MS renormalisation scheme, where the mass of the topquark is renormalised on-shell and the dimension-6 Wilson coefficients are renormalised in the MS scheme.The contribution of the chromomagnetic operator has been checked against a private version of GoSam [76,77].Moreover, we compared with the results of Ref. [78] for the total cross section evaluated at the central scale and, after adjustment to the described conventions, found complete agreement.The amplitude involving 4top operators has been checked in D dimensions against alibrary [79] in combination with Kira [80,81].The renormalised 4-top amplitudes were tested numerically in four dimensions by comparing the analytic implementation in the Powheg-Box-V2 [42][43][44] against the result obtained with alibrary and evaluated with pySecDec [82][83][84] for several phase-space points.The chiral structure of the 4-top couplings is treated in the Naive Dimensional Regularisation (NDR) scheme [85] assuming the cyclicity of traces of strings of gamma matrices.This is possible since (after reduction of loop integrals onto the integral basis of 't Hooft-Passarino-Veltman scalar integrals [86,87]) all appearing traces with an odd number of γ 5 matrices can be explicitly brought into the form γ µ 1 . . .γ µn γ 5 with n < 4 through anti-commutation and therefore vanish.In addition, the analytic calculation of the 4-top contributions in FeynCalc is repeated in the Breitenlohner-Maison-t'Hooft-Veltman (BMHV) scheme [88,89], with the symmetric definition for chiral vertices and the translation between the Lagrangian parameters obtained in Ref. [58] is verified.For convenience, the explicit form of the translation is also presented in Eq. (2.22).

Amplitude structure of chromomagnetic operator insertions
The contribution of the chromomagnetic operator to the amplitude leads to the diagram types shown in Fig. 2. At first sight, the diagrams are at one-loop order, such that, together with the explicit dimensional factor, the prefactor of the Wilson coefficient appears at O ((g 2 s L)Λ −2 ).However, the chromomagnetic operator belongs to the class of operators that, in renormalisable UV completions, can only be generated at loop level [69,70].Hence, the implicit loop factor of its Wilson coefficient promotes the order in power counting to , which is in that sense subleading with regards to the leading Born diagrams of Fig. 1.
The diagrams of type (a), (b) and (d) are UV divergent even though they constitute the leading order contribution of C tG to the gluon fusion process.However, this behaviour is well known [78,90,91] and leads to a renormalisation of C 0 HG = µ 2ϵ C HG + δ C i C HG (µ being the renormalisation scale) which in the MS scheme takes the form [17,78,91] With this renormalisation term the finiteness of the amplitude is restored, and it can be numerically evaluated using standard integral libraries.

Amplitude structure involving four-top operators
Four-top operators appear first at two-loop order in gluon-fusion Higgs-or di-Higgs production.Thus, their contribution is of the same order in the power counting as the one of the chromomagnetic operator, i.e.
. Following the reasoning of Ref. [92] in single Higgs production, we separate the contribution into different diagram classes, which are shown in Fig. 3.The ordering in columns is chosen in order to group in underlying Born topologies (i.e.triangles and boxes), the rows combine the type of one-loop correction (if applicable).The first column is thus analogous to single Higgs production as in Ref. [92], with one Higgs splitting into two, however we do not include bottom quark loops (and loops of other light quarks), since we apply a more restrictive flavour assumption in which the bottom quark remains massless and diagrams with bottom loops vanish in an explicit calculation, either due to the bottom-Yukawa coupling being zero or due to vanishing scaleless integrals.
The categories of diagrams in Fig. 3 can be structured in the following way: (a) and (b): loop corrections to top propagators, (c) and (d): loop corrections to the Yukawa interaction, (e): loop correction to the tthh vertex, (f) and (g): loop corrections to the gauge interaction (more precisely, a contraction of a one-loop subdiagram of (f) leads to the topologies of Fig. 2 (a) or (b)), and (h) without clear correspondence to a vertex correction of a Born structure (but related to type (d) diagrams of Fig. 2 after contraction of a one-loop subdiagram).
In the following we sketch the calculation of the contribution of those classes and then refer to the γ 5 -scheme dependence of the calculation, which first has been investigated in Ref. [58].We represent the results in terms of master integrals that are given by Passarino-Veltman scalar functions N 0 , N ∈ {A, B, C, . ..} in the convention of FeynCalc [73][74][75] (which is equivalent to the LoopTools [93] convention), such that loop factors are kept manifest in the formulas.
We begin with propagator corrections which have no momentum dependence and therefore contribute only proportional to a mass insertion Hence, after applying an on-shell renormalisation of the top quark mass m 0 t = m t + δm t with δm 4-top t = −m t C (1) the diagrams of class (a) and (b) are completely removed.Next, we consider loop corrections to Yukawa-type interactions.The explicit expression for h → tt for an off-shell Higgs is proportional to the SM Yukawa coupling where q denotes the momentum of the Higgs.The part involving the one-loop tadpole integral in Eq. (2.13) is expressed in terms of the on-shell mass counter term δm 4-top t such that the effect of on-shell m t renormalisation on the correction of the Yukawa interaction is made obvious.In order to derive the necessary counter term for C tH , it is sufficient to consider the case of the Higgs being on-shell.Renormalising in the MS scheme then leads to which coincides with δ 2 C j using the respective part of the anomalous dimension matrix γ C i ,C j of Refs.[15,16]. 6With the additional counter term diagrams of δm 4-top (1) (1) ) where and M gg→h SM and M gg→hh □,SM denote the SM gg → h amplitude and the SM box-type contribution to the gg → hh amplitude, respectively.Subsequently, we investigate contributions to the gauge interaction, as they appear in diagram classes (d) and (e) of Fig. 3.It is sufficient to consider the case of an on-shell external gluon.Thus, the vertex correction evaluates to (1) where we defined Since the Lorentz structure of the correction to the gauge vertex is similar to the insertion of a chromomagnetic operator, diagrams in class (d) of Fig. 3 acquire a UV divergence (class (e) remains finite) which, analogous to the case of the chromomagnetic operator, can be absorbed by a (now two-loop) counter term of C HG .In MS the explicit form is Schematically, we now have (1) where M (1) where M 4-top ∆QQ,tt,( 8) is a remaining amplitude piece for which we could not identify an expression in terms of a one-loop subamplitude.Note that Eq. (2.21) is the only appearance of a non-vanishing contribution of the operators in the class (LL)(LL) and (RR)(RR) of Ref. [5] with coefficients C (1) QQ and C tt .Evaluating the bubbles in Eqs.(2.13) and (2.17) (for on-shell gluons) without attaching the 4-top vertex only leads to scalar respective rank-2 tensor structures in Dirac space and therefore induces a chirality flip, which is incompatible with a 4-top interaction of the same chirality in both currents.Similarly, the tadpole in Eq. (2.11) and the triangle with two Higgs bosons attached in Eq. (2.15) have a scalar structure.The triangles of Eq. (2.21), each with one external gluon and one Higgs boson attached, are the only exception, since they also have contributions proportional to a single Dirac matrix.These parts lead to the combination C (1) QQ for the single trace contraction and in addition allow a contribution with two trace contractions involving the octet operators, which leads to the combination T F (C Qt ), both multiplying M 4-top ∆QQ,tt, (8) in Eq. (2.21).
A few comments about the difference between the NDR and BMHV schemes are in order.In our calculation, the treatment of γ 5 in the two schemes differs only by the 2ϵ-dimensional part of the Dirac algebra in D-dimensions.In the limit D → 4 the renormalised fixed order result between the two schemes therefore differs by terms stemming from the 2ϵ-dimensional parts of the Dirac algebra multiplying a pole of the loop integrals.In the 4-top calculation of this work, the BMHV results are obtained by removing the finite pieces in Eqs.(2.11), (2.12), (2.13) and (2.16) that do not multiply a Passarino-Veltman scalar function, i.e. the rational parts, and setting K tG = 0 in Eqs.(2.17), (2.20) and (2.21).These differences only affect the terms dependent on C (1) Qt and C (8) Qt .This scheme dependence has the same structure as the one in the process gg → h which was observed in Ref. [58], thus the scheme dependent amplitude structure of C (1)

Qt and C (8)
Qt is compensated by scheme dependent values for the other parameters of the Lagrangian, resulting in an overall scheme independence of the EFT prediction.The γ 5 schemes hence represent equivalent parameterisations of the new physics effects and a translation between the two schemes can be achieved by means of finite shifts of the Lagrangian parameters.The explicit form of the translation relation between the NDR and the BMHV scheme in terms of parameter shifts is as follows (1) which is equivalent to the relations presented in Eqs. ( 45)-( 47) of Ref. [58]. 7These relations are best understood in a top-down perspective: in an explicit matching calculation, different choices of γ 5 schemes naturally lead to relations like Eq. (2.22).Moreover, Eq. (2.22) describes a mutual relation.One could define parameter combinations in which the scheme dependence is absorbed, however this would require to define a 'canonical scheme' instead of using an intrinsically scheme independent form.In order to avoid such an arbitrary choice in physical predictions of the EFT, simultaneous contributions of several Wilson coefficients which allow to disentangle the scheme dependence at a given order, together with a documentation of the chosen scheme, would be necessary.
3 Implementation and usage of the code within the Powheg-Box The analytic formulas of the previous section are implemented as an extension to ggHH SMEFT [13] that already includes the combination of NLO QCD corrections with the leading operators and is publicly available in the framework of the POWHEG-BOX-V2 [42][43][44].Therefore, the calculation of the cross section at fixed order is extended by the subleading contributions in the form of Eqs.(2.6)-(2.8).
The subleading contributions enter the calculation as part of the Born contribution.Since the loop functions are expressed in terms of one-loop integrals, the evaluation time per phase-space point of the subleading contributions is of the order of the existing Born contribution, thus does not significantly change the run-time of the code.
The usage of the program ggHH SMEFT follows the existing version with the extension by a few parameters in the input card.An example is given in the folder testrun in the input card powheg.input-save.The new Wilson coefficients of the subleading operators in Eq. (2.1) can be set with: Qt , CQt8 : Wilson coefficient of 4-top operator C Qt , CQQtt : sum of Wilson coefficients of 4-top operators C (1) QQ .
The available options for the selection of cross section contributions from EFT operators are visualized in Table 1.
The structure of the code still allows the user truncation (a) (b) Table 1: Options to select EFT contributions for the calculation of the cross section.Columns denote the truncation options for the 1/Λ-expansion, rows show the selection of subleading operator contributions for the Born cross section in the upper part and the NLO cross section in the lower part which is untouched by the setting of includesubleading.The partial cross section contributions are understood to be added to the SM, a higher setting for the selection always includes the previous contributions as well.Note that includesubleading=2 requires the bornonly mode.
to choose all truncation options described in Ref. [13].However, including the subleading contributions, only options (a) (SM+linear dimension-6) and (b) (SM+linear dimension-6+quadratic dimension-6) are available, as the other options are not meaningful in combination with the subleading operators.The subleading contributions are activated through the keyword includesubleading which can be set to 0, 1 or 2. When includesubleading=0 the subleading contributions are not included and the program behaves as the previous ggHH SMEFT version, i.e. the values for CtG, CQt, CQt8, CQQtt and CQQ8 are ignored.With includesubleading=1 the subleading contributions enter -according to the power counting -only in the interference with the leading LO matrix elements.The setting includesubleading=2 is only available in bornonly mode.This allows the user to remain completely agnostic about possible UV extensions such that C tG is treated as if it was part of the leading operator contribution, i.e. allowing squared C tG -contributions to |M dim-6 | 2 in truncation option (b).However, no NLO QCD corrections to the squared C tG -part are available.
In addition, there is an option for 4-top contributions to choose between the NDR scheme (GAMMA5BMHV=0) and the BMHV scheme (GAMMA5BMHV=1) with the definition of chiral vertices according to Eq. (2.9).As described at the end of Section 2.2, this will only affect the dependence on CQt and CQt8.

Results
The results presented in the following were obtained for a centre-of-mass energy of √ s = 13.6 TeV using the PDF4LHC15 nlo 30 pdfas [94] parton distribution functions, interfaced to our code via LHAPDF [95], along with the corresponding value for α s .We used m h = 125 GeV for the mass of the Higgs boson; the top quark mass has been fixed to m t = 173 GeV to be coherent with the virtual two-loop amplitude calculated numerically, and the top quark and Higgs widths have been set to zero.Jets are clustered with the anti-k T algorithm [96] as implemented in the FastJet package [97,98], with jet radius R = 0.4 and a minimum transverse momentum p jet T,min = 20 GeV.We set the central renormalisation and factorisation scales to µ R = µ F = m hh /2.We use 3-point scale variations unless specified otherwise.

Total cross sections and heat maps
In this subsection we investigate the dependence of the total cross section on the contribution of subleading operators.Following the decomposition of the cross section in Eqs.(2.6), (2.7) and (2.8), these contributions only enter linearly in interference terms; we postpone the discussion of quadratic contributions from C tG to Section 4.3.The first part demonstrates the effect of variations of pairs of Wilson coefficients with respect to the SM configuration, where all contributions are included at LO QCD.In the second part, we present values for the total cross section of the SM and benchmark point 6 of Refs.[13,99] at NLO QCD and their dependence on variations of a single subleading Wilson coefficient.The definition of benchmark point 6 in terms of SMEFT Wilson coefficients is given in Table 2.The ranges for the variation of C H are oriented at a translation of the limits on κ λ from Ref. [101], the ranges for the other Wilson coefficients are taken from Ref. [62] based on O(Λ −2 ) individual bounds or O(Λ −2 ) marginalised fits over the other Wilson coefficients.Meanwhile, constraints on 4-fermion operators in the 3rd generation also have been derived from the measurement of 4-top-quark production [102,103], based on fits varying each Wilson coefficient individually, however we use the more conservative ranges here.Note that, besides a flavour assumption, no a priori assumptions on the Wilson coefficients were made for the derivation of those limits, such that their ranges include values where the truncation at O(Λ −2 ) and/or our power counting may not be valid, i.e. the value of C tG is not suppressed by a factor of (16π 2 ) −1 and the ranges for the 4-top Wilson coefficients, Table 2: Definition of benchmark scenarios considered here in terms of SMEFT Wilson coefficients.Benchmark point 6 refers to the set in Refs.[13,99], which is an updated version of Ref. [100].The benchmarks were originally derived in a non-linear theory (HEFT), where benchmark point 6 corresponds to c hhh = −0.684,c tth = 0.9, c tthh = − 1 6 , c ggh = 0.5, c gghh = 0.25.A value of Λ = 1 TeV is assumed for the translation between HEFT and SMEFT coefficients and C HG is determined using α s (m Z ) = 0.118.with values O(100), may be too large. 8The presented results using these ranges from marginalised fits should not be understood as predictions motivated by realistic UV effects, but rather investigate the potential for improvement in global fits, as the process gg → hh (and also gg → h) probe directions that are complementary to the data points included so far.
Nonetheless, for the ranges of the Wilson coefficients in the following heat maps we use the marginalised O(Λ −2 ) bounds of Ref. [62] in order to cover a conservative parameter range.In Fig. 4 we show heat maps illustrating the dependence of the LO QCD cross section on the variation of C tG at the level of linear dimension-6 truncation (option (a)), compared to the leading couplings C tH and C H , which corresponds to a comparison on equal footing.The allowed ranges of Wilson coefficients are still quite large, such that a sizeable fraction of the 2-dimensional parameter space leads to unphysical negative cross section values.As to be expected, the effect of a variation of C tG within the given range is less pronounced than the one from variations of the leading couplings C tH and C H within their range.From a power counting point of view, the allowed range for C tG should be much smaller, such that the difference of the impact on the cross section would be even more obvious.Nevertheless, it is reasonable to derive bounds while being agnostic about the size of Wilson coefficients as well as considering power counting arguments on the expected impact.The latter is the approach we follow.
In Fig. 5, heat maps for the dependence of the cross section on a variation of (independent) 4-top operator pairs C Qt , such that the contribution of the scheme translation in Eq. (2.22) can be by accident of the same order or even larger than the original coefficient, inserting the numbers naively.the right plot it is apparent that the (LL)(LL) and (RR)(RR) operators of Ref. [5] with coefficients C Qt , (left plot of Fig. 5) have a large impact on the cross section in the considered range of values, leading to modifications of more than 100% of the LO cross section.The effect on the total cross section of C (8) Qt is stronger than the effect of C (1) Qt (in NDR), which is due to a large impact following from a sign change of the interference with the SM, visible in the upper left diagram of Fig. 9.
Fig. 6 shows the dependence of the LO cross section on the variation of C tG and C (1) Qt , comparing the NDR and BMHV scheme choices for the chiral structure of the 4-top operator.We introduce C (1/8) Qt; BMHV as a short-hand notation to specify that the corresponding amplitude is calculated in the BMHV scheme.Hence, this does not mean that the value of C  teresting showcase, since in Ref. [58] it has been demonstrated that the two Wilson coefficients are closely related, because part of the translation between the schemes is achieved by shifting C tG , see Eq. (2.22).Supplementing SMEFT with a tree-loop classification of Wilson coefficients, these shifts are of equal order in the power counting as the original value of C tG .In Fig. 6, the gradient of the cross section in NDR (left) points in a completely different direction than the one in BMHV (right) and also the magnitude of the gradient changes significantly.The effect of the translation of C tG in Eq. ( 2 ) within the ellipsis and vice versa.Note that this does not describe the full scheme translation, as the shift in C tH of Eq. (2.22) is not considered, however it is not as relevant as the shift in C tG , as will become clear in the discussion of Fig. 11.In addition, the shift of C tG depends on a scale dependent coupling g s which was set to a constant, thus these areas should be only understood as a qualitative visualisation.This clearly highlights that predictions using just C Qt would suffer from significant ambiguities if they are not considered in combination with C tG , since the scheme differences can only be resolved if shifts of the form in Eq. (2.22) are considered.Moreover, when the tree-loop classification of Wilson coefficients is applied, C tG would similarly suffer from ambiguities if C Qt was neglected.In principle, this also holds for other processes where operators that are connected by similar relations enter at the same order.This demonstrates that bounds set on these operators individually, without considering cancellations of the scheme dependence between different operator contributions, may not be very meaningful.
In Table 3 we present values for the total cross section for the SM and benchmark point 6, using truncation options (a) and (b) at NLO QCD.We also demonstrate their dependence on the variation of a single subleading Wilson coefficient.In general, the relative difference due to the variation of these Wilson coefficients is more pronounced for the SM cross section than for benchmark point 6.
Due to the asymmetric range of C tG , its variation tends to a damping of the cross section, with up to −36% relative to the SM.For benchmark point 6, truncation (a) leads to a larger relative effect of C tG on the cross section than truncation (b).
The variation of single 4-top Wilson coefficients, on the other hand, is fairly symmetric for the marginalised limits and has larger relative impact for truncation option (b) than for truncation option (a).The cross section difference for a variation of Qt is larger when working in the BMHV scheme than in NDR, and the scheme difference is much more visible for C (1) Qt variation leads to up to ∼ 35% effects on the cross section in the NDR scheme and up to ∼ 100% in BMHV, whereas for C (8) Qt the maximum difference is in both schemes ≳ 100%.As already indicated by the heat map on the right of Fig. 5, the effect of C QQ variation is very small, with a relative difference of less than 4% and being only a fraction of the uncertainty due to 3-point scale variations.The effects of C Qt and C Qt is also shown for the BMHV scheme, which is denoted by C

Higgs boson pair invariant mass distributions
In this section we present differential distributions depending on the invariant mass of the Higgs boson pair, m hh , combining NLO QCD results and subleading operator contributions at LO QCD.Each plot demonstrates the variation of a single subleading Wilson coefficient w.r.t.either the SM or benchmark point 6 for truncations (a) (linear dimension-6 only) and (b) (linear+quadratic dimension-6).The ranges we used are oriented at the O(Λ −2 ) marginalised fits of Ref. [62].
In Fig. 7 the variation of the chromomagnetic operator coefficient C tG in the ranges specified in Table 3 is  the scale uncertainty band.Note that the C tG -variation range is asymmetric around zero and that the interference of the C tG -term with the SM contribution tends to decrease the cross section.
In Fig. 8 we present the variation of the 4-top operator coefficient C QQ and the combination C (1) QQ + C tt .As observed at the level of total cross sections in Section 4.1, the contribution of these operators remains within the scale uncertainties, except for small deviations in the tails for the case of C (1) QQ + C tt .Thus the process gg → hh is not sensitive to those operators even if the coefficients are varied in ranges as large as [−190, 190].The situation is different for the operators C Qt , as we will show below.However, the contribution of these Wilson coefficients depends on the chosen γ 5 -scheme in dimensional regularisation, as explained in Section 2.2.We begin with Fig. 9 which demonstrates the effect of varying C (1) Qt .We observe sizeable effects, differing from the baseline prediction (SM or benchmark 6) by more than 100% for some regions, which also leads to negative cross section values.In NDR, the low-and high m hh -regions exhibit large differences beyond the scale uncertainty, with unphysical cross sections at low m hh values and a sign change around m hh ∼ 460 TeV.This behaviour changes significantly in BMHV: there are visible, but weaker effects in the low m hh -region, the sign change occurs around m hh ∼ 360 TeV and the deviation in the high m hh -region begins for lower invariant masses and is also more pronounced.
The scheme dependent behaviour of C Qt is shown in Fig. 10.For both schemes we Qt -variations on m hh -distributions comparing γ 5 -schemes.Left: NDR scheme, right: BMHV scheme; upper: SM baseline scenario, lower: benchmark point 6 for truncation options (a) and (b).observe small effects in the low m hh -region, a sign change of the contribution around m hh ∼ 360 TeV and a pronounced effect in the high m hh -region.Overall, the difference between the schemes is not as significant as in the case of C (1) Qt .The contribution to the m hh -distribution in the BMHV scheme (right column of Fig. 10) is qualitatively very similar to the case of C variations in NDR, we investigate the effect of those rational terms contributing in NDR which are responsible for the scheme difference and eventually the translation relation Eq. (2.22).We distinguish in the following between the scheme dependent parts Qt -variations on m hh -distributions comparing γ 5 -schemes.Left: NDR scheme, right: BMHV scheme; upper: SM baseline scenario, lower: benchmark point 6 for truncation options (a) and (b).leading to the shift of C tH .In Fig. 11 we present the difference to the SM m hhdistribution originating from those scheme dependent terms, where we individually vary C (1) Qt , respectively.Considering all scheme dependent terms, there is a prominent contribution from C Qt , which is much larger than the scale uncertainty of the SM result for the whole m hh -range, especially apparent in the low to intermediate m hh -regime.Investigating the constituents, we notice that ∆C tG is much more relevant than ∆C tH when considering the contributions from C Qt to the shift.Comparing the change on the distribution related to ∆C tG and ∆C tH separately (middle and bottom left panels in Fig. 11) to the effect of the sum of both contributions (top left panel in Fig. 11), we observe that the range of the band in the top left panel is given by the sum

Linear versus linear+quadratic contributions from the chromomagnetic operator
So far, the results involving the chromomagnetic operator only include its linear contribution, as it is classified subleading in the scenario of weakly coupling and renormalisable new physics such that its square would be beyond the order we consider.In this subsection, however, we step back from this assumption and assess the effect of the quadratic chromomagnetic contribution at LO QCD.As in the previous subsections we vary C tG within the ranges from O(Λ −2 ) marginalised fits of Ref. [62]. 9n Fig.  than the linear when comparing the distance between the two lines.Beyond C tG ∼ 0.2 the quadratic piece gets relevant and has an effect of 5% of the SM cross section.However, this only reduces the destructive interference of the linear contribution due to the asymmetric range, making the overall difference to the SM smaller.Fig. 12 (right) demonstrates the SM distribution at LO QCD together with a variation of C tG comparing the linear and linear+quadratic insertions.Similar to the observations in Fig. 12 (left) on the level of the total cross section, the quadratic terms are most relevant for the largest values of C tG , which however leads to a reduction of the destructive interference with the SM, thus reducing the overall effect on the distribution.Note that the tails do not yet reach the energy range where, as predicted by the high energy expansion of the helicity amplitudes in Ref. [104], the quadratic contribution takes over.

Conclusions
We have calculated the matrix elements including the chromomagnetic operator and 4-top operators contributing to Higgs boson pair production in gluon fusion and demonstrated that these operators both appear at the same subleading order in a power counting scheme that takes into account a tree-loop classification of dimension-6 SMEFT operators.We emphasize again that this classification is based on the generic assumption of a renormalisable and weakly coupling new physics sector, so does not represent all potential UV effects in full generality.These subleading contributions, entering the cross section at LO QCD, have been combined with the leading SMEFT operators including NLO QCD corrections as described in Ref. [13], in the form of Eqs.(2.6)- (2.8).This combination is provided as an extension to the public ggHH SMEFT code as part of the POWHEG-Box-V2.We have also described the usage of the new features.
The matrix elements of the 4-top contributions have been decomposed analogous to the case of gg → h described in Refs.[58,92].In particular, the parts depending on the γ 5 -scheme in dimensional regularisation have been identified, such that we found a similar scheme dependence as in the gg → h case, which can be understood as a finite shift of Wilson coefficients, see Eq. (2.22) and Ref. [58].
The effect of the subleading operators on the total cross section and on the Higgs boson pair invariant mass distribution has been studied in detail, both with respect to the SM and for benchmark point 6.We observed that the operators O only marginally contribute, therefore gg → hh is not an adequate process to probe those coefficients.The cross section is noticeably affected by a variation of the Wilson coefficient C tG within current conservative bounds, which can lead to a damping of the invariant mass distribution in the low to intermediate m hh -region.However, the highest sensitivity is observed by a variation of C Qt within current bounds.Since the limits on the 4-top Wilson coefficients from marginalised fits are very loosely constrained so far, the inclusion of processes like gg → h and gg → hh, where the operators enter at higher orders, could potentially improve the global determination of bounds on C Qt .As has been investigated for single Higgs production in Ref. [58] and confirmed in this work, contributions of those Wilson coefficients are precisely the ones which, when considered individually, depend on the chosen γ 5 -scheme.Therefore, bounds for individual coefficients can turn out to be significantly different due to a (more or less arbitrary) calculational scheme choice, which makes their interpretation difficult.In general, this scheme dependence enters as soon as the calculation is performed at an order at which loop contributions of such chiral current-current operators are to be considered.This statement is of particular relevance when the effects of these chiral current-current operators are investigated in processes where they only enter at loop level, as in this case contributions of Wilson coefficients entering at lower loop order are necessary to resolve γ 5 -scheme ambiguities.Considering a tree-loop classification of Wilson coefficients, the requirement for γ 5 -scheme independence has even stronger impications: In case of a clear hierarchy, e.g. by the loop suppression of the shift translating C tH in Eq. (2.22), the shift would only be a higher order effect.For loopinduced Wilson coefficients, this would however not be the case, as the shifts can be of the same order in the power counting.This holds, for example, for the Wilson coefficient C tG , which, at the same order in the power counting, can contain a contribution from C Qt , depending on the scheme choice.Inserting numerical values for current bounds on these Wilson coefficients [62] into Eq.(2.22) illustrates that the shift induced by a scheme change can even be larger than the interval given by the original bounds.To obtain more meaningful results, it is therefore recommended to study those Wilson coefficients which are connected through the scheme translation relations together, such that their combination is a scheme independent parametrisation of BSM physics at the studied order in the power counting.
In the future it would be desirable to have QCD corrections to those subleading operators as well, in order to compare on equal footing with the leading operators, at NLO QCD.However, including NLO corrections to the 4-top operators would require a 3-loop calculation involving Higgs and top-quark masses.The two-loop contributions to the chromomagnetic operator would be more feasible, but also challenging due to the high tensor rank induced by this operator.Therefore these calculations would be clearly beyond the scope of this paper.Furthermore, operators of the class ψ 2 ϕ 2 D have not been considered in this work, even though they would enter at the same power counting order, because they are considered as electroweak-type.However, this indicates that the strict separation between QCD and electroweak contributions becomes ambiguous once SMEFT operators beyond the leading contributions are included and combined with higher order corrections.
Finally, we note that renormalisation group running effects have not been included in the present study, even though they may lead to sizeable effects.This is left to upcoming work.

8 )
. The third line in Eq. (2.1) contains the chromomagnetic operator and lines 4-6 show the relevant 4-top operators.The operator O (3) 3333 qq, Warsaw of the Warsaw basis [5] has been replaced by O (QQ where the relation in terms of the Wilson coefficients has the form [68] top operators are already present in the 3rd generation 4-fermion operators of the Warsaw basis.

Figure 1 :
Figure 1: Feynman diagrams of the leading SMEFT contributions to gg → hh (Born level).Black dots denote insertions of SM couplings, gray dots (potentially) treeinduced EFT operators, gray squares denote insertions of loop-induced couplings (here C HG ).

Figure 2 :
Figure 2: Feynman diagrams involving insertions of the chromomagnetic operator.The gray squares denote insertions of the (loop-suppressed) chromomagnetic operator.

Figure 3 :
Figure 3: Feynman diagrams involving insertions of 4-top operators.The gray dots denote insertions of 4-top operators.

t and δ 4
-top C tH the diagram classes (a), (b) and (c) of Fig. 3 are made finite, and we write schematically (a/b/c/d) tG denote the amplitude of diagram types (a), (b), (c) and (d) of Fig. 2, respectively.The remaining diagrams of class (h) of Fig. 3 are made UV finite by the gghh counter term vertex using precisely the same value of δ 4-top C HG which is an indication that eq.(2.19) is indeed the correct two-loop counter term.Finally, we obtain QQ are shown.Looking at8 Interestingly, the conservative limits from the marginalised fits have values below one for C tG and values of O(100) for C(1)

σσ 1 Figure 4 :σσ 6 6 × 10
Figure 4: Heat maps showing the dependence of the LO cross section on the pair of Wilson coefficients C tG , C tH (left) and C tG , C H (right), respectively, with Λ = 1 TeV for the linear dimension-6 truncation.The ranges for C H are oriented at a translation of the limits on κ λ from Ref. [101], the ranges for the other Wilson coefficients are obtained at O (Λ −2 ) constraints from Ref.[62] (marginalised over the other coefficients).The white areas denote regions in parameter space where the corresponding cross section would be negative.

Figure 5 :QQ + C tt and C ( 8 )
Figure 5: Heat maps showing the dependence of the cross section on the couplings C (1) Qt

( 1 ) 8 )
QQ , C tt and C (QQ hardly affect the cross section.This can be understood by the very limited contribution to the amplitude, given only by the residual structure M 4-top ∆QQ,tt,(8) in Eq. (2.21).On the other hand, the (LL)(RR) operators, with coefficients C

(1/ 8 )σσ
Qtitself is changed by the scheme choice.This selection is an in-SMEFT / σ SM at LO with linear dim-6 SMEFT / σ SM at LO with linear dim-6 using BMHV

Figure 6 : 1 )
Figure 6: Heat maps demonstrating the effect of γ 5 -scheme choice on the dependence of the cross section on the couplings C tG and C (1) Qt with = 1 TeV.Left plot NDR, right plot BMHV.The ranges are taken from Ref. [62], based on an O(Λ −2 ) fit marginalised over the other Wilson coefficients.The areas within the black circle (left) and within the ellipsis (right) demonstrate value pairs of Wilson coefficients that would be mapped into each other by using the relation for C tG in Eq. (2.22).

1 ) 1 )
.22) is visualised by the areas surrounded by the black circle (left) and black ellip-sis (right), respectively: The relation for the scheme translation would map coefficient value pairs (C (Qt , C tG ) from within the circle onto value pairs (C (Qt , C BMHV tG

Table 3 :
Qt on the difference ∆C tG := C BMHV tG − C tG and ∆C tH := C BMHV tH − C tH are illustrated later at distribution level in Fig. 11.Total cross sections for Higgs-boson pair production at NLO QCD for the SM and benchmark point 6 using truncation option (a) or (b) at 13.6 TeV.The modification of the cross section due to a variation of the subleading Wilson coefficients is given as relative change to the base value in the second row.The uncertainties in the second row are scale uncertainties based on 3-point scale variations.The ranges of the subleading Wilson coefficients are oriented at O (Λ −2 ) constraints from Ref. [62] (Upper values: individual bounds, lower values: marginalised over the other coefficients).The effect of the Wilson coefficients C

Figure 9 :
Figure 9: Effects of C

( 1 )
Qt shown in Fig.9.In order to better understand the qualitative difference between the C

Figure 10 :
Figure 10: Effects of C Difference to SM at NLO for contribution from C Difference to SM at NLO for contribution from C Difference to SM at NLO for contribution from C Difference to SM at NLO for contribution from C

Figure 11 :
Figure 11: Demonstration of the difference ∆σ = dσ dm hh − dσ SM dm hh to the SM invariant mass distribution only including contributions of the scheme dependent terms, ∆C tG := C BMHV tG

Figure 12 :
Figure 12: Diagrams comparing linear only with linear+quadratic contribution of C tG (with Λ = 1 TeV) using variations within the marginalised O (Λ −2 ) constraints from Ref. [62].(Left) Total cross section normalised to the SM at LO QCD as a function of C tG , (right) envelope of C tG -variations on the SM m hh -distributions at LO QCD.

( 1 )
QQ , O tt and O (8) QQ the dimension-6 amplitude with the SM amplitude and the terms inside {. . .} are the |M dim6 | 2 parts of the cross section, which can be switched on or off in the ggHH SMEFT code.The EFT contribution only based on leading operators is 12 (left)we present the total cross section normalised to the SM value as a function of C tG .For moderate values of C tG the quadratic contribution is less dominant