SMEFT at NNLO+PS: Vh production

In the context of the Standard Model effective field theory (SMEFT) the next-to-next-to-leading (NNLO) QCD corrections to the Higgsstrahlungs ($Vh$) processes in hadronic collisions are calculated and matched to a parton shower (PS). NNLO+PS precision is achieved for the complete sets of SMEFT operators that describe the interactions between the Higgs and two vector bosons and the couplings of the Higgs, a $W$ or a $Z$ boson, and light fermions. A POWHEG-BOX implementation of the computed NNLO SMEFT corrections is provided that allows for a realistic exclusive description of $Vh$ production at the level of hadronic events. This feature makes it an essential tool for future Higgs characterisation studies by the ATLAS and CMS collaborations. Utilising our new Monte Carlo code the numerical impact of NNLO$+$PS corrections on the kinematic distributions in $pp \to Zh \to \ell^+ \ell^- h$ production is explored, employing well-motivated SMEFT benchmark scenarios.


Introduction
The dominant decay of the Standard Model (SM) Higgs boson is into a pair of bottom quarks, with an expected branching ratio close to 60% for a Higgs mass of 125 GeV.Large QCD backgrounds from multi-jet production, however, make a search in the dominant gluon-gluon fusion (ggF) Higgs production channel very challenging at the LHC.The most sensitive production mode for detecting h → b b decays is the associated production of a Higgs boson and a massive gauge boson (V h), where the leptonic decay of the vector boson enables a clean selection, leading to a significant background reduction.The h → b b decay mode has been observed by both ATLAS and CMS in LHC Run II [1,2], and these measurements constrain the h → b b signal strength in the Higgsstrahlungs processes µ b b V h to be SM-like within about 25%.With LHC Run III ongoing and the high-luminosity upgrade (HL-LHC) on the horizon, the precision of the µ b b V h measurements is expected to improve significantly, resulting in an ultimate projected HL-LHC accuracy of 15% (5%) in the case of the W h (Zh) production channel [3,4].
Besides providing a probe of the dominant decay mode of the Higgs boson, precision V h measurements also play an important role in the Higgs characterisation programme which is commonly performed in the framework of the SM effective field theory (SMEFT) [5][6][7].In fact, radiative corrections in the SMEFT to both V h production [8][9][10][11][12][13] and the h → f f decays [14][15][16][17][18] have been calculated.The existing studies for V h production have mostly focused on the subset of higher-dimensional interactions that modify the couplings of the Higgs to two vector bosons achieving next-to-leading order (NLO) [8-10, 19, 20] and nextto-next-to-leading order (NNLO) [21] in QCD, respectively, while in the case of h → b b both NLO QCD and NLO electroweak (EW) corrections to the total decay width have been calculated for the full set of relevant dimension-six SMEFT operators [14][15][16].In the publications [11][12][13] special attention has finally been paid to the class of SMEFT operators that lead to interactions between a Higgs, W or Z boson, and light quarks.
The goal of the present work is to generalise and to extend the recent SMEFT calculation [18] which has achieved NNLO plus parton shower (NNLO+PS) accuracy for the dimension-six operators that contribute to the subprocesses pp → Zh and h → b b directly in QCD.This class of operators includes effective Yukawa-and chromomagnetic dipole-type interactions of the bottom quark that modify the h → b b decay but do not play a role in pp → Zh production.Purely EW effective interactions that alter the couplings of the Higgs to gauge bosons are instead not included in the NNLO+PS Monte Carlo (MC) generator presented in [18].Since these types of SMEFT contributions can lead to phenomenologically relevant effects in the Higgsstrahlungs processes [8][9][10][19][20][21], we include these type of interactions in the current article, extending the NLO SMEFT calculations [8][9][10] to the NNLO level.Likewise, we improve the precision of the calculations of SMEFT corrections to pp → V h production that are associated to couplings between a Higgs, W or Z boson, and light quarks [11][12][13] to NNLO in QCD.The obtained fixed-order SMEFT predictions are implemented into the POWHEG-BOX [22] and consistently matched to a parton shower (PS) using the MiNNLO PS method [23,24].In this way, NNLO QCD accuracy is retained for both production and decays, while the matching to the PS ensures a realistic exclusive description of the pp → Zh → ℓ + ℓ − h and the pp → W h → ℓνh process at the level of hadronic events.These features make our new NNLO+PS generator a precision tool for future LHC Higgs characterisation studies in the SMEFT framework.
This paper is structured as follows: in Section 2 we specify the dimension-six SMEFT operators that are relevant in the context of our work.A comprehensive description of the basic ingredients of the SMEFT calculation of pp → V h production and their implementation in the context of our NNLO+PS event generator is presented in Section 3. We motivate simple SMEFT benchmark scenarios in Section 4 by discussing the leading constraints on the Wilson coefficients of the relevant operators.The impact of the SMEFT corrections on kinematic distributions in pp → Zh → ℓ + ℓ − h production at the LHC is studied in Section 5 employing the simple benchmark scenarios for the Wilson coefficients identified earlier.In Section 6 we outline how the NNLO+PS calculations of pp → V h and h → b b can be combined, while Section 7 contains our conclusions and an outlook.The analytic formulae for the parameters and couplings that have been implemented into our MC code are relegated to Appendix A. Details on our implementation of the gg → Zh process can be found in Appendix B. In Appendix C we finally discuss the impact of NNLO QCD effects by comparing results for pp → Zh → ℓ + ℓ − h production obtained at NLO+PS and NNLO+PS, respectively.

SMEFT operators
Throughout this work we neglect all light fermion masses in both the SM and SMEFT corrections to the pp → Zh and pp → W h processes.The full set of dimension-six SMEFT operators has been presented in the so-called Warsaw basis in the article [6].This basis contains the following three independent operators µν W a,µν , Q HWB = H † σ a H W a µν B µν , (2.1) that modify the couplings between the Higgs and two vector bosons at tree level.The SM Higgs doublet is denoted by H, while B µν and W a µν are the U (1) Y and SU (2) L gauge field strength tensors and σ a are the Pauli matrices.In the case of the operators that result in couplings between the Higgs, a W or a Z boson, and light quarks, we consider the following four effective interactions where H with D µ the usual covariant derivative.The symbol q denotes left-handed quark doublets, while u and d are the right-handed quark singlets of up and down type, respectively.Illustrative diagrams that contribute to Zh production and involve an insertion of one of the operators in (2.1) or (2.2) are displayed in Figure 1.
Besides the two sets of operators (2.1) and (2.2) that alter the pp → V h production process, we also consider effective interactions that modify the Z → ℓ + ℓ − and W → ℓν decays at tree level.In the Warsaw basis there are three such operators, namely Here ℓ and e denote a left-handed lepton doublet and right-handed lepton singlet field, respectively.Notice that in writing (2.2) and (2.3) we have assumed that the full SMEFT Lagrangian respects an approximate U (3) 5 flavour symmetry which allows us to drop all flavour indices.The final type of SMEFT corrections that change the Higgsstrahlungs processes indirectly is provided by the Wilson coefficients of the operators that shift the Higgs kinetic term and/or the EW SM input parameters.In order to fully describe these shifts the following three additional operators are needed at tree level: Figure 1: Tree-level SMEFT contributions to q q → Zh production.The diagram on the left involves an insertion of one of the operators defined in (2.1), while the two graphs on the right stem from an insertion of one of the operators given in (2.2).The operator insertions are indicated by the blue squares.

