A fully differential SMEFT analysis of the golden channel using the method of moments

The Method of Moments is a powerful framework to disentangle the relative contributions of amplitudes of a specific process to its various phase space regions. We apply this method to carry out a fully differential analysis of the Higgs decay channel $h \to 4\ell$ and constrain gauge-Higgs coupling modifications parametrised by dimension-six effective operators. We find that this analysis approach provides very good constraints and minimises degeneracies in the parameter space of the effective theory. By combining the decay $h \to 4\ell$ with Higgs-associated production processes, $Wh$ and $Zh$, we obtain the strongest reported bounds on anomalous gauge-Higgs couplings.


Introduction
One of the main goals of the Large Hadron Collider (LHC) is to understand the precise nature of electroweak symmetry breaking. The most direct way to probe this is via a measurements of the couplings of the Higgs boson to the weak bosons. While the LHC measurements have already established that the Higgs boson couples to the gauge bosons [1][2][3][4][5], a full resolution of the tensor structure of these couplings is not possible without a thorough differential study of the relevant processes, namely-Higgs decays to gauge bosons, Higgsstrahlung and Higgs production in vector boson fusion. New tensor structures for these couplings are unavoidable once we include operators from the next order in the Standard Model Effective Field Theory (SMEFT) expansion (see Ref.  for other relevant SMEFT studies). As we await the arrival of large volumes of data from the high luminosity runs of the LHC, the development of differential strategies to probe the gauge-Higgs couplings is thus of great relevance.
In this work we perform a fully differential analysis of the golden Higgs decay channel, h → 4 [55][56][57][58][59] using the 'method of moments' [60,61]. This technique would utilise the In the SMEFT the h → 4 decay can, in principle, arise also from topologies with anomalously large h¯ couplings but such contributions will still not be comparable to the h → Z ( * ) Z * → 4 contribution if we impose current bounds on the h¯ couplings [67,68]. Another possibility is that the leptons arise from intermediate photons, however production of intermediate photons is loop suppressed and would require enhancement by at least an order of magnitude to have any impact; this would be easily ruled out by the bounds on the branching ratios for h → γγ [69,70] and Zγ [71,72] (see Ref. [73] for HL-LHC projections) and will not be considered further. In addition the contact interaction hZ¯ gives a new diagram not present in the SM. Finally, the SM diagram for this process gets EFT corrections at various vertices, i.e., ggh, hZZ, Z¯ . All these corrections are summarised by the following Lagrangian written in the broken phase [19,74] 2 , where, for brevity, we have just included the first generation leptons for the couplings with the Z-bosons, such that = e L , e R , ν e L . The Lagrangian is assumed to be extended to the remaining two generations, such that the couplings δg Z and g h Z are flavour diagonal and universal in the interaction basis. This allows us to impose strong constraints on these couplings [17,75]. This assumption is theoretically well-motivated and can be obtained by including the leading terms after imposing Minimal Flavour Violation (MFV) [76]. In the above Lagrangian we have omitted any EFT corrections related to the production of the Higgs boson as all these corrections cannot be parametrised by local Lagrangian terms. We discuss these at the end of the section.
The parameterisation in the above Lagrangian also holds for a non-linearly realised electroweak symmetry [77] and in this scenario, all the above couplings must be considered as independent. On the other hand, if electroweak symmetry is linearly realised, the aforementioned vertices, in the unitary gauge, arise from operators containing the Higgs doublet. The list of operators in the Warsaw basis [8] that contribute to this process, including those that affect the Higgs boson production are shown in Table 1; these contribute Table 1. List of dimension-six operators in the Warsaw basis which contribute to the anomalous hV V * /hVf f , the effective Higgs-gluon, Yukawa and chromomagnetic couplings in Eq. (2.1). Details about the notations can be found in Ref. [8].
to the said vertices as follows, where, (m W , m Z , α em ) are our input parameters. In the equation for δg Z above, the term explicitly shows the contribution of two of the aforementioned operators to the shift in m Z ; one of the input parameters. Of all the anomalous couplings in Eq. 2.1, only δĝ h ZZ , κ ZZ andκ ZZ would be eventually relevant for us. This is because the other couplings can be measured or constrained much more stringently in other processes. First, note that LEP1 [78] has put per-mille level constraints on δg Z from the partial Z-decay measurements of Γ(Z → ¯ ) parameters [17,75]; the corrections due to these couplings would thus be neglected.
As far as the corrections due to the couplings, g h Z , are concerned they can be ignored as these couplings can be very stringently constrained at HL-LHC using the D6 SMEFT correlations. This is because these couplings receive contributions from the same operators as δg Z , apart from c W B and c HD that only contribute to δg Z . The Wilson coefficients c W B and c HD can be constrained by their contribution to the anomalous charged Triple Gauge Couplings (TGCs) [79], (2.4) Using the expressions for g h Z , δg Z in Eq. 2.2 and the above expressions for the TGCs we obtain the following relationship, derived also in [19,74]. Whereas, the δg Z couplings are very stringently constrained by Z-pole measurements as discussed above, per-mille level bounds are also expected for the TGCs at the HL-LHC [80,81]. The contact term couplings, g h Z can, therefore, be tightly constrained at the HL-LHC using the correlation in Eq. 2.5. While we will consider the effect of the g h Z couplings in our theoretical discussions, these couplings will be eventually neglected in our final numerical analysis that would lead to bounds only at the percent level 3 .
We have still not considered corrections to Higgs production in the pp → h → Z ( * ) Z * → 4 process. These involve five other operators, namely the Yukawa operators, O y b ,yt , the chromomagnetic operator, O Htg,Hbg , and the Higgs gluon effective operator, O HG as discussed in Refs. [83,84]. In this work, we perform our analysis in the Higgs rest-frame and hence the effects of these operators only appear as an overall factor that affects the total rate but can be decoupled as far as the lepton distributions are concerned. These effects thus effectively redefine δĝ h ZZ , which also affects only the rate and has no differential signature, where f (c HG , c yb , c yt , c HtG , c HbG ), a linear combination of the aforementioned Wilson coefficients, can be obtained from the results of Refs. [83,84]. In order to constrain and disentangle these EFT corrections in the production sector, it is necessary to study the production of h+ jet and tth, in conjunction and include other Higgs decay channels. 3 The four-point contact interactions, g h Z can also be directly constrained at future e + e − colliders, running at TeV scale energies. It was shown in Ref. [82], that the analogous contact term couplings involving quarks can be constrained at the per-mille level or stronger at the HL-LHC; one can thus expect bounds at a similar level for g h Z as at lepton colliders.
As far as EFT contributions to the non-Higgs background are concerned, the main corrections to the dominant background from the qq → 4 process come from the TGCs and the Z-coupling deviations, δg Z , both of which would be strongly constrained in other processes as discussed above. As the gg → 4 background is much smaller (about a percent, see Sec. 5) we will not consider EFT modifications, for instance via anomalousttg andttZ vertices. Finally, one should also notice that a deviation from SM in the gg → h → 4 amplitude modifies its interference with the gg → ZZ → 4 continuum as well. However, we will neglect in our analysis this contribution to the total cross-section since at invariant masses around m h this effect is negligible with respect to the Higgs production channel [85].
3 Angular dependence of the h → 4l amplitude As we discussed in the previous section, the predominant contribution to the pp → h → 4 process is from gg → h → Z ( * ) Z * → 4 in the SM. In SMEFT at dimension-six level there can be two different kinds of processes that dominantly contribute to h → 4 . The first one is mediated, as in SM, by two intermediate Z bosons, taking into account the δĝ h ZZ , κ ZZ andκ ZZ (see Eq. 2.5) modifications of the hZZ vertex. The second one corresponds to an amplitude containing an effective hZ¯ contact interaction (induced by the O HL, SMEFT operators, see Sec. 2), followed by the decay of the single Z that is produced. In both cases, the 4 final state is made of two fermion currents, at least one of which is emitted via vector decay. Each current contains an outgoing fermion and an outgoing anti-fermion having opposite helicities, in the massless limit; thus, in the following, each one of these 2-fermion states will be denote as an + − system.
The goal here is the analysis of the differential distribution with respect to the angular variables described in the following. The three angles required to define the final state fermions are shown in Fig. 1, where the angle definitions assume that the leptons have a fixed charge. Theoretically, however, it is much more convenient to express the angular distribution for final-state fermions with definite helicity. For the purposes of the subsections, Sec. 3.1 and Sec. 3.2, the reader should assume that these angles refer to the positive helicity lepton and not the negatively charged lepton. We will provide the translation to the experimentally accessible distribution with respect to final state leptons with a fixed charge later in Sec. 3.3. The angles θ 1 and θ 2 are the polar angles that the momentum of a lepton + , with chosen positive helicity, forms with the direction of motion of the center of mass of its parent + − system. This is evaluated in the + − center of mass frame where the fermion and the antifermion are back to back. For example, when + and − come from a Z decay, θ i is measured in the Z rest frame, where θ i = 0 corresponds to the direction of motion of the parent Z in the Higgs rest frame, where i = 1, 2 refers to the two + − system. Furthermore, we consider, for each + − system, the azimuthal angle ϕ i that describes the orientation of the plane individuated by the two fermion momenta, evaluated in the Higgs (or equivalently 4 ) center of mass frame; we, then, define φ as the relative azimuthal angle between the two planes. For more details on definition of these angles see Sec. 5. 2 ). Note that two different reference frames are used to compute i) the azimuthal angle φ between the planes formed from the lepton pairs in the Higgs rest frame, and ii) the polar angle θ i (i = 1, 2) of the lepton in the rest frame of its parent Z boson. Each frame is characterised by the presence of back-to-back objects. We adopt the convention used in ref. [86].
As we have commented above, the structure of the interactions is such that each + − system has two opposite helicity fermions with λ = ±1/2. Then, the angular dependence of the amplitude is determined uniquely by the angular momentum quantum numbers of both the 2-fermion states, in the + − center of mass frame, namely by the total angular momentum J of each 2-fermion state and by the projection M along the + − direction of motion. More specifically, the total amplitude is proportional to the Wigner functions d J M,∆λ=1 (θ i , ϕ i ) for both the + − final states, where the angles are evaluated considering the momentum of the positive helicity fermions and ∆λ = 1 is the helicity difference between + and − . When + − come from a decay of a single particle, J and M are given respectively by the total spin and helicity of the intermediate unstable state. In our case, the decaying particle is a Z boson, implying J = 1 and M = 0, ±1 according to the helicity of the vector. Therefore, the angular dependence will be described by the following Wigner functions: As a consequence, the study of h → ZZ * and h → Z helicity amplitudes is crucial and it will be done in the next subsections. Even in the case of amplitude with insertion of a hZ¯ contact term, the angular modulation turns out to be completely described by d 1 M,∆λ=1 (θ, ϕ) functions, as discussed in details in Sec. 3.2.

