Higgs decay into four charged leptons in the presence of dimension-six operators

We study the indirect effects of New Physics in the Higgs decay into four charged leptons, using an Effective Field Theory (EFT) approach to Higgs interactions. We evaluate the deviations induced by the EFT dimension-six operators in observables like partial decay width and various kinematic distributions, including angular observables, and compare them with the contribution of the full SM electroweak corrections. The calculation is implemented in an improved version of the event generator Hto4l, which can provide predictions in terms of different EFT-bases and is available for data analysis at the LHC. We also perform a phenomenological study in order to assess the benefits coming from the inclusion of differential information in the future analyses of very precise data which will be collected during the high luminosity phase of the LHC.


Introduction
Now that a scalar particle, resembling the Higgs boson of the Standard Model (SM), has been discovered [1,2], the characterization of its properties represents one of the major tasks of the LHC physics programme. Besides the intrinsic importance of confirming that the new 125 GeV spin-0 boson is the Higgs particle of the SM, the precise measurement of its properties represents the opportunity to search for indirect hints of new physics (NP). Up to now experimental analyses have extracted bounds on NP parameters in the so called κ− framework, which considers modifications proportional to SM couplings [3,4]. However, the κ− framework does not provide a gauge invariant parametrization of NP and it cannot capture the effects of physics beyond SM (BSM) on kinematic distributions. The current experimental bounds allow a deviation of 10% in the Higgs to gauge bosons (HV V ) and about 20% deviation in Higgs fermion couplings.
Given the lack of a clear evidence of NP signals in the LHC data already analyzed, it is reasonable to assume that the scale Λ, where new particles would eventually appear, is well separated from the energy scale of the SM spectrum. If this is the case, physics at the electroweak (EW) scale can be adequately described by Effective Field Theory (EFT) methods. The building blocks of the EFT Lagrangian are the SM fields. The low-energy effects of new possible heavy degrees of freedom are captured by effective operators with mass dimension D larger than four. Since it provides a model-independent parametrization of possible NP effects, the EFT approach has become the phenomenological standard for the study of indirect signals of NP. Regarding the Higgs sector, the majority of these studies have interpreted the LHC data on Higgs production and decay modes to derive constraints on the D = 6 parameters. It should be noted that these constraints do depend on certain model dependent assumptions. Model independent approaches to Higgs physics have been also applied to differential cross sections in order to investigate their resolving power to extract information on the presence of anomalous couplings. In particular, due to its particularly clean signature and non-trivial kinematics, the Higgs decay into four leptons, i.e. H → ZZ * → 4 , has been considered in a number of works appeared in the literature [5][6][7][8][9][10][11][12]. The signal strength in H → ZZ * and in the gluon-gluon fusion (ggF) production channel after combining the CMS and ATLAS results is µ = 1.13 +0. 34 −0.31 [13]. In Refs. [5][6][7] NP effects in H → 4 decays are parametrized in terms of specific anomalous Higgs vertices, while in Refs. [9,11,12] the language of pseudo-observables (PO) is adopted. Finally, in Ref. [8], the observability of anomalous Higgs couplings in the H → Z(→ + − ) + − channel has been studied in a EFT framework by considering the differential decay width dΓ/dq 2 , as well as the relevant angular asymmetries.
While in EFT new gauge-invariant operators are added to the SM Lagrangian, the POs provide a parametrization of NP effects at the amplitude level and consequently they are process specific. On the other hand, the PO approach is more general, in the sense that it does not require any assumption on the underlying UV-complete theory. It is important to stress that the connections among POs of different observables become transparent once the mapping from EFT Wilson coefficients to POs has been set up.
In this paper, we study the H → 4 decay in the Standard Model Effective Field Theory (SMEFT) framework. In particular, we perform a reanalysis on the effects of the effective operators entering H → 4 decay channel both at the level of integrated partial width and on the relevant and experimentally accessible distributions. We compare the NP effects with the contributions of the full SM electroweak corrections. We also perform a phenomenological study in view of the outstanding integrated luminosity which is expected to be reached with the High Luminosity LHC (HL-LHC) project (3 ab −1 ), that will allow to test the SM validity at a precision level which has never been achieved before. With this study we aim at highlighting the importance of angular distributions in constraining D = 6 Wilson coefficients.
The rest of the paper is organized as follows. In Section 2, we introduce the phenomenological EFT Lagrangian in the so-called Higgs basis [14], which is advocated in the literature to study the NP signatures in the Higgs sector. In Section 3, we provide information on the H → 4 matrix elements implemented in the new version of the Hto4l 1 code [15]. Our numerical results and our study in the context of HL-LHC are presented in Section 4. We draw our conclusions in Section 5. Further details are provided in the Appendices. In Appendix A, we present our results for the H → 4e integrated partial width. In Appendix B we detail the computation of the H → 4 BSM matrix elements in the Warsaw basis [16] and the SILH basis [17,18]. In Appendix C we collect the formulae used in the analysis outlined in Section 4.

