Towards the ultimate differential SMEFT analysis

We obtain SMEFT bounds using an approach that utilises the complete multi-dimensional differential information of a process. This approach is based on the fact that at a given EFT order, the full angular distribution in the most important electroweak processes can be expressed as a sum of a fixed number of basis functions. The coefficients of these basis functions - the so-called angular moments - and their energy dependance, thus form an ideal set of experimental observables that encapsulates the complete multi-dimensional differential information of the process. This approach is generic and the observables constructed allow to avoid blind directions in the SMEFT parameter space. While this method is applicable to many of the important electroweak processes, as a first example we study the $pp \to V(\ell\ell)h(bb)$ process ($V \equiv Z/W^{\pm}, \; \ell\ell \equiv \ell^+\ell^-/\ell^\pm\nu$), including QCD NLO effects, differentially. We show that using the full differential data in this way plays a crucial role in simultaneously and maximally constraining the different vertex structures of the Higgs coupling to gauge bosons. In particular, our method yields bounds on the $h V_{\mu \nu}V^{\mu \nu}$, $h V_{\mu \nu}\tilde{V}^{\mu \nu}$ and $h Vff$ ($ff \equiv f\bar{f}/f\bar{f}'$) couplings, stronger than projected bounds reported in any other process. This matrix-element-based method can provide a transparent alternative to complement machine learning techniques that also aim to disentangle correlations in the SMEFT parameter space.


Introduction
The data being collected by the LHC is the first record of interactions of the Higgs and other Standard Model (SM) particles at the sub-attometre (multi TeV) scale. As long as beyond SM (BSM) physics is significantly heavier than the mass of electroweak particles, these interactions can be described in a model independent way by the the Standard Model Effective Field Theory (SMEFT) Lagrangian. The SMEFT Lagrangian is thus a statement of the laws of nature at the most fundamental scale ever probed. The measurement of (or constraints on) the SMEFT parameters  may well turn out to be the main legacy of the LHC after the Higgs discovery.
It is thus of great importance to maximally exploit all the data that the LHC would provide us. To constrain the SMEFT Lagrangian, it is especially important to extract the full multi-dimensional differential information available in a process. This is because the effect of new vertex structures arising at the dimension-6 (D6) level is often more pronounced in certain regions of the phase space, the most common example being the growth of EFT rates at high energies. A more subtle example is that of operators whose contributions do not interfere with the SM amplitude at the inclusive level [33]. These operators can generate large excesses differentially [32,[34][35][36][37] in certain regions of the phase space, which are cancelled by corresponding deficits in other regions. These effects can, therefore, get lost unless a sophisticated study is carried out to isolate these phase space regions. As discussed in Ref. [32], and as we will also see in this work, sometimes in order to resurrect these interference terms one has to go even beyond differential distributions with respect to a single variable and use multidimensional distributions. More generally, using the full differential information enlarges the list of observables and lifts flat directions in EFT space that can otherwise remain unconstrained. In order to optimally reconstruct the SMEFT lagrangian, it is thus essential to systematically and completely extract all the available differential information.
In the way experimental measurements are communicated, there is a large reduction in differential information, as often only a few intuitively chosen distributions are presented. To estimate this, consider a three body final state where the phase space in the center of mass frame can be completely described by four variables: an energy variable and three angles. For a given energy, taking for instance 10 bins for each of the angular variables results in 1000 units of data to capture the entire information contained in this process, at this level of experimental precision. However, often individual angles are analysed in isolation and the correlations contained in the full set of data are projected onto only 30 units of data, i.e. 10 for each angle, resulting in a loss of accessible information to search for new physics contributions.
Interestingly, for many important processes the 1000 units of data, contain redundant information. We argue, that with an understanding of the underlying theoretical structure of process the number of physical quantities required to completely characterise the full differential distribution can be drastically reduced. The main fact that we will utilise in this work is that, for some of the most important processes in Higgs and electroweak physics, the full angular distribution at a given energy can be expressed as a sum of a fixed number of basis functions as long as we limit ourselves to a certain order in the EFT expansion. The reason for this is that only a finite number of helicity amplitudes get corrections up to the given EFT order, see for instance Ref. [38,39]. The coefficients of these basis functions, the so called angular moments [40][41][42][43], and their energy dependance, thus, contain the full differential information available in a process. The effect of EFT operators on differential distributions can therefore be summarised by their contribution to these angular moments. As such angular moments can be used to construct any possible differential distribution, an analysis utilising them has the potential to reach maximal sensitivity in probing EFT coefficients.
These methods would complement other techniques that aim to employ a maximum-information approach, e.g. the matrix element method [49][50][51][52] or machine learning techniques that have recently gained popularity [53][54][55][56][57]. One advantage of this approach over other multivariate techniques is its more physical and transparent nature. The angular moments described above can be directly related to physical experimental quantities, e.g. they have well defined symmetry properties, than the abstract neural network outputs used in machine learning approaches. Another important distinction of the methods proposed here from some multivariate approaches like the matrix element method, is that the process of extraction of the angular moments is hypothesis-independent; for instance it would be independent of our assumptions about whether electroweak symmetry is linearly or non-linearly realised.
In this work we will show how these angular moments can be extracted and mapped back to the EFT lagrangian. While in this study we will focus on Higgs-strahlung at the LHC as a first example, this approach can be extended to all the important Higgs/electroweak production and decay processes, namely weak boson fusion, Higgs decay to weak bosons and diboson production. For the Higgs-strahlung process at the partonic level there are 9 angular moments, although a smaller number of these are measurable at the LHC for the final states that we are interested in. We will see that extracting all the experimentally available angular moments can simultaneously constrain all the possible hV V * /hV f f (V ≡ Z/W ± , f f ≡ ff /ff ) tensor structures. An essential prerequisite for our methods to be applicable is that the final angular distributions measured by the experiments should preserve, to a large extent, the initial theoretical form of EFT signal governed by the angular moments. To truly establish the usefulness of our methods, we therefore carry out a detailed and realistic collider study. In particular we include differentially QCD NLO effects that can potentially improve partonic contributions to the EFT signal reducing scale uncertainties. In our final results we find, despite these effects, a marked improvement in sensitivity compared to existing projections for most of the EFT couplings.
The paper is divided as follows. In Sec. 2, we write the most general Lagrangian for the pp → V ( )h(bb) at Dimension 6 in SMEFT and list the relevant operators in the Warsaw basis. Sec. 3 is dedicated in deriving the most general angular moments for the pp → V h processes in the SMEFT. In Sec. 4, we discuss the core idea of the method of moments which forms the backbone of this paper. In Sec. 5, we detail the collider studies that we undertake for the pp → V h processes. Sec. 6 is where we discuss the details of the angular analyses and obtain the bounds on the various couplings. We finally conclude in Sec. 7.

The pp → V ( )h(bb) process in the Dimension 6 SMEFT
We want to study the process pp → V ( )h(bb) where denotes + − ( + ν, −ν ) for V = Z (V = W ± ). The EFT corrections to pp → V ( )h(bb) are either due to corrections of the V f f , hbb and hV V /hZγ vertices or due to the new hV f f contact terms. In the unitary Table 1. Dimension-6 operators in the Warsaw basis that contribute to the anomalous hV V * /hVf f couplings in Eq. (2.1).Other details regarding the notation can be found in [3].
gauge all these corrections are contained in the following Lagrangian [9,58]), where for brevity we have only included the first generation for the couplings involving W ± , Z bosons, so that f = u L , d L , u R , d R , e L , e R , ν e L ; F = Q(L), the first generation quark (lepton) doublet. We assume that the above Lagrangian is extended to the other generations in a way such that the couplings δg Z,W f and g h Zf,W f are flavour diagonal and universal in the interaction basis, allowing us to impose strong constraints on them [59,60] (this is well motivated theoretically and can be obtained, for instance, by including the leading terms after imposing Minimal Flavour Violation [61]). If we limit ourselves to only universal corrections, the contact terms above must be replaced by hV µ ∂ ν V µν (note that ∂ µ hV ν V µν is equivalent to this vertex and the hV µν V µν vertices via integration by parts). The above parametrisation can be used even for non-linearly realised electroweak symmetry (see for eg. [62]) and in this case all the above couplings should be thought of as independent.
If electroweak symmetry is linearly realised, the above vertices arise in the unitary gauge from electroweak invariant operators containing the Higgs doublet. For instance, the operators of the Warsaw basis [3] in Table 1, give the following contributions to these vertices, where we have used (m W , m Z , α em , m b ) as our input parameters. In the equations for δg W,Z f above, the term, makes explicit the contribution to the shift in the input parameter, m Z , due to the above operators.
The pp → W ± ( ν)h(bb) process directly constrains the couplings δĝ h W W , κ W W and g h W Q , whereas the pp → Z(l + l − )h(bb) process constrains the couplings δĝ h ZZ , a linear combination of κ ZZ and κ Zγ , and the following linear combination of the contact terms [28], This linear combination arises by summing over the polarisations of the initial quarks as well as including the possibility of both up and down type initial-state quarks weighted by their respective PDF luminosities; the precise linear combination changes very little with energy.
For the case of linearly realised electroweak symmetry, the CP -even couplings involved in W ± h production can be correlated to those involved in Zh production using the fact that the same set of operators in Table 1 generate all the anomalous couplings as shown in Eq. (2.1). To derive these correlations we can trade the 13 CP -even Wilson coefficients above for the 13 independent (pseudo-)observables δĝ h bb , δg Z f (7 couplings), g h W Q , δĝ h W W , κ W W , κ Zγ and κ γγ , the coefficient of h 2v F µν F µν 1 . This can be done using the expressions in Eq. (2.1) and the corresponding expression for κ γγ , The rest of the anomalous couplings can then be expressed as functions of these independent ones; for example we obtain, Some of the couplings on the right-hand side of the above equations can be measured extremely precisely. For instance, the two couplings, κ Zγ and κ γγ , can be bounded very strongly (below per-mille level) by measuring the h → γγ/γZ branching ratios [6,59] 2 . In addition, the Z-coupling deviations, δg Z f , are constrained at the per-mille level by LEP data [60]. As we will see later, studying W ± h production at high energies would allow us to constrain g h W Q at the per-mille level. On the other hand, the couplings κ V V and δĝ h V V can be constrained at most at the 1-10% level. Thus, one can safely ignore the strongly-constrained couplings to obtain the direct relationships, which hold up to a very good approximation. We will utilise these relationships in order to combine our results from W ± h and Zh modes to obtain our final bounds on the CP -even vertices.
As far as the CP -odd couplings are concerned there are 4 of them including those corresponding to h 2v F µνF µν and h 2v F µνZ µν . The latter two couplings are, however, not precisely measurable as in the CP -even case. Thus an analog of the above procedure to correlateκ W W andκ ZZ is not possible.
1 This analysis is in the spirit of Ref. [9] but with a different choice of primary/independent observables.
Indeed, we include in our list the anomalous Higgs couplings, g h W Q and κZZ , rather than the anomalous triple gauge couplings (TGC) δκγ and δg Z 1 . As we will see, the bounds on the anomalous Higgs couplings are comparable or better than those expected for the TGCs. 2 This might seem surprising, as the branching ratios themselves are not constrained at this level. Recall, however, that the SM h → γγ/γZ rates are loop suppressed, so that even an O(10%) uncertainty in the branching ratios translate to per-mille level bounds on these couplings.

Beam Axis
Plane of pp-Vh Plane of V-ll In Vh CoM In ll CoM Finally we also have the correlation, which can also be translated to a correlation between the coupling g h Zp in Eq. (2.3) and those in the right hand side above.

Angular moments for the pp → V ( )h(bb) process in the Dimension-6 SMEFT
In this section we come to the central topic of this work and discuss how the full angular distributions in the pp → V ( )h(bb) processes, at a given energy, can be expressed in terms of a finite number of basis functions, both in the SM and D6 SMEFT. The corresponding coefficients of these functions are the so called angular moments for these processes. We start at the level of f f → V ( )h(bb) and then discuss the experimental subtleties that arise in the extraction of these angular moments for pp → W ± ( ν)h(bb) and pp → Z( + − )h(bb).
As we will require the two b-jets arising from the Higgs decay to form a fat jet in our analysis, we will effectively consider the three body final state of the fat jet and two leptons in this section.

Angular moments at the f f → V h level
The helicity amplitude formalism is the most convenient way to arrive at the full angular and energy dependance of the f f → V ( )h(bb) amplitude. Starting at the 2→2 level, f (σ)f (−σ) → V h, these helicity amplitudes are given by, and , λ = ±1 and σ = ±1 are, respectively, the helicities of the Z-boson and initial-state fermions, and √ŝ is the partonic centre-of-mass energy. The above expressions hold both for quark and leptonic initial states. In Eq. (3.1) above, we have kept only the terms with leading powers of √ŝ /m V both for the SM and D6 SMEFT (the subdominant terms are smaller by, at least, factors of m 2 V /ŝ). We have, however, retained the next-to-leading EFT contribution for the λ = 0 mode, as an exception, in order to keep the leading effect amongst the terms proportional to δĝ h V V . The full expressions for the helicity amplitudes including the SMEFT corrections can be found in Ref. [63]. The above expressions assume that the fermion momentum is in the positive z-direction of the lab frame. The expressions for the case where the anti-fermion has momentum in the positive z-direction can be obtained by making the replacement σ → −σ. Above, we have not included the effect of a V f f coupling deviation (δg V f in Eq. (2.1)) above which we will incorporate at the end of this section. It is worth emphasising that for both the SM and D6 SMEFT, only contributions up to the J = 1 helicity amplitude appear. For the SM this is because the f f → V h process is mediated by a spin-1 gauge boson. For the D6 SMEFT, in addition to diagrams with spin 1 exchange, there is also the contribution from the contact term in Eq. (2.1). As this contact term is exactly the vertex that would arise by integrating out a heavy spin-1 particle, even in the D6 SMEFT only contributions up to J = 1 exist. This fact will no longer be true at higher orders in the EFT expansion where higher-J amplitudes will also get contributions. At the 2 → 3 level, the amplitude can be most conveniently written in terms of ϕ and θ, the azimuthal and polar angle of the of the negatively charged lepton for V = W − , Z and the neutrino for V = W + , in the V rest frame in the coordinate system defined in Fig. 2, with τ being the lepton helicity. We have assumed a SM amplitude for the V -decay; modifications due to a V coupling deviation will be included at the end of this section. For V = W ± we always have τ = −1. We can now obtain the squared amplitude with the full angular dependence using Eq.
where we have summed over the final lepton helicity. The f i (Θ, θ, ϕ) are the 9 functions we obtained by squaring the sum of the 3 helicity amplitudes in the right-hand side of Eq. (3.3), see also [30,64,65]. Explicitly these are, where S α = sin α, C α = cos α. The subscripts of the above functions denote the Vpolarisation of the two interfering amplitudes, with T T denoting the interference of two transverse amplitudes with opposite polarisations. The corresponding coefficients a i are the so-called angular moments for this process, which completely characterise the multidimensional angular dependance of this process at a given energyŝ. The expressions for these angular moments in terms of the vertex couplings in Eq. 2.1 are given in Table 2.
Note the factor, in some of the moments, which arises from the sum over τ in Eq. (3.5). It is worth emphasising an important conceptual point here. The cross-helicity moment functions, i.e. the last six functions in Eq. (3.6), integrate to zero over the full phase space of the V -decay products. This is expected as the two amplitudes corresponding to different helicities at the level of the V -boson cannot interfere. If we look at the phase space of the decay products differentially, however, the corresponding angular moments carry very useful information. As one can verify from Table 2, for instance, the leading contribution of the κ ZZ (κ ZZ ) coupling is to to the moment a 2 LT (ã 2 LT ). As pointed out in Ref. [32], this effect can be recovered only if we study the triple differential with respect to all three angles, i.e. an integration over any of the three angles makes the basis functions f 2 LT and f 2 LT vanish. This is an example of an 'interference resurrection' study, see also Refs. [34][35][36][37], where interference terms absent at the inclusive level are 'recovered' by analysing the phase space of the decay products differentially.
It is possible that not all of these angular moments will be relevant or observable for a given initial and final state. Before considering in detail the case of the pp → V (ll)h process, our main focus, let us briefly comment on which of these angular moments are accessible to lepton colliders. For the e + e − → Z( + − )h process in lepton colliders, all nine angular moments can be measured. However, three of them, namely a 1 T T , a 1 LT andã 1 LT , are suppressed by the factor of | RL | = 0.16, which is accidentally small due to the numerical closeness of the couplings g Z l L and g Z l R . For the e + e − → W ± ( ± ν)h process, | RL | = 1 and thus this suppression is absent. There is, however, the new complication that the neutrino four-momentum is not completely accessible. As we will discuss in detail in Sec. 3.3, once conservation of energy and momentum are imposed, the neutrino four-momentum can be determined up to a two-fold ambiguity. While the other kinematical variables converge for the two solutions, the two corresponding values for ϕ, ϕ 1 and ϕ 2 , satisfy to a very good approximation ϕ 2 = π − ϕ 1 . Consequently, if one considers both solutions with equal weight, the angular moments a 1 LT , a 2 LT andã T T vanish, see Eq. (3.6).

Angular moments for the pp → Z( )h(bb) process
The first thing to note about the LHC is that the direction of the quark is not always in the same direction in the lab frame. The expressions in Table 2 are for the case where the quark moves in the positive z-direction. For the other case where the antiquark momentum is in the z-direction, as stated below Eq. (3.3), one can obtain the corresponding expressions for the angular moments by making the substitution σ → −σ. The angular moments a 1 T T , a 1 LT andã 1 LT thus vanish once we average over both these possibilities. We are thus left with the 6 moments. At high energy, a LL dominates over all other moments in the SM. The largest BSM contribution at high energies is also to a LL from the linear combination g h Zp , see Eq. (2.3), that arises from averaging over the initial state  ) are neglected, with the exception of the next-to-leading EFT contribution to a LL , which has been retained in order to keep the leading effect of the δĝ h V V term. The factor RL is defined in text and The SM part of our results can also be found in [66]. flavour and polarisation [28]. The contribution due to g h Zp grows quadratically with energy and this coupling can thus be measured very precisely as we will see in Sec. 6.2, this was also discussed in detail in Ref. [28].
Once g h Zp has been precisely measured we can use the remaining information in the angular moments to constrain the coupling δĝ h ZZ and the linear combinations, that enter, respectively, the CP -even and odd angular moments at the pp → Z( )h(bb) level. The coefficient of κ Zγ andκ Zγ above arise again by appropriately averaging Eq. (3.2) over the initial-state flavours and polarisations. Recall, however, that there is a very strong bound on κ Zγ , see Sec. 2, so that the above linear combination effectively reduces to only κ ZZ to a very good approximation. Consider now the angular moment a 2 T T and the contribution to a LL sub-dominant in γ, see Table 2, which is unconstrained even after the strong bound on g h Zp . First of all, the total rate of the pp → Z(l + l − )h(bb) process depends only on the two moments a LL and a 2 T T as all other non-vanishing moments are coefficients of cross-helicity terms that vanish upon integration over ϕ, see Eq. (3.6). The rate itself can constrain a linear combination of δĝ h ZZ and κ p ZZ . Additionally, these two moments also carry the information of the joint distribution of the events with respect to (θ, Θ), which, along with the total rate, can in  principle be used to constrain δĝ h ZZ and κ p ZZ simultaneously. We find in our final analysis, however, that the joint (θ, Θ) distribution in the events surviving our cuts is not very effective in simultaneously constraining these couplings. The main reason for this is that the Θ-distribution gets distorted with respect to the original theoretical form because of the experimental cuts necessary for our boosted Higgs analysis. In particular, we require p h T > 150 GeV, which eliminates forward events. Another effect that could further distort the distribution is radiation of hard jets. 3 As θ and Θ appear in a correlated way in the amplitude, these effects also deform the θ-distribution, but to a smaller extent. For this reason, as discussed in Sec. 4.2, we will isolate a LL and a 2 T T using only the θ-distribution in our final analysis, in order to obtain better bounds.
Much more reliable are the ϕ distributions, which preserve their original shape to a large extent. We show in Fig. 3(a), for instance, that the ϕ distributions corresponding to an enhanced a 2 LT andã 2 LT , for events that include the effect of jet radiation and pass all experimental cuts to be described in Sec. 5. We see the expected sinusoidal and cosinusoidal ϕ-dependances despite all these effects.
The information for the ϕ-dependance is carried by the angular moments a 2 LT and a T T in the CP -even case, which can be measured to constrain the linear combination κ p ZZ , assuming again that g h V f is already precisely constrained. Among these, as identified in Ref. [32], the leading contribution is from a 2 LT , as it is larger relative to a T T by a factor of γ, see Table 2. This moment provides the strongest bound on the above linear combination in our analysis but can be accessed only by looking at the joint distribution of (θ, Θ, ϕ). A standard analysis that integrates over any of these three angles would miss this effect completely.
Finally the CP -odd coupling,κ p ZZ , cannot be constrained without using ϕ information contained inã 2 LT andã T T . Again, the leading effect contained inã 2 LT is highly non-trivial and can only be accessed by utilising the triple differential distribution with respect to (θ, Θ, ϕ).
Before moving to the next subsection, we would like to comment that the distortion of the distribution due to experimental cuts and jet radiation does not invalidate our analysis. That is to say, while these effects perhaps reduce our sensitivity compared to the idealised case, as we will discuss later, these effects will already be factored into our uncertainty estimates. Moreover, our final analysis does not depend too much on the precise shape of the Θ-distribution, as we rely more on the θ and especially ϕ distributions.