Helicity amplitudes for ZZ * production
We consider a scalar particle which decays into 2 states; in the scalar rest frame, the system has zero total angular momentum. Therefore, in this frame the two final objects, with opposite momenta, must have same helicities, in order to guarantee M = 0 for the full final state. It is the case of h → ZZ * in the Higgs rest frame: the helicity configurations allowed for the final states are Z + Z + , Z − Z − and Z 0 Z 0 . We call the corresponding helicity amplitudes as A ++ , A −− and A 00 . The possible deviations from SM in the hZZ interactions are induced by the two CPeven effective couplings δĝ h ZZ and κ ZZ , that contribute to all the three helicity amplitudes, and the CP-odd couplingκ ZZ that contributes only to the transverse amplitudes. The three h → ZZ * amplitudes are: where q = q Z = − q Z * are the 3-momenta in the Higgs rest frame, m 2 Z = q 2 Z and m 2 Z * = q 2 Z * are the squared invariant masses of the on-shell and off-shell Z respectively. Above, the term independent from effective couplings δĝ h ZZ , κ ZZ andκ ZZ describes the SM amplitude. We have analysed above the h → ZZ * helicity amplitudes; however, as previously discussed, the helicities of the vector bosons determine the dependence of the amplitude on the final lepton angles. In fact, a h → ZZ * → + − + − helicity amplitude, which is obtained summing over the polarisations of the intermediate vector bosons, is given by: − ) are the amplitudes with the Z-lepton couplings factorised in front.
The parameter g Z,Z * is the SM coupling of the Z boson which depends on the chirality of the fermions, = L,R . For a fixed helicity the chirality is different for fermion and anti fermion fields: positive helicity corresponds to LH chirality in the case of positively charged leptons and to RH chirality in the case of negatively charged leptons. A more detailed discussion about relation between helicity and charge (and chirality) will be presented in Sec. 3.3. In Eq. 3.7,λ andλ are the helicities of Z and Z * and we have used in the second line the fact they must be equal in the Higgs rest frame. The angles are defined, for the two lepton pairs, as explained above and the minus for ϕ 2 arises due to the choice of the same reference frame for both the azimuthal angles as shown in Sec. 5. We see explicitly that the form of the angular dependence, encoded by the Wigner functions, is determined for each term in the sum by the polarisation of intermediate vectors. The production helicity amplitude A(h → ZλZ * λ ) is different for different polarisationsλ and receives BSM corrections (it is the only one under the assumptions applied here). Therefore, deviations from SM modify the differential distribution in the lepton angular variables. On the other hand, the Breit Wigner propagators do not depend on the helicity of the vector and thus can be factored out for all the terms with different angular modulations. Thus, the A ++ , A −− and A 00 amplitudes are respectively multiplied by the following angular functions, that are products of two d 1 λ,∆λ=1 : In the discussion above, we have assumed that one of the two Z bosons is produced on-shell, but both Z's could be off-shell, in which case the m Z m Z * terms in Eq. 3.2-3.6 become m Z * 1 m Z * 2 . However, the cross-section of a process with one on-shell Z is enhanced, as one can also see from Eq. 3.7; as a consequence, in the majority of the cases one pair of leptons will have invariant mass around the m Z mass (see Sec. 5 for a more detailed discussion). We should notice that one could consider the angles for the negative helicity fermion, in which case the d 1 λ,∆λ=−1 Wigner functions enter in the amplitude.