Theoretical framework
As mentioned above, the EFT approach is based on the hypothesis that the scale Λ of NP is much heavier than the EW scale. In this framework the decoupling of new particles is described by the Appelquist-Carazzone theorem [19]. Once the heavy degrees of freedom have been integrated out, the low-energy effects of new particles are captured by an arbitrary number of effective operators. The resulting effective Lagrangian takes the form Moreover, in linear EFT the spontaneous breaking of the SU (2) L ⊗ U (1) Y down to U (1) em arises from the non-vanishing vacuum expectation value (vev) of the complex Higgs doublet. Neglecting the D = 5 lepton flavor violating operator [20], the leading BSM effects are expected to be captured by D = 6 operators. Under the hypotheses of lepton and baryon number conservation, flavor universality and a linear realization of the EWSB, all possible BSM deviations can be parametrized by a basis of 59 D = 6 CP-even operators and 6 additional CP-odd bosonic operators. Different bases of D = 6 operators, which are related by equations-of-motion for fields, have been proposed in the literature. The most popular choices are the Warsaw basis [16] and the SILH (Strong Interacting Light Higgs) basis [17,18]. The choice of the basis is usually led by the convenience to minimize the number of operators that are necessary to parametrize the BSM effects on a given class of processes. However, since the operators of these two bases are manifestly invariant under the SU (2) L ⊗ U (1) Y symmetry, the connection between Wilson coefficients and phenomenology can be rather cumbersome. In this work we adopt the socalled Higgs basis [14], which has been designed to parametrize the effects of new physics in the Higgs sector in a more transparent way. As in the BSM primaries [21], the Higgs basis operators are written in terms of mass eigenstates. It has been argued that the coefficients of this parametrization of NP can be obtained as a linear transformation from any other basis. These transformations are chosen in order to map particular combinations of Wilson coefficients of a given basis into a subset of anomalous couplings of the mass-eigenstates Lagrangian extended to D = 6 operators. These are called independent couplings. The number of independent couplings is the same of any other basis. Once a maximal subset of independent couplings has been identified, the remaining dependent couplings can be written as linear combinations of the independent ones. We would like to point out that the Higgs basis is advocated in the literature to perform the leading order EFT analyses of the Higgs data. A complete picture of next-to-leading order EFT calculations in the Higgs basis is not yet clear [14].
In this section we limit ourselves to describe the parts of the effective Lagrangian which are relevant for the Higgs decay into four leptons. For the derivation of the complete effective Lagrangian in the Higgs basis framework we refer to [14]. All the kinetic terms are canonically normalized and there is no Z-γ kinetic mixing. The kinetic and mass terms for the gauge bosons are the same of the SM, except the W boson mass, which receives a correction of the form: Although the precision measurement of W mass gives the possibility to derive information on BSM physics and an EFT framework can be used in this context [22], it is important to stress that δm is presently very well constrained by experiments: δm = (2.6±1.9)·10 −4 [23], so that the effects proportional to δm would be irrelevant for Higgs physics. Moreover, if the underlying UV-complete theory preserves the custodial symmetry, δm = 0 by hypothesis. For these reasons δm = 0 is assumed in the following analysis. The operators giving rise to anomalous contributions entering the Higgs decay into four charged leptons can be divided in five classes. The first and most relevant class includes the effective operators affecting the Higgs couplings to gauge bosons. Regarding the neutral sector, the effective Lagrangian takes the form: A µνÃ µν +c Zγ e g 2 1 + g 2 where the convention to absorb the suppression factor 1/Λ 2 in the effective coefficients has been adopted. In the above, V µν = ∂ µ V ν − ∂ ν V µ and,Ṽ µν = 1 2 µνρσ V ρσ for both V = A, Z. g 1 and g 2 are coupling parameters of the U (1) Y and SU (2) L gauge groups, respectively. Of the six CP-even couplings in Eq. (2.3) only five are independent. We choose c γ as dependent coupling, which is then expressed as the following linear combination: For the sake of generality, we include CP-odd couplings parametrized byc V V in our calculation. Note that if one assumes that the Higgs particle is a pure CP-even eigenstate, CP-odd operators are not allowed 2 . The second class of operators is given by the anomalous contributions to Z vertex while the third class gives rise to HV contact interactions L HZ D=6 = 2 If a linear realization of the SU (2) L ⊗ U (1) Y symmetry is assumed, the contact terms in Eq. (2.6) are generated by the same operators which give rise to vertex corrections in Eq. (2.5). In the Higgs basis they are set to be the same 3 , In this scenario, the coefficients for the contact interactions are constrained by the EW precisions tests performed at LEP and their effects are expected to be rather small. However, there are scenarios in which the coefficients in Eqs. (2.5-2.6) can be independent (see for instance Refs. [14,28]). In this work we also assume flavor universality, so that g Zee L,R = g Zµµ L,R and g HZee L,R = g HZµµ L,R . This assumption is very much consistent with LEP data on Z couplings. Any violation of this assumption can be checked by comparing H → 2e2µ with H → 4e and H → 4µ [11].
The last two contributions involve dipole interactions between Z bosons and leptons and the dipole contact interactions of the Higgs boson. These terms are proportional to lepton masses and in the m l → 0 limit can be safely neglected. Moreover, as a consequence of the linearly realized electroweak symmetry in the D = 6 Lagrangian, the dipole parameters are proportional to the respective lepton dipole moments, which are tightly constrained by experimental data and usually neglected in the LHC analyses. Note that a contact term involving H and four leptons can only be generated by D = 8 operators. One would be sensitive to such a contact term in the kinematic region where the 4 invariant mass is much higher than the Higgs mass which is not the case in the on-shell decay of the Higgs boson.

