Electroweak Higgs boson production in the standard model effective field theory beyond leading order in QCD

We study the impact of dimension-six operators of the standard model effective field theory relevant for vector-boson fusion and associated Higgs boson production at the LHC. We present predictions at the next-to-leading order accuracy in QCD that include matching to parton showers and that rely on fully automated simulations. We show the importance of the subsequent reduction of the theoretical uncertainties in improving the possible discrimination between effective field theory and standard model results, and we demonstrate that the range of the Wilson coefficient values allowed by a global fit to LEP and LHC Run~I data can be further constrained by LHC Run~II future results.


Introduction
The LHC Run I and early Run II data have not yet put forward any strong evidence of physics beyond the standard model (SM) and limits on new states have instead been pushed to higher and higher energies. As a consequence, the effective field theory (EFT) extension of the SM (SMEFT) has become increasingly relevant. The SMEFT is built from the SM symmetries and degrees of freedom (including the Higgs sector) by adding new operators of dimension higher than four to the SM Lagrangian. Being a tool to parameterise the search for new anomalous interactions, it is fully complementary to direct searches for new particles. Interpreting data in the context of the SMEFT hence allows us to be sensitive to new physics beyond the current energy reach of the LHC in a modelindependent way.
The formulation of the effective Lagrangian restricted to operators of dimension of at most six relies on the definition of a complete and non-redundant operator basis [1][2][3] and should additionally include the translations among the possible choices [4]. This has been intensively a celine.degrande@cern.ch b fuks@lpthe.jussieu.fr c kentarou.mawatari@lpsc.in2p3.fr d ken.mimasu@uclouvain.be e v.sanz@sussex.ac.uk discussed and will be soon reported by the Higgs cross section working group [5]. Moreover, since we try to observe small deviations from the SM, precise theoretical predictions are required both in the SM and in the SMEFT framework. The accumulation of LHC data and the subsequent precision obtained indeed call for a similar accuracy on the theoretical side, which demands the inclusion of higher-order corrections.
What we present in this paper is a part of the current theoretical activities aiming for precision predictions for electroweak Higgs-boson production at the LHC, i.e. Higgs boson production in association with a weak boson (VH) and via vector-boson fusion (VBF). In this context, NLO+PS (next-to-leading order plus parton-shower) matched predictions for VH and VBF production in the SM have been released both in the MC@NLO [6][7][8] and Powheg [9][10][11] frameworks, and merged NLO samples describing VH production including up to one additional jet have been generated in both the Powheg [12] and Sherpa [13] platforms. NLO QCD corrections along with the inclusion of anomalous interactions have been further investigated for VH [14] and VBF [15] Higgs boson production, and matched to parton showers in the Higgs characterisation framework [16,17]. Finally, electroweak corrections as well as anomalous coupling effects for VH production have been included in the Hawk program [18] that also contains NLO QCD contributions. In contrast, fixed-order predictions are known to a higher accuracy for both VH and VBF SM Higgs production processes [19][20][21][22].
In the SMEFT framework (in contrast to the anomalous coupling approach), the VH process has been studied at the NLO+PS accuracy within the Powheg-Box framework [23], where a subset of the 59 independent dimension-six operators was taken into account. In this paper, similarly, we consider five operators which are relevant for VH production. Firstly, we independently provide NLO+PS predictions for the VH process by using a different framework via a joint use of FeynRules [24], NloCT [25] and MadGraph5 aMC@NLO (MG5 a-MC) [26] programs. This approach provides a fully automatic procedure linking the model Lagrangian to event generation matched to parton showers at NLO. Our work hence not only independently validates the previous results obtained with Powheg-Box but also includes the additional benefits stemming from the flexibility of the FeynRules program. As a result, one can exploit generators like MG5 aMC for simulating any desired process at the NLO+PS accuracy (i.e. the VBF process in our case) for which the same operators play a role. This is the second part of our paper, which presents the first SMEFT results for this process. Although we only present results for a couple of benchmark scenarios motivated by global fit results, our predictions can be straightforwardly generalised to any scenario by using our public Universal FeynRules Output (UFO) model [27] within MG5 aMC [28][29][30].
We emphasise that, following some recent results in the ttH channel [31], this work represents a step towards a complete SMEFT operator basis implementation for Higgs physics at the NLO QCD accuracy, which will be beneficial to both the theoretical and experimental communities.
In Sect. 2 we provide the necessary theoretical ingredients to calculate NLO-QCD corrections for VH and VBF Higgs production in the SMEFT. We also discuss current constraints on the Wilson coefficients originating from a LEP and LHC Run I global fit analysis with which we inform our benchmark point selection. In Sect. 3 we describe our setup for NLO computations matched to parton showers. We present our numerical results in Sects. 4 and 5, and also assess the validity of the EFT given the current constraints on the Wilson coefficients. We asses the reach of the future LHC reach in Sect. 6, before concluding in Sect. 7. Practical information for event simulation and model validations are provided in the appendix. In the SM of particle physics, the elementary particles and their interactions are described by a quantum field theory based on the SU (3) C × SU (2) L × U (1) Y gauge symmetry. The vector fields mediating the gauge interactions lie in the adjoint representation of the relevant gauge group, where the notations for the representation refer to the full SM symmetry group. The chiral content of the theory is defined by three generations of left-handed and righthanded quark (Q L , u R and d R ) and lepton (L L and e R ) fields whose representation under the SM gauge group is given by The Higgs sector contains a single SU (2) L doublet of fields that is responsible for the breaking of electroweak symmetry, where the components of the Φ doublet are given in terms of the physical Higgs field h shifted by its vacuum expectation value v and the Goldstone bosons G ± and G 0 that are eaten by the weak bosons to give them their longitudinal degree of freedom.
In the EFT framework, new physics is expected to appear at a scale Λ large enough so that the new degrees of freedom can be integrated out. As a result, the SM Lagrangian L SM is supplemented by higher-dimensional operators O i parameterising all effects beyond the SM, Restricting ourselves to operators of dimension six, the most general gauge-invariant Lagrangian L has been known for a long time [32][33][34] and can be expressed in a suitable form by choosing a convenient basis of independent operators O i [1][2][3]. In this work, we focus on five specific, bosonic operators, 1 which are relevant to the VH and VBF processes, taken from the strongly interacting light Higgs (SILH) basis [2,36,37], 2 The Wilson coefficientsc are free parameters, T 2k are the generators of SU (2) (with Tr(T 2k T 2l ) = δ kl /2) in the fundamental representation and the Hermitian derivative operators are defined by In our conventions, the gauge-covariant derivatives and the gauge field strength tensors read where ij k are the structure constants of SU (2). In addition, g and g denote the coupling constants of SU (2) L and U (1) Y respectively.
After the breaking of the electroweak symmetry down to electromagnetism, the weak and hypercharge gauge eigenstates mix to the physical W -boson, Z-boson and the photon A, We have introduced in this expression the sine and cosine of the weak mixing angleŝ W ≡ sinθ W andĉ W ≡ cosθ W which diagonalise the neutral electroweak gauge boson mass matrix. The higher-dimensional operators of Eq. (5) induce a modification of the gauge boson kinetic terms that become, in the mass basis and after integration by parts, where W ± µν , Z µν and A µν denote the W -boson, Z-boson and photon field strength tensors, respectively. Consequently, canonical normalisation has to be restored by redefining the electroweak boson fields, We have made use here of the freedom related to the removal of the photon and Z-boson mixing terms induced by the higher-order operators. This mixing can indeed be absorbed either in a photon field redefinition, or in a Zboson field redefinition, or in both (as in Eq. (10)). In order to minimise the modification of the weak interactions with respect to the SM, we additionally redefine the weak and hypercharge coupling constants As a result of this choice, the relations between the measured values for the electroweak inputs and all internal electroweak parameters are simplified. The Z-boson mass m Z is now given by while the photon stays massless and the expression of the W -boson mass m W is unchanged respect to the SM one. We define the electroweak sector of the theory in terms of the Fermi coupling constant G F as extracted from the muon decay data, the measured Z-boson mass m Z and the electromagnetic coupling constant α in the low-energy limit of the Compton scattering. The vacuum expectation value of the Higgs field can therefore be derived from the Fermi constant as in the SM, v 2 = 1/( √ 2G F ) . After the field redefinitions of Eq. (10), the electromagnetic interactions of the fermions to the photon field turn out to be solely modified by the O W operator, so that the electromagnetic coupling constant e is related to the input parameter α as Furthermore, the shift in the cosine of the Weinberg mixing angle cosθ W can be derived, at first order in 1/Λ 2 , from the Z-boson mass relation of Eq. (12) along with Eq. (13), withc 2W ≡ cos 2θ W ,s W ≡ sinθ W and As a consequence, thec W andc B parameters are constrained by the measurement of the W -boson mass and by the Z-boson decay data. Those constraints can nonetheless be modified if other dimension-six operators are added to the Lagrangian of Eq. (5).
In unitary gauge and rotating all field to the mass basis, all three-point interactions involving a single (physical) Higgs boson and a pair of electroweak gauge bosons are given by where integration by parts has been used to reduce the number of independent Lorentz structures. Table 1 shows the relation between the couplings in Eq. (16) and the Wilson coefficients in Eq. (5). As a reference, we also compare our conventions to those of the previous SILH Lagrangian implementation of Ref. [37] and of the Higgs characterisation Lagrangian of Ref. [16].