Helicity amplitudes for Z production through hZ¯ contact interactions
We now consider the g h Z contribution to the angular distributions derived in Sec. 3.1. Although we do not consider the effect of these interactions in our numerical analysis-as they can be probed more precisely in other processes as discussed in Sec. 2-we will keep them in all our analytical expressions.
We will show here how the g h Z contribution to the full h → 4 amplitude can be simply seen as a shift in the coupling of the leptons to the off-shell Z, g Z * 2 , in Eq. 3.7. In the considered hZ µ¯ γ µ operator, the fermion bilinear transforms as a vector under Lorentz transformations; in other words, the contact interaction has the form hZ µ V µ , where V µ is a spin 1 object. Therefore, the + − system involved in the g h Z vertex has intrinsic angular momentum J = 1. Furthermore, as previously explained, independently from the form of the interaction, in a h → ZV decay evaluated in the Higgs rest frame, Z and V must have same helicity. As a consequence, the M angular momentum quantum number for the + − system, in its center of mass frame, is fixed to be equal to the helicity of the emitted Z. Therefore, the dependence on lepton angular variables is given by d 1 M,∆λ=1 (θ, φ) and is determined by the helicityλ = M of the Z boson involved in the contact interaction; it is the same angular modulation that we would have if the leptons were coming from the decay of an additional intermediate Z. Furthermore, as the hZ µ V µ form of the contact term is analogous to the SM hZ µ Z µ vertex, we obtain the same amplitude for all the three helicity configurations. Thus, the g h Z contribution to h → Z2 → 4 corresponds to a shift in the SM M(h → ZZ * → 4 ) amplitude in Eq. 3.7 that can be express as a shift in the g Z * 2 coupling. Explicit computation confirms this heuristic reasoning and we obtain for the shift: So far, we have assumed that the Z is produced on-shell, but one could have a g h Z contact interaction with emission of an off-shell Z as well. However, it will be suppressed with respect to the on-shell case, due to the absence of resonance enhancement.