Computational details
In order to study the possible BSM deviations in the Higgs decay into four charged leptons we have considered the effective Lagrangian 4 where L SM is supplemented by the D = 6 contributions in Eqs. (2.3)-(2.6). The master formula for the LO decay width, in the presence of D = 6 operators reads In addition to the anomalous part in the HZZ and Z ¯ couplings, the presence of D = 6 operators gives rise to tree-level Hγγ and HZγ and HZ ¯ vertices which are not present in the SM Lagrangian. The Feynman rules for these anomalous vertices have been derived by implementing the effective Lagrangian of Eq. (3.1) in FeynRules 2.0 [30]. For massless leptons we get, Q . In the previous expressions p 1 and p 2 are the incoming momenta of gauge bosons and The calculation of new matrix elements for H → 2e2µ and H → 4e/4µ has been carried out by means of the symbolic manipulation program FORM [31], and they have been included in a new version of the code Hto4l, which is publicly available. As in other Monte Carlo tools for Higgs physics, such as HiGlu [32,33], Hawk [34][35][36][37] and HPair [38,39], the new version of Hto4l provides the possibility to compare present and future Higgs data with theoretical predictions derived in an EFT context.
Since we have neglected the lepton masses, the matrix elements for 4e and 4µ are the same. As a consistency check we have compared the value of the matrix elements implemented in Hto4l with the ones generated with MadGraph5@MC_NLO [40] for several phase-space points, finding excellent agreement.
Few important remarks are in order: first of all we note that the quadratic part |M D=6 | 2 of Eq. (3.2) is suppressed by a factor 1/Λ 4 . From the point of view of the EFT expansion, it contributes at the same level of D = 8 operators. Moreover, different bases of D = 6 operators are equivalent only at the order of 1/Λ 2 and they differ by terms which are of order 1/Λ 4 . It follows that predictions obtained by using only D = 6 operators are not complete at the O(1/Λ 4 ). There are different approaches in the literature regarding the treatment of quadratic contributions in the analyses. One approach consists in making linear approximation for the theoretical predictions and including the quadratic contributions in the estimation of the theoretical uncertainty. In this context, the constraints derived in one basis can be translated to other bases. Another approach keeps always the quadratic contributions in the calculations. The latter improves the accuracy of the calculation whenever the contribution of D=8 operators is subdominant [41]. Pragmatically, we have included the quadratic contributions in our calculation with the possibility of switching them on and off in the code.
In order to guarantee flexibility in the choice of the basis, a provision of calculating H → 4 matrix elements in SILH and Warsaw bases which are not affected by the basic assumptions of the Higgs basis, is also made. For that a separate dictionary between the anomalous coupling parameters appearing in the Feynman rules (Eqs. 3.3-3.7) and the Wilson coefficients of the SILH and Warsaw bases is implemented in the code and it is listed in Appendix B.