Constraints from global fits of LEP and LHC Run I data
In this section we summarise the current bounds on the Wilson coefficients associated with the effective operators under consideration. We start from the results of previous works [38,39], where a global fit to LEP and LHC Run I data has been performed. The results imply constraints on several linear combinations of thec coefficients appearing in Eq. (5) that we present in Table 2, each limit having been obtained by marginalising over all other coefficients. Leading-order (LO) theoretical predictions have been used and in addition, the modifications of the electroweak parameters computed in Eq. (13) and Eq. (14) have not been considered for LHC predictions. We have nevertheless checked that the corresponding effects are small compared with the LHC Run I sensitivity, as also noted by the ATLAS collaboration [40].
In many classes of SM extensions (featuring in particular an extended Higgs sector), certain relations among the coefficients appear. For instance, it is common that matching conditions such that g (2) hww ∝c HW +c W = 0 appear [41]. In this case, the global fit generates the more stringent constraintc HW = −c W = [0.0008, 0.04] when one sets the effective scale to Λ = m W [38].

Benchmark points
For both production processes of interest, we consider two benchmark scenarios in the Wilson coefficient parameter space. These two points are selected to be compatible with the global fit results discussed in Sect. 2.2.
We first make use of the fact that, as seen in Table 2, electroweak precision observables strongly constrain a particular linear combination of thec W andc B Wilson coefficients beyond a precision than can be hoped for at the LHC. We therefore imposec B = −c W /2, which in turn leads to an allowed range (setting Λ = m W ) forc W of [−0.035, 0.005], as obtained from the second constraint on these two parameters. In order to highlight the impact of the two new Lorentz structures appearing in the interaction vertices of the Lagrangian of Eq. (16), we allow for non-zero values for both thec HW andc W coefficients.
Our benchmark scenarios are defined in Table 3. In the first setup, we only switch on the O HW operator (which induces new physics contributions to both the g (1) hvv and g (2) hvv structures). With the second point, we additionally fixc W to an equal and opposite value relying on the constraint relation brought up in Sect. 2.2. This allows for turning on solely the g (1) hvv coupling (see Table 1).