Angular moments for the pp → W ( )h(bb) process
Much of the discussion in the previous section is also relevant here. Once again averaging over the initial quark-antiquark direction gets rid of the angular moments a 1 T T , a 1 LT and a 1 LT . The high energy amplitude is again dominated by a LL both in the SM and EFT. In To show the effect of the angular moments, a 1 LT andã 1 LT , we take the weight of each event to be the sign of sin(2θ) sin(2Θ). We then show the histogram with respect to ϕ and obtain the expected shapes for the two samples; (b) Regular ϕ-distributions for a Monte-Carlo sample for the W h mode with a non-zero value for the EFT coupling κ W W . We see the effect of the angular moment a T T , the only angular moment that survives after integrating over θ and Θ, and averaging over the two solutions. The events used are those passing all cuts. The angular momentã 1 LT can also be extracted in W h production but its effect can be seen only in a weighted distribution like in (a). the EFT case, the quadratically growing contribution due to g h W Q can be used to strongly constrain it. The discussion about the distortion of the Θ-distributions and its effect on extracting the moments a LL and a 2 T T also holds for this case. The main difference from pp → Z( )h(bb) arises in the ϕ-distributions. A complication arises from the fact that the neutrino four momentum is experimentally inaccessible. Imposing energy and momentum condition and assuming an on-shell W -boson yields two possible solutions for the neutrino four momentum, i.e. two solutions for the z-component of the neutrino momentum in the lab frame, the p T being equal for both solutions. While Θ, θ and the final-state invariant mass converge for the two solutions, especially at high energies [36], the values of ϕ for the two solutions do not converge, and in fact are related to each other as ϕ 2 = π − ϕ 1 to a very good approximation. In our analysis we average over Θ, θ and the final-state invariant mass, but keep both ϕ solutions with equal weight. This has the consequence that the functions cos ϕ and sin 2ϕ vanish when averaged over these two possibilities, resulting in the vanishing of the moments a 1 LT , a 2 LT andã T T , see Eq. (3.6).
In Fig. 2(a)-2(c) we show, for the three angles, a scatter plot between the truth and reconstructed values obtained after our collider analysis described in Sec. 5. For Θ and θ, we use for the reconstructed value the mean of the two solutions, whereas for ϕ, we populate the scatter plot with both solutions. It is clear from Fig. 2(c) that we have ϕ 1 +ϕ 2 = π to a very good approximation. While Fig. 2(a)-2(c) show that the angles can be reconstructed quite well, the procedure is not exact, as we have assumed that W is on-shell and did not properly take into account radiation of hard extra jets. In fact, for some rare events the virtuality of the W -boson is so high that no real solutions exist for the neutrino p z , if we assume an on-shell W -boson; we neglect such events in our analysis.
In Fig. 3(b) we show the ϕ-distribution for EFT events that finally survive the collider analysis discussed in Sec. 5. We again see the expected cos(2ϕ) shape corresponding to a T T , which is the only moment that survives integration over the other two angles and the averaging over the two solutions. The difference in the true and reconstructed distributions at ϕ = ±π/2 is related to fact that we discard events where the neutrino four momentum solutions are complex.
So far we have not considered the effect of V f f , V lland hbb coupling deviations due to D6 operators. All these coupling deviations are like δĝ h V V in that they simply rescaling the SM amplitude and thus all SM distributions. Their effect can thus be incorporated by making the replacement in Table 2 and elsewhere, (3.9) Of the above couplings, while the δg V f,l couplings are very precisely constrained to be close to zero, the effect of δĝ h bb cannot be ignored.