Numerical results
In this section we present some numerical results obtained with the new version of Hto4l for the H → 4 decay channel in the presence of D = 6 operators of the Higgs basis. The results have been obtained with the same SM input parameters as in Ref. [15]. In the Higgs basis, the {G F , α, M Z } input parameter scheme is assumed. A shift to the {G F , M W , M Z } input parameter scheme, which we have adopted, introduces corrections proportional to δm (see Eq. 2.2) in couplings and parameters dependent on the input parameters.
Since the anomalous vertices V Hγγ and V HZγ enter the calculation of the H → 4 partial decay width we expect an important BSM contribution coming from the kinematic configurations with one of the lepton pair invariant masses close to zero. In order to get rid of these contributions which would be hardly accessible by the experiments, we have implemented a lower cut of 15 GeV on the leading and subleading same-flavor opposite-sign (SFOS) lepton pair invariant masses.

BSM predictions for the partial decay width
The modification of the H → 2e2µ decay width in the presence of the CP-even and CP-odd parameters of the Higgs basis can be parametrized as, The absence of linear terms in CP-odd parameters is related to the fact that the partial decay width is a CP-even quantity. The coefficients of the linear and quadratic terms are given by 5 , The corresponding coefficients for H → 4e are given in Appendix A. Note that in the above we have intentionally kept δg HZ i and δg Z i independent of each other to cover the scenario in which new physics parametrization leads to additional contributions in δg Z i . In the Higgs basis we must set δg HZ i = δg Z i . The relative importance of various parameters of the Higgs basis in modifying the partial decay width can be inferred from the size of the coefficients derived above. To illustrate the relative effect of the parameters more clearly, in Fig. 1 we plot the ratio in Eq. (4.1) by scanning each parameter in the range between -1.0 and +1.0. Among CP-even parameters related to the HV V (V = γ, Z) couplings, the change in partial decay width due to c γγ is the smallest, while δc Z , which gives rise to a SM-like anomalous coupling, changes the width maximally. Due to different propagator effects the effect of c Zγ is larger than that of c γγ . The contact interaction parameters, however, modify the width the most because of no propagator suppression. In the CP-even case, these scan plots display the importance of the linear terms with respect to the quadratic terms. For instance, we find that for c γγ and c Zγ the quadratic contributions dominate over the linear ones in most of the parameter space, leading to an overall enhancement of the decay width. On the other hand, for δc Z , c ZZ and c Z , the linear terms play an important role and the decay width can become smaller than its SM value in certain regions of parameter space. Also, the effect of c ZZ and c Z on the partial width is opposite in nature. For contact interaction parameters the quadratic terms dominate over the linear ones, except for a small region of parameter space between 0 and 0.5 (-0.5) for δg HZ L (δg HZ R ) where the ratio goes below 1. As mentioned before, the CP-odd parameters contribute to the total Higgs decay rate only at the quadratic level leading to the ratio always greater than 1. Among the CPodd parameters, the change of the decay width due toc Zγ is the largest one while the corresponding change due toc ZZ is the smallest one. Information on CP-odd linear terms can be accessed from specific kinematic distributions which we discuss later.
It is important to stress that some of these parameters are already constrained by the available experimental data from LEP and LHC. For instance, by using LHC Run-I data [42], c γγ and c Zγ are constrained respectively at the 10 −3 and 10 −2 level. On the contrary δc Z , c ZZ and c Z are loosely constrained. An approximate degeneracy, which quad. corresponds to a strong correlation, is found between c ZZ and c Z (ρ ij = −0.997). Including the LEP data on W W production, δc Z and c Zγ become more constrained and the flat direction between c ZZ and c Z is also lifted to some extent (ρ ij = −0.96) [43]. These conclusions assume linear dependence of Higgs signal strength observables on parameters.
It has been argued that there is no model independent constraint on c ZZ and c Z because including contributions which are quadratic in these parameters would dramatically change the corresponding best-fit values and the relative uncertainties. To this end, more data are needed and the complementary information coming from kinematic observables will be helpful to improve the constraints on these coefficients [44,45].
Furthermore, the couplings of the Z boson to charged leptons are constrained by considering Electroweak precision data (see Refs. [46,47] for recent analyses where SMEFT theoretical errors are taken into account). In our framework these constraints are also applicable to the parameters of the ZH contact interactions. To obtain any constraint on the CP-odd parameters, it is necessary to go beyond the linear approximation for Higgs observables. Interpreting the results obtained on CP-odd parameters of the SILH basis in Ref. [48] using current Higgs data, we find thatc γγ is constrained at 1% level. However, the allowed values for |c Zγ | and |c ZZ | can be as large as 0.7 and 0.5 respectively. In the following we focus on the parameters which are loosely constrained by the data and have non-negligible effects on partial decay width.