Visible angular modulation
In the previous sections, we have considered h → 4 helicity amplitudes, evaluating the angles for leptons with λ = + 1 2 . However, one cannot experimentally have access to the helicities of the final state leptons. For this reason, in the angular analysis, the electric charge is fixed instead; in particular, in our analysis we consider angles that describe the emission direction of negatively charged fermions (see Fig. 1 and Sec. 5 for details). Then, in the physical h → 4 process that can be studied at colliders, each final lepton has definite charge but not helicity. The squared amplitude, if expressed in terms of helicity amplitudes, should thus be summed over the four possible helicity configurations¯ 1 ∓ , where¯ and are respectively positively and negatively charged leptons, namely anti-fermions and fermions, and ± stand for the helicities (recall that these are equal and opposite within each fermion pair): The differential distribution is computed with respect to the measured i emission angles, that coincide with the θ i and ϕ i entering in Eq. 3.9 in the cases in which i λ has positive helicity λ = + 1 2 . On the other hand, in the terms where the negatively charged fermion is i − the measured angles θ i and ϕ i refer to the negative helicity leptons, which implies that the angular distribution is described by the Wigner functions d 1 λ,∆λ=−1 (θ, ϕ) = d 1 λ,∆λ=1 (π − θ, π + ϕ). Then, remembering that positive helicity corresponds to positive charge in the case of LH fermions and to negative charge for RH fermions, if we take into account the θ and φ (or equivalently ϕ) variables for negatively charged leptons, the angular definition will be left unchanged for RH leptons, while in case of LH chirality in the angular functions f ++ , f −− , f 00 of Eqs. 3.10-3.12 we should apply the substitution (θ, φ) → (π − θ, φ + π).
Thus, different helicity h → 4 amplitudes correspond to different chirality configurations for the final fermions and therefore are associated to different Z-lepton couplings. Then, the total squared amplitude of the observed process can be expressed as: is the helicity amplitude in Eq. 3.7 in which the angles are evaluated for positive helicity fermions and where the Z-lepton couplings g Z ( * ) l L(R) have been factorised out.

The Method of Moments
In this section we use the results of the previous section to show that the h → 4 squared amplitude can be written as a sum of a set of angular functions both in the SM and D6 SMEFT. We will describe the method for extraction of the coefficients of these functions, the so-called angular moments, and discuss the associated uncertainty estimates.

Angular moments for h → 4
The cross-section of the process gg → h → 4 , induced by the two contributions studied in Secs. 3.1 and 3.2, is obtained summing the squared helicity amplitudes as in Eq. 3.15. Therefore, as one can see by considering Eqs. 3.10-3.12, it is a linear combination of the following 9 functions of the final lepton angles: where the last 3 functions appear in the CP-odd terms, linear in theκ ZZ coupling and the angles are defined for negatively charged leptons as in Fig. 1. Our results agree with Ref. [86]. The angular moments are the coefficients of the 9 angular functions above, in the differential cross-section evaluated in the Higgs rest frame. As a function of the δĝ h ZZ , κ ZZ , κ ZZ and g h Zf coefficients they are: where contain the effect of the contact terms, g h Z via Eq. 3.13 and we have chosen a normalisation such that a SM 1 = G 4 . Here, γ a and γ b are defined in Eqs. 3.5 and 3.6 and The angular moments, a 5 , a 6 and a 8 , are numerically suppressed in the SM as well as the EFT interference term. This is because, once the g h Z are stringently constrained by other processes as discussed in Sec. 2, the SM value acts as a suppression factor. Thus, these moments contribute only marginally as far as the final numerical bounds are concerned.
To understand this suppression better, consider for example how the cos(θ 1 ) cos(θ 2 ) dependence arises in f 6 in Eq. 4.1. The cos(θ 1 ) cos(θ 2 ) term in the first and second terms in Eq. 3.15 does not change sign, but the substitution cos(θ 1 ) cos(θ 2 ) → − cos(θ 1 ) cos(θ 2 ) is applied in the third and fourth terms where one of 1 or 2 is a LH and negative helicity fermion. This gives a factor of 2 in the total squared amplitude, making this contribution numerically small. In fact, there are cases, like this one, in which the helicity-charge interplay leads to a partial or almost complete cancellation of the angular differential distributions that we would have observed if we could have access to the helicities of the final state leptons. In general, by averaging over all the helicity configurations, we loose part of the information contained a priori in the angular modulations of Eqs. 3.10-3.12.