Basic idea
As we have seen in Sec. 3, the squared amplitude for our processes can be decomposed into a set of angular structures, f i (Θ, θ, ϕ), whose contributions are parameterised by the associated coefficients, the so called angular moments, a i . We would like to extract these coefficients in a way that best takes advantage of all the available angular information.
In principle, this can be done by a full likelihood fit, but here we use the method of moments [40,42,43]. This method has its advantages -especially if the number of events is not too large [43]. This method involves the use of an analog of Fourier analysis to extract the angular moments. Essentially, we look for weight functions, w i (Θ, θ, ϕ), that can uniquely extract the coefficients, a i , i.e., Assuming that the weight functions are linear combinations of the original basis functions, we can use Eq. (4.1) to show that the matrix λ ij = M −1 ij , where, For the set of basis functions in Eq. (3.6), the resulting matrix is given by, where we have organised the basis functions in the order in which they appear in Eq. (3.6).
It is convenient to go to a basis such that, M ij and thus its inverse λ ij , are diagonal. This can be achieved by an orthogonal rotation, by an angle, In the new fully orthogonal basis, with ξ ± = (53 ± 9 √ 29). This is the matrixλ −1 ij , so that the weight functions in the rotated basis are, We are now able to convolute our event distributions with these weight functions to extract values for the coefficients,â i in the new basis, which can then be rotated back if we are interested in the moments in the original basis.
4.2 Alternative weight functions for a LL and a 2