BSM predictions for kinematic distributions
In this section we use the new version of Hto4l to simulate the decay of the Higgs boson into four charged leptons in the presence of D = 6 operators at the differential level. The study of distributions can provide complementary information to the analyses of signal strengths and BRs. For the sake of simplicity we consider one parameter at a time while the remaining ones are artificially set to zero. More sophisticated analyses, where correlations among various coefficients are taken into account, are beyond the scope of this article.
The parameters of our interest are c ZZ , c Z andc Zγ . Moreover, since H → 4 decay can provide information on the contact ZH interaction, we will also consider the effect of δg HZ i independent of δg Z i . To emphasize the characteristic effects of these parameters on distributions, we consider a scenario in which the parameters lead to the same deviation in partial decay width. In particular, we choose the benchmark values for these parameters by considering an excess of 30% in Γ BSM (H → 2e2µ). In table 1, the benchmark values are reported by keeping only the interference terms and also by including the quadratic terms in the calculation.
Among the observables taken into account, the most sensitive ones to BSM kinematic effects turn out to be • the subleading lepton pair invariant mass M sub 6 ; • the angle φ between the decay planes of the two intermediate gauge-bosons in the Higgs rest-frame; • the angle ∆θ e − µ − between the electron and the muon in the Higgs rest-frame.
In Figs. 2-5, we compare the BSM predictions for the normalized distributions of these observables with the SM ones at Leading Order (LO) and at Next-to-Leading Order EW accuracy matched to a QED Parton Shower (NLOPS in the following), i.e. the highest SM theoretical accuracy achievable with Hto4l. In order to better highlight the kinematic effects we also plot the normalized ratios R norm. , defined as: where X is a generic observable, while i = c i or i = NLOPS. Note that to calculate the BSM excess in each bin this ratio has to be multiplied by 1.3. Continuous lines in the plots refer to distributions obtained by considering only the effects of interference, while for the dashed ones quadratics effects have been also taken into account. Several remarks are in order: • The angular variables turn out to be more sensitive to BSM kinematic effects than M sub .
• Among the CP-even parameters considered in the analysis, c ZZ and c Z have a larger impact on the normalized distributions than δg HZ L,R (Figs. 2-4). As far as φ and ∆θ e − µ − are concerned, the BSM effects are larger than the SM higher order corrections, while the effects of contact interactions seem to be of the same order of magnitude as EW corrections (Figs. 3-4).
• The effect of c ZZ on M sub monotonically increases as we go towards the tail of the distribution reaching an excess close to 40%. In case of c Z , the ratio grows mildly in the beginning and starts decreasing beyond 33 GeV (upper panel plot in Fig. 2).
• The effects of c ZZ and c Z on angular observables are opposite in nature. In the presence of c ZZ more events fall in the central φ region, while in the presence of c Z , 6 The leading lepton pair invariant mass is defined as the SFOS lepton pair invariant mass closest to the Z boson mass.
the edges get more populated (upper panel plot in Fig. 3). Similarly, looking at the ∆θ e − µ − distributions, we find that c ZZ , unlike c Z , puts more events in the region where the angle between electron and muon is greater than 90 degrees (upper panel plot in Fig. 4).
• For c ZZ and c Z the effects of quadratic terms depends on the considered observable and in general are not negligible. The difference between predictions obtained by including only 1/Λ 2 terms, with respect to those including also the quadratic contributions, turns out to be larger for the angular observables than on M sub (upper panel plots in Figs. 2-4). For instance, as far as c Z is concerned, the quadratic contributions can give up to a further 5% difference at the level of normalized ratios in some of the bins.
• On M sub , the effects of δg HZ L and δg HZ R are the same (lower panel plot in Fig. 2). However, the angular observables can be used to discriminate the two parameters (lower panel plots in Figs. [3][4]. Since the interference and quadratic values obtained for them are small and close to each other, the contribution of quadratic terms over the linear one is very minute.
• We find that for our choice of values forc Zγ , φ is the most sensitive observable. The angle φ is a CP-odd observable and, unlike the partial decay width, it is sensitive to the linear term inc Zγ . Also, for the same reason, it can provide information on the sign of the parameter. These features are clearly visible in Fig. 5.

Future prospects at HL-LHC
One of the main opportunities of the HL-LHC program is to enable precise measurements of the Higgs boson properties, such as the presence of anomalous couplings to bosons and fermions. It has been shown that kinematic distributions, such as the p T of the Higgs boson, can significantly improve the multi-dimensional parameter fit [49]. In this section we present the results of a χ 2 analysis carried out in the context of the High-Luminosity LHC (HL-LHC). The study has the illustrative purpose to assess how H → 2 2 angular observables can be exploited to constrain SMEFT coefficients in future analyses of LHC data (see Ref [44,45]). Due to the large limits and to the strong correlation arising from current constraints, the analysis is restricted to the c ZZ − c Z . At LHC, the H → 4 decay has been observed mainly in the gluon-gluon fusion channel. The observed signal strength using the 7 and 8 TeV LHC data is given by µ 4 ggF = 1.13 +0.34 −0.31 [13], while using 13 TeV LHC data the observed signal strength is µ 4 ggF = 1.20 +0.22 −0.21 [50]. The current data in H → 4 channel alone cannot be used to constrain the parameters c ZZ and c Z . Therefore, at present, any meaningful bounds on these parameters can be obtained by including data in other decay channels which have been observed in production modes sensitive to c ZZ and c Z , i.e. vector boson fusion (VBF) and associate production of Higgs and vector boson (VH) [42].

χ 2 fit with normalized distributions and asymmetries
In the first stage of the analysis we consider normalized distributions, and we look for the kinematic observables turning out to be particularly sensitive to c ZZ and c Z effects. The analysis is performed through a sample of pp → H → 4 pseudo-events. For the sake of simplicity the sample is restricted to the ggF production mode. The sample has been generated by interfacing POWHEG [51] to Hto4l, according to the procedure described in Ref. [15] and exploiting the Narrow Width Approximation (NWA). The expected number of SM events is derived by assuming 3 ab −1 of integrated luminosity. The adopted values for ggF cross section and H → 2 2 ( , = e, µ) branching ratios are taken from Ref. [14]. The events are then selected according to the experimental cuts adopted in ATLAS [52]. Eventually, the accepted events are scaled down by 20% to take into account the lepton reconstruction efficiency (95% for each lepton), leading to a sample of ∼ 6000 reconstructed events, in good agreement with the number found in Ref. [53]. Besides the distributions defined in the previous section, the two asymmetries

3)
A cθ 1 cθ 2 = 1 σ dΩ sgn {cos θ 1 cos θ 2 } dσ dΩ , (4.4) are sensitive to CP-even D = 6 coefficients (as already pointed out in Ref. [8]). In the above definitions φ is the angle between the decay planes of the two intermediate vector bosons, while θ 1 (θ 2 ) is the angle between the lepton produced in the decay of the nonresonant (resonant) Z boson and the direction opposite the Higgs boson, in the non-resonant (resonant) Z boson rest-frame. The χ 2 for distributions and asymmetries can be written as follows: N D is the number of bins of the D-th distribution. The quantity f SM i is the fraction of events, generated as described above, falling in the i-th bin of the SM distribution, while f BSM i is the fraction of expected events in the presence of a given combination of c ZZ and c Z . This last quantity is calculated by reweighting the events with a program in which the Hto4l BSM matrix elements have been implemented. As we deal with normalized quantities, we assume that the systematic and theoretical uncertainties are cancelled to a large extent in the ratio. Accordingly, σ i and σ A in Eqs. (4.5-4.5) are just the one-sigma statistical uncertainties where n i is the number of events falling in the i-th bin. The 68.3% Confidence Level (CL) contour plots for the aforementioned distributions and asymmetries are displayed in Fig. 6. The contour plot for the asymmetry A φ overlaps exactly the one for the φ angle and therefore is not shown. The contour plot for the combined χ 2 defined by the sum is also displayed. In the next section, we perform a global analysis using signal strengths where production channels other than ggF are also considered. Unlike ggF, the production channels VBF and V H depend upon c ZZ and c Z and this dependence is quite strong (C.1-C.3). This feature is taken into account in Fig. 6. The regions marked by green lines correspond to parameter-space points, i.e. to c ZZ and c Z values, driving any of these cross sections to negative values. We remark that these "unphysical" regions arise because in the linear approximation, the cross sections are not positive definite. Few remarks are in order: • the χ 2 analysis of single distributions and asymmetries is not sufficient to get closed contour plots in the range (-1,1) for c ZZ and c Z . This is mainly motivated by the fact that one can choose values for c ZZ and c Z whose effects cancel in the sum; • among the four angular observables taken into account, the asymmetry A cθ 1 cθ 2 turns out to be the least sensitive observable to c ZZ and c Z ; • the negative correlation between c ZZ and c Z resulting from the analysis of Run-I data arise in our analysis of distributions as well. However, the correlation for φ distribution and A φ asymmetry is larger than for ∆θ e − µ − and A cθ 1 cθ 2 . This feature allows to rule out a region of the parameter space, with c ZZ > 0 and c Z < 0.
• the contour plot derived from Eq. (4.9) lies inside the "physical region" of the parameter space. The constraints on c ZZ and c Z are tighter, while the correlation is not removed.