Calculation in a nutshell
In this section, we sketch the different ingredients of our NNLO+PS SMEFT calculation of pp → V h production.We begin by recalling the basic steps of the NNLO QCD computation within the SM and then detail the general method that we employ to calculate the relevant squared matrix elements in the SMEFT and their implementation into the POWHEG-BOX.It is then explained how the fixed-order NNLO SMEFT calculations of the pp → V h processes are consistently matched to a PS using the MiNNLO PS method.

SM calculation
A core input of the NNLO QCD calculation are the squared matrix elements up to O(α 2 s ) in the SMEFT.To better explain how the calculation of these objects is performed, we first revisit the structure of the NNLO computation in the SM, which we have repeated, and also implemented into the POWHEG-BOX.Before doing so, we note that we will generically refer to the process pp → V h (and its corresponding subprocesses) in both the text and corresponding figures in what follows, but it should be understood that V = W, Z refers to a final-state lepton pair, and that the calculation does include spin-correlation effects in the gauge-boson decays.
In the NNLO calculation of pp → V h, the contributing partonic channels can be classified according to the number of external quark lines (A = 0, B = 1, C = D = 2), external gluons, and also by the number of loops at the squared amplitude level -see also [25].Starting with the B-type corrections (i.e.those with a single external quark line), the required squared matrix elements are called B0g0V, B1g0V, B0g1V, B1g1V, B2g0V, B0g2V, where the number before g refers to the number of additional external gluons relative to the leading-order (LO) contribution for that type, and the number after the g refers to the number of loops at the squared level.For example, the left most diagram in Figure 2 contains one additional gluon (relative to the Born-level contribution in the quark-antiquark fusion or qqF channel) and is a one-loop graph, and would therefore contribute to B1g1V through interference with the corresponding tree-level amplitude.In the case of the SM, the analytic expressions for the corresponding spinor-helicity amplitudes can be found in [25][26][27][28].To obtain the desired squared matrix elements, the spinor-helicity amplitudes can be squared and then summed over all contributing helicities numerically -an explicit example of this procedure is given below.The C-type corrections feature two external quarks lines, or in other words a double real emission contribution with two final-state quarks, and arise from the interference of diagrams such as that represented in the centre of Figure 2. The D-type corrections account for the additional structures that can appear when same-flavour quarks are considered.These squared matrix elements are called C0g0V and D0g0V and within the SM the analytic expressions for the corresponding spinor-helicity amplitudes are provided in the work [25].Finally, the ggF contributions shown on the right in Figure 2 constitute the third type of correction which are considered (A-type).They are referred to as A0g2V and the corresponding SM spinor-helicity amplitudes are given in [29].Notice that due to charge conservation the third type of corrections only contributes to the pp → Zh but not the pp → W h process.We add that the corrections called V I,II and R I,II that are related to top-quark loops and involve one external quark line [30] are neglected in our SM calculation.Since in total the numerical effect of these contributions amounts to only around 1% [18,30,31], ignoring the V I,II and R I,II terms seems justified at present.
The corresponding calculation including the impact of SMEFT operators (which will be discussed in the following subsection) can also be performed using spinor-helicity techniques.That calculation requires new helicity amplitudes which can (in part) be obtained from knowledge of the SM amplitudes.For clarity of explanation, it will be useful to first consider an explicit example in the SM.To do that we consider the case of B1g0Z which involves a single external quark line and one external gluon at tree level.A corresponding SM Feynman diagram is displayed on the left-hand side in Figure 3.Note that we consider the leptons (quarks) to be outgoing (incoming).The corresponding spinor-helicity amplitude with left-handed fermion chiralities and a physical gluon with a negative helicity reads where ⟨ij⟩ and [ij] denote the usual spinor products -see for example [32] for a review of the spinor-helicity formalism.Notice that the semicolon in the expression on the left-hand side of (3.1) separates the particles with incoming and outgoing convention, respectively.
The amplitudes for the remaining helicity combinations can be obtained via the following parity and charge conjugation relations The resulting spin-averaged matrix element B1g0Z then takes the form where are the usual Mandelstam invariants with p i the four-momentum of particle i.We have furthermore introduced Zf and g hZZ represent the Zf f and hZZ coupling strengths, respectively.The explicit expressions for these quantities are given in Appendix A.
To compute the SMEFT contributions that involve modified couplings between the Higgs and two vector bosons, it is important to notice that by using the spinor identity the result (3.1) can be rewritten as Here the spinor-helicity amplitude corresponding to the q qg subprocess with the indicated helicities is given by where the structures A µ hZZ and A µ hγZ encode the modified hZZ and hγZ vertices, respectively, and p 123 denotes the four-momentum of the incoming vector boson.The symbols g hq γq are the γq q coupling strengths while δg hZZ , δg hγZ and δg (2) hγZ are anomalous couplings that describe the interactions between the Higgs boson and the relevant vector bosons as indicated by the subscript.The explicit expressions for all the couplings appearing in (3.9) can be found in Appendix A. We stress that although the anomalous couplings δg  hγZ do not receive corrections from the Wilson coefficients C HB , C HW and C HWB our POWHEG-BOX implementation contains the full generalised neutral currents (3.9).The presented MC code can therefore be used to extend the Higgsstrahlungs computations in the anomalous-coupling framework [19][20][21] to the NNLO+PS level.
By looking at (3.7) and (3.9) it is now readily seen that in order to obtain the spinaveraged squared matrix element B1g0Z that contains the contributions from the SM as well as the Wilson coefficients C HB , C HW and C HWB one just has to replace the coupling and propagator dressed helicity amplitude appearing in the modulus of (3.3) by the following spinor contraction A schematic depiction of (3.10) is given on the right in Figure 3. Notice that all helicity configurations of A µ qgq can be obtained from (3.7) and (3.8) using the relations (3.2) while in the case of A µ hZZ and A µ hγZ one just has to perform the replacements Insertions of the operators (2.2) and (2.3) lead to the Feynman diagrams shown on the right-hand side in Figure 1 at tree level.In order to capture this contribution in the case of the squared matrix element B1g0Z, one simply has to add the following term to the corresponding dressed SM amplitude appearing within the modulus of (3.3).The analytic expressions for the couplings δg Hq , Q Hd and Q Hu , while the second term is induced by Hℓ and Q He .Notice that compared to the corresponding SM contribution in (3.3) the SMEFT correction proportional to δg (1)hq hZq in (3.11) is missing the Z-boson propagator depending on s 123 .This feature explains the high-energy growth [11][12][13] of the SMEFT pp → V h amplitudes involving the Wilson coefficients C Hq , C Hd and C Hu .
The last type of SMEFT corrections to the squared matrix element B1g0Z is associated to the tree-level shifts of the SM parameters and couplings.EW input scheme corrections from C HD , C HWB , C Hℓ and C ℓℓ lead to the shifts δg 1 , δg 2 and δv of the U (1) Y and SU (2) L gauge coupling and the Higgs vacuum expectation value (VEV), respectively, that in turn induce the shifts δg (0) hZZ and δg (0)h f Zf in the respective couplings of the Z boson.The expressions for these shifts are listed in Appendix A. In practice, the input scheme corrections can be accounted for by applying the replacements g hZZ → g hZZ + δg (0) hZZ and g Hq , C Hd , C Hu , C Hℓ , C Hℓ and C He lead to direct shifts in the Zf f couplings that we include through the shifts g 3).The expressions for the latter shifts are again given in Appendix A.
While we have used the spinor-helicity amplitudes A B1g0Z in this section as examples to illustrate the general approach that we have employed in our SMEFT calculation of pp → V h at NNLO in QCD, it is important to realise that the computation of all other spinor-helicity amplitudes and squared matrix elements proceeds in an analogous manner.In the case of the SMEFT corrections arising from the choice of EW input scheme as well as those associated to insertions of the operators (2.2) and (2.3) this is clear in view of the factorisation properties of these contributions cf.(3.11) .Likewise, since the spinorhelicity structure of the partonic part of a given SMEFT amplitude remains the same as in the SM, it can simply be extracted from the SM expressions and contracted with the part of the SMEFT amplitude that does change.It follows that an explicit calculation of the full partonic structures for the higher-order corrections to pp → V h in the SMEFT can always be avoided, because the relevant amplitudes can be obtained from those in the SM by applying relations à la (3.7) and (3.10) which involves only spinor algebra.Still the calculation of all spinor contractions needed to achieve NNLO accuracy for the pp → V h processes in the SMEFT is a non-trivial task and the final expressions for the spinor-helicity amplitudes turn out to be too lengthy to be reported here.All algebraic manipulations of spinor products needed in the context of this work have been performed with the aid of the Mathematica package S@M [33].

