Measurements of inclusive and differential cross sections for the Higgs boson production and decay to four-leptons in proton-proton collisions at $\sqrt{s}$ = 13 TeV

Measurements of the inclusive and differential fiducial cross sections for the Higgs boson production in the H $\to$ ZZ $\to$ 4$\ell$ ($\ell$ = e,$\mu$) decay channel are presented. The results are obtained from the analysis of proton-proton collision data recorded by the CMS experiment at the CERN LHC at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 138 fb$^{-1}$. The measured inclusive fiducial cross section is 2.73 $\pm$ 0.26 fb, in agreement with the standard model expectation of 2.86 $\pm$ 0.1 fb. Differential cross sections are measured as a function of several kinematic observables sensitive to the Higgs boson production and decay to four leptons. A set of double-differential measurements is also performed, yielding a comprehensive characterization of the four leptons final state. Constraints on the Higgs boson trilinear coupling and on the bottom and charm quark coupling modifiers are derived from its transverse momentum distribution. All results are consistent with theoretical predictions from the standard model.


Introduction
The discovery of the Higgs (H) boson in 2012 by the ATLAS and CMS Collaborations [1][2][3] is a major confirmation of the correctness of the theoretical approach involving the electroweak (EW) symmetry breaking mechanism [4][5][6][7][8][9]. Subsequent measurements of the properties of this particle [10-13], including its mass, quantum numbers, and couplings, further confirmed the consistency of these measurements with the standard model (SM) predictions.
The H boson decay into four charged leptons (H → ZZ → 4ℓ, ℓ = e, µ), with its fully reconstructible final state and large signal-to-background ratio, has been one of the pillars for the characterization of this particle since its discovery. Several properties of the H boson were measured in this decay channel at the CERN LHC, based on the Run 1 data set at center-of-mass energies of 7 and 8 TeV and on the Run 2 data set at 13 TeV. These include the determination of its mass, spin and parity [14][15][16][17][18][19], width [20-23], inclusive and differential fiducial cross sections [18,[24][25][26][27][28], and tensor structure for interactions with a pair of gauge bosons [17,19,21,[29][30][31]. The most precise value of the H boson mass to date, measured by the CMS Collaboration, is m H = 125.38 ± 0.14 GeV, obtained from the combination of the H → ZZ → 4ℓ and H → γγ decay channels from the analysis of the Run 1 and 2016 Run 2 data sets [32]. From the analysis of the full Run 2 data set, the CMS Collaboration reported the first evidence for the off-shell H boson production in events with a final state of two Z bosons decaying into either four charged leptons, or two charged leptons and two neutrinos, with a measured value of the H boson width of Γ H = 3.2 +2.4 −1.7 MeV [33]. The H boson production is often experimentally characterized via the so called simplified template cross section (STXS) framework, which defines mutually exclusive phase space regions designed to maximize the experimental sensitivity to physics beyond the SM (BSM) effects and reduce, at the same time, the theoretical model dependence in the measurements [34]. The ATLAS and CMS Collaborations published results of cross section measurements in the STXS framework using the full Run 2 data set in the H → ZZ → 4ℓ [27,35] and other decay channels [36][37][38][39][40][41][42]. The ATLAS and CMS Collaborations recently published results obtained from the combination of all the decay channels focusing on measurements of simplified template cross sections (STXS) [43] and H boson couplings [44], respectively.
Fiducial cross section measurements constitute a complementary approach for the characterization of the H boson production and decay that provide a set of less model-dependent results by unfolding detector effects from the data, thus allowing a direct comparison with state-ofthe-art theoretical predictions. The ATLAS and CMS Collaborations published fiducial cross section measurements in the H → γγ [45,46], H → WW [47][48][49], and H → ZZ → 4ℓ [28,35] decay channels using the full Run 2 data set. The CMS Collaboration also published results in the H → ττ [50] decay channel, while the ATLAS Collaboration presented results from the combination of the H → γγ and H → ZZ → 4ℓ decay channels [51]. This paper presents measurements of inclusive and differential cross sections for the H boson production in the H → ZZ → 4ℓ decay channel using data from proton-proton (pp) collisions recorded with the CMS detector at the LHC in 2016-2018 and corresponding to an integrated luminosity of 138 fb −1 . To reduce the model dependence, all the measurements are performed within a fiducial phase space region defined to closely reproduce the experimental acceptance and reconstruction-level selection criteria.
Differential cross sections are measured for several kinematic observables sensitive to the H boson production and its decay into four leptons, providing a complete characterization of this channel and coverage of the entire fiducial phase space. This includes the measurement of six double-differential cross sections. Fiducial cross sections are also measured in bins of matrix element kinematic discriminants sensitive to possible anomalous couplings of the H boson to vector bosons. This provides a valuable test of the SM predictions and may reveal possible BSM physics.
The analysis builds upon the methods used in previous measurements of H boson properties in the four-lepton decay channel [25,35], featuring the latest CMS Run 2 calibrations and a reduction by ∼ 40% for the leading systematic uncertainty in lepton reconstruction and selection efficiencies.
The measurement of the fiducial cross section in bins of transverse momentum p T of the H boson (p H T ) is also used to set constraints on the H boson trilinear self-coupling and on the coupling modifiers of the H boson to b and c quarks. The CMS Collaboration recently reported constraints on the H boson self-coupling from the combination of several decay channels using the full Run 2 data set [44]. These constraints are obtained from the interpretation of the STXS results, while in this paper an alternative and complementary approach using differential cross section measurements is explored. The ATLAS Collaboration set limits on the coupling modifiers of the H to b and c quarks from the combination of the H → ZZ → 4ℓ and H → γγ decay channels using the 2016-2018 Run 2 data set [51]. The CMS Collaboration reported similar results using the 2016 Run 2 data set and combining the H → ZZ → 4ℓ and H → γγ decay channels [52]. The constraints derived from the analysis presented in this paper supersede the ones reported for the H → ZZ → 4ℓ channel in Ref. [52]. This paper is organized as follows. The CMS detector is briefly described in Section 2. The data set used is presented in Section 3, along with a description of the simulated signal and background samples. The event reconstruction techniques and the selection criteria used to identify H boson candidates are outlined in Section 4. The definition of the restricted phase space region where the differential cross sections are measured is given in Section 5. A complete description of all the kinematic observables, with a particular emphasis on matrix element discriminants, is presented in Section 6. The background modeling is presented in Section 7. The signal modeling and the statistical procedure adopted in the extraction of the inclusive and differential cross sections are presented in Section 8. The systematic uncertainties that affect the measurement are described in Section 9. The results of the analysis and their comparison to the SM expectations are outlined in Section 10. In Section 11 the measurement of the fiducial cross section in differential bins of p H T are used to set constraints to the trilinear self-coupling of the H boson and to its couplings with charm and bottom quarks. A summary highlighting the main findings of the analysis is given in Section 12.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. The electromagnetic calorimeter consists of 75 848 lead tungstate crystals, which provide coverage in pseudorapidity |η| < 1.48 in a barrel region (EB) and 1.48 < |η| < 3.0 in two endcap regions (EE). Preshower detectors consisting of two planes of silicon sensors interleaved with a total of three radiation lengths of lead are located in front of Table 1: Thresholds applied on the p T of the leading/subleading leptons in each data-taking period for the main dielectron (e/e), dimuon (µ/µ), and electron-muon (e/µ, µ/e) HLT algorithms.
e/e (GeV) µ/µ (GeV) e/µ, µ/e (GeV) 2016 17 with a pair of top quarks (ttH) [70]. Events produced via the ggH mechanism are simulated at NLO with POWHEG 2.0 and reweighted to match the predictions at next-to-next-to-leading order in the strong coupling, including matching to a parton shower (NNLOPS) [71] as a function of the p H T and of the number of jets in the event. The gg → ZH contribution to the ZH production mode is simulated at leading order (LO) using JHUGEN 7.3.0 [72][73][74][75][76]. The ggH production mechanism is simulated at NLO also with MADGRAPH5 aMC@NLO using the HC NLO X0 UFO-HEFT model in the 5 flavor scheme [77]. The H boson is produced in association with 0, 1, or 2 jets in the final states, merged with the FxFx scheme. The top quark mass is set to 173.0 GeV in the simulation, but finite top mass effects in loops are filtered out. The H boson production in association with b quarks is not considered in this analysis as its impact on the unfolded distributions is expected to be negligible with respect to all the other production modes. The decay of the H boson to four leptons is modeled with JHUGEN 7.0.2. The simulation of the various production and decay modes is based on the theoretical predictions from Refs. , which are summarized in Ref. [34].
The main background processes originate from ZZ production from quark-antiquark annihilation and gluon fusion. The former is simulated at NLO in pQCD with POWHEG 2.0 [101], while the latter is generated at LO with MCFM 7.0.1 [102][103][104][105]. The reducible background contribution arising from the production of Z bosons with associated jets (Z+jets) is estimated with the data-driven technique already used in Ref. [35] and described in Section 7.2.
An additional sample of Drell-Yan plus jets (DY+jets) events is produced with MADGRAPH5 aMC@NLO 2.4.2 for validation studies and for the training of the boosted decision tree (BDT) used for the identification and isolation requirements on electrons, as described in Section 4. All other simulated samples are used to model the signal shape, estimate backgrounds, optimize the analysis strategy, and evaluate the systematic uncertainties.
All Monte Carlo (MC) generators are interfaced with PYTHIA to simulate the parton showering and hadronization effects. Version 8.230 [106] is used for the three data-taking years with the CUETP8M1 tune [107] for 2016 and the CP5 tune [108] for 2017 and 2018. Parton distribution functions (PDFs) are taken from the NNPDF3.0 set [109] for the three data taking periods.
The response of the CMS detector is modeled using the GEANT4 [110,111] package. The simulated events are reconstructed with the same algorithms used for data and the distribution of the number of pileup events per bunch crossing is reweighted to match that observed in the data.

Event reconstruction and selection
The particle-flow (PF) algorithm [112] aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The energy of photons is obtained from the ECAL measurement. The energy of electrons is determined from a combination of the electron momentum at the PV, as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The momentum of muons is obtained from the combined information of the tracker and the muon chambers. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energy deposits.
The information from the ECAL and the tracker is combined to reconstruct electrons [56] with p e T > 7 GeV within the geometrical acceptance of the detector, defined by the pseudorapidity region |η e | < 2.5. The identification of electrons is performed with a BDT algorithm sensitive to the presence of bremsstrahlung along the electron trajectory, the geometrical and momentumenergy matching with the corresponding cluster in the ECAL, the features of the electromagnetic shower in the ECAL, and observables that discriminate against electrons originating from photon conversions. The isolation sums for electrons, defined similarly as for muons, are included in the BDT discriminant. This choice is proven to enhance the suppression of nonprompt electrons originating from hadron decays and from overlap of neutral and charged hadrons within jets [56] and has a better performance than a cutoff-based approach using the relative isolation. The BDT for the electron identification and isolation is implemented using the XGBOOST library [113]. The training is performed on a dedicated sample of DY+jets simulated events. Electron samples are divided into six mutually exclusive categories defined by two p T ranges (7 < p e T < 10 GeV and p e T > 10 GeV) and three η selections corresponding to the central barrel (|η e | < 0.8), outer barrel (0.8 < |η e | < 1.479), and endcaps (1.479 < |η e | < 2.5). The BDT is trained separately for the three data-taking periods and the selection requirements are defined to achieve the same signal efficiency for the three data taking periods (97% for p e T > 10 GeV; 80% for p e T < 10 GeV in the barrel; 74% for p e T < 10 GeV in the endcap). The information from the silicon tracker and the muon system [58] is combined to reconstruct muons with p µ T > 5 GeV and |η µ | < 2.4. The matching between inner and outer tracks is performed starting either from the tracks in the silicon trackers or from those reconstructed in the muon system. Cases where inner tracks are matched to segments in only one or two muon detector layers are also considered, to cope with very-low-p T muons that do not traverse the entire detector. Muon objects are selected from the muon track candidates by applying loose requirements on the track in the muon system and the inner tracker, taking into account also their compatibility with small energy deposits in the ECAL and HCAL.
A requirement on the relative isolation, I µ < 0.35, is introduced to discriminate between muons from Z boson decays and those originating from hadron decays within jets, where I µ is defined as: and where ∑ p charged T is the scalar sum of the transverse momenta of charged hadrons originating from the PV, whereas ∑ p neutral T and ∑ p γ T are the scalar sums for neutral hadrons and photons, respectively. The isolation requirement is defined using a cone of radius ∆R = 0.3 around the muon direction at the PV, with the angular distance between two particles i and j defined as ∆R(i, j) = √ (∆η i,j ) 2 + (∆ϕ i,j ) 2 . The quantity p µ,PU T in Eq. (1) is defined from the p T sum of all the charged hadrons i not originating from the PV as p µ,PU T ≡ 0.5 ∑ i p µ,PU T,i , where the factor of 0.5 corrects for using only the charged particles in the isolation cone [114]. The p µ,PU T contribution is subtracted in the definition of I µ to correct for energy deposits arising from pileup interactions.
Final-state radiation (FSR) photons arising from Z boson decays are recovered as follows. The PF photon candidates with |η γ | < 2.4 are considered as FSR objects if they have p γ T > 2 GeV and a relative isolation I γ < 1.8, where I γ is defined similarly as for muons in Eq. (1). These FSR candidates are associated with the closest lepton in the event and are not retained if ∆R(γ, ℓ)/(p γ T ) 2 > 0.012 GeV −2 and ∆R(γ, ℓ) > 0.5. For each lepton, the FSR candidate with the lowest value of ∆R(γ, ℓ)/(p γ T ) 2 , if any, is selected. The photon candidates identified from the FSR recovery algorithm are excluded from the computation of the muon isolation.
Nonprompt leptons from decays of hadrons or photon conversions are suppressed based on the impact parameter significance. This variable is defined as the ratio of the 3-dimensional impact parameter, computed with respect to the position of the PV, to its uncertainty, and leptons are rejected if the value of this quantity is greater than 4.
The leptonic decays of known dilepton resonances are used to calibrate the momentum scale and resolution of electrons and muons in bins of p ℓ T and η ℓ , as described in Refs. [56,58]. Efficiencies for the lepton reconstruction and selection are measured in several bins of p ℓ T and η ℓ by means of the tag-and-probe technique using samples of Z boson events both in data and simulation. Simulated yields are corrected by the measured efficiency ratio between data and simulation.
For each event, hadronic jets are clustered from reconstructed particles using the infrared-and collinear-safe anti-k T algorithm [115] with a distance parameter of 0.4 [116]. The jet momentum is computed from the vectorial sum of all particle momenta in the jet, and is found in simulation to be, on average, within 5 to 10% from the true momentum over the whole p T spectrum and detector acceptance. Additional pp interactions within the same or nearby bunch crossings can contribute with additional tracks and calorimetric energy deposits, increasing the apparent jet momentum. To mitigate this effect, tracks identified as originating from pileup vertices are discarded and an offset correction is applied to correct for remaining contributions. Jet energy corrections are derived from simulation studies so that the average measured energy of jets becomes identical to that of particle-level jets. In situ measurements of the momentum balance in dijet, photon + jet, Z + jet, and multijet events are used to account for any residual differences in the jet energy scale between data and simulation [117]. Additional selection criteria are applied to remove jets potentially dominated by instrumental effects or reconstruction failures. Only jets with p jet T > 30 GeV, |η jet | < 4.7, and a distance parameter of ∆R(ℓ/γ, jet) > 0.4 from all selected leptons and FSR photons, are considered. Jets not satisfying the tight identification criteria and the criteria corresponding to the tight working point of the pileup jet identification algorithm described in Ref. [114] are also discarded.
The PF objects mentioned above serve as input to the event selection, which targets events containing at least four well-identified and isolated leptons originating from the PV and possibly accompanied by a FSR photon. The FSR photons are included in the invariant mass computations. The event selection, which closely follows that employed in Ref. [35], is detailed below.
The Z boson candidates are formed from pairs of same-flavor and opposite-charge leptons (e + e − , µ + µ − ) with an invariant mass within 12 < m ℓ + ℓ − < 120 GeV. Two such pairs are required to create ZZ candidates, where the Z boson candidate with invariant mass closest to the world-average Z boson mass [118] is referred to as Z 1 , whereas Z 2 denotes the other Z boson candidate. Three mutually exclusive subchannels are defined from the flavors of the four leptons in the event: 4e, 4µ, and 2e2µ.
The ZZ candidates must fulfill additional requirements designed to improve the sensitivity to H boson decays. The Z 1 candidates are required to have an invariant mass larger than 40 GeV. All lepton pairs (ℓ i , ℓ j ) must be separated by an angular distance of ∆R(ℓ i , ℓ j ) > 0.02. Events must contain two leptons with p T > 10 GeV and at least one with p T > 20 GeV. In the 4e and 4µ channels, where the same four leptons can be used to build an alternative Z a Z b candidate, candidates with m Z b < 12 GeV are not considered if Z a is closer to the world-average Z boson mass than Z 1 and the event is rejected. This suppresses events with an on-shell Z accompanied by a low-mass dilepton resonance (e.g., J/ψ or Υ). The invariant mass of the four possible opposite-charge lepton pairs (irrespective of flavor), computed without FSR photons, must satisfy m ℓ + ℓ ′− > 4 GeV in order to further suppress events with leptons originating from hadron decays in jet fragmentation or from leptonic decays of low-mass resonances. The ZZ candidates are retained if the invariant mass of the four-lepton system m 4ℓ is larger than 70 GeV.
In events where more than one ZZ candidate satisfies the selection requirements above, the one with the largest scalar sum of transverse momenta of the two leptons defining the Z 2 is retained.
Finally, only events with 105 < m 4ℓ < 160 GeV are considered for the statistical analysis.

Fiducial phase space definition
Cross sections are measured in a fiducial phase space defined to match closely the experimental acceptance of the reconstruction-level selections. The fiducial phase space is defined at generator-level, following the strategy adopted in previous H → ZZ → 4ℓ analyses [18,35]. It relies on requirements on the lepton kinematics and isolation, and on the event topology, in order to minimize the model dependence of the results.
The definition of the fiducial phase space is summarized in Table 2. The events are retained if the leading (subleading) lepton has p T > 20 (10) GeV. Additional electrons (muons) that may be present in the event are required to have p T > 7 (5) GeV and |η| < 2.5 (2.4). Lepton isolation is ensured by requiring the scalar sum of the p T of all stable particles, i.e., those particles not decaying in the detector volume, within a cone of radius ∆R = 0.3 to be less than 0.35 times the p T of the lepton. Neutrinos, FSR photons, and leptons (electrons and muons) are not included in the computation of the isolation sum to enhance the model independence of the measurements, following the findings of Ref. [25]. Events passing these requirements are retained if they have at least two same-flavor, opposite-sign lepton pairs. The pair with invariant mass closest to the world-average Z boson mass [118] is labeled as Z 1 and it must have 40 < m Z 1 < 120 GeV. The second Z boson candidate is referred to as Z 2 and it must have 12 < m Z 2 < 120 GeV. Each lepton pair ℓ i , ℓ j must be separated by ∆R(ℓ i , ℓ j ) > 0.02, while any opposite-sign lepton pair must satisfy m ℓ + ℓ ′− > 4 GeV, reflecting the selection criteria used at reconstruction level.
Leptons at the fiducial level are considered as dressed, i.e., FSR photons are collected within a cone of radius 0.3. Jets do not enter in the definition of the fiducial phase space, but they are used when dealing with jet observables. Jets at the fiducial level are built with the anti-k T clustering algorithm with a distance parameter of 0.4 out of stable particles, excluding neutrinos. Jets are retained if they satisfy p jet T > 30 GeV and |η jet | < 4.7, similarly to the condition used at reconstruction level. Only jets with no leptons inside a cone of radius 0.4 are kept. Requirements for the H → ZZ → 4ℓ fiducial phase space Lepton kinematics and isolation Leading lepton p T p T > 20 GeV Sub-leading lepton p T p T > 10 GeV Additional electrons (muons) p T p T > 7(5) GeV Pseudorapidity of electrons (muons) |η| < 2.5 (2.4) Sum of scalar p T of all stable particles within ∆R < 0.3 from lepton < 0.35p T Event topology Existence of at least two same-flavor OS lepton pairs, where leptons satisfy criteria above Inv. mass of the Z 1 candidate 40 < m Z 1 < 120 GeV Inv. mass of the Z 2 candidate 12 < m Z 2 < 120 GeV Distance between selected four leptons ∆R(ℓ i , ℓ j ) > 0.02 for any i ̸ = j Inv. mass of any opposite sign lepton pair m ℓ + ℓ ′− > 4 GeV Inv. mass of the selected four leptons 105 < m 4ℓ < 160 GeV

Observables
Fiducial cross sections are measured in bins of several kinematic observables sensitive to the H boson production and decay pp → H → ZZ → 4ℓ, of which a schematic representation is given in Fig. 1.
The decay of the H boson to four leptons is fully described by the invariant mass of the two Z boson candidates, three angles describing the Z boson decays (Φ, θ 1 , θ 2 ), and two angles connecting production to decay (Φ 1 , θ * ). The angle θ * is defined in the H rest frame as the angle between the beam axis and the direction of the Z 1 candidate. Φ and Φ 1 are the azimuthal angles between the three planes constructed from the H decay products and the decay products of the two Z bosons in the H rest frame. The θ 1 and θ 2 angles are defined in the Z 1 and Z 2 rest frames, respectively, as the angles between the Z boson direction in the H boson rest frame and the direction of the negative decay lepton. The set comprising these seven observables is hereafter referred to as ⃗ Ω H→ZZ→4ℓ (θ * , θ 1 , θ 2 , Φ, Φ 1 , m Z1 , m Z2 |m 4ℓ ) and can be used to build matrix element discriminants sensitive to the 4ℓ decay, as detailed in Section 6.1. Fiducial cross sections are measured in differential bins of these observables, except for the θ angles for which their cosine is used. The distributions as a function of p H T and pseudorapidity of the reconstructed H boson are also measured.
Fiducial cross sections are also measured in differential bins of the number of associated jets (N jets ) and p T of the leading (p j 1 T ) and subleading (p j 2 T ) jet in the event. For events with two or more jets the properties of the dijet system constituted by the two leading jets are assessed by measuring differential cross sections in bins of its invariant mass (m jj ), of the difference in pseudorapidity (∆η jj ), and of the difference in azimuthal angle (∆ϕ jj ) between the two jets. The angle ∆ϕ jj is defined to be invariant under the exchange of the two jets as follow: where the vectors ⃗ j 1,2 represent the direction of the leading and subleading jet in the laboratory frame, and the unit vectorsĵ T1,2 the corresponding transverse component. This definition is also independent of the choice of the positive z axis direction,ẑ.
The rapidity-weighted jet vetoes T max C and T max B are also studied. These are defined, following Ref. [119], as: where y j and m j T are the rapidity and transverse mass of the jet, defined from its mass m and momentum p as m j T = m 2 + p 2 x + p 2 y , while y H is the rapidity of the H boson. The value of each observable is computed for each jet in the event and its maximum value is taken for each event. Since their resummation structure is different from the canonical p j T , they give complementary information on the properties of jets in an event and can be used as a test of quantum chromodynamics (QCD) resummation. The 0-jet phase space can be redefined using these observables. The events with no jets are defined as the ones with T max C < 15 GeV and T max B < 30 GeV, where the values of these cuts are chosen accordingly to the findings of Ref. [119]. In the following, these events will be defined as 0-jet|T max C and 0-jet|T max B , respectively.
The properties of the H+jet(s) system are also studied by measuring differential cross sections in bins of the transverse momentum and invariant mass of the H plus leading jet system for events with at least one jet, or of the H plus leading and subleading jet system for events with at least two jets. The observables characteristic of the H + j(j) system can be defined only in events with at least one (two) jets. In all other cases, an underflow bin is introduced to consider all events for which the observable is undefined.

Matrix element discriminants
The JHUGEN and MCFM generators are used to compute the matrix element probability P i for an event to arise from a physical process i, given the value of the reconstructed invariant mass of the four-lepton system m 4ℓ . These probabilities are defined as a function of ⃗ Ω H→ZZ→4ℓ and retain the maximal information on the underlying physics content of each event. Hence, the P i ( ⃗ Ω H→ZZ→4ℓ ) probabilities are used to construct likelihood-ratio-like matrix element discriminants sensitive to the difference between two physical processes a and b, when considering two production mechanisms, or to be used to test a BSM hypothesis against the SM scenario. These matrix element discriminants have been widely used in the context of H → ZZ → 4ℓ analyses, from the measurement of the H boson properties [35] to the constraints on possible anomalous couplings [31]. The general structure of these discriminants is an adaptation of the more classic likelihood ratio, properly rescaled to ensure that the discriminants are always bounded between 0 and 1. Two types of kinematic discriminants can be built to test the compatibility between signal ("sig") and alternative ("alt") hypotheses and their interference ("int"): where P sig and P alt are the probabilities of an event under the two considered hypotheses, given their kinematic properties ⃗ Ω, and P int is the probability for the interference between the two model contributions ("sig" and "alt"). This definition of D int is bounded between −1 and 1 for any value of D alt .
A total of six matrix element discriminants sensitive to different values of possible anomalous couplings of the H boson to vector bosons are considered. The general scattering amplitude describing the interaction between a spin-zero H boson and two spin-one gauge bosons V 1 and V 2 can be written, following the conventions of Ref. [31], as: where v is the vacuum expectation value of the H potential, The a i and κ i terms correspond to the strengths of vector boson couplings, following the notation adopted in Ref. [31]. In particular, the a 3 CP-odd term is expected to be null in the SM and is sensitive to possible BSM effects that would result in CP violation. The a 2 term corresponds to the CP-even contribution to the HVV coupling and is sensitive to possible BSM contributions from heavy H bosons. The κ 1,2 /(Λ 1 ) 2 and κ 3 /(Λ Q ) 2 terms are sensitive to possible physics at a Table 3: Matrix element kinematic discriminants considered in the analysis. Some discriminants have a special label to identify the targeted Higgs boson property rather than the name of the coupling. D dec 0-is sensitive to a CP-odd Higgs boson, D dec CP is the observable sensitive to the CP-mixing, and D dec 0h+ is sensitive to heavy CP-even Higgs boson. coupling in any process involving an on-shell photon. However, the H → 4ℓ channel contains events featuring an off-shell photon, i.e., H → Zγ * → 4ℓ, that can be used to study the Λ Zγ 1 coupling. Table 3 details the set of kinematic discriminants considered and the couplings to which they are sensitive. The index "dec" indicates that only decay information is used to build these discriminants.
Differential cross sections are measured in bins of these six matrix element discriminants under the SM hypothesis. The compatibility of the measurements with the SM predictions is assessed by comparing the results with the discriminants built for alternative BSM scenarios, where HVV anomalous couplings are introduced by modifying the a i and κ i values in Eq. (6) with respect to their SM values.
Tables 4 and 5 summarize the bin boundaries for all the observables considered in this analysis that target the H boson production and the H → ZZ → 4ℓ decay, respectively.

Irreducible backgrounds
Irreducible ZZ background contributions arising from qq annihilation or gluon fusion are estimated from simulation. The former is simulated at NLO in pQCD with POWHEG 2.0 and reweighted to NNLO using a K factor computed as a function of m ZZ exploiting the NNLO computation of the qq → ZZ fully differential cross section [120]. The K factor ranges 1.0-1.2 and is 1.1 at m ZZ = 125 GeV. The NLO EW corrections are applied as a function of m ZZ according to the computation presented in Ref. [121].
The soft collinear approximation has been shown to describe accurately the cross section and the interference term for the gluon fusion ZZ production at NNLO in pQCD [122]. Additional calculations demonstrate that the K factors are very similar at NLO for signal and background [123] and at NNLO for the signal and interference terms [124]. Hence, the same K factor is used for the signal and the background [125]. The HNNLO v2 program [126][127][128] is used to obtain the signal NNLO K factor as a function of m ZZ from the ratio of the NNLO and LO gg → H → 2ℓ2ℓ ′ cross sections for the predicted SM H boson decay width of 4.07 MeV [118]. The NNLO/LO K factor for gg → ZZ varies from ≈2.0 to 2.6 and is 2.27 at m ZZ = 125 GeV, and a 10% systematic uncertainty is used when it is applied to the background.
The irreducible background contributions are included as binned templates in the likelihood function separately for the three considered final states (4e, 4µ, and 2e2µ). The templates are normalized to the most accurate theoretical calculations for the qq → ZZ → 4ℓ and gg → ZZ → 4ℓ cross sections [120][121][122][123][124][125]. A second method for the measurement of the inclusive fiducial cross section is presented in Section 10, where the normalization of these processes is treated as an unconstrained parameter in the fit to assess the constraint that can be derived from sidebands in data.

Reducible background
The reducible background contribution to the H boson signal in the 4ℓ channel mainly comes from the Z+jets, tt+jets, Zγ+jets, WW+jets, and WZ+jets production, hereafter collectively referred to as "Z+X" since the Z+jets contribution is the dominant one.
The contribution from the reducible background is estimated with the technique explained in Ref. [35]. The method is based on lepton misidentification rates, which are defined as the fraction of non-signal leptons that satisfy the selection criteria, computed in a control region in data that includes a Z boson and exactly one additional "loose" lepton (Z + ℓ), i.e., leptons with p T , η, and PV cuts but without identification nor isolation cuts applied. The lepton misidentification rates are then applied to another control region, comprised of a Z boson candidate and two opposite-sign or same-sign "loose" leptons (Z + ℓℓ), to reweigh the number of events to the signal region.
The distributions as functions of m 4ℓ of the Z+jets reducible background are derived for the three final states (4µ, 4e, and 2e2µ) separately and are included as binned templates in the likelihood function.

Measurement methodology
The pp → H → ZZ → 4ℓ fiducial cross sections are extracted from a maximum likelihood fit of the signal and background expected distributions to the observed 4ℓ mass distribution, N obs (m 4ℓ ), parametrized for each final state f , in each kinematic bin i of a given observable, and year of data taking y as: The The H resonant signal distribution is parametrized with a double-sided Crystal Ball (DCB) function [129] around m H = 125 GeV. The corresponding probability density function, P res (m 4ℓ ), is scaled by the fiducial cross section, σ fid , and the integrated luminosity L. The DCB function parameters are obtained from a simultaneous fit of the m 4ℓ distributions corresponding to the various mass points in the m H range of 105-160 GeV, which allows to express the dependency of the fitted parameters in m H directly in the fit, following the same strategy of Ref. [35].
A Landau distribution is introduced to empirically model the shape of the nonresonant signal contribution, P nonres (m 4ℓ ), for the WH, ZH, and ttH processes where one of the leptons from the H boson decay is either not selected or falls outside the acceptance. The fraction of such events in the mass range considered is about 5, 22, and 17% for WH, ZH, and ttH, respectively. These nonresonant events are treated as a background in the measurements. The reducible and irreducible backgrounds are included in the fit as normalized binned templates, P bkg , of the mass distribution of these processes.
An additional contribution ( f nonfid ) is introduced to take into account the presence of events not originating from the fiducial volume but satisfying the selections and is treated as background in the measurements. This contribution is referred to as the "nonfiducial signal" and is estimated from simulation for each signal model. The values of f nonfid are found to be consistent across the different observables considered, for the same production mechanism. Dedicated simulations have shown that the m 4ℓ distribution of these events is identical to that of the resonant fiducial signal. To minimize the model dependence of the measurement, the value of f nonfid is fixed to be a fraction of the fiducial signal component. The values of this fraction are reported in Table 7 and range between 4% for the VBF production mechanism and up to 18% for the ttH mode. The acceptance of the events originating from VH or ttH is lower than ggH and VBF events, reflecting the possible presence of leptons in the final states not produced by the H boson decay and resulting in larger values of f nonfid for these production mechanisms.
Generator-level observables used in the definition of the fiducial phase space are smeared by detector effects at reconstruction level. The ϵ f i,j response matrix is obtained from simulation, for each final state f , and is used to unfold the number of expected events in bin i at the reconstruction level to the number of expected events of a given observable in bin j at the fiducial level. For the measurement of the inclusive fiducial cross section, ϵ f i,j corresponds to a single number, the efficiency, listed in the second column of Table 7 for the various production mechanisms. The table also shows the acceptance A fid , defined as the fraction of signal events that fall within the fiducial phase space. Table 7: Summary of the inputs to the maximum likelihood based unfolding. The fraction of signal events within the fiducial phase space (acceptance A fid ), the reconstruction efficiency (ϵ) in the fiducial phase space, and the ratio of the number of reconstructed events outside the fiducial phase space to that of the ones inside the fiducial phase space ( f nonfid ) are quoted for each production mechanism for m H = 125.38 GeV. The last column shows the value of (1 + f nonfid )ϵ, which regulates the signal yield for a given fiducial cross section. All values are shown with their statistical uncertainty. The values for the ggH production mode are obtained using the POWHEG generator.

Signal process
A Systematic uncertainties are included in the form of nuisance parameters and the fiducial cross section measurements are obtained using an asymptotic approach [130] with a test statistic based on the profile likelihood ratio [131]. A maximum likelihood fit is performed simultaneously in all final states and bins of each observable, assuming m H = 125.38 GeV. The branching fractions of the H boson to the different final states (4e, 4µ, 2e2µ) are unconstrained parameters in the fit to increase the model independence of the measurements, following the strategy adopted in Ref. [35]. A likelihood-based unfolding is performed to resolve the detector effects from the observed distributions to the fiducial phase space. This approach is the same as in Refs. [35,132] and allows to simultaneously unfold detector effects and perform the fit to extract the fiducial cross section. The analysis strategy of Ref.
[35] is extended by measuring separately the fiducial cross sections in 4e + 4µ and 2e2µ final states for observables targeting the H → ZZ → 4ℓ decay. This choice is driven by the different physics in the final states containing different-and same-flavor leptons arising from the destructive interference between the two alternative methods of constructing the H → ZZ → 4ℓ diagrams in the same-helicity states in the case of identical leptons.

Systematic uncertainties
The integrated luminosities of the 2016, 2017, and 2018 data-taking periods are individually known with uncertainties in the range 1.2-2.5% [60][61][62], while the 2016-2018 integrated luminosity has an uncertainty of 1.6%. The partial correlation scheme considered for this systematic uncertainty is summarized in Table 8.
Experimental systematic uncertainties in trigger and lepton reconstruction and selection efficiencies are estimated from data for different final states. These uncertainties are derived from a tag-and-probe technique using J/ψ and Z decays into a pair of leptons and range 4.3-10.9% in the 4e channel and 0.6-1.9% in the 4µ channel, depending on the p T region. In this paper, a  new way to estimate the systematic uncertainty in the electron efficiency measurements is introduced. Alternative variations in the tag-and-probe fit used to derive the systematic uncertainty are combined and the RMS of the values is used for the systematic uncertainty evaluation. This makes the evaluation of the systematic uncertainty for the low-p T bins more solid, and leads to their reduction of a factor approximately 40% in the electron reconstruction and selection efficiency with respect to the values reported in Ref. [35].
The systematic uncertainties in the lepton momentum scale and resolution are estimated from dedicated studies of the Z → ℓ + ℓ − mass distribution in data and simulation. The momentum scale uncertainty is 0.06% in the 4e channel and 0.01% in the 4µ channel, while the resolution uncertainty is 10% in the 4e channel and 3% in the 4µ channel. The effect of these uncertainties is evaluated by allowing the corresponding parameters of the DCB function used to model the resonant signal to remain unconstrained in the fit.
Jet-related observables are affected by systematic uncertainties in the estimation of the jet energy scale. These uncertainties affect the normalization of the processes and are modeled with a set of nuisance parameters representing the various sources and accounting for the partial correlation among the various final states and years. Their values depend on the kinematic bin and range 0.1-33%, with an average value of 3%.
Furthermore, experimental systematic uncertainties in the reducible background estimation are considered. These uncertainties originate from the evaluation of the lepton misidentification rates and vary between 23 and 43%, depending on the final state. The impact of these uncertainties is negligible. Theoretical uncertainties in the renormalization and factorization scales, and in the choice of the PDF set affect both the signal and background rates. The scale uncertainty is determined by varying renormalization and factorization scales between 0.5 and 2 times their nominal value, while keeping their ratio between 0.5 and 2. The uncertainty in the PDF set is determined following the PDF4LHC recommendations taking the root mean square of the variation of the results when using different replicas of the default NNPDF set [133]. An additional 10% uncertainty in the K factor used for the gg → ZZ background prediction is applied. A systematic uncertainty of 2% [34] in the branching fraction of H → ZZ → 4ℓ is considered and only affects the signal yield. The theoretical uncertainties affecting the signal are not included in the fit but evaluated and indicated in Figs. 4-24.

Results
This section reports the measurements of the fiducial cross sections in differential bins of the kinematic observables introduced in Section 6. All fiducial cross section measurements are in agreement with the SM predictions within uncertainties. The compatibility of the results with the theoretical predictions is quantified by reporting the p-value for each observable, computed using the negative log-likelihood ratio as test statistic evaluated at the SM point. The p-value is calculated using a χ 2 probability density function with the number of bins used in the measurements taken as the number of degrees of freedom. The observed p-values range from 0.05 to 0.99. The inclusive cross section is measured with an overall precision of 10%, with statistical and systematic contributions of 8% and 6%, respectively. All the differential measurements are limited by the statistical uncertainty. The differential cross sections as functions of p H T and |y H | are measured with an average precision of 35%, while the most precise cross sections are measured with a precision of 20%. The measurements are compared to the theoretical predictions from various generators. The uncertainties in these predictions come from the uncertainty in the fiducial acceptance, the H → ZZ → 4ℓ branching fraction and variations of the PDF replicas, α s value, and renormalization and factorization scales. Figure 2 depicts the distributions of these two observables comparing data to predictions from simulation. With respect to Ref.
[35], the data set used in this analysis benefits from an improved object calibration. This leads to a better precision in the final results and permits measurements of jet-related observables in a phase space region that extends up to |η| = 4.7. The inclusive fiducial result features a reduction of 15% of the uncertainty, particularly evident in the 40% reduction of the systematic component obtained using a root-mean-square approach to compute the uncertainties in the electrons selection efficiency [134], which are the leading source of systematic uncertainty on the measurements performed in this analysis. Tabulated results are provided in the HEPData record for this analysis [135]. The kinematic properties of the four-lepton system are studied by measuring differential cross sections in bins of angular observables sensitive to the HVV decay. These results are reported for the inclusive four-lepton final state and for the same-flavor and different-flavor final states separately to account for interference effects in the case of identical helicity states. In all cases the results agree with the distributions predicted by the SM. The largest deviations with respect to the expected values are observed in the central bins of cos θ 2 and Φ and are compatible with statistical fluctuations in the observed data. The p-values of these two measurements are found to be 0.23 and 0.24, respectively, thus corroborating the compatibility with the SM predictions. For the first time, differential cross sections are measured in bins of kinematic discriminants sensitive to the presence of possible HVV anomalous couplings. These measurements are compared to the distributions of these discriminants computed under the SM hypothesis and under various anomalous coupling hypotheses. The former is always favored in the comparison with data.
A total of six double-differential measurements are also reported.

Inclusive cross section
The measured inclusive fiducial cross section for the H → ZZ → 4ℓ process is for a H boson mass of m H = 125.38 GeV, in agreement with the SM expectation of σ SM fid = 2.86 ± 0.15 fb. Figure 3 shows the corresponding log-likelihood scan. The systematic component of the uncertainty is dominated by electron-related nuisance parameters (electrons), especially the electron selection efficiency that is the leading nuisance parameter in the four-lepton decay channel. The muon-related nuisance parameters (muons) and the uncertainties on the luminosity measurement (lumi) and on the background predictions (bkg) play a minor role on the overall systematic uncertainty on σ fid . The inclusive fiducial cross section measured in the three final states (4e, 4µ, and 2e2µ) is shown in the left panel of Fig. 4, while the right panel depicts the evolution of the H → ZZ → 4ℓ fiducial cross section as a function of the centerof-mass energy. The results are compared with the cross sections predicted by the POWHEG, MADGRAPH5 aMC@NLO, and NNLOPS generators for the H boson production and parton showering, while the decay is always modeled by JHUGEN. The measurement of the inclusive fiducial cross section is repeated treating the normalization of the ZZ irreducible background processes as an unconstrained parameter in the fit. The results are presented in the left panel of Fig. 5 for the inclusive H → ZZ → 4ℓ measurement and the three final states considered. The correlation coefficient (ρ) between the inclusive fiducial cross section measurement and the ZZ normalization in the 4ℓ final state is found to be ρ = −0.03, while the correlations between the ZZ normalization in each final state and the corresponding fiducial cross section measurements are shown in the right panel of Fig. 5. The positive correlations observed between the ZZ background normalizations and the cross sections measured in the 4e and 2e2µ final states are driven by the systematic uncertainties on the electrons reconstruction and identification, which are the leading nuisance parameters in the analysis.
The results are summarized in Table 9. The measured cross sections, obtained with the ZZ normalization treated as an unconstrained parameter in the fit, are in agreement with the results obtained when the irreducible background normalization is constrained to the theoretical expectation. The uncertainty in this parameter when extracted from sidebands in data (7.5%) is larger than the theoretical uncertainty in its predictions (6.3%). For these reasons, in the following the ZZ normalization is fixed to the SM prediction.  In the left panel the acceptance and theoretical uncertainties are calculated using POWHEG (blue), NNLOPS (orange), and MAD-GRAPH5 aMC@NLO (pink). The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. In the right panel the acceptance is calculated using MINLOHJ at √ s = 13 TeV and NNLOPS [128,136] at √ s = 7 and 8 TeV.

Differential cross sections: production
Fiducial cross sections are measured in differential bins of observables sensitive to the H boson production. The results for the p H T and |y H | are shown in Fig. 6. Figure 7 shows the measurements of the fiducial cross sections in bins of the number of associated jets and of the p T of the leading and subleading jet in the event. The fiducial cross section is also measured in bins of the invariant mass and difference in pseudorapidity of the dijet system, as shown in Fig. 8. These measurements enhance the sensitivity to phase space regions where VBF and ttH production   mechanisms dominate and where a larger jet multiplicity is expected.
Cross sections in bins of observables of the H + j and H + jj systems are also measured. The results in differential bins of the invariant mass and p T of the H + j system are presented in Fig. 9 together with the results in differential bins of the p T of the H + jj system. Cross sections are also measured in differential bins of the rapidity-weighted jet vetoes introduced in Section 6, to enhance the sensitivity to phase space regions that probe directly the departure from LO kinematics and the QCD emission pattern. Figure 10  The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-toleading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions.  T is undefined. Lower: the fiducial cross section in the last bin is measured for events with p j2 T > 90 GeV and normalized to a bin width of 150 GeV. The first bin comprises all events with less than two jet, for which p j2 T is undefined. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions. , the difference in azimuthal angle ∆ϕ jj (upper right) the difference in pseudorapidity |∆η jj | (lower) of the dijet system. Upper Left: the fiducial cross section in the last bin is measured for events with m jj > 300 GeV and normalized to a bin width of 225 GeV. The first bin comprises all events with less than two jets, for which m jj is undefined. Upper right: the first bin comprises all events with less than two jet, for which |∆ϕ jj | is undefined. Lower: the first bin comprises all events with less than two jet, for which |∆η jj | is undefined. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions.  Figure 9: Upper left: differential cross sections as functions of the invariant mass of the H + j system m H j , where j is the leading jet in the event. The fiducial cross section in the last bin is measured for events with m H j > 600 GeV and normalized to a bin width of 280 GeV. The first bin comprises all events with less than one jet, for which m H j is undefined. Upper right: differential cross sections as functions of the transverse momentum of the H + j system p H j T . The fiducial cross section in the last bin is measured for events with p H j T > 110 GeV and normalized to a bin width of 90 GeV. The first bin comprises all events with less than one jet, for which p H j T is undefined. Lower: differential cross sections as functions of the transverse momentum of the H + jj system p H jj T . The fiducial cross section in the last bin is measured for events with p H jj T > 60 GeV and normalized to a bin width of 40 GeV. The first bin comprises all events with less than two jet, for which p H jj T is undefined. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions. , i.e., events with less than one jet, for which T max B is undefined, and events with T max B < 30 GeV. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions.

Differential cross sections: decay
In this section the measurements of fiducial cross sections in differential bins of observables sensitive to the H → ZZ → 4ℓ decay are presented. Since the final state is sensitive to interference effects in the case of identical particles, the results for decay observables are also presented separately for same-and different-flavor final states. This ensures a complete coverage of the whole phase space and a more model-independent set of results.
The cross sections measured in bins of the invariant mass of the two Z boson candidates are shown in Figs. 11 and 12. The additional degrees of freedom that characterize the H → ZZ → 4ℓ decay are the five angles introduced in Section 6. The cross sections in differential bins of the cosine of the θ angles are presented in Figs. 13, 14, and 15, respectively. Figures 16 and 17 show the results for the measurements in bins of the Φ and Φ 1 angles. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to nextto-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions.  Figure 13: Differential cross sections as functions of cos θ * in the 4ℓ (upper) and in the sameflavor (lower left) and different-flavor (lower right) final states. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions.  Figure 15: Differential cross sections as functions of cos θ 2 in the 4ℓ (upper) and in the sameflavor (lower left) and different-flavor (lower right) final states. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions. Ratio to NNLOPS Figure 16: Differential cross sections as functions of the Φ angle in the 4ℓ (upper) and in the same-flavor (lower left) and different-flavor (lower right) final states. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions. The acceptance and theoretical uncertainties in the differential bins are calculated using the ggH predictions from three different generators normalized to next-to-next-to-next-to-leading order (N 3 LO) [34]. The subdominant component of the signal (VBF + VH + ttH) is denoted as XH and is fixed to the SM prediction. The measured cross sections are compared with the ggH predictions from POWHEG (blue), NNLOPS (orange), and MADGRAPH5 aMC@NLO (pink). The hatched areas correspond to the systematic uncertainties in the theoretical predictions. Black points represent the measured fiducial cross sections in each bin, black error bars the total uncertainty in each measurement, red boxes the systematic uncertainties. The lower panels display the ratios of the measured cross sections and of the predictions from POWHEG and MADGRAPH5 aMC@NLO to the NNLOPS theoretical predictions.
These observables can be used to compute matrix element discriminants sensitive to the presence of possible BSM physics effects as described in Section 6.1 and Ref. [31].
Cross sections are measured in differential bins of six kinematic discriminants sensitive to various HVV anomalous couplings and the interference between two model contributions (SM and a BSM scenario): D dec 0-, D dec 0h+ , D dec CP , D dec int , D dec Λ1 , and D Zγ,dec Λ1 . The results of these measurements, shown in Figs. 18-23, are compared to distributions for the matrix element discriminants corresponding to various anomalous coupling hypotheses. Following the conventions adopted in Ref. [31], rather than using the value of the coupling to identify the type of the anomalous coupling sample, the cross sections fractions f ai are used: where σ i is the cross section for the process corresponding to a i = 1, a j̸ =i = 0 in Eq. (6). The term for Λ 1 isσ Λ1 /(Λ 1 ) 4 instead of |a i | 2 σ i , whereσ Λ1 is the effective cross section for the process corresponding to Λ 1 = 1 TeV, given in units of fb·TeV 4 . To study the a 2 and a 3 couplings, discriminants of the form D alt and D int are built. The former are compared to the distributions obtained for pure anomalous coupling scenarios corresponding to f a3 = 1 and f a2 = 1, while the latter are compared to the interference scenario where f a3 = 0.5 and f a2 = 0.5. A value of f ai = 0.5 corresponds to a maximal mixing between the SM and the BSM scenarios. To inspect the couplings κ 1 and κ Zγ 2 , the interference discriminant is not built since it does not provide additional information and the corresponding D alt can also be used to study the interference.

Double-differential cross sections
The differential cross section measurements presented so far ensure a good coverage of the production and decay phase spaces in the H → ZZ → 4ℓ channel, together with a separation of possible interference effects present in the same-and different-flavor final states. To improve the characterization of this decay channel and to maximize the coverage and separation of the different phase space regions, a set of double-differential measurements is also performed. The results are shown in Fig. 24 for |y H | vs. p         Table 6. The content of each plot is described in the caption of     Table 6. The content of each plot is described in the caption of Fig. 6.

Constraints on the H boson self-coupling
The differential cross section for the H boson production as a function of p H T can be used to extract limits on the H boson self-coupling, following the approach described in Refs. [137][138][139]. At NLO in pQCD the H boson production includes processes sensitive to the trilinear self-coupling (λ 3 ). The production modes ttH and VH introduce sizeable contributions to the H boson self-coupling due to the large vector boson and the top quark masses, whereas ggH and VBF production lead to much smaller contributions to the loop correction and are therefore less sensitive to possible modifications of λ 3 .
The cross sections for the various production mechanisms of the H boson are parametrized as functions of a coupling modifier κ λ = λ 3 /λ SM 3 in order to account for NLO terms arising from the H boson trilinear self-coupling. The signal model defined in Section 8 is modified by fixing the cross sections and branching fractions to their SM expectation values and by introducing scaling functions µ i,j (κ λ ) in each bin i of p H T , for each production mechanism j. The dominant production mechanism is ggH, for which a differential parametrization of the cross section as a function of κ λ is not available yet, as discussed in Refs. [137][138][139]. The inclusive value is used for the parametrization of the H boson cross section for this production mechanism, taking into account an inclusive O(λ 3 ) correction factor.
In order to compute the scaling functions µ i,j (κ λ ) for the other production modes, LO partonlevel events are generated using MADGRAPH5 aMC@NLO 2.5.5 and are reweighted on an eventby-event basis using a dedicated EW reweighting tool, which computes the corresponding NLO λ 3 -corrections (O(λ 3 )). The ratio of the O(λ 3 ) to the LO distributions in bins of p H T is used to derive the scaling functions µ i,j (κ λ ) as detailed in Ref. [138].
Constraints on κ λ are extracted from the maximum likelihood scan in the range −10 < κ λ < 20, outside which the model is no longer valid as NLO effects start to dominate, while the other H couplings are fixed to their SM value. The likelihood scan as a function of κ λ is shown in Fig. 26. The minimum of the negative log-likelihood ratio corresponds to a measured value of: for an expected value of: The corresponding observed (expected) excluded κ λ range at the 95% confidence level (CL) is: −5.4 (−7.6) < κ λ < 14.9 (17.7).
The current best available constraints on κ λ are obtained from the combination of measurements of H boson pair production performed with the Run 2 data-set. The limits set by the ATLAS and CMS Collaborations correspond to observed limits at the 95 % CL of −0.6 < κ λ < 6.6 [140] and −1.24 < κ λ < 6.49 [44], respectively.

Constraints on the charm and bottom quark couplings
The p H T differential cross section of the ggH production mechanism is used to set constraints on the H boson coupling modifiers to b and c quarks in the context of the κ-framework [141]. In fact, because of the presence of b and c quarks in the ggH loop [142], anomalous values of these couplings can result in modifications of the ggH cross section. The other production mechanisms are set to their SM expectation and no dependence from these couplings is assumed. The H boson coupling to the top quark is fixed to the SM value. The effects of the associated production with b quarks, whose contribution increases with increasing values of the coupling of the H boson to the b quark, are taken into account in the theoretical inputs to compute the parametrization.
The results are extracted from a maximum likelihood fit where the approach described in Section 8 is modified by separating the ggH production from the other mechanisms, which are considered as background and constrained to the SM predictions with their respective uncertainties. The combined effect of the H boson couplings to b (κ b ) and c quarks (κ c ) is modeled independently in each bin of the p H T spectrum by means of a quadratic polynomial, following the strategy of Ref. [52]. A simultaneous constraint on κ b and κ c is also derived by treating the H → ZZ branching fraction as an unconstrained parameter in the fit. The constraint from the total width and the overall normalization is removed in this way, and what remains is purely the constraint obtained from the shape of the p H T spectrum. The result is shown in Fig. 27 (right). Confidence intervals on κ b and κ c are obtained from a maximum likelihood fit leaving one of the two parameters unconstrained in the fit and scanning the other. The observed (expected) exclusion limits at the 95% CL are: assuming a dependence of the branching fraction on κ b and κ c , and: −5.6 (−5.5) < κ b < 8.9 (7.4) −20 (−19) < κ c < 23 (20), treating the branching fraction as an unconstrained parameter in the fit.

Summary
This paper presents a comprehensive characterization of the H → ZZ → 4ℓ decay channel via the measurement of fiducial differential cross sections as functions of several kinematic observables. The H boson production is characterized via measurements of differential cross sections in bins of p H T and |y H |, the p T of the leading and subleading jets and observables of the dijet system, when associated with jets. Fiducial cross sections are measured in bins of the seven kinematic observables that completely define the four-lepton decay: the invariant mass of the two Z bosons and the five angles that describe the fermions kinematical properties and the production and decay planes. Differential cross sections are also measured in bins of six matrix element kinematic discriminants sensitive to various anomalous couplings of the H boson to vector bosons. The dynamical evolution of the renormalization and factorization scales, and resummation effects are probed by measuring cross sections in bins of rapidity-weighted jet vetoes, and in bins of observables of the H plus jets system. An extensive set of double-differential measurements is presented, providing a complete coverage of the phase space under study. The H → ZZ → 4ℓ inclusive fiducial cross section is σ fid = 2.73 ± 0.26 fb = 2.73 ± 0.22 (stat) ± 0.15 (syst) fb, in agreement with the SM expectation of 2.86 ± 0.15 fb. The measurement of the fiducial cross section in differential bins of p H T is used to set constraints on the trilinear self-coupling of the H boson, with an observed (expected) limit of −5.4 (−7.6) < κ λ < 14.9 (17.7) at the 95% CL. Finally, constraints on the modifiers of H boson couplings to b and c quarks (κ b and κ c ) are also determined with an observed (expected) limit of −1.1 (−1.3) < κ b < 1.1 (1.2) and −5.3 (−5.7) < κ c < 5.2 (5.7) at the 95% CL, obtained assuming a dependence of the branching fraction on κ b and κ c . All results are consistent with the SM predictions for the H → ZZ → 4ℓ decay channel in the considered fiducial phase space.
[36] ATLAS Collaboration, "Measurement of the properties of Higgs boson production at √ s = 13 TeV in the H → γγ channel using 139 fb −1 of pp collision data with the ATLAS experiment", 2022. arXiv:2207.00348. Submitted to JHEP.
[39] CMS Collaboration, "Measurements of the Higgs boson production cross section and couplings in the W boson pair decay channel in proton-proton collisions at √ s = 13 TeV", 2022. arXiv:2206.09466. Submitted to Eur. Phys. J. C.
[46] CMS Collaboration, "Measurement of the Higgs boson inclusive and differential fiducial production cross sections in the diphoton decay channel with pp collisions at √ s = 13 TeV", 2022. arXiv:2208.12279. Submitted to JHEP.
[48] ATLAS Collaboration, "Fiducial and differential cross-section measurements for the vector-boson-fusion production of the Higgs boson in the H → WW * → eνµν decay channel at 13 TeV with the ATLAS detector", 2023. arXiv:2304.03053. Submitted to Phys. Rev. D.