χ 2 fit using signal strengths
At the high-luminosity run of the LHC, the H → 4 decay channel is likely to be observable in production channels other than ggF. Moreover, we will also have access to the kinematic distributions in H → 4 . Therefore, in the context of HL-LHC, it would be interesting to study the sensitivity of the future data in constraining parameters c ZZ and c Z mainly using the H → 4 decay channel. Our main motivation is to highlight the effect of angular distributions φ and ∆θ e − µ − in the fit. Our analysis is based on minimizing the χ 2 function built using signal strengths as observables. The signal strength in a given production channel i and decay channel f is defined as The expressions for µ i and µ f , in presence of parameters c ZZ and c Z are taken from [44] where linear approximation in parameters is adopted. The µ 4 both at inclusive and differential levels is calculated using the Hto4l code with input parameter choice and kinematic cuts mentioned above. We have summarized all these expressions in Appendix C. The χ 2 function in terms of signal strengths is given by, where the one-sigma uncertainties σ f i are taken from Ref. [54]. For future data we assume . We find that the effect of ∆θ e − µ − distribution in the fit is similar but less important than that of φ distribution. All the results in the following are presented using φ distribution.
In Fig. 7 we display the region plots with 68.3% CL in c ZZ − c Z plane. In each plot we compare the fits obtained by using the decay signal strength µ 4 at the inclusive (1 bin) and differential (3 bins) levels. To understand the effect better we divide the fit in many categories depending upon which production channels are used in the fit. The differential effects are the largest when only ggF and ttH production channels are included in the fit (Fig. 7a). This is not surprising given the fact that both these production channels do not depend on parameters c ZZ and c Z . The positive correlation between the parameters is governed by the decay signal strength µ 4 (see Eqs. (C.6-C.9)). When we use ggF and VBF channels (Fig. 7b), the correlation becomes negative due to a stronger dependence of µ VBF (C.1) on the parameters which is opposite in nature to that of µ 4 . Once again, the φ distribution in the fit improves the bounds significantly. Using ggF and V H in the fit (Fig. 7c) the constraints on the parameters become tighter and the effect of including distributions is mostly visible at the edges. Note that the constraints become stronger because the dependence of V H channel on parameters is very strong (C.2, C.3). Thus, when all the production channels are combined (Fig. 7d) the constraints on parameters are governed by the production channels rather than the H → 4 decay channel and the distribution still leads to a noticeable improvement in the fit. In this case, we also derive 1σ constraints on each parameter when the other parameter is ignored in the fit. The 1σ errors on c ZZ and c Z resulting from the inclusive (differential) fit are ±0.032 (±0.026) and ±0.014 (±0.013) respectively. In Fig. 8a, we relax theoretical and experimental systematic uncertainties in the future data, and perform the fit using only µ 4 at the differential level in all production channels. Clearly, improvements in precision calculations would reduce theoretical uncertainties, allowing to put tighter constraints on the parameters. We also study the relevance of φ distribution in the fit when other decay channels like 2 2ν and γγ, which depend on the parameters c ZZ and c Z , are included (Fig. 8b). The γγ decay channel is included in all the production channels while, the 2 2ν decay channel is included only in ggF and VBF production channels. Notice that the γγ partial decay width does not depend on the parameters at the LO and in our analysis we have used the one-loop expression derived in [55] and quoted in [44] for the Higgs basis parameters. At inclusive level, the 1σ errors on c ZZ and c Z become ±0.02 and ±0.01 respectively. These errors on individual parameters are consistent with those obtained in Ref. [44]. We find that the improvements in bounds due to distribution are marginal but still visible.