Squared matrix element library
We provide all squared matrix elements discussed in the previous sections in a self-contained Fortran library [34].Our library includes the spinor-helicity amplitudes for the dimensionfour SM and dimension-six SMEFT contributions as well as the definitions for the couplings and the propagators depending on the EW input scheme, which are combined and evaluated numerically in the squared matrix elements up to the desired SMEFT power counting order.Since several of these parts may be relevant in a broader context, we briefly outline the structure of the library.The file Bridge contains the general routines that are required to evaluate the squared matrix elements for a phase-space point, which is represented internally by an object of type Event_t.These routines allow to set up the numerical expressions for the spinor-helicity brackets, pass the input parameters to the Event_t object and calculate the dependent parameters for a chosen EW input scheme.The file squaredamps contains the squared matrix elements.The squared matrix element B1g0Z discussed above, for example, has the form B1g0Z(i1,i2,i3,i4,i5,K,f1,f2), where the integers i1, ..., i5 ∈ {1, ..., 5} allow to specify the crossing of the external legs, K is the Event_t object of the event, and f1 and f2 indicate the flavours of the quark and lepton lines present in the relevant topology, respectively.Our implementation employs the Monte Carlo Particle Numbering Scheme conventions of the PDG [35].The dimension-four, -six and -eight contributions to the squared matrix elements are calculated individually and their inclusion can be controlled via the flags SM, Linear and Quadratic of the Event_t object, respectively.Notice that the dimension-six or linear (dimension-eight or quadratic) SMEFT contributions arise from the interference of the SMEFT and SM amplitudes (self-interference of the SMEFT amplitudes).The spinor-helicity amplitudes, the loop coefficients and the functions implementing the parity and charge conjugation relations are collected in the amplib file.We point out that the required spinor-helicity amplitudes for pp → Z → ℓ + ℓ − and pp → W → ℓν production at NNLO in the SMEFT are just special cases of the amplitudes in our library.They can be obtained by setting up the Event_t object in the exact same way (for example giving it the four physical momenta of the external fermions in the case of B0g0Z), only that in this case the momenta will be linearly dependent because of momentum conservation.
Two further comments seem to be in order.First, besides including the squared matrix elements described above, we also provide the corresponding colour-and spin-correlated squared matrix elements that are required to build the infrared (IR) subtraction terms in the NNLO+PS implementation of pp → V h production.In the case of B1g0V for instance the colour-and spin-correlated squared matrix elements are called B1g0V_colour and B1g0V_spin, respectively.The definition of these squared matrix elements follows the POWHEG conventions specified in (2.6) and (2.8) of the publication [22].While the elements of B1g0V_colour are simply equal to B1g0V times colour factors, calculating B1g0V_spin requires a bit more care.In our notation, it takes the form where the ϵ µ ± are polarisation vectors normalised as with where the form of the four-vectors ϵ µ 1 and ϵ µ 2 can be found in (A.12) of [39].In order to obtain A B1g0Z as given in (3.1), we have however employed and therefore we have to take into account a complex phase See the discussion around (3.22) of [40].We compute the phase φ numerically and by taking it into account in our definition of the polarisation vectors allows us to compare our results for the squared matrix elements with previous implementations that are based on (3.14).
Our implementation is further validated by the fact that the IR poles are properly cancelled in the POWHEG procedure.Second, for the squared matrix element B1g1Z, we follow the discussion presented in the article [28].There, the spinor-helicity amplitude is split into three primitive amplitudes A Ω with Ω = α, β, γ and their corresponding loop coefficients Ω: (3.17) The loop coefficients Ω can be further divided into an IR divergent and an IR finite part as follows with I (1) (ϵ) the usual IR singularity operator as given for example in (C.9) of [28].We include the O(1/ϵ 2 ), O(1/ϵ) and O(1) pieces of I (1) (ϵ) Ω (0) in the array entries 1, 2 and 3 of the loop coefficients Ω, respectively.Note that α (0) = β (0) = 1 and γ (0) = 0.The finite parts Ω (1), finite are instead decomposed as and include the leading colour, the subleading colour and the β 0 = (11C A − 4 T F N f )/6 pieces individually.Here T F = 1/2 and N f = 4 denotes the number of active quark flavours.Analytic expressions for the coefficients Ω (1), finite i with i = 1, 2, 3 were provided as FORM output in the arXiv submission of [28] for the three kinematical regions (i.e.s 13 > 0, s 12 > 0 and s 23 > 0) relevant at hadron colliders.They are expressed as one-and twodimensional harmonic polylogarithms (HPLs) and we translate the appearing HPLs to logarithms and dilogarithms using the relevant formulae in [36,37], computing the latter numerically with the help of LoopTools [38].The finite contributions Ω (1), finite are added to the array entry 3 of the loop coefficients Ω in our code.