The basic idea behind the method of moments
As seen in Sec. 4.1, the squared amplitudes for our present process, can be written as a set of angular structures, f i (θ 1 , θ 2 , φ), which are parameterised by the corresponding coefficients, the angular moments, a i . In this section, we explain how to extract these coefficients by taking the best-possible advantage of all the available angular information. Even though a full likelihood fit can be appropriate, here we consider the method of moments [60,87,88]. This method is transparent and advantageous, especially when the number of events is not very large [88]. For this method, an analog of Fourier analysis is used in order to extract the angular moments. Essentially, we seek weight functions, w i (θ 1 , θ 2 , φ), that can extract all the coefficients, a i uniquely, i.e., Upon assuming that these weight functions are linear combinations of the functions in the original basis, we can write We can use Eq. (4.6) to show that λ ij = M −1 ij , with, For the set of basis functions listed in Eq. (4.1), the corresponding matrix is given by, where we organise the basis functions according to their order in Eq. (4.1).
It is advantageous to translate to a basis such that M ij , and thus correspondingly its inverse λ ij , are diagonal. We can achieve this by the following orthogonal rotation, with ξ ± = (53 ± 9 √ 29). Thus, the weight functions in the rotated basis can be expressed as We can now convolute the various distributions from our events with these weight functions and extract the coefficients in the new diagonalised basis, {â 1 ,â 2 , a 3 , a 4 , a 5 , a 6 , a 7 , a 8 , a 9 }. (4.14) These coefficients can be rotated back in case we are seeking the moments in the original basis.

Moments estimates and estimation of uncertainties
In the last sections we have defined the angular moments a i and their extraction through weight functions defined over a continuous phase space; we want now to show how to estimate them starting from the observed experimental dataset, or in order to obtain projections, from Monte Carlo samples. We will also discuss how to estimate the uncertainty in this procedure.
One can notice, from Eq. 4.6, that the angular moments indicated there are the expectation values of the weights w j for a probability distribution i a i f i normalised to the squared amplitude. Changing the normalisation to the total number N of observed events, we can consider the random variableã for each event in the experimental dataset. Averaging over all events the observed value for the angular moments can be obtained, which is indeed the discretised version of Eq. 4.6. This is the procedure that must be used by the experiments to extract the angular moments.
In the absence of the true experimental dataset we have used, for our projections, separate Monte Carlo samples both with and without the EFT couplings turned on. The SM sample includes both Higgs and non-Higgs backgrounds. These Monte Carlo samples have a much larger number of events, N SM,EF T M C , than the number of events expected at 3 ab −1 , that we denote asN SM,EF T . As we will discuss in Sec. 5, for our statistical analysis to estimate the final bounds, we will take the SM expectation to be the null-hypothesis and assume that the experiments observe an excess over the SM because of the presence of EFT terms.
The weight functions w i , for a sufficiently large number of events, converge to a multivariate Gaussian distribution for which the estimates of the expectation values and of the covariance matrix are given bȳ In the above equation N SM,EF T is assumed to be a variable following the Poisson distribution with both mean and variance given byN SM,EF T . Note that thew i SM,EF T is supposed to be evaluated in an expected/observed dataset with a number of events equal toN SM,EF T , which is much smaller than the number of Monte-Carlo events. Thew i SM,EF T in Eq. 4.17, which gives a very good estimate of the w i expectation value, is a random variable whose covariance is given bȳ both for the SM and EFT cases. Theσ ij values decrease for increasingN , which is related to the fact that thew i estimates become more precise for larger number of events. The statistical uncertainties of the estimated mean values are computed as covariances of functions of the random variables N andw i : As we take SM to be our null-hypothesis, the uncertainties that will enter our final χ 2 function in Sec. 5 are the ones related to the SM expectation. We also take into account a flat systematic covariance on the SM prediction given by κ 2 syst a SM i a SM j where we will take κ syst = 0.02 following Ref. [89] 4 . a SM i includes all SM contributions, including the Higgs process. Then, the total covariance of the set of estimated SM angular moments is (4.23)