T T
The above algorithm to extract the moments systematically generates the set of weight functions, but this set is not unique. For instance, a function proportional to cos 2ϕ can also be the weight function for f T T . As we already mentioned, the Θ distribution suffers distortions to its original shape due to experimental cuts and other effects. For this reason the extraction of a LL and a 2 T T using the weight functions derived above does not give optimal results. To avoid this, we can use weight functions only involving θ to extract these two moments.
Consider for instance, let us integrate Eq. (3.5) over the (Θ, ϕ) to keep only the θ dependance, where a LL and a 2 T T are related to the original moments a LL and a 2 T T as follows, to make diagonal the matrix in Eq. (4.3). In this case the angle of rotation is given by tan β = 1. In this basis the weight functions for are proportional tof 1 andf 3 , given by, Convoluting the observed distribution with above with these weight functions yieldsâ 1 and a 3 which can be rotated back to give a LL andâ 2 T T and finally a LL andâ 2 T T using Eq. (4.10). Using these alternative weight functions is equivalent to using only the information in the θ-distribution to extract these two moments and ignoring the distorted Θ distribution. This will improve the final bounds we obtain in Sec. 6.2.

Extraction of angular moments and uncertainty estimate
For our simulated samples, which are generated following the procedure detailed in Sec. 5, the events are already distributed according to the squared amplitude so that the convolution in Eq. (4.1) becomes a simple summation over all the events in our sample, where in order to also take into account energy dependance we have split the events in bins of the final state invariant mass, M being the central value of a given bin. Here N is the number of Monte-Carlo events in the sample andN the actual number of events expected in the particular bin for a given integrated luminosity. Note that we have changed the normalisation of theâ i in Eq. (4.13); now iâ ifi yields the distribution of the actual number of events expected at a certain integrated luminosity and not the squared amplitude integrated over the full phase space as in Eq. (3.5). For a sufficiently large number of events N , theŵ i , converge to a multivariate Gaussian distribution with a mean and covariance matrix given byw (4.14) We find that if we keep increasing N , as soon as it is large enough (say 100), thew i and σ ij approach fixed values. In the orthonormal basis, involving the functionsf i , we find a covariance matrix that is nearly diagonal. Assuming a diagonal covariance matrix, the angular moments in the orthonormal basis converge to Gaussians with mean and standard deviation given by, see Eq. (4.13), As a cross-check we also computed the second term above by splitting our Monte-Carlo sample into parts withN events each and computingâ i in each case; the standard deviation of theâ i thus obtained matches very closely the second term above. This way of estimating the error also shows that any deformation of the original angular distribution due to experimental or QCD effects (see Sec. 3.2), has been already factored into our uncertainty estimate. To estimate the uncertainty on theâ i one must also consider that fact that,N , the expected number of events in the given bin, itself fluctuates statistically. Moreover, there are systematic uncertainties in the value ofâ i we obtain in this way. These two effects result in additional uncertainties on the mean value above. Adding all these errors in quadrature we obtain for the uncertainty in each moment, a i , where κ syst represents the percentage systematic error that we take to be 0.05 in this work.