Conclusions
In the present work we have investigated possible NP effects in the Higgs decay into four charged leptons using an EFT approach to Higgs interactions. We have adopted the Higgs basis for the computation of the BSM matrix elements for the H → 2e2µ and H → 4e/4µ channels. We considered both CP-even and CP-odd operators and we mostly focused on those parameters, which are weakly constrained by LHC Run-I data. Since the H → 4 channel can provide information about the presence of ZH contact interactions, we have also considered the scenario in which they are independent of anomalous Z interactions and therefore unconstrained by electroweak precision data.
For the sake of illustration, we have presented numerical results for the H → 2e2µ channel. As a first step, we have analyzed the impact of D = 6 operators on the partial decay width. As the information on partial decay width is not sufficient in discriminating different parameters, we have also studied some kinematic distributions of particular experimental interest. We have found that with the help of angular observables (φ and ∆θ e − µ − ) it is possible to distinguish different parameters with values which would lead to the same modification in the partial decay width. In the case ofc Zγ , the angle φ would be quite useful in deriving a stronger constraint on the parameter as it captures the information on the linear piece in the parameter.
Aiming to assess the impact of differential information in future analyses, we have performed a global analyses in the context of HL-LHC. From the preliminary study on differential distributions and asymmetries we derived that the angle φ is the most sensitive observable to c ZZ and c Z effects. In our global analysis based on signal strengths we find that the impact of the angular information is quite dependent on the production channels that are taken into account. The largest improvements are observed in the ggF and ttH channels. When also VBF and V H channels are included, the benefits coming from the inclusion of angular information are moderate but are still noticeable. More sophisticated analyses, where other coefficients and differential information coming from other production channels are also considered, are beyond the purposes of this work and will be considered in future extensions of the present study.
The above phenomenological study has been carried out through a new version of the Hto4l event generator, which allows to study the effects of D = 6 operators in the H → 2e2µ and H → 4e/4µ channels. The BSM matrix elements are calculated in the Higgs basis. The code also allows independent calculations in SILH and Warsaw bases. As an option, the possibility of including quadratic D = 6 contributions (of the order of 1/Λ 4 ) on top of pure 1/Λ 2 interference contributions is given. Since it can be easily interfaced with any event generator for the Higgs production, Hto4l can be used in association with other MC tools for the full simulation of Higgs events in an EFT framework.
A H → 4e partial decay width