NNLO+PS implementation
In the following, we briefly describe how we consistently match the fixed-order NNLO SMEFT calculations of pp → Zh → ℓ + ℓ − h and pp → W h → ℓνh production to a PS.We first recall that within the SM computations of the pp → V h processes have recently reached NNLO+PS accuracy [31,[41][42][43].In fact, here we follow the approach presented in [31] which employs the MiNNLO PS method to match fixed-order QCD and PS effects.The interested reader is referred to the latter work and the articles [23,24] for additional technical details not covered below.
In the case of Higgsstrahlung from qqF, the starting point of our NNLO+PS implementation is the computation of the q q → V h channel in association with one light QCD parton at NLO according to the POWHEG method [44,45], inclusive in the radiation of a second light QCD parton.The computation of the relevant matrix elements has been outlined in Section 3.1 and relies on the SM spinor-helicity amplitudes calculated in [25][26][27][28].In a second step, an appropriate Sudakov form factor and higher-order corrections are applied such that the calculation remains finite in the unresolved limit of the light partons and NNLO accurate for inclusive q q → V h production.In the third step, the kinematics of the second radiated parton (accounted for inclusively in the first step) is generated following the POWHEG method to preserve the NLO accuracy of the V h plus jet cross section, including subsequent radiation through Pythia 8.2 [46].We stress that since all emissions are ordered in transverse momentum (p T ) and the used Sudakov form factor matches the leading logarithms generated by Pythia 8.2, the MiNNLO PS approach maintains the leading-logarithmic accuracy of the PS.The ggF one-loop contributions to Higgsstrahlung, i.e. the gg → Zh process, can instead be computed independently at LO+PS and simply added to the qqF results.Following the general discussion in Section 3.2 the relevant SMEFT spinor-helicity amplitudes can be obtained from those in the SM.We take the SM amplitudes from the MCFM implementation [47] of the results obtained in [29].Further details on our implementation of the gg → Zh process are given in Appendix B.
Our MiNNLO PS generator for q q → V h production has been implemented into the POWHEG-BOX.While it uses the infrastructure of the NNLO+PS SM Higgsstrahlungs generator [31] the parts of the code that calculate the squared matrix elements are entirely new and independent.To validate our implementation of the SM computation we have performed extensive numerical checks against [31].The individual spinor-helicity amplitudes were furthermore compared to a private implementation of the results presented in [25] (which entered the calculation [48]), and results for the squared matrix elements were numerically validated with OpenLoops 2 [49].We exploit the Frixione-Kunszt-Signer subtraction [50,51] to deal with soft and collinear singularities of the real contributions and to cancel the IR poles of the virtual corrections.In fact, the full POWHEG-BOX machinery is used that automatically builds the soft and collinear counterterms and remnants, and also verifies the behaviour in the soft and collinear limits of the real squared matrix elements against their soft and collinear approximations.These checks provide non-trivial tests of our SMEFT calculation of Higgsstrahlung in the qqF channel.In the case of the gg → Zh channel we have instead written a simple LO+PS generator using the POWHEG-BOX framework.Our implementation of the corresponding SM spinor-helicity amplitudes has been validated numerically against OpenLoops 2 at the level of the squared matrix elements.
Besides the direct SMEFT contributions described in Section 3.2 the pp → V h processes also receive corrections from the propagators of the Z and W boson, because SMEFT operators generically modify the masses and the total decay widths of all unstable particles.Our POWHEG-BOX implementation contains the complete tree-level shifts of the relevant masses and total decay widths that are induced by the Wilson coefficients of the operators (2.1) to (2.4) for the α, the α µ as well as the LEP scheme -cf.Appendix A for a comprehensive discussion of the different EW input schemes in the SMEFT.For instance, in the case of the total decay width of the Z and W boson in the LEP scheme we find the relevant shifts Hq − 3.89C (3) Hℓ + 4.98C (3) Hq + 5.31C (3) These results agree with those presented in the literature cf. for instance [52] .We stress that we do not perform an expansion in the SMEFT corrections to the propagators of the Z and W boson.In consequence, the dependence on the Wilson coefficients of our numerical results is generally non-linear.The resulting non-linear effects are however always very small.

Anatomy of SMEFT effects
In this section, we motivate simple SMEFT benchmark scenarios by discussing the leading experimental constraints on the Wilson coefficients of the dimension-six operators intro-duced in (2.1) to (2.4).In this discussion the choice of an EW input scheme is an important ingredient.While in our POWHEG-BOX implementation the user can choose between the α, α µ , and LEP schemes (see Appendix A for a brief discussion of EW input schemes in the SMEFT) the following discussion is based on the LEP scheme which uses as inputs {α, G F , m Z }.Here α is the fine-structure constant, G F is the Fermi constant as extracted from muon decay and m Z is the mass of the Z boson in the on-shell scheme.
In the LEP scheme, the weak mixing angle θ w , the U (1) Y and SU (2) L couplings g 1 and g 2 and the VEV of the Higgs boson v can all be written in terms of the EW input parameters {α, G F , m Z }.One finds the following relations and where the given numerical values correspond to α = 1/127.951,G F = 1.1663788•10 −5 GeV −2 and m Z = 91.1876GeV [35].Notice that in (4.1) and (4.2) we have used the abbreviations s w and c w to denote the sine and cosine of the weak mixing angle, respectively.

W -boson mass
The on-shell mass of the W boson is a predicted quantity in the LEP scheme as well.In terms of the results (4.1) and (4.2) one has 3) The corresponding relative tree-level shift of m W induced in the U (3) 5 symmetric SMEFT is given by where we have employed (4.1), (4.2) and assumed Λ = 1 TeV for the energy scale that suppresses the dimension-six operators to obtain the numerical results presented in the second line.Employing the latest world average of the measured W -boson mass obtained by the PDG [35] and the state-of-the-art SM prediction [55], we obtain the following 95% confidence level (CL) limit on the allowed relative shift of m W .This result in general sets stringent constraints on the Wilson coefficients entering (4.4).For instance, in the case of the operator Q HWB we find the bound if all other Wilson coefficients in (4.4) are set to zero.Similar though weaker limits also apply in the case of C HD , C Hℓ and C ℓℓ if each of the Wilson coefficients is treated as the only non-zero contribution.

Z-boson couplings
The linear combinations of Wilson coefficients that modify the couplings of the W and Z boson to fermions at tree level are strongly constrained by the EW precision measurements performed at LEP and SLD.In the case of the Z boson the shifts of the left-and righthanded couplings to a fermion ψ = {e, ν, d, u} take the following form (1) with Here the coupling g 2 is defined via (4.2), ψ L = {ℓ, q} and ψ R = {e, d, u}, while T 3 ψ = ±1/2 denotes the weak isospin eigenvalue, Q ψ is the electric charge of the fermion ψ in units of e = √ 4πα.Notice that the expressions given in (4.8) correspond to the LEP scheme and that we have assumed that the Wilson coefficients of the operators that involve fermionic fields are flavour universal.Below we will assume that such a flavour universality holds for all SM fermions apart from the charm, bottom and top quark.
In the case of the electron the formulae (4.7) and (4.8) lead to the following left-and right-handed coupling shifts δg e L ≃ 0.036 C HWB − 0.022C Hℓ + 0.020C (3) when (4.1), (4.2) and Λ = 1 TeV are used as numerical inputs.The EW precision measurements performed by LEP and SLD [35,56] imply that at 95% CL one has Ignoring cancellations these limits again put severe constraints on the Wilson coefficients that appear in (4.9).Numerically, we obtain for instance the bounds from the constraint on δg e L and δg e R , respectively, if we allow only the considered Wilson coefficient to take a non-vanishing value.We add that the Wilson coefficients C HD , C Hℓ and C ℓℓ can also be constrained by the LEP and SLD measurements of the Z-boson coupling to neutrinos.The obtained bounds, however, turn out to be weaker than those that derive from (4.9) and (4.10).
Employing again (4.1), (4.2), (4.7) and (4.8) the left-and right-handed coupling shifts in the case of the down and up quark take the following form Hq + C (3) Hq − C (3) for a common operator suppression scale of Λ = 1 TeV.The measurements by LEP and SLD performed at the Z-pole lead to the following limits at 95% CL.These bounds translate into if each Wilson coefficient is treated independently and the limit on δg d L as given in (4.13) is used to constrain C (1)

Hq and C
(3) Hq .We emphasise that notably more stringent limits on the Wilson coefficients in (4.14) would be obtained if one were to perform a complete EW fit including heavy flavour measurements and assume a SMEFT Lagrangian with an approximate U (3) 5 symmetry.This is a simple consequence of the fact that in such a case one has ) and that the bottom (charm) couplings were significantly better determined than the down (up) couplings at LEP and SLD -see for example Figure F.3 in [56].The limits (4.14) are relaxed because they assume flavour universality only for the down and strange quark as far as quarks are concerned.