Collider Simulation
In this study, we take into account NLO QCD effects. We work under the MG5 aMC@NLO [68] environment to generate NLO events showerd with Pythia8 [69,70]. Inside this framework, real emission corrections are perfomed following the FKS substraction method [71], whereas virtual corrections are done using the OPP reduction technique [72]. The MC@NLO formalism [73] takes care of the matching between the NLO matrix element and parton shower, thus avoiding double counting. Decay of heavy bosons has been carried out with the help of MadSpin [74], which retains spin information at tree-level accuracy. We construct our NLO model using FeynRules [75] and then employ NLOCT [76] to compute the U V and R 2 counterterms, which are required for the one-loop calculation. U V counterterms are essential to remove ultraviolet divergences that appear at the loop level, whereas R 2 terms originate from the one-loop integrands that carry (n − 4)-dimensional pieces in the numerators and n-dimensional terms in the denominators. As and when required, we manually insert the R 2 terms in the NLO model as the usage of publicly available NLOCT version is restricted to renormalisable interactions only.
In this work, we focus on three different processes, i.e. pp → Zh and pp → W ± h, with the Higgs decaying to a pair of b-quarks and the Z/W decaying leptonically. Thus, for the Zh (W h) process, we study the + − bb ( νbb) final states, where = e, µ, τ . The qq → Zh and qq → W ± h processes are generated at NLO QCD, whereas the gg → Zh channel is generated at LO (which is at one loop). The following analyses are performed at 14 TeV centre-of-mass energy and the predictions are shown for the HL-LHC for an integrated luminosity of 3 ab −1 .
First we outline the generations of the signal and background samples for the pp → Zh → bb + − analysis. While generating the signal samples, i.e. qq → Zh, we use the aforementioned NLO model file and interface it with Pythia8. We choose dynamic renormalisation and factorisation scales, µ F = µ R = m Zh . We choose NNPDF2.3@NLO as our parton distribution function (PDF) for the NLO signal samples. As mentioned above, for the NLO signal samples we use MadSpin [74] to decay the heavy bosons. This step is done at LO and hence we correct for the branching ratios following the Higgs working group recommendations. We follow Refs. [28,32] while generating the background samples. All background samples are generated at LO with NNPDF2.3@LO as the PDF. The dominant backgrounds comprises the Zbb and the irreducible SM Zh production. For the Zbb production, we consider the tree-level mode as well as the gg → ZZ mode at one-loop. Furthermore, we consider reducible backgrounds like Z+ jets and the light jets are misindentified as b-tagged jets (c-jet misidentification is not considered separately), and the fully leptonic decay of tt. Rather than performing a standard resolved analysis, where one would consider two separate narrow b-tagged jets, here we require a fat jet with its jet parameter R = 1.2. We utilise a modified version of the BDRS algorithm [77] in order to maximise sensitivity. This procedure helps us in maximising the signal by retaining extra radiations and in discriminating electroweak-scale resonant signals from strong QCD backgrounds, see also [78,79].
To briefly review the BDRS approach, the jets are recombined upon using the Cambridge-Aachen (CA) algorithm [80,81] with a considerably large cone radius in order to contain the maximum number of decay products ensuing from a resonance. The jet clustering process is then read through backwards and one stops when the mass of a subjet, m j 1 < µm j with µ = 0.66, where m j is the mass of the fatjet. This step is called the mass drop and is required to occur without a significant asymmetric splitting, where y cut = 0.09. When this condition is not satisfied, the softer subjet, j 2 , is removed from the list and the subjets of j 1 are subjected to the aforementioned criteria. This procedure is repeated iteratively until the aforementioned condition is met. This algorithm terminates when one obtains two subjets, j 1,2 which abide by the mass drop condition. However, the mass drop algorithm does not improve the resonance reconstruction significantly and more fine-tuning is necessary to segregate the signal from the background. A further step is performed: filtering. In this algorithm, the constituents of the subjets j 1 and j 2 are further recombined using the CA algorithm but with a cone radius R filt = min(0.3, R bb /2). This algorithm chooses only the hardest three filtered subjets in order to reconstruct the resonance. In the original paper [77], the resonance in question is the SM-like Higgs boson and thus the hardest two filtered subjets are required to be b-tagged. In the present work, we find that the filtered cone radius R filt = max(0.2, R bb /2) performs better in reducing the backgrounds. As shown in Ref. [77], the filtering step significantly reduces the active area of the initial fatjet. Finally, we require the hardest two filtered subjets to be b-tagged with tagging efficiencies of 70%. Moreover, the misidentification rate of light subjets faking as b-subjets is taken as 2%.
One of our goals is to look for new physics effects in high-energy bins and hence it is imperative to generate the signal and background samples with certain generation-level cuts in order to improve statistics. For the qq → Zh samples generated at NLO, we require a cut on the p T of the Higgs boson, p T,h > 150 GeV. The Zbb and tt samples are generated with the following cuts: p T,(j,b) > 15 GeV, p T, > 5 GeV, |y j | < 4, |y b/ | < 3, ∆R bb/bj/b > 0.2, ∆R > 0.15, 70 GeV < m < 110 GeV, 75 GeV < m bb < 155 GeV and p T, > 150 GeV. The Zbb sample is generated upon merging with an additional matrix element (ME) parton upon using the MLM merging scheme [82]. For the Z+ jets samples, we do not impose any invariant mass cuts in the jets. Furthermore, the sample is merged with three additional partons. Since the backgrounds are generated at LO, we use flat K-factors to bring them at a similar footing to the signal. For the tree-level Zbb, one loop gg → ZZ, one loop gg → Zh and Z+ jets, we respectively use K-factor values of 1.4 (computed within MG5 aMC@NLO), 1.8 [83], 2 [84] and 1.13, computed within MCFM [85][86][87].
A cut-based analysis has been done in Ref. [28] and it has been shown that the prowess of a multivariate analysis exceeds that of a simple cut-and-count analysis. Thus, in this work we do not revisit the cut-and-count analysis and delve directly into the multivariate formulation. We start by constructing fatjets with cone radii of R = 1.2. Furthermore, we require these fatjets to have p T > 80 GeV and to lie within a rapidity, |y| < 2.5. We employ FastJet [88] in constructing the jets. Moreover, we isolate the leptons (e, µ) upon demanding that the total hadronic activity deposited around a cone radius of R = 0.3 can at most be 10% of its transverse momentum. The leptons are also required to have p T > 20 GeV and have rapidity, |y| < 2.5. In our setup, every non-isolated object is considered to be part of the fatjet construction. Before performing the multivariate analysis, we require each event to have exactly two oppositely charged same flavour (OSSF) isolated leptons. Moreover, we apply loose cuts on certain kinematic variables. We require the invariant mass of the leptons to be in the range 70 GeV < m < 110 GeV, the transverse momentum of the di-lepton system, p T, > 160 GeV. We also require ∆R > 0.2 4 , p T,fatjet > 60 GeV, the reconstructed Higgs mass, 95 GeV < m h < 155 GeV, ∆R b i , j > 0.4 (i = 1, 2) and / E T < 30 GeV. We also require that there is at least one fat jet with at least two B-meson tracks, there are exactly two mass-drop subjets and at least three filtered subjets. We also require that the hardest two filtered subjets are b-tagged. Owing to the smallness of the Z+ jets and tt backgrounds compared to Zbb, we train our boosted decision tree (BDT) upon only considering the NLO Zh and the tree-level Zbb samples. We use the following variables to train the BDT: p T of both isolated leptons, ∆R between the b-subjets and the isolated leptons (four combinations), between the isolated leptons and also between the two b-subjets in the fatjet, the reconstructed dilepton mass and its p T , the ∆φ separation between the fatjet and the reconstructed dilepton system, the missing transverse energy, / E T , the mass of the Higgs fatjet and its transverse momentum, p T of the two b-tagged filtered subjets, the ratio of the p T of these b-tagged subjets and finally the rapidity of the reconstructed Higgs fatjet. During our training process, we do not require variables that are 100% correlated but retain every other variable. Given that one of our final variables of interest is the reconstructed Zh invariant mass, we refrain from using it as an input variable. For the BDT analysis, we use the TMVA [89] pacakge in the root framework. During the analysis, we use 50% of the samples for training and always ensure that there is no overtraining by requiring that the Kolmogorov-Smirnov statistic is at least O(0.1) [90]. After optimising the cut on the BDT variable, one finds that there are around 463 qq → Zh (SM) and 820 Zbb events at 3 ab −1 , which amounts to the SM qq → Zh (SM ) over rest of the background (B) ratio, SM/B ∼ 0.56. Using the same training, we have respectively 44, 7 and 57 Z+ jets, gg → ZZ and gg → Zbb backgrounds after the BDT cut. This yields SM/B ∼ 0.5.
For the W ± h → bb ν analysis, we follow a very similar framework as before. The dominant backgrounds are the irreducible SM W ± h and the reducible W ± bb channels. We also consider the fully and semi-leptonic tt events, W ± + jets and Z+ jets, where Z → + − . The W ± samples are generated at NLO QCD using the aforementioned method. The W ± bb samples are generated upon merging with an additional parton as described above. Unlike the Zh channel, the W ± h channel only has quark-initiated production mode. For the Zh channel, it was quite simple to reduce the tt background by imposing a lower cut on / E T . For the W ± study, the signal itself contains a final state with a neutrino and hence demanding a cut on / E T will not only reduce the tt backgrounds but also a significant fraction of the signal. The signal samples are generated with p T,h > 150 GeV and the invariant mass of the W h system, m W h > 500 GeV (we clarify this choice later). We use the same PDF choice as for the Zh samples and the scales are chosen to be µ F = µ R = m W h . The backgrounds are generated with the same PDF choice at LO. The scales chosen for the background generation are m W for the W bb and W + jets samples and 2m t for the tt samples. Moreover, weak cuts are imposed on the background samples at the generation level. These include, p T,(j,b) > 15 GeV, p T, > 5 GeV, |y b/ | < 3, |y j | < 5, ∆R bb > 0.1, ∆R b > 0.2 and 70 GeV m bb < 155 GeV. We separate the W h analysis into two parts depending on the charge of the isolated lepton. For the analysis, we require one isolated charged lepton. In contrast to the Zh analysis, the W ± h has a known ambiguity in the form of the p z component of the neutrino momentum. We deal with this by requiring that the invariant mass of the neutrino and the isolated lepton peaks around the W -boson mass. This gives us two solutions to p z,ν and we demand that the solutions are always real. We discard events where complex solutions are encountered. We construct two invariant masses for the W h system for the two neutrino p z solutions, m fatjet ν 1,2 . Before implementing the BDT analysis, we employ certain loose cuts like p T,fatjet > 150 GeV, 95 GeV < m h < 155 GeV, m fatjet ν 1,2 > 500 GeV and ∆R b i , > 0.4. On top of this we require certain number of fatjets, mass-drop and filtered subjets as discussed for the Zh scenario. For the BDT analyses (one for W + h and another for W − h), we train the samples upon considering the SM W h sample as the signal and the W bb, semi-leptonic and fully leptonic tt samples as backgrounds. Owing to multiple backgrounds, we impose relative weight factors to these backgrounds which are defined as 1/L gen , where L gen is the generated luminosity that depends on the production cross-section, including the K-factors, and the number of Monte Carlo generated events. Besides, NLO samples also contain negative weights for certain events, which we include while training the BDT samples. We also find that the effect of including the weight factor in our training is small, owing to the very small number of signal events having negative weights (less than 4% percent). We optimise the BDT analysis for W + h (W − h) and find 1326 (901) events for the signal and 4473 (3476) W + bb (W − bb) events at 3 ab −1 . The number of surviving events for tt, W + jets and Z+ jets are much smaller. Ultimately, we find SM/B ∼ 0.28 (0.24) for W + h (W − h).