Setup for NLO+PS simulations
Our numerical results are derived at the NLO accuracy in QCD thanks to a joint use of the FeynRules/NloCT and MG5 aMC packages. The EFT Lagrangian of Eq. (5) has been implemented in FeynRules [24], while the computation of the ultraviolet counterterms and the rational R 2 terms necessary for numerical loop-integral evaluation has been done by NloCT [25] that relies on Fey-nArts [42]. The model information is then provided to MG5 aMC [26] in the UFO format [29]. Within MG5 aMC, loop-diagram contributions are numerically evaluated [43] and combined with the real emission pieces within the Eq. (16) Our conventions Ref. [37] (HEL) Ref. [16] (HC) gmW cακSMgHWW Table 1. New physics effects in three-point interactions involving a Higgs boson and electroweak gauge bosons. The loopinduced SM contributions to the Higgs-boson couplings to two photons aH and to one Z-boson and one photon a H have been explicitly indicated. Table 2. Current 95% confidence level constraints on the considered effective coefficients marginalised in a global fit to LEP and LHC Run I data [38]. Table 3. EFT benchmark points under consideration, with Λ = mW . We additionally setcBB =cHB = 0 for simplicity.

Coefficients Bounds
FKS subtraction scheme [44,45]. Short-distance events are finally matched to parton showers according to the MC@NLO prescription [6]. We generate events for 13 TeV LHC collisions using the LO and NLO NNPDF2.3 set of parton densities [46] for LO and NLO simulations, respectively. Events are then showered and hadronised within the Pythia8 infrastructure [47], which is also used for handling Higgs-boson decays. This latter step relies on eHdecay [48] that computes all branching fractions of the Higgs boson into the relevant final states to the first order in the Wilson coefficients. This procedure has the advantage of providing a correct normalisation for the production rates that includes all effects originating from the EFT operators. For the two adopted benchmark points, the deviations from the SM branching ratios are found to be very small.
Event reconstruction and analysis are performed using the MadAnalysis5 [49] framework, which makes use of all jet algorithms implemented in the FastJet program [50]. Jets are defined using the anti-k T algorithm [51] with a radius parameter of 0.4.
Theoretical uncertainties due to renormalisation (µ R ) and factorisation (µ F ) scale variations are accounted for thanks to the reweighting features of MG5 aMC [52]. At the event generation stage, nine alternative weights are stored for each event, corresponding to the independent variation of the two scales by a factor of two up and down with respect to a central scale µ 0 . Since the parton shower is unitary, these could be used to reweight the events after showering and reconstruction, saving a great deal of computational time and storage. The scale variation uncertainty is taken to be the largest difference between the central scale and the alternative scale choice predictions. We use as a central scale µ 0 = H T /2 and m W for VH and VBF processes, respectively, where H T is defined at the parton-level as the scalar sum of the transverse momentum of all visible final-state particles and the missing transverse energy. We refer the reader to the appendix for further technical details on event generation.