Monte Carlo samples and analysis setup
In this section, we discuss the collider analysis that helps us in obtaining bounds on the relevant operator combinations. We implement our UFO [90] model with the help of FeynRules [91]. This model is required to generate the signal samples, including the interference and the squared terms ensuing from the dimension-six interactions. We consider the gg → h → 4 process, where the irreducible backgrounds are composed of the quark-initiated qq → 4 and gluon-initiated gg → 4 processes, with = {e, µ, τ }. Reducible backgrounds arise from processes where jets can be misidentified as charged leptons in the fiducial region of the detector; these are dominated by Z/γ * + jets, which 4 One can deduce a total fractional systematic error of 0.039 = √ 0.035 2 + 0.016 2 + 0.006 2 , from Table  2 of Ref. [89] where the three values-respectively corresponding to the experimental error, the theory error for the SM Higgs process and the theory error for the non-Higgs background-have been added in quadrature. Note, however, that this is normalised with respect to the SM gg → h → ZZ * process whereas κsyst above is normalised with respect to the full SM cross-section including non-Higgs backgrounds. Our value κsyst = 0.02 is obtained using κsystNSM = 0.039N h SM , whereN h SM is the number of SM Higgs events. Note that, while Ref. [89] expresses this systematic error as a fraction of the SM Higgs rate, it actually includes, as a separate contribution, the error from the predominantly qq initiated non-Higgs background.
we generate as pp → + − + 2 jets. Negligible sources of fake backgrounds include tt, W W + jets, and W Z + jets 5 .
We generate Monte Carlo events considering a centre-of-mass energy of √ s = 14 TeV. The SM-and EFT-driven gg → h → 4 samples, as well as the reducible background Z/γ * + jets, are generated at leading order (LO) with MadGraph [92], including the full decay chain. The quark-initiated qq → 4 background samples are generated at nextto-leading order (NLO) with POWHEG BOX V2 [93][94][95]. The NNPDF31 lo as 0130 and NNPDF31 nlo hessian pdfs [96] PDF sets are used to generate MadGraph samples and qq → 4 events, respectively. The gluon-initiated gg → 4 background samples are generated at LO with MCFM 7 [97] using the CTEQ6L [98] PDF set. All events are further passed on to Pythia 8 [99] for parton shower and hadronisation. For the quarkinitiated background events, we apply a generator-level invariant-mass cut for each pair of opposite-sign same-flavour (OSSF) leptons of M + , − ≥ 4 GeV. For the remaining samples we require η ≤ 3, as well as ∆R( i , j ) ≥ 0.015, where ∆R = ∆η 2 + ∆φ 2 is the separation in the η − φ plane. In MadGraph, we impose an additional set of cuts, namely p 1 T ≥ 15 GeV, p 2,3 T ≥ 8 GeV and p 4 T ≥ 3 GeV. For the Z/γ * + jets samples we further apply p j T > 20 GeV, y j ≤ 3, ∆R(j, ) ≥ 0.015, ∆R(j m , j n ) ≥ 0.015, and M 2 ,2j ∈ [95,155] GeV. In the case of events generated using MCFM 7, we require p T ≥ 3 GeV, M + , − ≥ 2.5 GeV, as well as M 4 ≥ 70 GeV. Here, the indices on the leptons indicate their p T ordering, with 1 being the hardest lepton.
Following the recommendations of the LHC Higgs cross-section working group [100] (LHC HXSWG, CERN Report 4), we scale the LO gg → h production cross-section of the Monte Carlo sample to the N 3 LO-accurate prediction for M h = 125 GeV, obtaining an overall K-factor of 3.155, which we apply to the SM-, as well as to the EFT-driven samples 6 . We assume a flat NNLO/NLO K-factor of 1.1 for the qq → 4 background, given the differential cross-section for the quark-initiated process shown in Ref. [101]. For the NNLO/LO scaling of the gg → 4 samples, we apply a flat K-factor of 2.27, as considered in the experimental search described in Ref. [102]. We further adopt a conservative approach by using a flat K-factor of 0.91 [103] for the Z/γ * + jets events. After reweighting, the signal-to-irreducible background ratio S/B irr is found to be 0.00734. Here, by signal we mean the SM production of gg → h → 4 .
We base our analysis strategy on the experimental search 7 described in Ref. [2]. A simplified detector analysis is performed on the stable final-state particles using HepMC [104] and FastJet [105]. Visible objects are selected if they fulfill |η| < 4.7 and p T > 0.5 GeV. Electrons (muons) are preselected within the geometrical acceptance |η| < 2.5 (2.4), with p T > 7 (5) GeV, and are in turn isolated by demanding that the total hadronic activity 5 We take into account those processes that yield exactly four parton-level visible objects (charged leptons + jets) in the final state, where we consider channels with up to three jets. We find that the only non-negligible contribution arises from Z/γ * + jets. 6 Within our simulation framework we further set the width of the Higgs boson Γ h = 4.088 MeV, for consistency with the LHC HXSWG. 7 We validate our analysis against the experimental search, and find that our results are compatible with the experimental numbers within 96%.
Invariant mass distribution M (4 ) of the 4-lepton system after reconstruction of Z 1 Z 2 pairs. Stacked histograms for SM gg → h (red) signal, quark-initiated qq → 4 (light blue), gluoninitiated gg → 4 (dark blue), and Z/γ * + jets (green) backgrounds, are normalised to the expected number of events at the HL-LHC. Except for the quark-initiated background, distributions are scaled (2×) for visualisation purposes.
around a cone radius of R = 0.3 centred in the lepton's direction must be less than 35% of its p T . The overall missing transverse momentum / E T is the magnitude of the total transverse momentum calculated from all preselected particles, and its direction is opposite to these transverse momenta. Jets are clustered using the anti-k t algorithm [106] with a radius parameter 0.4 and p T > 30 GeV. We simulate the detector response in the reconstruction of electrons, muons, and jets, by applying a Gaussian smearing [107], as implemented in Rivet [108], to their energy, p T , and 3-momentum components, respectively. For leptons, these mass-and direction-preserving smearing functions are applied before selection, whereas for jets the mass-preserving smearing is applied after clustering. We assume a flat leptonic reconstruction efficiency of 0.92 and consider a rapidity-dependent jet-to-electron fake rate of 0.016 (0.044) for jets with y j < 1.48 1.48 < y j < 2.5 [109].
The selection is designed to extract signal candidates from events with no jet activity, / E T < 25 GeV, and exactly 2 pairs of OSSF leptons. The experimental treatment of preselected leptons is mimicked by requiring ∆R( i , j ) > 0.02, as well as M + , − > 4 GeV (irrespective of flavour), in order to suppress events with leptons originating from the decay of low-mass resonances. We further impose a cut on the leading lepton's p T > 20 GeV, and require that at least two of the sub-leading leptons have p T > 10 GeV. Pairs of OSSF leptons are combined into Z 1 Z 2 candidates, where Z 1 corresponds to the Z candidate with an invariant mass closest to the nominal Z-boson mass (91.1876 GeV) [110], and Z 2 is the   remaining one. Low-mass dilepton resonances produced along with an on-shell Z-boson are rejected by requiring M (Z 1 ) ∈ [40, 120] GeV, as well as M (Z 2 ) ∈ [12, 120] GeV. Finally, the mass range that characterises our signal region is defined as M (4 ) ∈ [118, 130] GeV, as shown in Fig. 2, which results in a S/B irr ratio of 1.37. Upon including the yield of reducible-background events into account, the signal-to-background ratio S/B gets reduced to 1.09. The effect of each cut on the fraction of retained events of the SM signal and irreducible backgrounds is shown in Table 2, and the invariant mass distributions of our surviving Z i candidates are shown in Fig. 3.