Analysis and Results
In this section we describe how we obtain our final sensitivity estimates and present our main results. We will consider only the interference contribution in this study which in any case is expected to be dominant piece below the EFT cut-off. There is no conceptual hurdle in including also the squared terms, as Eq. (3.5) is still equally valid, and the reasons for neglecting them are only practical. We first consider the contact terms g h V f which can be very precisely constrained in the high energy bins. Once these couplings are very precisely constrained we will turn to the lower energy bins where there are a sufficient number of events to carry out an angular moment analysis to constrain the other couplings. All the results we will present in this section will be for an integrated luminosity of 3 ab −1 .

Bounds on contact terms
As already discussed, at high energies the EFT deviations are dominated by the contribution of the contact interactions, g h V f , to a LL . Because this contribution grows quadratically with energy relative to the SM V h contribution, it can be very precisely constrained by probing high energy bins. Unfortunately some of the bins providing maximum sensitivity have too few events for an angular moment analysis. We thus constrain these couplings simply using the final state invariant mass distribution. Following Ref. [28], where this procedure was carried out for the Zh mode, we construct a bin-by-bin χ 2 function assuming the expected number of events is given by the SM and the observed by the SMEFT. To ensure that we do not violate EFT validity we neglect any event with a final state invariant mass above the cut-off, which is evaluated for a given value of the anomalous couplings, by setting the Wilson coefficients in Eq. (2.1) to unity. For an integrated luminosity of 3 ab −1 , we obtain the sub-per-mille level bounds at the one sigma level,