Higgs-boson observables
The Wilson coefficients of the operators (2.1) modify the Higgs signal strengths in final states with two vector bosons.In the LEP scheme and including only the tree-level SMEFT corrections due to C HB , C HW and C HWB the coupling modifiers relevant for h → W W and h → ZZ can be written in the κ framework as follows [57] (4.15) In the tree-level approximation the corresponding expressions for the h → γγ and h → γZ decays do not depend on the choice of the EW input scheme.One finds where g hγγ ≃ −2.02 • 10 −3 and g hγZ ≃ −7.23 • 10 −3 parameterise the loop-induced hγγ and hγZ couplings within the SM.Explicit analytic formulae for g hγγ and g hγZ can be found for instance in [58].
Assuming that the ggF Higgs production cross section is SM-like and taking the SM predictions for the total decay width of the Higgs and its branching ratios from [59], we obtain for Λ = 1 TeV the following linearised expressions for the relevant signal strenghts Here the first three results represent unofficial weighted averages of the ATLAS [60] and CMS [61] measurements, while the fourth result stems from an official combination performed by the ATLAS and CMS collaborations in [62].By comparing (4.17) and (4.18) it is evident that under the assumption of C HWB ≃ 0 the constraints on the Wilson coefficients C HB and C HW from h → γγ and h → γZ are in general more stringent than those that arise from h → W W and h → ZZ.In fact, since the signal strength µ γγ ggF is at present significantly better measured than µ γZ ggF , all combinations of Wilson coefficients that obey are most weakly constrained by (4.18).Notice that for Wilson coefficients C HB and C W B satisfying (4.19) and |C HW | sufficiently small the signal strengths µ W W ggF , µ ZZ ggF and µ γγ ggF are all predicted to be SM-like, while µ γZ ggF can at the same time be notably different from one.

Discussion
We are now in a position to identify simple benchmark scenarios for the full set of Wilson coefficients of the operators defined in (2.1) to (2.4).As a common operator suppression scale we hereafter take Λ = 1 TeV.In view of the stringent constraints arising from (4.4), (4.5), (4.9) and (4.10) we always employ

C
(1) when studying the numerical impact of the SMEFT dimension-six operators (2.1) to (2.4) on Zh production at the LHC in Section 5.
In the case of the dimension-six operators that modify the couplings between the Higgs and two vector bosons at tree level, we choose the following values for the Wilson coefficients In the case of the operators that result in couplings between the Higgs, a W or a Z boson, and light quarks, we consider the four benchmark scenarios C Hu = 0.1 , C Hq = C (3) C (1) where one of the four relevant Wilson coefficients is non-zero while the remaining ones vanish identically.As in the case of (4.21) we set all Wilson coefficients not present in (4.23) to (4.26) to zero when considering the class of operators introduced in (2.2).

Phenomenological analysis
In the following, we present NNLO+PS accurate results for pp → Zh → ℓ + ℓ − h production with a stable Higgs boson including the SMEFT effects discussed in Section 3.All SM input parameters are taken from the most recent PDG review [35].In particular, we use The values of the weak mixing angle, the U (1) Y and SU (2) L gauge coupling and the Higgs VEV are calculated from (4.1) and (4.2), respectively, meaning that we work in the LEP scheme.The NNPDF31_nnlo_as_01180 parton distribution functions (PDFs) [63] are employed in our MC simulations and events are showered with Pythia 8.2 [46] utilising the Monash tune [64].To target the Z → ℓ + ℓ − decay, we select events with two charged leptons (electrons or muons).The leptons are defined at the dressed level, meaning the lepton four-momentum is combined with the four-momenta   [23,24] and the grey bands in the lower panels represent the corresponding perturbative uncertainties in the SM.These uncertainties have been obtained from seven-point scale variations enforcing the constraint 1/2 ≤ µ R /µ F ≤ 2. Within the SM they do not exceed 5% for what concerns the considered distributions, and relative scale uncertainties of very similar size are also found in the case of the SMEFT spectra.The SM predictions are indicated by the solid black lines while the solid (dotted) orange curves represent the SMEFT contributions that are linear (quadratic) in the Wilson coefficients.Notice that the linear (quadratic) SMEFT contributions arise from the interference of the SMEFT and SM amplitudes (self-interference of the SMEFT amplitudes).The solid dark orange lines finally correspond to the sums of the linear and quadratic SMEFT contributions.We observe that in the case of the |η Z | distribution the SMEFT effects related to C HB and C HW to first approximation simply shift the spectra by a constant amount.In the cases of the p T,Z , the |η Z − η h | and the m Zh spectra the relative sizes of the SMEFT corrections instead grow with increasing p T,Z , |η Z − η h | and m Zh , respectively.It is also evident from all panels that the quadratic SMEFT contributions are negligibly small compared to the linear terms.Numerically, we find that the benchmark scenario (4.21) leads to relative corrections of around +2% to +5% in the studied distributions -the fiducial cross section is enhanced by +3.5% compared to the SM.Notice that while the predicted deviations are sometimes larger than the corresponding SM scale uncertainties, they are typically smaller than the ultimate projected HL-LHC accuracy in Zh production channel that amounts to 5% [3,4].From this numerical exercise one can conclude that future constraints on the Wilson coefficients C HB , C HW and C HWB from V h production are unlikely to be as stringent as the limits that future determinations of the Higgs signal strengths in h → γγ and h → γZ will allow to set.This will in particular be the case if both Higgs signal strengths turn out to be SM-like see (4.17) and (4.18) .
Figures 5 to 7 contain our NNLO+PS predictions for pp → Zh → ℓ + ℓ − h production in the SMEFT benchmark scenarios (4.23) to (4.25).In all cases we have employed an operator suppression scale of Λ = 1 TeV.The results depicted in the first three figures show very similar features.In all cases the relative SMEFT corrections are rather flat in the |η Z | and |η Z − η h | distributions, while in the case of the p T,Z and m Zh spectra they get larger with increasing p T,Z and m Zh .The observed high-energy growth is expected from (3.11) and in all three cases most pronounced in the p T,Z distribution.One also sees that the Hq and δg u L ∝ C (3) Hq -see (4.12).Notice also that the size of the quadratic SMEFT corrections is relatively small in the case of (4.23) while these effects are comparable to or even larger than the linear terms for the benchmark scenarios specified  in (4.24) and (4.25).In the case of the SMEFT benchmark scenario (4.26) the pattern of SMEFT deviations turns out to be more complicated.This is illustrated in the four panels of Figure 8.One observes that the linear SMEFT effects are typically small and even change sign in some distributions as indicated by the transitions from solid to dashed lines.As Figure 4 but for benchmark scenario (4.26) assuming Λ = 1 TeV.The SMEFT predictions are coloured in magenta.In the case of the linear SMEFT contributions the solid (dashed) lines correspond to positive (negative) corrections to the relevant distribution.
To understand these features it is important to realise that g Hq and to keep in mind that the down-quark luminosity in a proton is smaller than the up-quark luminosity at large x while the two luminosities are of similar size at small x.For the choice C (1) Hq = 0.05 the down-and up-quark contributions thus tend to cancel leading to a numerical suppression of the full linear SMEFT effects compared to naive expectation.It is also evident from the shown results that the quadratic SMEFT corrections are as large in magnitude as the linear terms.In fact, in the case of the p T,Z and m Zh spectra the two types of SMEFT effects have opposite relative signs resulting in very small combined BSM contributions not exceeding the level of +2% in the case of the p T,Z distribution.Notice finally that due to the energy growth most of the sensitivity to the SMEFT effects considered in Figures 5 to 8 comes from the high-energy tails of kinematic distributions such as the p T,Z and m Zh spectra.In such a case the Higgs-boson decay products are significantly boosted, giving rise to very specific kinematic features and providing additional handles to distinguish signal from background events.The articles [11][12][13] have exploited this feature to obtain stringent constraints on the dimension-six operators (2.2) using future hypothetical hadron collider measurements of V h production.The MiNNLO PS generator presented in this work would allow to improve the accuracy of these studies from the NLO+PS to the NNLO+PS level.We leave such an analysis for future study.