B EFT Dictionaries
In the latest version of the Hto4l code, the calculation of H → 4 BSM matrix elements can be performed independently in the Higgs basis, the SILH basis and the Warsaw basis. The BSM matrix elements are implemented in the Higgs basis. For predictions in SILH and Warsaw bases the Higgs basis parameters are seen just as the coefficients corresponding to specific Lorentz structures in the Feynman rules of section 3, and following dictionaries between the Higgs basis parameters and the Wilson coefficients of SILH and Warsaw bases are used. δg Z R = − 1 2 K He + g 2 1 g 2 2 (g 2 1 + g 2 2 ) 2 (K W + K B ). (B.13) The corrections to weak boson masses are given by, δm 2 W = 0, (B.14) The dependence of s W and e on g 1 and g 2 becomes, c ZZ = 4v 2 1 (g 2 1 + g 2 2 ) 2 g 2 1 C HB + g 2 2 C HW + g 1 g 2 C HW B (B.24) c ZZ = 4v 2 1 (g 2 1 + g 2 2 ) 2 g 2 1 C HB + g 2 2 C HW + g 1 g 2 C HW B (B.25) H + C In addition, s W and e have a modified dependence on g 1 and g 2 given by,

B.3 Input parameter scheme
In the input parameter scheme {G F , M Z , M W }, the parameters in the Feynman rules of the Warsaw basis are given by, where, and, corresponding changes in s W and e should also be taken into account [56]. These values have to be used in the SM Feynman rules. In the Feynman rules proportional to the Wilson coefficients, the parameters should be simply replaced by their SM definitions. For SILH basis, one needs to replace C

C Production and decay signal strengths
The signal strengths for production channels at √ŝ = 13 TeV are given by, All the above expressions are taken from Refs. [42,44]. The partial decay widths for the H → 4 decay at the inclusive and differential levels are calculated using the Hto4l code. The signal strengths in this channel are given by,