Gluon-fusion Higgs production in the Standard Model Effective Field Theory

We provide the complete set of predictions needed to achieve NLO accuracy in the Standard Model Effective Field Theory at dimension six for Higgs production in gluon fusion. In particular, we compute for the first time the contribution of the chromomagnetic operator $ \bar Q_L \Phi \sigma q_R G$ at NLO in QCD, which entails two-loop virtual and one-loop real contributions, as well as renormalisation and mixing with the Yukawa operator $\Phi^\dagger \Phi\, \bar Q_L \Phi q_R$ and the gluon-fusion operator $\Phi^\dagger \Phi\, GG$. Focusing on the top-quark-Higgs couplings, we consider the phenomenological impact of the NLO corrections in constraining the three relevant operators by implementing the results into the MadGraph5_aMC@NLO framework. This allows us to compute total cross sections as well as to perform event generation at NLO that can be directly employed in experimental analyses.


Introduction
Five years into its discovery at the LHC, the Higgs boson is still the centre of attention of the high-energy physics community.A wealth of information has been collected on its properties by the ATLAS and CMS experiments [1][2][3][4][5], all of which so far support the predictions of the Standard Model (SM).In particular, the size of the couplings to the weak vector bosons and to the electrically charged third generation fermions has been confirmed, and the first evidence of the coupling to second generation fermions (either charm quark or muon) could arrive in the coming years, if SM-like.
The steady improvement in the precision of the current and forthcoming Higgs measurements invites to explore physics beyond the SM not only via the search of new resonances, as widely pursued at the LHC, but also via indirect effects on the couplings of the Higgs boson to the known SM particles.The most appealing aspect of such an approach is that, despite being much more challenging than direct searches both experimentally and theoretically, it has the potential to probe new physics scales that are beyond the kinematical reach of the LHC.A powerful and predictive framework to analyse possible deviations in the absence of resonant BSM production is provided by the SM Effective Field Theory (SMEFT) [6][7][8], i.e., the SM augmented by higher-dimensional operators.Among the most interesting features of this framework is the possibility to compute radiative corrections in the gauge couplings, thus allowing for systematic improvements of the predictions and a reduction of the theoretical uncertainties [9].In particular, higher-order corrections in the strong coupling constant typically entail large effects at the LHC both in the accuracy and the precision.They are therefore being calculated for a continuously growing set of processes involving operators of dimension six featuring the Higgs boson, the bottom and top quarks and the vector bosons.Currently, predictions for the most important associated production channels for the Higgs boson are available in this framework, e.g., VH, VBF and t tH [10][11][12].For top-quark production, NLO results for EW and QCD inclusive production, i.e., tj and t t, and for top-quark associated production t tZ, t tγ have also appeared [13][14][15][16][17][18].The effect of dimension-six operators has also become available recently for top-quark and Higgs decays [19][20][21][22][23].
The situation is somewhat less satisfactory for gluon fusion, which, despite being a loop-induced process in the SM, is highly enhanced by the gluon density in the proton and provides the most important Higgs-production channel at the LHC.In the SM, the QCD corrections are now known up to N 3 LO in the limit of a heavy top quark [24][25][26].The full quark-mass dependence is known up to NLO [27][28][29][30], while at NNLO only subleading terms in the heavy top-mass expansion [31][32][33][34] and leading contributions to the top/bottom interference [35,36] are known.Beyond inclusive production, the only available NNLO result is the production of a Higgs boson in association with a jet in the infinite top-mass limit [37][38][39], while cross sections for H + n-jets, n = 2, 3, are known only at NLO in the heavy top-mass expansion [40,41].
In the SMEFT, most studies have been performed at LO, typically using approximate rescaling factors obtained from SM calculations.Higher-order results have only been considered when existing SM calculations could be readily used within the SMEFT.The simplest examples are the inclusion of higher orders in the strong coupling to the contribution of two specific dimension-six operators, namely the Yukawa operator (Φ † Φ) QL Φq R and the gluon-fusion operator (Φ † Φ)GG.The former can be accounted for by a straightforward modification of the Yukawa coupling of the corresponding heavy quark, b or t, while the latter involves the computation of contributions identical to SM calculations in the limit of an infinitely-heavy top quark.Results for the inclusive production cross section including modified top and bottom Yukawa couplings and an additional direct Higgs-gluons interaction are available at NNLO [42] and at N 3 LO [43,44].At the differential level, phenomenological studies at LO have shown the relevance of the high transverse momentum region of the Higgs boson in order to resolve degeneracies among operators present at the inclusive level [12,[45][46][47].Recently, the calculation of the Higgs spectrum at NLO+NNLL level for the Yukawa (both b and t) and Higgs-gluons operator has appeared [48,49].
The purpose of this work is to provide the contribution of the chromomagnetic operator QL Φσq R G to inclusive Higgs production at NLO in QCD, thereby completing the set of predictions (involving only CP -even interactions) needed to achieve NLO accuracy in the SMEFT for this process.The first correct computation at one-loop of the contribution of chromomagnetic operator of the top quark to gg → H has appeared in the erratum of ref. [50] and later confirmed in refs.[12,49].The LO contribution of the chromomagnetic operator of the top-quark to H+jet was computed in ref. [12].An important conclusion drawn in ref. [12] was that even when the most stringent (and still approximate) constraints from t t production are considered [14], this operator sizably affects Higgs production, both in gluon fusion (single and double Higgs) and t tH production.
At LO the chromomagnetic operator enters Higgs production in gluon fusion at one loop.Therefore NLO corrections in QCD entail two-loop virtual and one-loop real contributions.The latter can nowadays easily be computed using an automated approach.The former, however, involve a non-trivial two-loop computation that requires analytic multiloop techniques and a careful treatment of the renormalisation and mixing in the SMEFT, both of which are presented in this work for the first time.In particular, while the full mixing pattern of the SMEFT at one loop is known [51][52][53], a new two-loop counterterm enters our computation, and we provide its value for the first time here.Moreover, we present very compact analytic results for all the relevant amplitudes up to two loop order.Focusing on possibly anomalous contributions in top-quark-Higgs interactions, we then consider the phenomenological impact of the NLO corrections, including also the Yukawa operator and the gluon-fusion operator at NLO by implementing the respective virtual two-loop matrix elements into the MadGraph5_aMC@NLO framework [54].This allows us to compute total cross sections as well as to perform event generation at NLO plus parton shower (NLO+PS) that can be directly employed in experimental analyses.
The paper is organised as follows.In section 2 we establish our notations and set up the calculation by identifying the terms in the perturbative expansion that are unknown and need to be calculated.In section 3 we describe in detail the computation of the two-loop virtual contributions and the renormalisation procedure and we provide compact analytic expressions for the finite parts of the two-loop amplitudes.We also briefly discuss the leading logarithmic renormalisation group running of the Wilson coefficients.In section 4 we perform a phenomenological study at NLO, in particular of the behaviour of the QCD and EFT expansion at the total inclusive level and provide predictions for the p T spectrum of the Higgs via a NLO+PS approach.

Gluon fusion in the SM Effective Field Theory
The goal of this paper is to study the production of a Higgs boson in hadron collisions in the SMEFT, i.e., the SM supplemented by a complete set of operators of dimension six, The sum in eq.(2.1) runs over a basis of operators O i of dimension six, Λ is the scale of new physics and C b i are the (bare) Wilson coefficients, multiplying the effective operators.A complete and independent set of operators of dimension six is known [7,55].In this paper, we are only interested in those operators that modify the contribution of the heavy quarks, bottom and top quarks, to Higgs production in gluon fusion.Focusing on the top quark, Representative diagrams contributing to gluon-fusion amplitudes with one insertion of the three relevant operators.Heavy quarks, b or t, provide the leading contributions to the first and third amplitudes.Note that for chromomagnetic operator, QL Φσq R G, a diagram featuring the four point gluon-quark-quark-Higgs interaction is also present (not shown).
there are three operators of dimension six that contribute to the gluon-fusion process, where g s is the (bare) strong coupling constant and v denotes the vacuum expectation value (vev) of the Higgs field Φ ( Φ = iσ 2 Φ).Q L is the left-handed quark SU (2)-doublet containing the top quark, t R is the right-handed SU (2)-singlet top quark, and G a µν is the gluon field strength tensor.Finally, T a is the generator of the fundamental representation of SU (3) (with [T a , T b ] =1 2 δ ab ) and σ µν = i 2 [γ µ , γ ν ], with γ µ the Dirac gamma matrices.Two comments are in order.First, the corresponding operators O 1 and O 3 for the b quark can be obtained by simply making the substitutions { Φ → Φ, t R → b R }.Second, while O 2 is hermitian O 1 and O 3 are not. 1 In this work, we focus on the CP -even contributions of O 1 and O 3 .For this reason, all the Wilson coefficients C i with i = 1, 2, 3 are real.Representative Feynman diagrams contributing at LO are shown in fig. 1.
In the SM and at leading order (LO) in the strong coupling the gluon-fusion process is mediated only by quark loops.This contribution is proportional to the mass of the corresponding quark and therefore heavy quarks dominate.While we comment on the b (and possibly c) contributions later, let us focus on the leading contributions coming from the top quark, i.e., the contributions from the operators of dimension six shown in eqs.(2.2 -2.4).The (unrenormalised) amplitude can be cast in the form where α b s = g2 s /(4π) denotes the bare QCD coupling constant and m H and m b t are the bare masses of the Higgs boson and the top quark.The factor S = e −γ E (4π) is the usual MS factor, with γ E = −Γ (1) the Euler-Mascheroni constant and µ is the scale introduced by dimensional regularisation.For i = 0, the form factor A b,i denotes the unrenormalised SM contribution to gluon fusion [56], while for i > 0 it denotes the form factor with a single 2 operator O i inserted [48,50,57].The normalisation of the amplitudes is chosen such that all coupling constants, as well as all powers of the vev v, are explicitly factored out.Each form factor admits a perturbative expansion in the strong coupling, Some comments about these amplitudes are in order.First, after electroweak symmetry breaking, the operator O 1 only amounts to a rescaling of the Yukawa coupling, i.e., A b,1 is simply proportional to the bare SM amplitude.Second, at LO the operator O 2 contributes at tree level, while the SM amplitude and the contributions from O 1 and O 3 are loop-induced.Finally, this process has the unusual feature that the amplitude involving the chromomagnetic operator O 3 is ultraviolet (UV) divergent, and thus requires renormalisation, already at LO [12,49,50].The UV divergence is absorbed into the effective coupling that multiplies the operator O 2 , which only enters at tree level at LO.The renormalisation at NLO will be discussed in detail in section 3. The goal of this paper is to compute the NLO corrections to the gluon-fusion process with an insertion of one of the dimension six operators in eqs.(2.2 -2.4).We emphasise that a complete NLO computation requires one to consider the set of all three operators in eq.(2.2 -2.4), because they mix under renormalisation [51][52][53].At NLO, we need to consider both virtual corrections to the LO process g g → H as well as real corrections due to the emission of an additional parton in the final state.Starting from NLO, also partonic channels with a quark in the initial state contribute.Since the contribution from O 1 is proportional to the SM amplitude, the corresponding NLO corrections can be obtained from the NLO corrections to gluon-fusion in the SM including the full top-mass dependence [27,28,30,58].The NLO contributions from O 2 are also known, because they are proportional to the NLO corrections to gluon-fusion in the SM in the limit where the top quark is infinitely heavy [59] (without the higher-order corrections to the matching coefficient).In particular, the virtual corrections to the insertion of O 2 are related to the QCD form factor, which is known through three loops in the strong coupling [60][61][62][63][64][65][66][67][68][69].Hence, the only missing ingredient is the NLO contributions to the process where the chromomagnetic operator O 3 is inserted.The computation of this ingredient, which is one of the main results of this paper, will be presented in detail in the next section.
As a final comment, we note that starting at two loops other operators of EW and QCD nature will affect gg → H.In the case of EW interactions, by just looking at the SM EW contributions [70,71], it is easy to see that many operators featuring the Higgs field will enter, which in a few cases could also lead to constraints, see, e.g., the trilinear Higgs self coupling [72,73].In the case of QCD interactions, operators not featuring the Higgs field will enter, which, in general, can be more efficiently bounded from other observables.For example, the operator g s f abc G ν aµ G λ bν G µ cλ contributes at two loops in gg → H and at one loop in gg → Hg.The latter process has been considered in ref. [74], where effects on the transverse momentum of the Higgs were studied.For the sake of completeness, we have reproduced these results in our framework, and by considering the recent constraints on this operator from multi-jet observables [75], we have confirmed that the Higgs p T cannot be significantly affected.For this reason we do not discuss further this operator in this paper.Four-fermion operators also contribute starting at two loops to gluon fusion but as these modify observables related to top quark physics at leading order [76,77] we expect them to be independently constrained and work under the assumption that they cannot significantly affect gluon fusion.

Computation of the two-loop amplitudes
In this section we describe the virtual corrections to the LO amplitudes in eq.(2.5).For the sake of the presentation we focus here on the calculation involving a top quark and discuss later on how to obtain the corresponding results for the bottom quark.With the exception of the contributions from O 2 , all processes are loop-induced, and so the virtual corrections require the computation of two-loop form factor integrals with a closed heavyquark loop and two external gluons.We have implemented the operators in eqs.(2.2 -2.4) into QGraf [78], and we use the latter to generate all the relevant Feynman diagrams.The QGraf output is translated into FORM [79,80] and Mathematica using a custom-made code.The tensor structure of the amplitude is fixed by gauge-invariance to all loop orders, cf.eq.(2.5), and we can simply project each Feynman diagram onto the transverse polarisation tensor.The resulting scalar amplitudes are then classified into distinct integral topologies, which are reduced to master integrals using FIRE and LiteRed [81][82][83][84][85].After reduction, we can express all LO and NLO amplitudes as a linear combination of one and two-loop master integrals.
The complete set of one-and two-loop master integrals is available in the literature [58,[86][87][88] in terms of harmonic polylogarithms (HPLs) [89], In the case where all the a i 's are zero, we define, H(0, . . ., 0 3) The number of integrations w is called the weight of the HPL.The only non-trivial functional dependence of the master integrals is through the ratio of the Higgs and the top masses, and it is useful to introduce the following variable, or equivalently The change of variables in eq. ( 3.4) has the advantage that the master integrals can be written as a linear combination of HPLs in x.In the kinematic range that we are interested in, 0 < m 2 H < 4m 2 t , the variable x is a unimodular complex number, |x| = 1, and so it can be conveniently parametrised in this kinematics range by an angle θ, In terms of this angle, the master integrals can be expressed in terms of (generalisations of) Clausen functions (cf.ref. [58,[90][91][92][93] and references therein), where we used the notation , σ 1 , . . ., 0, . . ., 0 The number k of non-zero indices is called the depth of the HPL.
Inserting the analytic expressions for the master integrals into the amplitudes, we can express each amplitude as a Laurent expansion in whose coefficients are linear combinations of the special functions we have just described.The amplitudes have poles in which are of both ultraviolet (UV) and infrared (IR) nature, whose structure is discussed in the next section.

UV & IR pole structure
In this section we discuss the UV renormalisation and the IR pole structure of the LO and NLO amplitudes.We start by discussing the UV singularities.We work in the MS scheme, and we write the bare amplitudes as a function of the renormalised amplitudes as, where Z g is the field renormalisation constant of the gluon field and α s (µ 2 ), C i (µ 2 ) and m t (µ 2 ) are the renormalised strong coupling constant, Wilson coefficients and top mass in the MS scheme, and µ denotes the renormalisation scale.The renormalised parameters are related to their bare analogues through with (a 1 , a 2 , a 3 ) = (3, 0, 1).Unless stated otherwise, all renormalised quantities are assumed to be evaluated at the arbitrary scale µ 2 throughout this section.We can decompose the renormalised amplitude into the contributions from the SM and the effective operators, similar to the decomposition of the bare amplitude in eq.(2.5) and each renormalised amplitude admits a perturbative expansion in the renormalised strong coupling constant, The presence of the effective operators alters the renormalisation of the SM parameters.Throughout this section we closely follow the approach of ref. [12], where the renormalisation of the operators at one loop was described.The one-loop UV counterterms for the strong coupling constant and the gluon field are given by where δZ g,SM and δZ αs,SM denote the one-loop UV counterterms in the SM, and β 0 is the one-loop QCD β function, where N c = 3 is the number of colours and N f = 5 is the number of massless flavours.We work in a decoupling scheme and we include a factor µ 2 /m 2 t into the counterterm.As a result only massless flavours contribute to the running of the strong coupling, while the top quark effectively decouples [59].The renormalisation of the strong coupling and the gluon field are modified by the presence of the dimension six operators, but the effects cancel each other out [50].Similarly, the renormalisation of the top mass is modified by the presence of the effective operators, where the SM contribution is In eq.(3.16) we again include the factor µ 2 /m 2 t into the counterterm in order to decouple the effects from operators of dimension six from the running of the top mass in the MS scheme.
The renormalisation of the effective couplings C b i is more involved, because the operators in eqs.(2.2 -2.4) mix under renormalisation.The matrix Z C of counterterms can be written in the form We have already mentioned that the amplitude A b,3 requires renormalisation at LO in the strong coupling, and the UV divergence is proportional to the LO amplitude C is non-trivial at LO in the strong coupling, At NLO, we also need the contribution δZ C to eq. (3.18).We have δZ where, apart from z 23 , all the entries are known [51][52][53].z 23 corresponds to the counterterm that absorbs the two-loop UV divergence of the operator O 3 , which is proportional to the tree-level amplitude A (0) b,2 in our case.This counterterm is not available in the literature, yet we can extract it from our computation.NLO amplitudes have both UV and IR poles, and so we need to disentangle the two types of divergences if we want to isolate the counterterm z 23 .We therefore first review the structure of the IR divergences of NLO amplitudes, and we will return to the determination of the counterterm z 23 at the end of this section.
A one-loop amplitude with massless gauge bosons has IR divergences, arising from regions in the loop integration where the loop momentum is soft or collinear to an external massless leg.The structure of the IR divergences is universal in the sense that it factorises from the underlying hard scattering process.More precisely, if A (1) denotes a renormalised one-loop amplitude describing the production of a colourless state from the scattering of two massless gauge bosons, then we can write [94] A (1) = I (1) ( ) where A (0) is the tree-level amplitude for the process and R is a process-dependent remainder that is finite in the limit → 0. The quantity I (1) ( ) is universal (in the sense that it does not depend on the details of the hard scattering) and is given by where s 12 = 2p 1 p 2 denotes the center-of-mass energy squared of the incoming gluons.Since in our case most amplitudes are at one loop already at LO, we have to deal with two-loop amplitudes at NLO.However, since the structure of the IR singularities is independent of the details of the underlying hard scattering, eq.(3.21) remains valid for two-loop amplitudes describing loop-induced processes, and we can write We have checked that our results for amplitudes which do not involve the operator O 3 have the correct IR pole structure at NLO.For A 3 , instead, we can use eq.(3.23) as a constraint on the singularities of the amplitude.This allows us to extract the two-loop UV counterterm z 23 .We find Note that the coefficient of the double pole is in fact fixed by requiring the anomalous dimension of the effective couplings to be finite.We have checked that eq.(3.24) satisfies this criterion, which is a strong consistency check on our computation.Let us conclude our discussion of the renormalisation with a comment on the relationship between the renormalised amplitudes in the SM and the insertion of the operator O 1 .We know that the corresponding unrenormalised amplitudes are related by a simple rescaling, and the constant of proportionality is proportional to the ratio C b 1 /m b t .There is a priori no reason why such a simple relationship should be preserved by the renormalisation procedure.In (the variant of) the MS-scheme that we use, the renormalised amplitudes are still related by this simple scaling.This can be traced back to the fact that the MS counterterms are related by 1 will in general not hold after renormalisation.

Analytic results for the two-loop amplitudes
In this section we present the analytic results for the renormalised amplitudes that enter the computation of the gluon-fusion cross section at NLO with the operators in eqs.(2.2 -2.4) included.We show explicitly the one-loop amplitudes up to O( 2 ) in dimensional regularisation, as well as the finite two-loop remainders R i defined in eq.(3.21).The amplitudes have been renormalised using the scheme described in the previous section and all scales are fixed to the mass of the Higgs boson, µ 2 = m 2 H .The operator O 2 only contributes at one loop at NLO, and agrees (up to normalisation) with the one-loop corrections to Higgs production via gluon-fusion [59].The amplitude is independent of the top mass through one loop, and so it evaluates to a pure number, where β 0 is defined in eq.(3.15).The remaining amplitudes have a non-trivial functional dependence on the top mass through the variables τ and θ defined in eq.(3.4) and (3.6).
We have argued in the previous section that in the MS-scheme the renormalised amplitudes A (1) 0 and A 1 are related by a simple rescaling, A We therefore only present results for the SM contribution and the contribution from O 3 .We have checked that our result for the two-loop amplitude in the SM agrees with the results of ref. [27,28,30,58].The two-loop amplitude A 3 is genuinely new and is presented here for the time.
The one-loop amplitude in the SM can be cast in the form where the coefficients a i are given by The finite remainder of the two-loop SM amplitude is where we have defined the function (3.31) The one-loop amplitude involving the operator O 3 is where the coefficients b i are given by b The finite remainder of the two-loop amplitude A (1) 3 is Although the main focus of this paper is to include effects from dimension six operators that affect the gluon-fusion cross section through the top quark, let us conclude this section by making a comment about effects from the bottom, and to a lesser extent, the charm quark.The amplitudes presented in this section are only valid if the Higgs boson is lighter than the quark-pair threshold, τ < 4. It is, however, not difficult to analytically continue our results to the region above threshold where τ > 4. Above threshold, the variable x defined in eq.(3.5) is no longer a phase, but instead we have −1 < x < 0. As a consequence, the Clausen functions may develop an imaginary part.In the following we describe how one can extract the correct imaginary part of the amplitudes in the region above threshold (see also ref. [27,28,30,58]).
We start from eq. (3.7) and express all Clausen functions in terms of HPLs in x and its inverse, e.g., HPLs evaluated at 1/x can always be expressed in terms of HPLs in x.For example, one finds (3.36) Similar relations can be derived for all other HPLs in an algorithmic way [89,95,96].The previous equation, however, is not yet valid above threshold, because the logarithms H(0; x) = log x may develop an imaginary part.Indeed, when crossing the threshold x approaches the negative real axis from above, x → x + i0, and so the correct analytic continuation of the logarithms is The previous rule is sufficient to perform the analytic continuation of all HPLs appearing in our results.Indeed, it is known that an HPL of the form H(a 1 , . . ., a k ; x) has a branch point at x = 0 only if a k = 0, and, using the shuffle algebra properties of HPLs [89], any HPL of the form H(a 1 , . . ., a k , 0; x) can be expressed as a linear combination of products of HPLs such that if their last entry is zero, then all of its entries are zero.The amplitudes can therefore be expressed in terms of two categories of HPLs: those whose last entry is non-zero and so do not have a branch point at x = 0, and those of the form H(0, . . ., 0 which are continued according to eq. (3.37).
Using the procedure outlined above, it is possible to easily perform the analytic continuation of our amplitudes above threshold.The resulting amplitudes contribute to the gluon-fusion process when light quarks, e.g., massive bottom and/or charm quarks, are taken into account.Hence, although we focus primarily on the effects from the top quark in this paper, our results can be easily extended to include effects from bottom and charm quarks as well.

Renormalisation group running of the effective couplings
After renormalisation, our amplitudes depend explicitly on the scale µ, which in the following we identify with the factorisation scale µ F .It can, however, be desirable to choose different scales for the strong coupling constant, the top mass and the effective couplings.In this section we derive and solve the renormalisation group equations (RGEs) for these parameters.
Since we are working in a decoupling scheme for the top mass, the RGEs for the strong coupling constant and the top mass are identical to the SM with N f = 5 massless flavours.We have checked that we correctly reproduce the evolution of α s and m t in the MS scheme, and we do not discuss them here any further.For the RGEs satisfied by the effective couplings, we find where C = (C 1 , C 2 , C 3 ) T , and the anomalous dimension matrix is given by As already mentioned in the previous section, the double pole from the two-loop counterterm in eq.(3.24) cancels.We can solve the RGEs in eq.(3.38) to one loop, and we find We show in fig. 2 the quantitative impact of running and mixing by varying the renormalisation scale from 10 TeV to m H /2 in two scenarios: one where all Wilson coefficients are equal at 10 TeV and another where only C 3 is non-zero.This latter example serves as a reminder of the need to always consider the effect of all the relevant operators in phenomenological analyses as choosing a single operator to be non-zero is a scale-dependent choice.

Cross-section results
In this section we perform a phenomenological study of Higgs production in the SMEFT, focusing on anomalous contributions coming from the top quark.Results are obtained within the MadGraph5_aMC@NLO framework [54].The computation builds on the implementation of the dimension-six operators presented in ref. [12].Starting from the SMEFT Lagrangian, all tree-level and one-loop amplitudes can be obtained automatically using a series of packages [97][98][99][100][101][102].The two-loop amplitudes for the virtual corrections are implemented in the code through a reweighting method [103,104].Within the Mad-Graph5_aMC@NLO framework NLO results can be matched to parton shower programs, such as PYTHIA8 [105] and HERWIG++ [106], through the MC@NLO [107] formalism.
Results are obtained for the LHC at 13 TeV with MMHT2014 LO/NLO PDFs [108], for LO and NLO results respectively.The values of the input parameters are The values for the central scales for µ R , µ F and µ EF T are chosen as m H /2, and we work with the top mass in the on-shell scheme.We parametrise the contribution to the cross section from dimension-six operators as Within our setup we can obtain results for σ SM , σ i , and σ ij .We note here that results for single Higgs and H + j production in the SMEFT were presented at LO in QCD in ref. [12].
The normalisation of the operators used here differs from the one in ref. [12], but we have found full agreement between the LO results presented here and those of ref. [12] when this difference is taken into account.Furthermore, the SM top-quark results obtained here have been cross-checked with the NLO+PS implementation of aMCSusHi [109].
Our results for the total cross section at the LHC at 13 TeV at LO and NLO are shown in table 1.We include effects from bottom-quark loops (top-bottom interference and pure bottom contributions) into the SM prediction by using aMCSusHi.However, in this first study, we neglect bottom-quark effects from dimension-six operators in σ i and σ ij as we assume them to be subleading.As mentioned above, our analytic results and MC implementation can be extended to also include these effects.We see that the contributions from effective operators have K-factors that are slightly smaller then their SM counterpart, with a residual scale dependence that is almost identical to the SM.In the following we present an argument which explains this observation.We can describe the total cross section for Higgs boson production to a good accuracy by taking the limit of an infinitely heavy top quark, because most of the production happens near threshold.In this effective theory where the top quark is integrated out, all contributions from SMEFT operators can be described by the same contact interaction κG a µν G µν a H.The Wilson coefficient κ can be written as where κ 0 denotes the SM contribution and κ i those corresponding to each operator O i in the SMEFT.As a result each σ i is generated by the same Feynman diagrams both at LO and NLO in the infinite top-mass EFT.The effect of radiative corrections is, however, not entirely universal as NLO corrections to the infinite top-mass EFT amplitudes come both from diagrammatic corrections and corrections to the Wilson coefficients κ i , which can be obtained by matching the SMEFT amplitude to the infinite top mass amplitude, as illustrated in fig. 3. Indeed, each κ i can be expressed in terms of SMEFT parameters as a perturbative series i + O(α 2 s ).In the infinite top mass EFT, each K-factor K i can be decomposed as  where K U is the universal part of the K-factor, which is exactly equal to K 2 .By subtracting K 2 to each K i in the infinite top mass limit numerically (setting m t = 10TeV), we could extract the ratios α s κ (1) i κ i and check explicitly that these non-universal corrections are subdominant compared to the universal diagrammatic corrections, which explains the similarity of the effects of radiative corrections for each contribution.
Our results can be used to put bounds on the Wilson coefficients from measurements of the gluon-fusion signal strength µ ggF at the LHC.Whilst here we do not attempt to perform a rigorous fit of the Wilson coefficients, useful information can be extracted by a simple fit.For illustration purposes, we use the recent measurement of the gluon-fusion signal strength in the diphoton channel by the CMS experiment [110] µ ggF = 1.1 ± 0.19, ( which we compare to our predictions for this signal strength under the assumption that the experimental selection efficiency is not changed by BSM effects where we set Λ = 1TeV and kept only the O(1/Λ 2 ) terms.We therefore find that we can put the following constraint on the Wilson coefficients with 95% confidence level: While the correct method for putting bounds on the parameter space of the SMEFT is to consider the combined contribution of all relevant operators to a given observable, the presence of unconstrained linear combinations makes it interesting to consider how each operator would be bounded if the others were absent in order to obtain an estimate of the size of each individual Wilson coefficient.Of course such estimates must not be taken as actual bounds on the Wilson coefficients and should only be considered of illustrative value.We obtain For these individual operator constraints, the impact of the σ ii terms on the limits is at most 10%.
For reference we note that if one includes the O(1/Λ 4 ) contributions the linear combination in the bound becomes a quadratic one:

Differential distributions
In the light of differential Higgs measurements at the LHC, it is important to examine the impact of the dimension-six operators on the Higgs p T spectrum.It is known that measurements of the Higgs p T spectrum can be used to lift the degeneracy between O 1 and O 2 [12,45,111].For a realistic description of the p T spectrum, we match our NLO predictions to the parton shower with the MC@NLO method [107], and we use PYTHIA8 [105] for the parton shower.Note that we have kept the shower scale at its default value in MC@NLO, which gives results that are in good agreement with the optimised scale choice of ref. [112], as discussed in ref. [109].
The normalised distributions for the transverse momentum and rapidity of the Higgs boson are shown in figs.4 for the interference contributions.The impact of the O(1/Λ 4 ) terms is demonstrated in fig. 5 for the transverse momentum distribution.We find that the operators O 3 and O 2 give rise to harder transverse momentum tails, while for O 1 the    Finally we show the transverse momentum distributions for several benchmark points which respect the total cross-section bounds in fig.6.The operator coefficients are chosen such that eq.(4.10) is satisfied.We find that larger deviations can be seen in the tails of the distributions for coefficient values which respect the total cross-section bounds.

Renormalisation group effects
The impact of running and mixing between the operators is demonstrated in fig.7, where we show the individual (O(1/Λ 2 )) contributions from the three operators in gluon-fusion Higgs production at LO and NLO, as a function of µ EF T , assuming that C 3 = 1, C 1 = C 2 = 0 at µ EF T = m H /2 and Λ = 1 TeV.While at µ = m H /2 the only contribution is coming from the chromomagnetic operator, this contribution changes rapidly with the scale.While the effect of the running of C 3 is only at the percent level, σ 3 has a strong dependence on the scale.At the same time non-zero values of C 1 and C 2 are induced through renormalisation group running, which gives rise to large contributions from O 2 .We find that the dependence on the EFT scale is tamed when the sum of the three contributions is considered.This is the physical cross section coming from C 3 (m H /2) = 1 which has a weaker dependence on the EFT scale.The dependence of this quantity on the scale gives an estimation of the higher order corrections to the effective operators and should be reported as an additional uncertainty of the predictions.By comparing the total contributions at LO and NLO we find that the relative uncertainty is reduced at NLO.

Conclusion and outlook
A precise determination of the properties of the Higgs boson and, in particular, of its couplings to the other SM particles is one of the main goals of the LHC programme of the coming years.The interpretation of such measurements, and of possible deviations in the context of an EFT, allows one to put constraints on the type and strength of hypothetical new interactions, and therefore on the scale of new physics, in a model-independent way.The success of this endeavour will critically depend on having theoretical predictions that at least match the precision of the experimental measurements, both in the SM and in the SMEFT.
In this work we have computed for the first time the contribution of the (CP -even part of the) QL Φσq R G operator to the inclusive Higgs production at NLO in QCD.Since the NLO corrections for the other two (CP -even) operators entering the same process are available in the literature, this calculation completes the SMEFT predictions for this process at the NLO accuracy.Even though our results can be easily extended to include anomalous couplings of the bottom quark, we have considered in the detail the case where new physics mostly affects the top-quark couplings.Our results confirm the expectations based on previous calculations and on the general features of gluon-fusion Higgs production: at the inclusive level the K-factor is of the same order as that of the SM and of the other two operators.The residual uncertainties estimated by renormalisation and factorisation scale dependence also match extremely well.The result of the NLO calculation confirms that the chromomagnetic operator cannot be neglected for at least two reasons.The first is of purely theoretical nature: the individual effects of QL Φσt R G and Φ † ΦGG are very much dependent on the EFT scale, while their sum is stable and only mildly affected by the scale choice.The second draws from the present status of the constraints.Considering the uncertainties in inclusive Higgs production cross section measurements and the constraints from t t production, the impact of the chromomagnetic operator cannot be neglected in global fits of the Higgs couplings.As a result, a two-fold degeneracy is left unresolved by a three-operator fit using the total Higgs cross section and one is forced to look for other observables or processes to constrain all three of the operators.
The implementation of the finite part of the two-loop virtual corrections into Mad-Graph5_aMC@NLO has also allowed us to study the process at a fully differential level, including the effects of the parton shower resummation and in particular to compare the transverse momentum distributions of the SM and the three operators in the region of the parameter space where the total cross section bound is respected.Once again, we have found that the contributions from QL Φσt R G and Φ † ΦGG are similar and produce a shape with a harder tail substantially different from that of the SM and the Yukawa operator (which are the same).While QL Φσt R G and Φ † ΦGG cannot really be distinguished in gluon-fusion Higgs production, they do contribute in a very different way to t tH where the effect of Φ † ΦGG is extremely weak.Therefore, we expect that H, H+jet, and t tH (and possibly t t) can effectively constrain the set of the three operators.
In this work we have mostly focused our attention on the top-quark-Higgs boson interactions and only considered CP -even operators.As mentioned above and explained in section 3, extending it to include anomalous couplings for lighter quarks, the bottom and possibly the charm, is straightforward.On the other hand, extending it to include CP -odd operators requires a new independent calculation.We reckon both developments worth pursuing.

. 25 )
If the top mass and the Wilson coefficient C b 1 are renormalised using a different scheme which breaks this relation between the counterterms, the simple relation between the amplitudes A

Figure 3 .
Figure 3. Diagrammatic description of the matching between the SMEFT and the infinite top mass EFT at LO (left) and at NLO (right).The NLO amplitude in the infinite top-mass EFT contains two elements: diagrammatic corrections, which contribute universally to the K-factors and Wilson coefficient corrections, which are non-universal.

Figure 4 .
Figure 4.Higgs distributions, normalised for the interference contributions from σ i .Left: Higgs transverse momentum.Right: Higgs rapidity.SM contributions and individual operator contributions are displayed.Lower panels give the ratio over the SM.

Figure 5 .
Figure 5. Higgs transverse momentum distributions, normalised.Left: squared contributions σ ii .Right: interference between operators, σ ij .SM contributions and operator contributions are displayed.Lower panels give the ratio over the SM.

Figure 6 .
Figure 6.Transverse momentum distributions of the Higgs for different values of the Wilson coefficients.The lower panel shows the ratio over the SM prediction for the various benchmarks and the SM scale variation band.

σ(total) c 3 σ 3 c 2 σ 2 c 1 σ 1 µFigure 7 .
Figure 7. Contributions of the three operators to the inclusive Higgs production cross section at the LHC at 13 TeV as a function of the EFT scale.Starting from one non-zero coefficient at µ EF T = m H /2 we compute the EFT contributions at different scales, taking into account the running and mixing of the operators.LO and NLO predictions are shown in dashed and solid lines respectively.

Table 1 .
Total cross section in pb for pp → H at 13 TeV, as parametrised in eq.(4.3).