Combining production and decay at NNLO+PS
Let us finally outline how one can combine the NNLO+PS calculation of V h production with the decay of the Higgs boson to bottom quarks at NNLO+PS.The starting point is to replace the Higgs boson in each V h event by the possible decay products (i.e.b bg, b bq q or b bgg) taken from a given decay event.Working in the narrow width approximation (NWA), the full weight w full of each event is then calculated as where w SMEFT prod denotes the weight of the production event obtained at NNLO+PS using the MiNNLO PS pp → V h generator described earlier in this section, while w SMEFT dec is the weight of the decay event computed at NNLO+PS by means of the MiNLO ′ method [53,54] employing the procedure detailed in Section 3.2 of the paper [43].Notice that insertions of the operators (2.1) to (2.4) do not modify the kinematics of any Higgs decay channel involving bottom-quark pairs.In fact, the SMEFT decay weights are simply related to the SM decay weights by where the term 1 + 2c kin arises from the canonical normalisation of the Higgs kinetic term.The factor Γ SMEFT h in (6.1) takes into account that the total decay width of the Higgs boson is modified by SMEFT effects.For example, one can employ the result where we have factorised the contribution 1 + 2c kin due to the canonical normalisation of the Higgs kinetic term and Γ SM h denotes the total decay width of the 125 GeV Higgs boson within the SM.Notice that the factor 1+2c kin drops out in the ratio of (6.2) and (6.3) and therefore in (6.1) when production and decay events are combined.Based on (6.1) for each event the Les Houches event (LHE) file in addition to the weight also contains the value of the hardest radiation allowed by the shower (i.e. the scale scalup).In the combined LHE file the value of scalup prod for the production process is recorded and, once the event is passed to the PS, the value of scalup dec for the specific decay kinematics is recomputed.Emissions of the Higgs decay in all the available phase space is then generated and after the shower is complete it is checked whether the hardness of the splittings is below the veto scale scalup dec .If this is not the case the event is again showered until the latter condition is met.
Notice that while above we have considered the SMEFT corrections to h → b b that arise from (2.1) to (2.4) the modifications required to achieve NNLO+PS accuracy for other operators is relatively straightforward as long as one works in the NWA for the Higgs propagator.A non-trivial example for the combination of Higgsstrahlung and h → b b at NNLO+PS has been provided in the recent article [18].In fact, this work achieved NNLO+PS precision for the dimension-six operators that contribute to the subprocesses pp → Zh and h → b b directly in QCD.This class of operators includes effective Yukawa-and chromomagnetic dipole-type interactions of the bottom quark that modify the h → b b decay but do not play a role in pp → Zh production.Combining the results obtained here and in [18] via (6.1) it is now possible to obtain NNLO+PS accurate predictions for the full pp and pp → W h → ℓνb b processes taking into account 17 different dimension-six SMEFT operators.Since in the NWA the Higgsstrahlungs production processes including the leptonic decays of the massive gauge bosons factorises from the subsequent decay of the Higgs boson, extending the general method of combining production and decay (6.1) to incorporate additional operators modifying h → b b or considering further Higgs-boson decay modes is relatively simple.