Angular extraction
Given the experimental limitations to determine the helicity of the final-state leptons, in what follows we restrict ourselves to define the various scattering angles with respect to the lepton with negative charge coming from the decay of the parent Z i boson (i = 1, 2), and which we refer to as − i as depicted in Fig. 1. In order to implement the analysis of the  . Azimuthal angle φ S differential distribution of the negatively-charged lepton (defined in the Higgs rest frame S ) for the SM-driven gg → h → 4 process (solid black), as well as the SM + interference terms for the CP -even κ ZZ = 0.5 (dotted red) and CP -oddκ ZZ = 0.5 (dashed blue) operators. differential distributions with respect to the angular variables described in Secs. 3-4, we Lorentz-boost the 4-lepton system to the centre-of-momentum frame S , where the Higgs boson is at rest, and we have back-to-back Z-boson momenta. In frame S , we construct a Cartesian coordinate system {x,ŷ,ẑ} as follows: theẑ axis points in the direction of motion of Z 1 ;ŷ is normal to the plane generated byẑ andB, whereB = (0, 0, 1) corresponds to the unit vector defining the beam direction in the laboratory frame; finally,x completes the right-handed set. Hence, we calculate the azimuthal angle ϕ i,S formed in frame S between − i and the scattering plane as where − i,ê corresponds to the projection of − i onto theê axis in frame S , and ϕ ∈ [0, 2π). The azimuthal angle φ S between the planes formed by the lepton pairs in frame S , used in our angular analysis, is then defined as tan φ S ≡ tan(ϕ 2,S − ϕ 1,S ), (5.1) with φ ∈ [0, 2π). The azimuthal distribution for the process gg → h → 4 is depicted in Fig. 4. Finally, each pair of OSSF leptons is further boosted to the rest frame S of its parent Z-boson, where we have back-to-back leptons momenta. The polar angle θ i,S

Total Rate
All moments Figure 5. Bounds at 68% CL on the CP -even anomalous couplings. The green band shows the bound from the total rate which keeps a flat direction, κ ZZ ≈ 3.7 δĝ h ZZ , unconstrained. The blue ellipse shows our final bounds including all the angular moments. The red ellipse also shows the results of an angular moment analysis but considering only the interference between the EFT and SM terms. between − i (in frame S ) and the direction of motion of its parent Z i boson 8 (in frame S ) is defined as where k M i ,N corresponds to the 3-momentum of particle M i in frame N , and θ ∈ [0, π].  Figure 6. Bounds at 68% CL on the CP -even anomalous couplings after combination with the results of the angular moment analysis of the pp → W h/Zh processes carried out in Ref. [51]. The blue band shows the results of this work, the green ellipse the bounds from pp → W h/Zh processes. The red ellipse is the combination of the present work and moments analysis presented in Ref. [51] for the pp → V h processes. Finally, the dashed yellow ellipse shows the final bound after combination with pp → h → W W → 2 2ν process in Ref. [73].