Higgs production in association with a vector boson
Higgs-boson production in association with a vector boson is an excellent probe for new physics, as the momentum transfer in the process is directly sensitive to the Lorentz structure appearing in the interaction vertices [53]. The use of differential information at the LHC Run I has therefore enhanced the sensitivity of Higgs data to possible new physics effects [39]. Those Run I studies have, however, relied on predictions evaluated at the LO accuracy in QCD. With the improved capabilities of the LHC Run II, NLO QCD effects become more relevant and more precise predictions are in order.
To showcase our NLO simulation setup for associated VH production, we study various differential distributions in the channel, where / E stands for the final-state missing energy. We impose the requirement that both b-jets and leptons have a pseudorapidity, η, and a transverse momentum, p T , satisfying |η| < 2.5 and p T > 25 GeV, respectively, while non-b-tagged jets are instead allowed to be more forward, with |η| < 4, for the same p T requirement. We select events by demanding the presence of one lepton and two b-jets based on truth-level hadronic information, a b-tagged jet being defined by the presence of a b-hadron within a cone of radius R = 0.4 centred on the jet direction.
In Fig. 1, we present the transverse momentum spectrum of the bb system (upper left), of the leading (upper centre) and next-to-leading (upper right) b-jets, of the lepton (middle left) and of the leading jet (middle centre). We then focus on the distribution in pseudorapidity for the bb system (middle right), in the transverse mass of the W -boson and Higgs boson (lower left) and of the Wboson, Higgs boson and leading-jet system (lower centre) and in the total transverse energy (lower right). In each subfigure, the results are shown both at the LO+PS and NLO+PS accuracies, together with uncertainties related to scale variation.
For each studied observable, we investigate in the first two lower bands of each subfigure the relative difference between the predictions in the SM and in the case of both considered benchmark points A and B, In the last band, we additionally show differential Kfactors defined as the binned ratio of NLO to LO predictions taking only the total NLO uncertainty into account. The predictions are found to be stable under radiative corrections, as expected for any process with a Drell-Yanlike topology. The obtained K-factors are indeed relatively flat and independent of the EFT parameters, with the exception of the observables that rely on the leading-jet kinematics which turn out to be much harder at NLO. Hard QCD radiation contributions originating from the matrix element are in this case included, in contrast to the LO setup where QCD radiation is only described by the parton shower and thus modelled in the soft-collinear kinematical limit.
LO predictions are found to be inaccurate and do not overlap with the NLO results even after considering scale variation uncertainties. This is particularly true at high transverse momentum p T , transverse mass M T and total transverse energy E T . This behaviour is once again expected for a Drell-Yan-like process that does not depend on α S at fixed LO. If one were to use the difference between the LO and NLO results as an error estimate for the LO predictions and the scale variation only for the NLO, then the reduction of the theory error would be better reflected by Fig. 1. The δ i SM ratios also remain stable with respect to QCD corrections, except at very high energies for the benchmark point A where small differences appear between the LO and NLO predictions. These would, by construction, be covered by the aforementioned improved definition of the LO theoretical uncertainties. All distributions strongly depend on the value of the EFT Wilson coefficients. For the adopted scenario A, large enhancements are observed in the tails of the p T , M T and E T distributions, which correspond to a centrally produced bb system (with a small pseudorapidity). In contrast, event rates are only rescaled by about 15-20% with respect to the SM for the scenario B. This originates from the g (2) hvv coupling that vanishes in this scenario, so that only the g (1) hvv coupling drives the EFT behaviour in the high-energy tails. However, this latter coupling is known to yield a smaller impact than the g (2) hvv coupling [17,23], and it is therefore the presence of the g (2) hvv interaction vertex in scenario A that leads to the large observed deviations. This constitutes a very promising avenue for setting limits on EFT parameters from W H studies and similar behaviour can be observed for ZH production, where thē c HB andc BB coefficients additionally play a role. In this case, the gluon fusion initiated contribution should however also be considered, as discussed in Refs. [23,54].
While such large enhancements can be exploited to obtain powerful constraints on the SMEFT Wilson coefficients, they do raise the question of the validity of the EFT approach at large momentum transfer [39,[55][56][57]. This question could be addressed with the use of dedicated benchmark models to compare the breakdown of the EFT framework against well-motivated ultravioletcomplete models [41,58]. At a more simplistic level we can also make use of the MG5 aMC ability to select only interference contributions (at LO) to assess the impact of the squared EFT terms given our benchmark choices (technical details are described in Appendix A). Fig. 2 shows a selection of distributions, overlaying predictions with and without this squared term. We observe significant differences between the two choices that are greater than the scale uncertainty of the predictions. Depending on the observable, these can range from 40 to 100% on the interference-only prediction for the benchmark scenario A, while they are much milder for the benchmark scenario B. This suggests that current sensitivities on this region of the Wilson coefficient parameter space may not yet lend themselves to an EFT interpretation within the validity of the framework. A reduction of the production rate from the SM value, as seen for benchmark scenario B, moreover indicates the dominance of the interference term between the SM and EFT contributions given that the squared terms are positive-definite (Fig. 4).

Higgs production via vector boson fusion
Another powerful probe of anomalous higher-derivative interactions between weak and Higgs bosons consists of the VBF Higgs production mode where it is produced in as-   sociation with two forward jets, Our event selection requires the presence of at least two jets with a pseudorapidity |η| < 4.5 and a transverse momentum p T > 25 GeV, and we additionally impose the requirement that the Higgs boson decays into a pair of photons with a pseudorapidity satisfying |η| < 2.5 and a transverse momentum p T > 20 GeV. We moreover include a standard VBF selection on the invariant mass M jj and pseudorapidity separation ∆η jj of the pair of forward jets, M jj > 500 GeV and ∆η jj > 3 .
Several kinematical observables are sensitive to the momentum flow in the VBF process, for which EFT contributions deviate from the SM prediction. We consider in Fig. 4 the distribution in the transverse momentum of the diphoton system (upper left), in the p T of the leading (upper centre) and subleading (upper right) jets, in the invariant mass of the dijet system (lower left), as well as in its pseudorapidity (lower centre) and azimuthal angular (lower right) separations. The consistent definition of scale uncertainties that are possible with the NLO predictions helps to quantify the discriminatory power between the new physics benchmarks and the SM. Similarly to the VH process, the NLO corrections are independent of the EFT parameters and cannot be completely described by an overall K-factor.
In contrast to the VH process, we observe a depletion of the production rate for both benchmark scenarios, which mainly impacts the high-energy tails of the differential distributions. This indicates that our predictions may be safer with respect to the validity of the EFT, as it implies that the interference term dominates over the EFT squared one. In particular, the effects for the benchmark scenario B are more pronounced with respect to the SM compared to the V H production case and show some different shape deformations. This illustrates the complementarity between the VH and VBF processes in disentangling the possible EFT sources for any potential deviation. Although the correlations between the forward jets as well as between the jets and the Higgs boson are known to be sensitive to new physics effects [59,60], those are less sensitive than the individual Higgs and jet p T distributions for our two benchmark scenarios.
We also repeat the simple EFT validity analysis performed for the V H case and assess the impact of the EFT squared terms at LO, as shown in Fig. 4. As suggested by the depletion effect of the EFT operators in the high energy bins of the differential distributions, the squared terms appear much more under control in this process compared to the V H case. Within the ranges of our predictions the impact of the squared term is again most pronounced for the benchmark A, reaching at most 5-12% while for benchmark B their effect is much smaller.

Future LHC reach
Before concluding, we attempt to estimate the reach of the LHC Run II with respect to the Wilson coefficients considered in our benchmark scenarios. Our results so far suggest to make use of the high energy tails of differential distributions as handles for new physics. For concreteness, we focus on the associated production process p p → H W ± → bb ± + / E T which has already been searched for at both LHC Run I [61] and II [62]. In both analyses, a large number of signal and control regions are defined according to the lepton and additional jet multiplicities, as well as to the vector boson transverse momentum p V T . These are combined in a global fit to obtain the corresponding SM Higgs signal strength. In this fitting procedure, several dominant components of the background, namely tt and W -boson production in association with heavy-flavour jets, are left free to float. We consider as signal regions the p V T overflow bin in the single lepton channel for both the zero-jet and one-jet categories. The Run II study, however makes use of multivariate methods in the event selection process that are difficult to reproduce -a task that definitely lies beyond the scope of the simple estimate we are intending to derive. We therefore choose to consider only the cut-based signal selection procedure employed in the Run I analysis and then project the results for various Run II integrated luminosities.

Signal prediction and background estimation
In order to estimate the number of background events in the single lepton signal regions of a possible cut-based, LHC Run II analysis, we extrapolate the results of the corresponding Run I analysis. We first consider the dominating tt contribution which arises from semi-leptonic topantitop decays and which makes up 54 and 85% of the total background in the 0-jet and 1-jet categories, respectively. As a crude estimate for the corresponding 13 TeV yields, we compute a transfer factor f i tr (with i = 0, 1 for the 0-jet and 1-jet categories) by generating large statistics of SM semi-leptonic tt events at centre-of-mass energies of 8 and 13 TeV on which we apply the kinematic selection of the Run I analysis summarised in Appendix B.
The transfer factor f i tr is defined as the ratio of the two fiducial cross sections and we deduce the Run II analysis background contributions by multiplying the 8 TeV SM expectation σ i bkg inferred from the Run I background event counts N i bkg assuming 25 fb −1 of 8 TeV data. These should not depend much on the actual composition of 7 and 8 TeV data analysed, particularly in the high transverse momentum overflow bin which is dominated by 8 TeV data. Table 4 summarises the information obtained and used in the subsequent analysis. Our theoretical predictions for the tt contributions to the 0-and 1-jet signal regions at 8 TeV, σ overf. 8 , lie within a factor 2 of the cross-sections  Table 4. Information necessary to estimate the 13 TeV projections of the fitted background yields, N bkg , in the p V T > 200 GeV overflow bins quoted in the analysis of Ref. [61]. We show the tt fiducial cross sections after the kinematic selection of Appendix B for 8 and 13 TeV collisions and assuming a 29.2% semi-leptonic branching fraction, along with the derived transfer factor f i tr .
inferred from the post-selection, fitted background decomposition presented in Table 5 of Ref. [61]. Due to the multivariate nature of the recent Run II analysis, its fitted background yields cannot be used to validate the results of our projection, which rather represents the scenario in which a cut-based analysis similar to the Run I counterpart were performed at 13 TeV.
The signal predictions have been generated using the previously described UFO implementation, and both the tt and W H contributions have been simulated at the NLO accuracy in QCD as described in previous sections, the fixed-order results being matched with Pythia 8 for both handling the top decays and the parton showering. A grid of points spanning the allowed region of the (c HW ,c W ) parameter space, including the SM prediction, was simulated and the generated events were passed through the same kinematic selection of Appendix B. Following the previous discussion on the existing constraints, we assume the simplificationc W = −c B /2. Such a relation would not be retained in a complete global fit including, e.g., LEP data. However, it is instructive to follow this simplified path as it assesses the sensitivity of the LHC to the direction in the (c W ,c B ) plane that is orthogonal, and thus complementary, to the one tightly constrained by precision measurements at the Z-pole. We have derived least-squaresfitted quadratic polynomial forms for the zero and one-jet overflow bin cross sections in the two-dimensional parameter plane σ i W H (c HW ,c W ). Our results, moreover, embed a b-tagging efficiency of 70%, and more information (in particular on the explicit coefficients of the fits) is given in Appendix B.

Results
Our results have been derived from the fitted functional forms for the signal cross sections in combination with the projected background yields. We have performed a ∆χ 2 analysis to extract 95% confidence intervals assuming L =30, 300 and 3000 fb −1 of integrated luminosities of 13 TeV proton-proton collisions. Denoting by B i and S i the event counts in the signal region in the backgroundonly and signal-plus-background hypotheses, respectively, we have The 95% confidence intervals are obtained at the boundary of ∆χ 2 = 5.99 which equates to the corresponding p-value for a χ 2 distribution with two degrees of freedom. Fig. 5 depicts these confidence intervals for the zero-and one-jet bin separately as well as their combination, for the three integrated luminosity points. For comparison, the marginalised single parameter exclusion regions established in Table 2 and the benchmark discussion of Sect. 2.3 are indicated. This simplified projection shows that this type of analysis is likely to substantially improve the existing limits on these Wilson coefficients in combination with existing data. Since the one-jet category suffers both from a larger background and smaller signal contribution, its relative impact on the overall reach is small. The blind direction associated with this measurement lies very close to thē c HW = −c W line, corresponding to the benchmark choice B of the earlier sections. This is consistent with the very mild expected impact of this particular new physics scenario in the high energy tails of the differential distributions for the pp → W + H process (see Sect. 4). Nevertheless, our results suggest that in the general case, taking the full integrated luminosity of LHC Run II will individually . 95% confidence intervals in the (cHW ,cW ) plane depicting the projected reach at the LHC Run II extracted from data in the p W T overflow bin of the corresponding Run I analysis performed in Ref. [61]. We consider three different choices for the integrated luminosities, and the dashed lines indicate the previously obtained marginalised limits on the Wilson coefficients from the global fit of Refs. [38,39].
allow to constrain Wilson coefficients with a precision of a few per-mille and the results presented in Sec. 5 indicate that combining VBF and WH studies will break this degeneracy.

Conclusions
We have presented FeynRules and UFO implementations of dimension-six SMEFT operators affecting electroweak Higgs-boson production, which can be used for NLO(QCD)+PS accurate Monte Carlo event generation within the MG5 aMC framework. We have considered five SILH basis operators and have accounted for all field redefinitions that are necessary to canonically normalise the theory. Moreover, the ensuing modifications of both the gauge couplings and the relationships between the electroweak inputs and the derived parameters have also been included. We have showcased the strength of our approach by simulating both associated VH and VBF Higgsboson production at the 13 TeV LHC, selecting a pair of benchmark scenarios informed both by recent limits from global fits to the LEP and LHC Run I data and by theoretical motivations originating from integrating out certain popular ultraviolet realisations.
We have found that EFT predictions and deviations from the SM are stable under higher-order corrections. Overall, we have also observed a significant reduction of the theoretical errors, which would have an impact on the future measurements aiming to unravel dimension-six operator contributions.
Furthermore, as a test for the validity of the EFT approach, we have proposed to compare distributions that either include the full matrix element (embedding all SM and new physics contributions) or account solely for the interference of the SM component with the new physics component. For our benchmark choices that saturate current experimental limits, the differences were observed to be large in the kinematic extremes of some of our distri-butions, particularly for V H production. This points to the possibility that the EFT description is breaking down in these regions of the parameter space and that the most precise measurements undertaken at the LHC Run II may be required to probe the EFT (while staying in its region of validity).
When comparing results for the VH and VBF channels, we have found that both Higgs-boson production modes are sensitive to new physics, but the VH one seems to have a better handle on g (1) -type (V µν V µν h) structure, since several key distributions display deviations that may be more easily distinguished from the background. Although the QCD K-factors have been observed not to depend on the EFT parameters, the reduced theoretical uncertainties are crucial for disentangling a non-vanishing contribution of the g (1) -type structure to the predictions from the SM. This has been singled out in our study of the benchmark point B. Moreover, our results exhibit an interesting complementarity of the two Higgs production channels, since the interference pattern between the SM and the SMEFT contributions is quite different and benchmark-dependent.
In order to estimate the reach that might be possible at LHC Run II, we have performed a simplified analysis projecting the Run I SM background expectations in a search for W H associated production and combining this information with LHC Run II signal predictions obtained using our implementation. Using the overflow bin of the reconstructed W -boson transverse momentum distribution in the single lepton channel as a probe for EFT effects suggests that the LHC Run II will significantly improve the current limits obtained from global fits. Clearly both the V H and VBF processes deserve further investigation including detector effects and an analysis strategy to reject the SM backgrounds. In this case, the new physics contributions to the SM background processes should also correctly be accounted for, since effective operators affecting electroweak Higgs-boson production also impact the normalisation of the triple gauge-boson interactions both directly and indirectly via the aforementioned field redefinitions [23,38,41,63,64].
Finally, our work has demonstrated a proof-of-concept for automated NLO+PS simulations in the SMEFT framework. To this aim, we have limited ourselves to a small set of dimension-six operators and a pair of benchmark points. This is characterised as a first step towards a complete operator basis implementation, in which the renormalisation group running of the Wilson coefficients [31,[65][66][67][68][69] could also be supplemented in the future. Since the usual decay syntax of MG5 aMC is not available for NLO event generation, we directly request the presence of the W -boson decay products in the final state. An alternative way would require one to simulate the production of a Higgs boson in association with an on-shell W -boson that is subsequently decayed within the Mad-Spin infrastructure [70,71] before invoking the parton showering. On the other hand, VBF Higgs-boson production is achieved by typing in the generation command generate p p > h j j $$ w+ w-z a QCD=0 [QCD]

A Simulation in MadGraph5 aMC@NLO
The removal of the s-channel gauge boson contributions (via the $$ syntax) avoids the generation of VH topologies where the gauge boson decays hadronically.
In the parameter card, users may set values for the five Wilson coefficients defined in Eq. (5) as well as for the cutoff scale Λ. The modifications to the electroweak parameters in terms of the inputs in the (G F , m Z , α) scheme are taken into account, as well as the shifts induced by the field redefinitions discussed in Sect. 2.
At the FeynRules level, the final Lagrangian involving all of the redefined fields and parameters is expanded up to O(1/Λ 2 ). From the point at which this truncation occurs, all subsequent performed calculations will necessarily induce some O(1/Λ n ) (with n > 2) dependence from, e.g., higher powers in the electroweak couplings. Furthermore, MG5 aMC constructs its matrix elements by squaring helicity amplitudes, so that the squared EFT contribution terms that are formally of O(1/Λ 4 ) are by default retained on top of the leading interference terms with the SM of O(1/Λ 2 ). A positive definite cross section is ensured at the price of including higher-order terms. It is therefore important to remain in the regime where the EFT expansion can be trusted in that higher-order contributions are sufficiently suppressed. Departure from this safe zone may be reflected by a rapid growth of amplitudes with energy leading to extreme deviations in the tails of the distributions, as well as by discrepancies between independent Monte Carlo setups performing their truncation in different ways. Alternatively, this issue could be investigated thanks to some recent developments in MG5 aMC that allow the user to specify an order for the squared matrix element calculation, like QED^2 or QCD^2. However, this feature is currently only available for LO computations. For including the EFT effects in the simulation, the coupling order parameter NP can hence be set either to NP=1 to retain the full amplitude squared or to NP^2<=1 to throw away the aforementioned EFT squared terms, e.g.

generate p p > h ve e+ NP^2<=1
Finally, due to the details of the implementation, we advise users seeking to recover the SM limit to avoid setting the Wilson coefficients to zero. It is preferable to either set them to very small non-zero values or to use restrictions. The cutoff scale parameter NPl may also be alternatively fixed to a very large value.

A.2 Comparison with the MCFM/POWHEG-BOX implementation
We have verified that our results in the W H channel are compatible with those stemming from an alternative existing implementation, which is based on MCFM [72][73][74] and the Powheg-Box framework [75,76] and that has been introduced in Ref. [23]. We have scrutinised several differential distributions for both our benchmark points in the Powheg and MG5 aMC frameworks, as presented in Fig. 6. A good consistency has been found up to statistical uncertainties and despite the difference in the dynamical scale choice which is taken as H T /2 in MG5 aMC and the invariant mass of the Higgs-vector boson system in the Powheg-Box implementation. We summarise here the kinematic selection through which the tt and signal W H events have been passed in order to determine the background transfer factor and signal efficiencies of the analysis performed in Ref. [61]. The signal region used is the p W T > 200 GeV overflow bin in the 0 and 1-jet categories of the single-lepton channel. The kinematic selection for this channel applied to our event samples after parton shower is as follows: • We require the final state to contain exactly one lepton with |η| < 2.47 and E T > 25 GeV. • Jets reconstructed by means of the anti-k T jet algorithm with a radius parameter R = 0.4, and we discard the jet candidates for which the conditions |η| < 4.5 and p T > 20 GeV are not satisfied. • We require exactly two b-jets with a pseudorapidity |η| < 2.5, the hardest one being further imposed tho have a transverse momentum p T > 45 GeV.
Not more than one additional jet in the |η| < 2.5 region is allowed and any event with a jet with p T > 30 GeV in the |η| > 2.5 region is also rejected. Events are then split into the zero-and one-jet categories based on the presence of an extra, non-forward jet softer than the two b-jets. The vector boson transverse momentum p W T is defined as the vector sum of the transverse momentum of the lepton and the missing transverse energy. Additionally, in the p W T > 200 overflow bin, a cut of 50 GeV on the missing energy is imposed as well as a cut on the distance between the two b-jets, ∆R bb ≡ ∆η 2 bb + ∆φ 2 bb < 1.4. A flat btagging efficiency of 70% is in addition assumed. 34 signal samples have been simulated with the Wilson coefficients being taken in the ranges −0.02 <c HW < 0.03 and −0.03 <c W < 0.01, the total rates being rescaled by the Higgs branching fraction to bb obtained with eHDE-CAY. These events have then been passed through the above selection in order to determine the functional forms for the zero-and one-jet signal cross sections in p W T over-flow bin, σ 0 and σ 1 . A least-squares fit yields The fit coefficients are within 5-10% of one another, which is to be expected given that our signal process receives contributions from an electroweak vertex and should not be directly sensitive to additional QCD radiation.