Conclusions
In this article, we have presented novel SMEFT predictions for Higgsstrahlung in hadronic collisions.Specifically, we have calculated the NNLO QCD corrections for the complete sets of dimension-six operators that describe the interactions between the Higgs and two vector bosons and the couplings of the Higgs, a W or a Z boson, and light fermions.These fixedorder predictions have been consistently matched to a PS using the MiNNLO PS method and the matching has been implemented into the POWHEG-BOX.Our new MC implementation allows for a realistic exclusive description of V h production at the level of hadronic events including SMEFT effects.This feature makes it an essential tool for future Higgs characterisation studies by the ATLAS and CMS collaborations, and we therefore make the relevant code available for download on the official POWHEG-BOX web page [68].Notice that together with the MC code presented in [18] one can now simulate the pp → Zh → ℓ + ℓ − b b and pp → W h → ℓνb b processes at NNLO+PS including a total number of 17 dimension-six SMEFT operators.
To motivate simple SMEFT benchmark scenarios we have discussed the leading constraints on the Wilson coefficients of the 13 dimension-six operators listed in (2.1) to (2.4).In the case of the effective interactions between the Higgs boson and two vector bosons described by Q HB , Q HW and Q HWB , we found that the combination of the measurements of the W -boson mass and the Higgs signal strength in h → γγ and h → γZ allow to set stringent constraints on the respective Wilson coefficients.In order to avoid these bounds we have introduced the benchmark scenario (4.21) which as far as the Higgs signal strengths are concerned predicts h → W W , h → ZZ and h → γγ to be SM-like, while enhancing the h → γZ rate by a factor of around 2 with respect to the SM.Such a pattern of deviations is presently favoured by LHC data cf.(4.18) .In the case of the dimension-six terms that give rise to couplings between the Higgs, a W or a Z boson, and light fermions, we have pointed out that while the leptonic operators (2.3) are in general tightly constrained by LEP and SLD measurements the resulting bounds on the quark operators (2.2) depend sensitively on the flavour assumption in the SMEFT.While in the case of an approximate U (3) 5 flavour symmetry the constraints turn out to be stringent, the limits on the Wilson coefficients of (2.2) are notably relaxed if the assumption of flavour universality only applies to the light but not the heavy down-and up-quark flavours.The remaining operators (2.4) shift the Higgs kinetic term and/or the EW SM input parameters at tree level.The corresponding Wilson coefficients are therefore in general better probed by a global SMEFT fit than by pp → V h production alone.
We have then performed an NNLO+PS study of the impact of SMEFT contributions on several kinematic distributions in pp → Zh → ℓ + ℓ − h production for a stable Higgs boson considering the simple benchmark scenarios identified earlier.While in our POWHEG-BOX implementation the user can choose between the α, α µ , and LEP schemes, for concreteness our discussion was based on the LEP scheme which uses {α, G F , m Z } as EW input parameters.Another feature of our MC code worth highlighting is that it is able to compute separately both the SMEFT corrections that are linear and quadratic in the Wilson coefficients.Our numerical analysis showed that once the stringent constraints from m W , h → γγ and h → γZ are imposed the numerical impact of C HB , C HW and C HWB on the kinematic distributions in pp → Zh → ℓ + ℓ − h are rather limited, amounting to relative deviations of no more than 5%.Future limits on the Wilson coefficients of the effective interaction in (2.1) from V h production are therefore unlikely to be competitive with the limits that future determinations of the Higgs signal strengths in h → γγ and h → γZ will allow to set.This will in particular be the case if the latter measurements turn out to be SM-like.The situation turns out to be more promising in the case of the operators (2.2) that induce couplings of the Higgs, a W or a Z boson, and light quarks.The sensitivity of the process pp → Zh → ℓ + ℓ − h to the Wilson coefficients C Hq , C Hd and C Hu arises from the energy growth of the respective amplitudes which results in enhanced highenergy tails of kinematic distributions such as the p T,Z and m Zh spectra.Numerically, we found that these enhancements can reach 50% in the p T,Z spectrum in the region where the Higgs-boson decay products are significantly boosted.As shown in the papers [11][12][13], future hypothetical HL-LHC measurements of V h production can therefore provide constraints on the Wilson coefficients C (1) Hq , C Hd and C Hu that are competitive with the bounds obtained from projected global SMEFT fits.Utilising the MiNNLO PS generator presented in this work would allow to improve the accuracy of the studies [11][12][13] from the NLO+PS to the NNLO+PS level.
Notice that in our phenomenological analysis we have focused on the 0-jet categories of the stage 1.2 simplified template cross sections (STXS) framework [69][70][71] for the V h production processes.However, it is important to realise that our POWHEG-BOX implementation of the Higgsstrahlungs processes also allows to simulate the different 1-jet STXS categories with NLO+PS accuracy.This represents an important improvement compared to the MC code presented in [10] or the SMEFT@NLO package [72] which are only LO+PS accurate for 1-jet observables in V h production.Another novel feature of our implementation of Higgsstrahlung is that it is able to combine production and decay consistently at NNLO+PS including SMEFT effects.To our knowledge, this is presently not possible with any other publicly available tool even if one only aims at NLO+PS precision.Finally, the presented squared matrix element library contains all spinor-helicity amplitudes that are needed to obtain NNLO+PS predictions for Drell-Yan production taking into account the effects of the dimension-six operators (2.2) and (2.3).Modifying the code such that one can calculate the SMEFT effects in diboson production at NNLO+PS due to operators that induce anomalous triple gauge couplings is also relatively straightforward.
Table 1: The parameters g 1 , g 2 and v expressed in terms of the input parameters for the three EW input schemes implemented in the POWHEG-BOX code.
where α is the fine-structure constant, G F is the Fermi constant as extracted from muon decay and m Z (m W ) is the mass of the Z (W ) boson in the on-shell scheme.The relevant expressions for the U (1) Y and SU (2) L gauge couplings g 1 and g 2 and the Higgs VEV v in terms of the EW input parameters are given in Table 1 for the α, the α µ , and the LEP scheme.
In terms of the parameters g 1 , g 2 and v the Zf f , γf f and hZZ coupling strengths take the following form in the SM Notice that these relations are independent of the employed EW input scheme.Here the symbol Y f represents the weak hypercharge, T 3 f is the third component of the weak isospin and Q f denotes the electric charge.The fermions are f = q, ℓ with q = d, u and ℓ = e, ν, and the helicity states f + and f − are identical to the chirality states f R and f L in the massless limit.
The relations among the EW input parameters and g 1 , g 2 and v are modified at tree level by the presence of some of the dimension-six SMEFT operators listed in (2.1) to (2.4), leading to so-called input scheme corrections.These can be accounted for via the shifts x → x + δx for x = g 1 , g 2 , v. We summarise the relevant shifts in Table 2.The input scheme corrections δg 1 δg 2 and δv themselves lead to the shifts δg (0)± Zf and δg (0) hZZ of the Zf f and hZZ couplings, respectively.We find the following scheme-independent results At the same time, the SMEFT operators listed in (2.1) to (2.4) give direct contributions to the Z-boson couplings to two gauge bosons.We find the following analytic expressions Furthermore, we obtain δg (1) Hq , δg Hq − C Hq , δg (1) Hℓ , δg (1) Hℓ , δg are the relevant direct SMEFT corrections to the Zf f couplings.

B SMEFT corrections to gg → Zh process
Our calculation of gg → Zh production is based on the spinor-helicity amplitudes for the SM derived in [29] and implemented into MCFM [47].In unitary gauge, the expression for the triangle contributions with positive gluon helicities and left-handed fermion chiralities reads Notice that we have followed the convention of [29] and written the amplitude for all momenta outgoing.In (B.1) the two terms in the last factor in the first line stem from the transversal and longitudinal part of the Z-boson propagator in unitary gauge, respectively, q is the quark running in the loop with mass m q and C 0 is the scalar Passarino-Veltman (PV) triangle integral defined as in [38,74].The corresponding SM Feynman diagram is displayed on the right-hand side in Figure 2. Similarly, we have implemented the box amplitudes which are, however, too lengthy to be reported here but may be inspected in our squared matrix element library.The remaining non-zero helicity combinations may be obtained via parity and charge conjugation relations.In the case of the triangle contributions, these relations take the form where the overline means that the brackets should be exchanged, i.e. [. ..] ↔ ⟨. ..⟩. Analog relations hold for the box contributions including the cases where the gluons have opposite helicities which are only present for A q A0g2Z .Including the triangle and box contributions the resulting spin-averaged matrix element takes the form Here D Z (s) has been defined in (3.5) while the expressions for the couplings g ± Zf and g hZZ can be found in (A.2).The coupling g hZZ appearing in (B.6) requires some explanation.In fact, the box contributions do not involve a hZZ vertex but instead the Higgs boson couples directly to the quarks.However, since with a factor m q /v coming from the hq q vertex and another m q stemming from the mass insertion in the box diagram the expected mass dependence for A q A0g2Z is recovered.It is important to realise that as a result of the generalised Furry theorem the vectorcurrent coupling of the Z boson, which is proportional to the combination (g − Zq + g + Zq ) of couplings, does not contribute to the spin-averaged matrix element A0g2Z as given in (B.4).However, the axial-current part contributes, as signalled by the factor (g − Zq − g + Zq ) in both (B.5) and (B.6), and this contribution is directly connected to the U (1) A × SU (3) c gauge anomaly.In fact, a regulator and a loop routing scheme must be introduced to properly define the amplitude A q A0g2Z△ , rendering its expression scheme-dependent -for a detailed explanation of this point see for example [75].Within the SM, the axial parts of the top-and bottom-quark couplings obey and as a result all gauge anomalies are cancelled.It follows that the sum over q that appears in (B.4) evaluates to q=t,b and in consequence any scheme-dependent constant shift in the amplitude A q A0g2Z△ drops out in the combination A t A0g2Z△ −A b A0g2Z△ .Notice that in the degenerate or zero mass case the sum (B.9) vanishes identically.Since we treat the light-quark generations as massless, down-, up-, strange-and charm-quark loops hence do not need to be included in the spinaveraged matrix element (B.4).
The amplitudes including the SMEFT contributions to (B.4) were computed with the procedure outlined in Section 3.2.Since the SM amplitudes were derived in unitary gauge, only SMEFT contributions to vertices involving the Z boson have to be considered.We have checked explicitly that in Feynman gauge, the SMEFT effects in the Goldstone diagrams are equivalent to the effects in the longitudinal part of the amplitude in unitary gauge.Another important point is that in the SMEFT, the coupling shifts δg (0)± Zq and δg (1)± Zq as given in (A.3) and (A.6), respectively, in general do not obey (B.8) and the scheme-dependent constant shift in A q A0g2Z△ thus does not automatically cancel in the sum (B.9).One may incorrectly jump to the conclusion that additional cancellation conditions analogous to (B.8) have to be imposed on the SMEFT Wilson coefficients in order to ensure the cancellation of gauge anomalies.However, one can show [76][77][78][79][80] that the anomalous contributions depending on the Wilson coefficients are "irrelevant", in the sense that they can be removed by adding a local counterterm, i.e. a Wess-Zumino term [81].As a result, the condition for the cancellation of "relevant" gauge anomalies in the SMEFT is the same as in the SM and only dependent on the gauge quantum numbers of the fermionic sector, as one would naively expect from an effective field theory point of view.We cancel the unphysical terms in the amplitudes depicted in Figure 9 via the replacement for what concerns the SMEFT contribution to (B.4).Notice that for m 2 q /s 12 → ∞ the PV triangle integral appearing in (B.1) obeys m 2 q C 0 (s 12 , 0, 0, m q , m q , m q ) → −1/2 and hence the subtraction term in (B.10) is local and independent of the mass of the quark in the loop.Since the replacement (B.10) is the same for all q in the sum appearing in A0g2Z, it drops out for an anomaly-free UV theory.At the same time the replacement ensures that A q A0g2Z△ decouples as s 12 /m 2 q in the limit m 2 q ≫ s 12 , recovering the behaviour obtained in the covariant anomaly scheme [72,75,82,83] where the need for a local counterterm is circumvented by an adequate choice of the loop routing scheme.Notice that in contrast to the triangle contributions, the box corrections do not require a subtraction à la (B.10).An example SMEFT box graph is shown in the middle of Figure 9.It follows that, up to an overall factor of 2/v, the counterterms necessary to remove the anomalous SMEFT contributions in the ggZ and gghZ vertices are identical.
Some additional care has to be taken regarding the SMEFT contributions to (B.4) from Feynman diagrams that involve an effective hZq q or hZℓ + ℓ − vertex.A graph of the former type is depicted on the right in Figure 9.In these cases only the transversal part of (B.1) contributes -the longitudinal part vanishes because the Z boson couples directly to the leptons that are treated as massless -and therefore in addition to dropping the factor D Z (s 12 ) in (B.As Figure 10 but for benchmark scenario (4.21).The curves called SM+SMEFT correspond to the full squared matrix elements including the sum of both the SM and SMEFT contributions.The lower panels show the ratios between the SM+SMEFT and the SM predictions at the same order in QCD.For more details consult main text.the case in the SMEFT, where effects stemming from the gg → Zh channel are suppressed compared to the SM as a result of the additional subtraction (B.10).We add that the comparisons of SM+SMEFT predictions present this appendix represent a non-trivial validation of our new NNLO+PS MC code for Higgsstrahlung.

Figure 2 :
Figure 2: Examples of higher-order QCD corrections to pp → Zh production within the SM.The diagram on the left features additional virtual and real gluon lines (B-type), the diagram in the middle involves a second quark line (C-and D-type) and the diagram on the right is a ggF contribution (A-type).Consult the main text for further details.
) with Γ Z denoting the total decay width of the Z boson.In (3.3) the variable α s denotes the strong coupling constant while C F = 4/3 and C A = 3 are the relevant colour factors.The symbols g h f