Results
To obtain our bounds we define a χ 2 function as follows, where the covariance Σ ij is defined in Eq. 4.23. As explained in Sec. 2 and 3.2 we do not include the g h Zf parameters of the hZ¯ contact terms as they can be constrained very stringently using other processes. Using the above χ 2 function, we obtain the 68% CL bounds shown in Fig. 5. The green band shows the bound obtained if we include only 8 The definition of the aforementioned angles is taken with respect to the SM-driven gg → h → 4 process, where both pairs of leptons are produced from the decay of a Z boson. However, it is important to note that a g h Zf insertion, where a pair of leptons has no parent Z boson, is also possible, and hence the direction of motion of the dilepton system in the Higgs rest frame S needs to be taken into account.
information about the total rate of the process. It is clear that there is a flat direction, κ ZZ ≈ 3.7 δĝ h ZZ , which can be constrained only by introducing the differential information of the angular moments. We also show in red the bound obtained if we only include the interference term between the SM and EFT which almost coincides with the bounds that include the EFT squared term in blue. This implies that, if Λ is the UV cutoff of the EFT, a truncation of the cross-section at the 1/Λ 2 order includes already most of the BSM effects encoded by dimension-six operators; thus, there is no need of taking into account the 1/Λ 4 level and an analysis without dimension-eight operators can be considered consistent. If we assume no systematic uncertainties we obtain the bound |κ ZZ | < 0.05 for δĝ h ZZ = 0 which is comparable to the MELA bound |κ ZZ | < 0.04 [64].
To obtain the bounds in Fig. 5 we have taken,κ ZZ = 0. We have, however, checked that a non-zeroκ ZZ hardly changes the plot, i.e. to a very good approximation we will obtain the same bounds ifκ ZZ is marginalised over. This is because there is hardly any contribution to χ 2 from a 7 -a 9 , the only moments that contain a term linear inκ ZZ . The moments a 1 -a 6 that give the dominant contribution, on the other hand, depend on theκ 2 ZZ which results in the χ 2 function being only mildly dependent onκ ZZ . For the same reason we obtain a very weak bound onκ ZZ , |κ ZZ | 0.5 (6.2) after marginalising over δĝ h ZZ and κ ZZ . This is unfortunately not competitive with the projection, |κ ZZ | < 0.05 from the pp → W h/Zhh(bb) processes in Ref. [51].
We can combine the above bounds with the projections from the pp → W h(bb)/Zh(bb) and pp → h → W W → 2 2ν processes in Ref. [51] and Ref. [73]. In order to combine these different processes, we need to utilise the stringent constraints on the Zf f couplings, assume that the function f in Eq. 2.6 and EFT deformations that rescale the h → bb branching ratio can be independently constrained in a global fit including all the relevant Higgs physics processes. We also need to use EFT correlations, derived in Ref. [51], between anomalous gauge-Higgs couplings involving the W -boson and those involving the Z boson. We can then combine the results of the angular moment analysis of pp → W h(bb)/Zh(bb) with leptonic decays of the W/Z in Ref. [51] and the bound on the total rate for pp → h → W W → 2 2ν in Ref. [73]. 9 The final results are shown in Fig. 6. We see that the complementarity of the h → V V * → V and pp → V h processes, i.e., the fact that these processes probe very different linear combination of δĝ h ZZ and κ ZZ , results in strong percent level bounds on these couplings.
A possible criticism towards our approach is that it is based on leading order matrix elements and is not optimised to include the effects of parton shower, detector effects and selection cuts. Before concluding this section we want to compare our results with a BDT analysis which does not have this shortcoming. We perform a simple BDT analysis with three variables, i.e., θ 1 , θ 2 and φ. We choose the Gradient Boosted Decision Tree (BDTG) algorithm as it can properly handle negative-weight events which arise upon including NLO samples. As δĝ h ZZ only rescales the SM matrix element and leaves no differential signatures, we only vary κ ZZ to obtain the bound |κ ZZ | < 0.052, with a zero systematic uncertainty hypothesis. This can be compared with our earlier derived bound, |κ ZZ | < 0.051. As we can see, the numbers are very similar and the nominal difference is perhaps owing to statistical reasons. Thus, this validation shows that the presence of the experimental effects mentioned above fortunately do not affect the sensitivity of our method for this particular final state.

Conclusions
We have carried out a fully differential study of the golden Higgs decay channel h → 4 . The leptonic final state can be accurately reconstructed to give a wealth of differential information. In the Higgs rest frame three angles completely determine the direction of the final state leptons. A thorough differential study of the resulting three-dimensional space is one of the main experimental tools to probe the tensor structure of the Higgs coupling to gauge bosons which includes new contributions in the dimension-six SMEFT.
In this work we show that the full angular distribution can be written as a sum of a set of basis functions both in the SM as well as in the dimension-six SMEFT. The coefficients of these functions, the so-called angular moments, therefore encapsulate the full angular information of the process. We derive the analytical expressions for these angular moments including dimension-six SMEFT deformations. We then use the method of moments to extract these angular moments from our Monte Carlo sample, which is simulated and analysed using a strategy that closely follows the LHC experiments.
We finally use the extracted angular moments in the SM and dimension-six SMEFT to obtain projections for bounds on all the relevant gauge-Higgs coupling deformations parametrised by δĝ h ZZ , κ ZZ andκ ZZ in Eq. 2.1. Our final results in Fig. 5 show that the angular moment analysis is crucial in eliminating flat directions that arise if one takes into account only the total rate of the process in the SM and SMEFT. Finally we combine our results with those of Ref. [51] where a similar analysis using angular moments was carried out for the pp → V h process as shown in Fig. 6. This combination allows us to obtain the strongest reported bounds on the anomalous gauge-Higgs couplings.
acknowledges the grant received from the IPPP, where the major part of this work was done.