Angular Moment analysis
Now that g h W Q and g h Zp are strongly constrained from the higher energy bins, we turn to the lower energy bins with enough events to perform an angular moment analysis to constrain the other couplings. Ideally we should marginalise over the effect of contact terms also in the lower bins, but as we will see the expected bounds on the contact terms are almost two orders of magnitude smaller than that of the other couplings, and thus their effect is negligible in the lower energy bins. Therefore we will ignore them in further analysis. We first split our simulated events into 200 GeV bins of the final state invariant mass. To obtain the angular moments we first convolute the events in each energy bin with the weight functions using Eq. (4.13). As the CP -even and odd couplings contribute to a mutually exclusive set of angular moments we construct two separate bin-by-bin χ 2 functions as follows, where κ p V V ,κ p V V are same as κ W W ,κ W W for V = W and defined in Eq. (3.8) for V = Z. In the above equation, we include only the CP -even (CP -odd) angular moments in χ 2 (χ 2 ), the index i indicates the different moments and M j labels the invariant mass bins.
The squared error in the denominator is computed using Eq. (4.16) on the background sample (which includes the SM V h contribution) whereN in this case is the total number of background events in the j-th bin.
Once again the contributions due to κ p V V andκ p V V grow with energy and one must be careful about EFT validity. For a given value of the coupling we estimate the cut-off Λ using Eq. (2.1) putting the all the Wilson coefficients to unity. We ignore any event that has final state invariant mass above 1500 GeV, a value smaller than the cut-off corresponding to the size the couplings we will eventually constrain. The most sensitive bins for the analysis of the contact term, on the other hand, are bins higher than 1500 GeV. The contribution due toĝ h V V does not grow with energy with respect to the SM and thus the bounds on this coupling are in any case dominated by the contribution from the lowest energy bins in our analysis.
We now discuss the results for the Zh and W ± h modes separately before presenting our combined bounds. The individual bounds are important as they do not assume Eq. (2.6) whih has bee derived assuming that electroweak symmetry is linearly realised. In fact, the independent measurement of couplings involving the Z and W can be used to verify Eq. (2.6) as a prediction of linearly realised electroweak symmetry.

Zh mode
Following the discussion in Sec. 3.2 we include, in Eq. (6.2), the momentsâ 1 ,â 3 , a 1 LT and a T T in χ 2 . Recall thatâ 1 andâ 3 are linear combinations of the original angular moments a LL and a 2 T T defined in Sec. 4.2. The bound obtained for the two CP -even couplings is shown in Fig. 4(a). To show the power of our method we show the progression of the bounds obtained as the differential information used is gradually increased. The bound obtained, if one uses only total rate to constrain a linear combination of the two couplings, δĝ h ZZ and κ p ZZ is shown by the two dashed lines. Next we include distributions of the final state invariant mass and other differential information at the level of Z-boson four momentum, i.e. the decay products of the Z-boson are treated inclusively, and obtain the excluded region shown in purple; for this we include only the angular momentsâ 1 andâ 3 , extracted using the weights in Sec. 4.1, thus using information of the Θ-distribution. The analysis at this stage is comparable to a regular SMEFT analysis that includes a few standard differential distributions. Finally we include the effects of the the angular moments a 1 LT and a T T and obtain our final bound shown in red. The main improvement in sensitivity in the final bounds comes from a 1 LT the effect of which can be captured only by a careful study of the joint (Θ, θ, ϕ) distribution as pointed out in Ref. [32]. While this is clearly something beyond the scope of a regular cut-based analysis, as one would need to take into account all the correlations of the final state phase space, the angular moment approach captures it effortlessly.
We show also the projected bounds from the h → ZZ → 4l process in Fig. 4(a). The blue band shows the bound from the h → ZZ → 4l rate whereas the green bar is the bound on obtained using the Matrix Element Likelihood Analysis (MELA) framework [91]. As far as κ p ZZ is concerned, we see that the bound obtained from Zh production using our  Figure 4. (a) Bounds at 65% CL on the CP -even anomalous couplings from Zh production with 3 ab −1 integrated luminosity, assuming that the contact term has been very precisely constrained (see Eq. (6.1)). We show the improvement of the bounds as more and more differential information is included in the fit. The dashed lines show the bound just from the total rate. The purple region includes differential information at the level of the Z-boson four momentum such as the final state invariant mass distribution and Θ-distribution. Finally the red region includes information from all the angular moments including the cross-helicity interference terms. The blue band shows the bound from h → ZZ → 4l rate. The bars show the bounds on one of the couplings when the other coupling is 0. The green bar shows the bound obtained using the Matrix Element Likelihood Analysis (MELA) in Ref. [91] and assuming δĝ h ZZ = 0. (b) Same as in (a) but for the W ± h mode where there is no bound from MELA. methods surpass the other existing projections shown in Fig. 4(a) 5 . In the horizontal direction our bounds might seem redundant once the h → ZZ → 4l process is taken into account, but if one allows for hbb coupling deviations our bounds become the measurement of a truly independent effect, see Eq. (3.9).
The CP odd coupling,κ p ZZ is constrained using the functionχ 2 in Eq. (6.2) which includes the momentsã 1 LT andã T T . We finally obtain the one sigma level bound,