Figure 3 :
Figure 3: Example graphs that contribute to the B1g0Z matrix elements.The diagram on the left shows a SM contribution.On the right we instead depict a SMEFT correction that receives contributions from the generalised hZZ and hγZ currents introduced in (3.9).The four-momentum flow is indicated by the grey arrows and labels.See the main text for additional explanations.

( 1
)h f hZf are given in Appendix A. In (3.11) the first term in the brackets describes the contribution from Q (1) Hq , Q

Figure 4 :
Figure 4: NNLO+PS predictions for pp → Zh → ℓ + ℓ − h production in the SMEFT benchmark scenario (4.21) assuming a common operator suppression scale of Λ = 1 TeV.The four panels show the fiducial cross section differential in |η Z | (upper left), p T,Z (upper right), |η Z − η h | (lower left) and m Zh (lower right) for proton-proton (pp) collisions at 13 TeV.The SM predictions are indicated by the solid black lines while the solid (dotted) orange curves represent the SMEFT contribution linear (quadratic) in the Wilson coefficients.The solid dark orange lines correspond to the sums of the linear and quadratic SMEFT contributions.The lower panels depict the ratios between the BSM and the SM distributions with the grey band representing the SM scale uncertainties.See main text for further details.

Figure 7 :
Figure 7: As Figure 4 but for benchmark scenario (4.25) with Λ = 1 TeV.The blue lines are the BSM results.

Figure 8 :
Figure 8:As Figure4but for benchmark scenario (4.26) assuming Λ = 1 TeV.The SMEFT predictions are coloured in magenta.In the case of the linear SMEFT contributions the solid (dashed) lines correspond to positive (negative) corrections to the relevant distribution.
hγZ = 0 meaning that the corresponding Dirac structures are not generated at the dimension-six level in the SMEFT.The expressions for the hZf f couplings can finally be written as δg

Figure 9 :
Figure 9: Examples of contributions to gg → Zh production within the SMEFT.All graphs involve an insertion of one of the operators given in (2.2) as indicated by the blue squares.Further details can be found in the main text.

m 2 Z,
(B.11)    from SMEFT diagrams with a hZq q or hZℓ + ℓ − vertex.This contribution can be included by simply adding the expression (B.11) to the sum over q in (B.4).In order to cancel anomalous terms in the amplitude (B.11) one in addition has to include an appropriate counterterm contribution by applying (B.10) to (B.11).

Figure 10 :
Figure 10: SM NLO+PS and NNLO+PS results for pp → Zh → ℓ + ℓ − h production.The |η Z | (upper left), p T,Z (upper right), |η Z − η h | (lower left) and m Zh (lower right) spectra are shown.The dashed (solid) lines illustrate the NLO+PS (NNLO+PS) results, while the dotted curves are the ggF NNLO+PS corrections.The solid (dotted) lines in the lower panels depict the ratios between the full NNLO+PS (NLO+PS plus ggF NNLO+PS) and the NLO+PS results.See main text for further explanations.

Figure 11 :
Figure 11:As Figure10but for benchmark scenario (4.21).The curves called SM+SMEFT correspond to the full squared matrix elements including the sum of both the SM and SMEFT contributions.The lower panels show the ratios between the SM+SMEFT and the SM predictions at the same order in QCD.For more details consult main text.
[8], Q HW and Q HWB generate modified hV V vertices with helicity structures different from the one present in the SM, i.e. the spinor chain ⟨4|γ µ |5] in (3.7).These modifications can be included at the level of (3.1) by means of generalised currents that describe the splitting of the initial vector boson V 1 into the outgoing vector boson V 2 and the Higgs boson h[8].If the initial-state quarks and finalstate leptons are left-handed the relevant generalised neutral currents are given by that the ggF Higgs production cross section in the SMEFT is equal to that of the SM and sets to zero the Wilson coefficients of all the operators appearing in (2.2) to(2.4).Notice that the predictions (4.22) are all SM-like apart from µ γZ ggF which is close to the central value of the measured signal strength given in(4.18).When considering the benchmark scenario (4.21) we will, besides employing(4.20),also set all the Wilson coefficients of the operators (2.2) to zero.

Table 2 :
SMEFT input scheme corrections for the three EW input schemes implemented in the MC code.