W h mode
As discussed in Sec. 3.3 the relevant angular moments for this case are a LL , a 2 T T and a T T in the CP -even case. Instead of the first two moments we use the linear combinationâ 1 and a 3 described in Sec. 4.2. Again we show the progression of the bounds at different stages of 5 A bound using the matrix element method for pp → Zh may potentially match our bounds but the results in Ref. [91] are unfortunately not comparable to ours as these studies include high energy phase space regions where the EFT contribution is many times that of the SM. The methodology iused to obtain these bounds, thus, violate our assumption of O(1) Wilson coefficients.
inclusion of differential information. The dashed lines show bounds from the total rate and the purple region shows the bound obtained by including only the angular moments, a LL and a 2 T T that encapsulate the differential information at the level of the Z-boson treating its decay products inclusively. Our final bound, also including the effect of a 1 LT and a T T , is shown in red. We show also the projected bounds from the h → W W → 2l2ν decay rate in blue to which our bounds are complementary (recall again that, what aour bounds actually probe is a linear combination also involving hbb coupling deviations, see Eq. (3.9)). In this case there is no competing bound on κ W W from the h → W W mode presumably because the neutrinos in the final state make much of the differential information inaccessible in this case. Thus our bounds on κ W W from the pp → W ± h process is likely to be the best bound on this coupling possible.
Again the CP odd coupling,κ W W is constrained by including the momentã LT 1 in the functionχ 2 in Eq. (6.2). We finally obtain the one sigma level bound, |κ W W | < 0.04. (6.4) We see that we obtain bounds of similar size from the pp → W h and pp → Zh processes on the respective anomalous couplings. The fact that the couplings and can be independently measured is very important as we can then use these measurements to test the correlations in Eq. (2.6) which in turn tests whether electroweak symmetry is linearly realised or not. An alternative approach would be to use the correlation in to combine the bounds from W h and Zh production as we show in the next subsection.

Combination of Zh and W h modes
In Fig. 5 we show the bounds obtained after combining the results of using Eq. (2.6), thus assuming electroweak symmetry is linearly realised. Again, we show the bound obtained at various levels of inclusion of differential data. The dashed lines show the bound just from the total rate, the purple region includes differential information at the level of the Z/W -boson four momentum and the red region is our final bound including all angular moments. The blue band shows the bound from a combination of h → W W → 2l2ν and h → ZZ → 4l rate. The green bar shows the MELA bound from Ref. [91] on κ ZZ assuming δĝ h ZZ = 0, translated to this plane.

Conclusions
The precise measurement of Higgs boson properties will be one of the legacies of the LHC's scientific achievements. Potential deformations of the Higgs boson's couplings to other particles compared to Standard Model predictions can be cast into limits on Wilson coefficients of effective operators originating in the SMEFT framework. To obtain predictive limits on the highly complex system of SMEFT operators, it is necessary to measure Higgs interactions in various ways production and decay channels. One of the most important ones to establish the nature of the Higgs boson and its embedding into the scalar sector are its couplings to massive gauge bosons, i.e. the W and Z bosons. . Bounds at 65% CL on the CP -even anomalous couplings, with 3 ab −1 integrated luminosity, after combining results from Zh and W h production using Eq. (2.6) and assuming that the contact terms have been very precisely constrained (see Eq. (6.1)). Again, we show the progression of the bounds as more and more differential information is included in the fit. The dashed lines show the bound just from the total rate in both processes. The purple region includes differential information at the level of the Z/W -boson four momentum. The red region is our final bound and includes information from all the angular moments. The blue band shows the bound from a combination of h → W W → 2l2ν and h → ZZ → 4l rate. The bars show the bounds on one of the couplings when the other coupling is 0. The green bar shows the bound implied by the bound on κ ZZ using the Matrix Element Likelihood Analysis (MELA) in Ref. [91] and assuming δĝ h ZZ = 0.
We proposed a novel method to probe the full structure of the Higgs-gauge boson interactions in Higgs-associated production. Using the helicity amplitude formalism and expanding the squared matrix elements into angular moments the whole process can be expressed in terms of nine trigonometric functions. This is true not only in the SM but also in the D6 SMEFT. Extracting the coefficients of these functions, the so called angular moments, is a powerful and predictive way of encapsulating the full differential information of this process. As differential information can encode signatures of EFT operators in subtle ways, maximally mining the differential information is essential to obtain the best possible bounds on the EFT operators. As the actual interpretation of the measurement relies now on a shape analysis of a small number of trigonometric functions, strong constraints can be obtained, provided experiments are going to publicise measurements of these functions. Thus, we encourage the experimental collaborations to provide such measurements for various Higgs production processes 6 .
The efficacy of this method relies crucially on whether the theoretical form of the original angular distribution can be preserved despite effects like experimental cuts, showering and hadronisation. In this article, we carried out a detailed collider simulation of the Higgs-strahlung process, including these effects, before applying the method of angular moments. The results we find are encouraging, indicating that a shape analysis using the trigonometric basis functions can set the most sensitive limits on effective operators within the SMEFT framework. While the high energy behaviour of the process results in the strongest possible bounds on the hV f f contact terms (see Eq. (6.1)), the full angular moment analysis leads to the strongest reported bounds on the hV µν V µν (see Figs. 4(a), 4(b) and 5) and hV µνṼ µν (see Eq. (6.3) and Eq. (6.4)).
We aim to extend this method to various other Higgs/electroweak production and decay processes such as weak boson fusion [93], the h → ZZ → 4l decay [94] and diboson production [95]. One can then perform a full global fit including this enlarged set of observables to obtain the best possible bounds on the SMEFT lagrangian.