Higgs boson production cross-section measurements and their EFT interpretation in the $4\ell$ decay channel at $\sqrt{s}=$13 TeV with the ATLAS detector

Abstract

Higgs boson properties are studied in the four-lepton decay channel (where lepton = e, $\mu$) using 139 $\hbox {fb}^{-1}$ of proton–proton collision data recorded at $\sqrt{s}=$13 TeV by the ATLAS experiment at the Large Hadron Collider. The inclusive cross-section times branching ratio for $H\rightarrow ZZ^*$ decay is measured to be $1.34 \pm 0.12$ pb for a Higgs boson with absolute rapidity below 2.5, in good agreement with the Standard Model prediction of $1.33 \pm 0.08$ pb. Cross-sections times branching ratio are measured for the main Higgs boson production modes in several exclusive phase-space regions. The measurements are interpreted in terms of coupling modifiers and of the tensor structure of Higgs boson interactions using an effective field theory approach. Exclusion limits are set on the CP-even and CP-odd ‘beyond the Standard Model’ couplings of the Higgs boson to vector bosons, gluons and top quarks.

1 Introduction

The observation of the Higgs boson by the ATLAS and CMS experiments [1, 2] with the Large Hadron Collider (LHC) Run 1 data set at centre-of-mass energies of $\sqrt{s}=$ 7 $\text {Te} \text {V}$ and 8 $\text {Te} \text {V}$was a major step towards an understanding of the electroweak (EW) symmetry breaking mechanism [3,4,5]. Tests of its spin and CP quantum numbers strongly indicate that the observed particle is of scalar nature and that the dominant coupling structure is CP-even, consistent with the Standard Model (SM) expectation [6,7,8]. The measurements of the Higgs boson production and differential cross-sections, branching ratios, and the derived constraints on coupling-strength modifiers, assuming the SM coupling structure, have also shown no significant deviation from the predictions for the SM Higgs boson with a mass of 125 $\text {Ge} \text {V}$ [9,10,11,12]. Furthermore, constraints have been set on various coupling parameters beyond the SM (BSM) that modify the tensor structure of the Higgs boson couplings to SM particles [8, 13,14,15,16,17,18,19,20].

Motivated by a clear Higgs boson signature and a high signal-to-background ratio in the $H \rightarrow ZZ^{*} \rightarrow 4\ell$ decay channel (where $\ell = e\text { or }\mu$), the updated measurements of the Higgs boson coupling properties in this channel are presented using the entire Run 2 data set with $139$ fb$^{-1}$of proton–proton (pp) collision data collected at $\sqrt{s}=13$ $\text {Te} \text {V}$ by the ATLAS detector between 2015 and 2018. Three types of results are presented in this paper: (i) measurements of the Higgs boson production cross-sections times branching ratio, hereafter referred to as cross-sections, for the main production modes in several exclusive phase-space bins in dedicated fiducial regions; (ii) interpretation of the measurements in terms of constraints on the Higgs boson coupling-strength modifiers within the $\kappa$-framework [21]; and (iii) interpretation of the measurements in terms of modifications to the tensor structure of Higgs boson couplings using an effective field theory (EFT) approach.

In addition to a nearly four times higher integrated luminosity, there are several other important differences compared to the previous results in this analysis channel [17]:

• an improved lepton isolation to mitigate the impact of additional pp interactions in the same or neighbouring bunch crossings (pile-up),

• an improved jet reconstruction using a particle flow algorithm [22],

• additional event categories for the classification of Higgs boson candidates,

• new discriminants to enhance the sensitivity to distinguish the various production modes of the SM Higgs boson,

• the use of data sidebands to constrain the dominant $ZZ^*$ background process,

• a dedicated control region to constrain the background in the reconstructed event categories probing $ttH$ production,

• improved estimates of $Z+\mathrm {jets}$, $tt$, and WZ backgrounds, and

• an EFT interpretation, based on a parameterisation of the cross-sections rather than a direct parameterisation of the reconstructed event yields.

1.1 Simplified template cross-sections

In the framework of Simplified Template Cross Sections (STXS) [23,24,25], exclusive regions of phase space are defined for each Higgs boson production mechanism. These phase-space regions, referred to as production bins, are defined to reduce the dependence on theoretical uncertainties that directly fold into the measurements and at the same time maximise the experimental sensitivity to measure the bins, enhance the contribution from possible BSM effects, and allow measurements from different Higgs boson decay modes to be combined. The number of production bins is limited to avoid loss of measurement sensitivity for a given amount of integrated luminosity.

The definitions of the production bins used for this measurement are shown in the left panel of Fig. 1 (shaded area). All production bins are defined for Higgs bosons with rapidity $|y_H|<$ 2.5 and no requirement is placed on the particle-level leptons. Two sets of production bins with different granularity are considered, as a trade-off between statistical and theoretical uncertainties.

Fig. 1

Two sets (Production Mode Stage and Reduced Stage 1.1) of exclusive phase-space regions (production bins) defined at particle-level for the measurement of the Higgs boson production cross-sections (left and middle-left shaded panels), and the corresponding reconstructed event categories for signal (middle-right panel) and sidebands (right panel). The description of the production bins is given in Sect. 1.1, while the reconstructed signal region and sideband event categories are described in Sects. 5 and 6, respectively. The $bbH$ (tH) contribution is included in the $\mathrm {ggF}$ ($ttH$) production bins. The colours of each reconstructed event category box indicates the contributions from the relevant production processes

The first set of production bins (Production Mode Stage) [24] is defined according to the Higgs boson production modes: gluon–gluon fusion ($\mathrm {ggF}$), vector-boson fusion ($\mathrm {VBF}$) and associated production with vector bosons (VH, where $V = W$ or Z) or top quark pairs ($ttH$). Since b-jets from $bbH$ associated production are emitted at small angles relative to the beam axis and usually outside of the detector acceptance, the $bbH$ and $\mathrm {ggF}$ Higgs boson production modes have similar signatures and acceptances. Their contributions are considered together with their relative ratio fixed to the SM prediction. In the following, the sum of their contributions is referred to as ggF. Similarly, single top production (tH) is considered together with $ttH$, with their relative ratio fixed to the SM prediction. In contrast to the Stage-0 production bins described in Ref. [24], the VH events with hadronic decays of the vector boson V are included in the VH production bin rather than in the $\mathrm {ggF}$ or $\mathrm {VBF}$ bins. In this way, each of the four main Higgs boson production modes can be measured separately.

The second set of production bins (Reduced Stage 1.1) is more exclusive than the first one. Starting from the production bins of a more granular Stage 1.1 set [25], several production bins are merged as the full set of bins cannot be measured separately in the $H\rightarrow ZZ^* \rightarrow 4\ell$ channel with the current data sample. The definitions of the bins are based on the multiplicity of particle-level jets, the Higgs boson transverse momentum $p_{\text {T}} ^H$ and the invariant mass $m_{jj}$ of the two jets with the highest transverse momentum. Particle-level jets are built from all stable particles (particles with lifetime c${\tau }{>}10$ mm) including neutrinos, photons, and leptons from hadron decays or those produced in the parton shower. The anti-$k_t$ jet reconstruction algorithm [26, 27] with a radius parameter $R=0.4$ is used. All Higgs boson decay products, as well as the leptons and neutrinos from the decays of the associated V bosons are excluded from the jet building, while the decay products from hadronically decaying associated V bosons, are included. The jets are required to have $p_{\text {T}}>$ 30 $\text {Ge} \text {V}$, with no restrictions on rapidity.

Events from $\mathrm {ggF}$ production and $gg\rightarrow ZH$ production with a hadronically decaying Z boson are split into seven common production bins. Six bins have a Higgs boson transverse momentum below 200 $\text {Ge} \text {V}$, while the seventh bin with Higgs boson transverse momentum above 200 $\text {Ge} \text {V}$ ($\mathrm {gg2H}$-$p_{\mathrm {T}}^H$-$\mathrm {High}$) is sensitive to contributions from BSM physics. For $p_{\text {T}} ^H$ below 200 $\text {Ge} \text {V}$, further splits are made according to the jet multiplicity and $p_{\text {T}} ^H$. Events with no jets are split into two bins with $p_{\text {T}} ^H$ below and above 10 $\text {Ge} \text {V}$. Events with one jet are split into three bins with $p_{\text {T}} ^H$ below 60 $\text {Ge} \text {V}$, between 60 and 120 $\text {Ge} \text {V}$, and above 120 $\text {Ge} \text {V}$. Finally, Higgs boson events with two or more jets are combined into one bin. The bins are respectively denoted by $\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {Low}$, $\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {High}$, $\mathrm {gg2H}$-$1j$-$p_{\mathrm {T}}^H$-$\mathrm {Low}$, $\mathrm {gg2H}$-$1j$-$p_{\mathrm {T}}^H$-$\mathrm {Med}$, $\mathrm {gg2H}$-$1j$-$p_{\mathrm {T}}^H$-$\mathrm {High}$ and $\mathrm {gg2H}$-$2j$.

As described in Ref. [25], $\mathrm {VBF}$ and VH production with hadronically decaying associated V bosons represent the t-channel and s-channel contributions to the same electroweak qqH production process and are therefore considered together for further splitting. Three bins are defined: one bin, sensitive to BSM contributions ($\mathrm {qq2Hqq}$-$\mathrm {BSM}$), with $p_{\text {T}} ^H$ above 200 $\text {Ge} \text {V}$ and $m_{jj}$ above 350 $\text {Ge} \text {V}$; one bin ($\mathrm {qq2Hqq}$-$\textit{VH}$) with $m_{jj}$ between 60 and 120 $\text {Ge} \text {V}$ to target the VH production mode; and one bin ($\mathrm {qq2Hqq}$-$\mathrm {VBF}$) with the Higgs boson not satisfying these criteria to ensure sensitivity to the $\mathrm {VBF}$ process. qqH events in which one or both jets have transverse momenta below the 30 $\text {Ge} \text {V}$ threshold are treated as a part of the $\mathrm {qq2Hqq}$-$\mathrm {VBF}$ bin.

The VH process with the associated V boson decaying leptonically is considered separately (VH-Lep). The leptonic decay includes the decays into $\tau$-leptons and neutrino pairs. The $ttH$ production bin remains the same as in the Production Mode Stage.

The middle-right and right panels of Fig. 1 summarise the corresponding categories of reconstructed events in which the cross-section measurements and background estimations are performed. These are described in detail in Sect. 5.

1.2 Higgs boson couplings in the $\kappa$-framework

To probe physics beyond the SM, the measured production cross-sections are interpreted within a leading-order-motivated $\kappa$-framework [21], in which a set of coupling modifiers $\vec {\kappa }$ is introduced to parameterise deviations from the SM predictions of the Higgs boson couplings to SM bosons and fermions. The framework assumes that the data originate from a single CP-even Higgs boson state with a mass of 125 $\text {Ge} \text {V}$ and the tensor coupling structure of the SM for its interactions. Only the coupling strengths are allowed to be modified by the BSM processes. The Higgs boson width is assumed to be small enough such that the narrow-width approximation is valid, allowing the Higgs boson production and decay to be factorised:

$$\begin{aligned} \sigma \cdot {\mathcal {B}} \ (i\rightarrow H \rightarrow f) = \sigma _{i}(\vec {\kappa })\cdot \frac{\Gamma _{f}(\vec {\kappa })}{ \Gamma _{H}(\vec {\kappa })}, \end{aligned}$$

where $\sigma _{i}$ is the production cross-section via the initial state i, ${\mathcal {B}}$ and $\Gamma _{f}$ are the branching ratio and partial decay width for the decay into the final state f, respectively, and $\Gamma _{H}$ is the total width of the Higgs boson. For a Higgs boson production and decay process via couplings i and f, respectively, coupling-strength modifiers are defined as

$$\begin{aligned} \kappa ^2_i = \frac{\sigma _i}{\sigma ^{\text {SM}}_i} \ \ \ \text {and} \ \ \ \kappa ^2_f = \frac{\Gamma _f}{\Gamma ^{\text {SM}}_f}, \end{aligned}$$

so that

$$\begin{aligned} \sigma \cdot {\mathcal {B}} \ (i\rightarrow H \rightarrow f) = \kappa _i^2 \cdot \kappa _f^2 \cdot \sigma _i^{\text {SM}} \cdot \frac{\Gamma _f^{\text {SM}}}{\Gamma _H(\kappa _i^2,\kappa _f^2)}. \end{aligned}$$

1.3 Tensor structure of Higgs boson couplings in the effective field theory approach

The $\kappa$-framework assumes that the tensor structure of the Higgs boson couplings is the same as in the SM. In order to probe for possible non-SM contributions to the tensor structure of the Higgs boson couplings, the measured simplified template cross-sections are interpreted using an EFT approach. In this approach, which exploits exclusive kinematical regions of the Higgs boson production and decay phase space, the BSM interactions are introduced via additional higher-dimensional operators ${\mathcal {O}}_i^{(d)}$ of dimension d, supplementing the SM Lagrangian ${{\mathcal {L}}}_{\text {SM}}$,

$$\begin{aligned} {\mathcal {L}}_{\text {EFT}} = {\mathcal {L}}_{\text {SM}} + \sum _{i}^{} \frac{C_i^{(d)}}{\Lambda ^{(d-4)}} {\mathcal {O}}_i^{(d)} \ \ \ \text {for} \ d>4. \end{aligned}$$

The parameters $C^{(d)}_i$ specify the strength of new interactions and are known as the Wilson coefficients, and $\Lambda$ is the scale of new physics. Only dimension-six operators are considered for this paper, since the dimension-five and dimension-seven operators violate lepton and baryon number conservation and the impact of higher-dimensional operators is expected to be suppressed by more powers of the cutoff scale $\Lambda$ [28]. For energies less than the scale of new physics, only the ratio $c_i = C_i^{(d=6)}/\Lambda ^2$ can be constrained by the data.

Constraints are set on the Wilson coefficients defined within the Standard Model Effective Field Theory (SMEFT) formalism [29] in the Warsaw basis [30]. The measurements in the $H\rightarrow ZZ^* \rightarrow 4\ell$ channel do not provide sensitivity for simultaneous constraints on the full set of these coefficients. To reduce the number of relevant parameters, a minimal flavour-violating scenario is assumed and only operators affecting the Higgs boson cross-section at tree level are considered. Operators affecting only double Higgs boson production and those affecting the Higgs boson couplings to down-type quarks and leptons are neglected due to limited sensitivity. The impact of these operators on the total Higgs boson decay width is also neglected.

Table 1 Summary of EFT operators in the SMEFT formalism that are probed in the $H\rightarrow ZZ^* \rightarrow 4\ell$ channel. The corresponding tensor structure in terms of the SM fields from Ref. [29] is shown together with the associated Wilson coefficients, the affected production vertices and the impact on the $H\rightarrow ZZ^*$ decay vertex. The Higgs doublet field and its complex conjugate are denoted as H and ${\widetilde{H}}$, respectively. The left-handed quark doublets of flavour p (the right-handed up-type quarks) are denoted $q_p$ ($u_r$). $V_{\mu \nu }$ (${\widetilde{V}}_{\mu \nu }=\epsilon ^{\mu \nu \rho \sigma }V_{\rho \sigma }$) is the (dual) field strength tensor for a given gauge field $V = G, W, B$. The bosonic operators with (without) a dual field strength tensor are CP-odd (CP-even). For the remaining operator with fermions (${\mathcal {O}}_{uH}$), the CP-odd contribution is introduced through the non-vanishing imaginary part of the corresponding Wilson coefficient, denoted as $c_{{\widetilde{u}}H}$

The remaining ten operators (see Table 1) comprise five CP-even and five CP-odd ones. The CP-even operators describing interactions between the Higgs boson and gluons and the top-Yukawa interactions are associated with the Wilson coefficients $c_{HG}$ and $c_{uH}$ from Ref. [29], respectively. Similarly, the CP-even Higgs boson interactions with vector bosons are related to $c_{HW}$, $c_{HB}$, and $c_{HWB}$ that impact the $\mathrm {VBF}$ and VH production and the Higgs boson decay into Z bosons. The Wilson coefficients for the corresponding CP-odd operators are $c_{{\widetilde{u}}H}$, $c_{H{\widetilde{G}}}$, $c_{H{\widetilde{W}}}$, $c_{H{\widetilde{B}}}$ and $c_{H{\widetilde{W}}B}$.

The constraints on the Wilson coefficients can be derived by comparing the expected with the measured simplified template cross-sections. For that purpose, the corresponding expected signal production cross-sections, the branching ratio and the signal acceptances are parameterised in terms of the Wilson coefficients. The dependence of signal production cross-sections on the EFT parameters can be obtained from its separation into three components:

$$\begin{aligned} \sigma \propto \left| {\mathcal {M}}_{\text {SMEFT}} \right| ^2= & {} \left| {\mathcal {M}}_{\text {SM}} + \sum _i \frac{C_i}{\Lambda ^2} {\mathcal {M}}_{i} \right| ^2 \\= & {} \left| {\mathcal {M}}_\text {SM} \right| ^2 + \sum _i 2Re \left( {\mathcal {M}}_\text {SM}^* {\mathcal {M}}_i \right) \frac{C_i}{\Lambda ^2} \\&+ \sum _{ij} 2Re \left( {\mathcal {M}}_i^* {\mathcal {M}}_j \right) \frac{C_i C_j}{\Lambda ^4}, \end{aligned}$$

where the first term on the right-hand side is the squared matrix element for the SM, the second term represents the interference between the SM and dimension-six EFT amplitudes and the third term comprises the pure BSM contribution from dimension-six EFT operators alone. Following this expression, the dependence of the Higgs boson cross-section $\sigma ^p (\vec {c})$ in a given production bin p on a set of Wilson coefficients $\vec {c}$ is parameterised relative to the SM prediction $\sigma ^p_{\text {SM}}$ as

$$\begin{aligned} \frac{\sigma ^{p} (\vec {c})}{\sigma ^{p}_\text {SM}} = 1 + \sum _i A^{p}_i c_i + \sum _{ij} B^{p}_{ij} c_i c_j, \end{aligned}$$

(1)

where the coefficients $A^p_i$ and $B^p_{ij}$ are independent of $\vec {c}$ and are determined from simulation. A similar procedure is applied to obtain from simulation the EFT parameterisation of the branching ratio $\mathcal {B} ^{4\ell }$ for the $H\rightarrow ZZ^*\rightarrow 4\ell$ decay from the partial $(\Gamma ^{4\ell })$ and total decay width $(\Gamma ^{\text {tot}})$ parameterisations,

$$\begin{aligned} \mathcal {B} ^{4\ell } (\vec {c})= & {} \frac{\Gamma ^{4\ell } (\vec {c})}{\Gamma ^{\text {tot}} (\vec {c})} \nonumber \\= & {} \mathcal {B} ^{4\ell }_{\text {SM}} \cdot \frac{1 + \sum _i A^{4\ell }_i c_i + \sum _{ij} B^{4\ell }_{ij} c_i c_j}{1 + \sum _f \left( \sum _i A^{f}_i c_i + \sum _{ij} B^{f}_{ij} c_i c_j \right) }, \end{aligned}$$

(2)

where the total decay width is the sum of all partial decay widths $\Gamma ^f$ related to the decay mode f. The procedure for the parameterisation of the cross-sections and the branching ratios is described in more detail in Ref. [31]. The criteria employed in the selection of four-lepton candidates introduce an additional dependence of the signal acceptance on the EFT parameters. This is taken into account in the interpretation, as discussed in Sect. 10.

2 ATLAS detector

The ATLAS detector [32,33,34] at the LHC is a multipurpose particle detector with a forward–backward symmetric cylindrical geometry¹ and a nearly $4\pi$ coverage in solid angle. It consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid, which provides a ${2}\,\hbox {T}$ axial magnetic field, electromagnetic (EM) and hadron calorimeters, and a muon spectrometer (MS). The inner tracking detector covers the pseudorapidity range $|\eta | < 2.5$. It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors. A lead/liquid-argon (LAr) sampling calorimeter provides electromagnetic energy measurements in the pseudorapidity range $|\eta | < 3.2$ with high granularity. A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range ($|\eta | < 1.7$). The endcap and forward regions are instrumented up to $|\eta | = 4.9$ with LAr calorimeters for both the EM and hadronic energy measurements. The calorimeters are surrounded by the MS and three large air-core toroidal superconducting magnets with eight coils each. The field integral of the toroid magnets ranges between 2.0 and 6.0 Tm across most of the detector. The MS includes a system of precision tracking chambers and fast detectors for triggering, covering the region $|\eta | < 2.7$. Events are selected using a first-level trigger implemented in custom electronics, which reduces the event rate to a maximum of 100 kHz using a subset of detector information. Software algorithms with access to the full detector information are then used in the high-level trigger to yield a recorded event rate of about 1 kHz [35].

3 Data set and event simulation

The full ATLAS Run 2 data set, consisting of pp collision data at $\sqrt{s}$ = 13 $\text {Te} \text {V}$ taken between 2015 and 2018, is used for this analysis. The total integrated luminosity after imposing data quality requirements [36] is 139 fb$^{-1}$.

The production of the SM Higgs boson via gluon–gluon fusion, via vector-boson fusion, with an associated vector boson and with a top quark pair was modelled with the Powheg-Box v2 Monte Carlo (MC) event generator [37,38,39]. For $\mathrm {ggF}$, the PDF4LHC next-to-next-to-leading-order (NNLO) set of parton distribution functions (PDF) was used, while for all other production modes, the PDF4LHC next-to-leading-order (NLO) set was used [40].

The simulation of $\mathrm {ggF}$ Higgs boson production used the Powheg method for merging the NLO Higgs boson + jet cross-section with the parton shower and the multi-scale improved NLO (MINLO) method [41,42,43,44] to simultaneously achieve NLO accuracy for the inclusive Higgs boson production. In a second step, a reweighting procedure (NNLOPS) [45, 46], exploiting the Higgs boson rapidity distribution, was applied using the HNNLO program [47, 48] to achieve NNLO accuracy in the strong coupling constant $\alpha _{\text {S}}$. The transverse momentum spectrum of the Higgs boson obtained with this sample is compatible with the fixed-order calculation from HNNLO and the resummed calculation at next-to-next-to-leading-logarithm accuracy matched to NNLO fixed-order with Hres2.3 [49, 50].

The matrix elements of the $\mathrm {VBF}$, $qq\rightarrow VH$, and $ttH$ production mechanisms were calculated up to NLO in QCD. For VH production, the MINLO method was used to merge 0-jet and 1-jet events [41, 43, 51,52,53,54]. The $gg\rightarrow ZH$ contribution was modelled at leading order (LO) in QCD.

The production of a Higgs boson in association with a bottom quark pair ($bbH$) was simulated at NLO with MadGraph5_aMC@NLO v2.3.3 [55, 56], using the CT10 NLO PDF [57]. The production in association with a single top quark (tH+X where X is either jb or W, defined in the following as tH) [58, 59] was simulated at NLO with MadGraph5_aMC@NLO v2.6.0 using the NNPDF3.0nlo PDF set [60].

For all production mechanisms, the Pythia 8 [61] generator was used for the $H\rightarrow ZZ^{*}\rightarrow 4\ell$ decay with $\ell = (e, \mu )$ as well as for parton showering, hadronisation and the underlying event. The contribution of the $Z\rightarrow \tau \tau$ decays is shown to have a negligible impact on the final result. The event generator was interfaced to EvtGen v1.2.0 [62] for simulation of the bottom and charm hadron decays. For the $\mathrm {ggF}$, $\mathrm {VBF}$ and VH processes, the AZNLO [63] set of tuned parameters was used, while the A14 [64] set was used for $ttH$, $bbH$ and tH processes. All signal samples were simulated for a Higgs boson mass $m_H=$ 125 $\text {Ge} \text {V}$.

For additional cross-checks, the $\mathrm {ggF}$ sample was also generated with MadGraph5_aMC@NLO. This simulation is accurate at NLO QCD accuracy for zero, one and two additional partons merged with the FxFx merging scheme [55, 65]. The events were showered using the Pythia 8 generator with the A14 set of tuned parameters.

The Higgs boson production cross-sections and decay branching ratios, as well as their uncertainties, are taken from Refs. [21, 24, 60, 66,67,68,69,70,71]. The $\mathrm {ggF}$ production is calculated with next-to-next-to-next-to-leading order ($\hbox {N}^{3}$LO) accuracy in QCD and has NLO electroweak (EW) corrections applied [72,73,74,75,76,77,78,79,80,81,82]. For $\mathrm {VBF}$ production, full NLO QCD and EW calculations are used with approximate NNLO QCD corrections [83,84,85]. The qq- and qg-initiated VH production is calculated at NNLO in QCD and NLO EW corrections are applied [86,87,88,89,90,91,92,93,94], while gg-initiated VH production is calculated at NLO in QCD. The $ttH$ [95,96,97,98], $bbH$ [99,100,101] and tH [58, 59] processes are calculated to NLO accuracy in QCD. The total branching ratio is calculated in the SM for the $H\rightarrow ZZ^{*}\rightarrow 4\ell$ decay with $m_H=$ 125 $\text {Ge} \text {V}$ and $\ell$ = (e, $\mu$) using PROPHECY4F [102, 103], which includes the complete NLO EW corrections, and the interference effects between identical final-state fermions. Due to the latter, the expected branching ratios of the 4e and 4$\mu$ final states are about 10% higher than the branching ratios to $2e2\mu$ and $2\mu 2e$ final states. Table 2 summarises the predicted SM production cross-sections and branching ratios for the $H\rightarrow ZZ^{*}\rightarrow 4\ell$ decay for $m_H=125$ $\text {Ge} \text {V}$.

Table 2 The predicted SM Higgs boson production cross-sections ($\sigma$) for $\mathrm {ggF}$, $\mathrm {VBF}$ and five associated production modes in pp collisions for $m_H=125$ $\text {Ge} \text {V}$ at $\sqrt{\mathrm {s}}=13~\text {Te} \text {V}$ [21, 24, 58,59,60, 66,67,68, 68,69,105]. The quoted uncertainties correspond to the total theoretical systematic uncertainties calculated by adding in quadrature the uncertainties due to missing higher-order corrections and PDF$+\alpha _{\text {S}}$. The decay branching ratios ($\mathcal {B}$) with the associated uncertainty for $H\rightarrow ZZ^{*}$ and $H\rightarrow ZZ^{*}\rightarrow 4\ell$, with $\ell = e, \mu$, are also given

For the study of the tensor structure of Higgs boson couplings within an effective field theory approach, several samples with different values of EFT parameters were simulated at LO in QCD separately for the $\mathrm {ggF} +bbH$, $\mathrm {VBF} + V(\rightarrow qq)H$, $qq\rightarrow Z(\rightarrow \ell \ell )H$, $qq\rightarrow W(\rightarrow \ell \nu )H$, $ttH$, tHW and tHjb production modes using MadGraph5_aMC@NLO and the NNPDF23lo PDF. The BSM signal is defined by the flavour symmetric SMEFTsim_A_U35_MwScheme_UFO_v2.1 model [29, 106], which incorporates the SMEFT dimension-six operators in the standard Universal FeynRules Output format created using the FeynRules framework [107, 108]. The light quarks (u, d, s and c) and leptons are assumed to be massless in the model. The generated events were showered with Pythia 8, using the CKKW-L matching scheme to match matrix element and parton shower computations with different jet multiplicities [61]. The A14 set of tuned parameters was used. All processes were simulated in the four-flavour scheme, apart from the tHW production, for which the five-flavour scheme was used [55].

The $ZZ^{*}$ continuum background from quark–antiquark annihilation was modelled using Sherpa v2.2.2 [109,110,111,112], which provides a matrix element calculation accurate to NLO in $\alpha _{\text {S}}$ for 0-jet and 1-jet final states and LO accuracy for 2-jets and 3-jets final states. The merging with the Sherpa parton shower [113] was performed using the ME+PS@NLO prescription [114]. The NLO EW corrections were applied as a function of the invariant mass $m_{ZZ^{*}}$ of the $ZZ^*$ system [115, 116].

The gluon-induced $ZZ^{*}$ production was modelled by Sherpa v2.2.2 [109,110,111] at LO in QCD for 0-jet and 1-jet final states. The higher-order QCD effects for the $gg\rightarrow ZZ^{*}$ continuum production cross-section were calculated for massless quark loops [117,118,119] in the heavy top-quark approximation [120], including the interference with $gg\rightarrow H^{*} \rightarrow ZZ$ processes [121, 122]. The $gg \rightarrow ZZ$ simulation was scaled by a K-factor of 1.7 ± 1.0, which is defined as the ratio of the higher-order to the leading-order cross-section predictions.

Production of $ZZ^{*}$ via vector-boson scattering was simulated with the Sherpa v2.2.2 [112] generator. The LO-accurate matrix elements were matched to a parton shower using the MEPS@LO prescription.

For all $ZZ^{*}$ processes modelled using Sherpa, the NNPDF3.0nnlo PDF set [60] was used, along with a dedicated set of tuned parton-shower parameters.

For additional checks, the $q{\bar{q}}$-initiated $ZZ^{*}$ continuum background was also modelled using Powheg-Box v2 and MadGraph5_aMC@NLO, using the CT10 [57] and the PDF4LHC NLO PDF set, respectively. For the former, the matrix element was generated at NLO accuracy in QCD and effects of singly resonant amplitudes and interference effects due to $Z/\gamma ^*$ were included. For the latter, the simulations are accurate to NLO in QCD for zero and one additional parton merged with the FxFx merging scheme. For both, the Pythia 8 generator was used for the modelling of parton showering, hadronisation, and the underlying event. The AZNLO and A14 sets of tuned parameters were used for the simulations performed with Powheg-Box v2 and MadGraph5_aMC@NLO generators, respectively.

The WZ background [123] was modelled at NLO accuracy in QCD using Powheg-Box v2 with the CT10 PDF set and was interfaced to Pythia 8, using the AZNLO set of tuned parameters for modelling of parton showering, hadronisation, and the underlying event and to EvtGen v1.2.0 for the simulation of bottom and charm hadron decays. The triboson backgrounds ZZZ, WZZ, and WWZ with four or more prompt leptons (VVV) were modelled at NLO accuracy for the inclusive process and at LO for up to two additional parton emissions using Sherpa v2.2.2.

The simulation of $ttZ$ events with both top quarks decaying semileptonically and the Z boson decaying leptonically was performed with MadGraph5_aMC@NLO using the NNPDF3.0nlo [60] PDF set interfaced to Pythia 8 using the A14 set of tuned parameters, and the total cross-section was normalised to a prediction computed at NLO in the QCD and EW couplings [98]. For modelling comparisons, Sherpa v2.2.1 was used to simulate $ttZ$ events at LO. The tWZ, $ttWW$, $ttWZ$, $ttZ\gamma$, $ttZZ$, $ttt$, $tttt$ and tZ background processes were simulated with MadGraph5_aMC@NLO interfaced to Pythia 8, using the A14 set of tuned parameters. These processes are collectively referred to as the tXX process.

The modelling of events containing Z bosons with associated jets ($Z+\,$jets) was performed using the Sherpa v2.2.1 generator. Matrix elements were calculated for up to two partons at NLO and four partons at LO using Comix [110] and OpenLoops [111], and merged with the Sherpa parton shower [113] using the ME+PS@NLO prescription [114]. The NNPDF3.0nnlo PDF set is used in conjunction with dedicated set of tuned parton-shower parameters.

The $tt$ background was modelled using Powheg-Box v2 with the NNPDF3.0nlo PDF set. This simulation was interfaced to Pythia 8, using the A14 set of tuned parameters, for parton showering, hadronisation, and the underlying event, and to EvtGen v1.2.0 for heavy-flavour hadron decays. Simulated $Z+\,$jets and $tt$ background samples were normalised to the data-driven estimates described in Sect. 6.

Generated events were processed through the ATLAS detector simulation [124] within the $\textsc {Geant} 4$ framework [125] and reconstructed in the same way as collision data. Additional pp interactions in the same and nearby bunch crossings were included in the simulation. Pile-up events were generated using Pythia 8 with the A2 set of tuned parameters [126] and the MSTW2008LO PDF set [127]. The simulation samples were weighted to reproduce the distribution of the number of interactions per bunch crossing observed in data.

4 Event selection

4.1 Event reconstruction

The selection and categorisation of the Higgs boson candidate events rely on the reconstruction and identification of electrons, muons, and jets, closely following the analyses reported in Refs. [17, 128].

Proton–proton collision vertices are constructed from reconstructed trajectories of charged particles in the ID with transverse momentum $p_{\text {T}}>$ 500 $\text {Me} \text {V}$. Events are required to have at least one collision vertex with at least two associated tracks. The vertex with the highest $\sum {p_{\text {T}} ^2}$ of reconstructed tracks is selected as the primary vertex of the hard interaction. The data are subjected to quality requirements to reject events in which detector components were not operating correctly.

Electron candidates are reconstructed from energy clusters in the electromagnetic calorimeter that are matched to ID tracks [129]. A Gaussian-sum filter algorithm [130] is used to compensate for radiative energy losses in the ID for the track reconstruction, while a dynamical, topological cell-based approach for cluster building is used to improve the energy resolution relative to the previous measurements in Refs. [17, 128], in particular for the case of bremsstrahlung photons. Electron identification is based on a likelihood discriminant combining the measured track properties, transition radiation response, electromagnetic shower shapes and the quality of the track–cluster matching. The ‘loose’ likelihood criteria, applied in combination with track hit requirements, provide an electron reconstruction and identification efficiency of at least 90% for isolated electrons with $p_{\text {T}} > 30$ $\text {Ge} \text {V}$ and 85%–90% below [129]. Electrons are required to have $E_{\text {T}}>$ 7 $\text {Ge} \text {V}$ and pseudorapidity $|\eta |<$ 2.47, with their energy calibrated as described in Ref. [129].

Muon candidate reconstruction [131] within the range $|\eta |<2.5$ is primarily performed by a global fit to fully reconstructed tracks in the ID and the MS, with a ‘loose’ [131] identification criterion applied. This criterion has an efficiency of at least 98% for isolated muons with $p_{\text {T}} = 5$ $\text {Ge} \text {V}$ and rises to 99.5% at higher $p_{\text {T}}$. At the centre of the detector ($|\eta |<$ 0.1), which has a reduced MS geometrical coverage, muons are also identified by matching a fully reconstructed ID track to either an MS track segment or a calorimeter energy deposit consistent with a minimum-ionising particle (calorimeter-tagged muons). For these two cases, the muon momentum is measured from the ID track alone. In the forward MS region (2.5 $<|\eta |<$ 2.7), outside the full ID coverage, MS tracks with hits in the three MS layers are accepted and combined with forward ID tracklets, if they exist (stand-alone muons). Calorimeter-tagged muons are required to have $p_{\text {T}}>$ 15 $\text {Ge} \text {V}$. For all other muon candidates, the transverse momentum is required to be greater than 5 $\text {Ge} \text {V}$. The muon momentum is calibrated using the procedure described in Ref. [131]. Muons with transverse impact parameter greater than 1 mm are rejected.² Additionally, muons and electrons are required to have a longitudinal impact parameter ($|z_0\sin \theta |$) less than 0.5 mm.

Jets are reconstructed using a particle flow algorithm [22] from noise-suppressed positive-energy topological clusters [132] in the calorimeter using the anti-$k_t$ algorithm [26, 27] with a radius parameter $R = 0.4$. Energy deposited in the calorimeter by charged particles is subtracted and replaced by the momenta of tracks that are matched to those topological clusters. Compared to only using topological clusters, jets reconstructed with the particle flow algorithm with $p_{\text {T}} >30~\text {Ge} \text {V}$ have approximately 10% better transverse momentum resolution. The two different algorithms have similar resolution for $p_{\text {T}}$ above 100 $\text {Ge} \text {V}$. The jet four-momentum is corrected for the calorimeter’s non-compensating response, signal losses due to noise threshold effects, energy lost in non-instrumented regions, and contributions from pile-up [22, 133, 134]. Jets are required to have $p_{\text {T}}>$ 30 $\text {Ge} \text {V}$ and $|\eta |<$ 4.5. Jets from pile-up with $|\eta |<$ 2.5 are suppressed using a jet-vertex-tagger multivariate discriminant [135, 136]. Jets with $|\eta |<$ 2.5 containing b-hadrons are identified using the MV2c10 b-tagging algorithm [137, 138], and its 60%, 70%, 77% and 85% efficiency working points are combined into a pseudo-continuous b-tagging weight [139] that is assigned to each jet.

Ambiguities are resolved if electron, muon, or jet candidates overlap in geometry or share the same detector information. If the two calorimeter energy clusters from the two electron candidates overlap, the electron with the higher $E_{\text {T}}$ is retained. If a reconstructed electron and muon share the same ID track, the muon is rejected if it is calorimeter-tagged; otherwise the electron is rejected. Reconstructed jets geometrically overlapping in a cone of radial size $\Delta R$ = 0.1 (0.2) with a muon (an electron) are also removed.

The missing transverse momentum vector, ${\mathop {{\textit{E}}}\limits ^\rightharpoonup {}_{\text {T}}^{\text {miss}}}$, is defined as the negative vector sum of the transverse momenta of all the identified and calibrated leptons, photons and jets and the remaining unclustered energy, where the latter is estimated from low-$p_{\text {T}}$ tracks associated with the primary vertex but not assigned to any lepton, photon, hadronically decaying $\tau$-lepton or jet candidate [140, 141]. The missing transverse momentum ($E_{\text {T}}^{\text {miss}}$) is defined as the magnitude of ${\mathop {{\textit{E}}}\limits ^\rightharpoonup {}_{\text {T}}^{\text {miss}}}$.

4.2 Selection of the Higgs boson candidates

A summary of the event selection criteria is given in Table 3. Events were triggered by a combination of single-lepton, dilepton and trilepton triggers with different transverse momentum thresholds. Single-lepton triggers with the lowest thresholds had strict identification and isolation requirements. Both the high-threshold single-lepton triggers and the multilepton triggers had looser selection criteria. Due to an increasing peak luminosity, these thresholds increased slightly during the data-taking periods [142, 143]. For single-muon triggers, the $p_{\text {T}}$ threshold ranged from between 20 and 26 $\text {Ge} \text {V}$, while for single-electron triggers, the $p_{\text {T}}$ threshold ranged from 24 to 26 $\text {Ge} \text {V}$. The global trigger efficiency for signal events passing the final selection is about 98%.

Table 3 Summary of the criteria applied to the selected Higgs boson candidate in each event. The mass threshold $m_{\text {min}}$ is defined in Sect. 4.1

In the analysis, at least two same-flavour and opposite-charge lepton pairs (hereafter referred to as lepton pairs) are required in the final state, resulting in one or more possible lepton quadruplets in each event. The three highest-$p_{\text {T}}$ leptons in each quadruplet are required to have transverse momenta above 20 $\text {Ge} \text {V}$, 15 $\text {Ge} \text {V}$ and 10 $\text {Ge} \text {V}$, respectively. To minimise the background contribution from non-prompt muons, at most one calorimeter-tagged or stand-alone muon is allowed per quadruplet.

The lepton pair with the invariant mass $m_{12}$ ($m_{34}$) closest (second closest) to the Z boson mass [144] in each quadruplet is referred to as the leading (subleading) lepton pair. Based on the lepton flavour, each quadruplet is classified into one of the following decay final states: 4$\mu$, 2e2$\mu$, 2$\mu$2e and 4e, with the first two leptons always representing the leading lepton pair. In each of these final states, the quadruplet with $m_{12}$ closest to the Z boson mass has priority to be considered for the selection of the final Higgs boson candidate. In case additional prompt leptons are present in the event, the priority may change due to the matrix-element based pairing as described later on. All quadruplets are therefore required to pass the following selection criteria.

To ensure that the leading lepton pair from the signal originates from a Z boson decay, the leading lepton pair is required to satisfy 50 $\text {Ge} \text {V}$ $<m_{12}<$ 106 $\text {Ge} \text {V}$. The subleading lepton pair is required to have a mass $m_{\mathrm {min}}<m_{34}<$ 115 $\text {Ge} \text {V}$, where $m_{\mathrm {min}}$ is 12 $\text {Ge} \text {V}$ for the four-lepton invariant mass $m_{4\ell }$ below 140 $\text {Ge} \text {V}$, rising linearly to 50 $\text {Ge} \text {V}$ at $m_{4\ell }=$ 190 $\text {Ge} \text {V}$ and then remaining at 50 $\text {Ge} \text {V}$ for all higher $m_{4\ell }$ values. This criterion suppresses the contributions from processes in which an on-shell Z boson is produced in association with a leptonically decaying meson or virtual photon. In the 4e and 4$\mu$ final states, the two alternative opposite-charge lepton pairings within a quadruplet are required to have a dilepton mass above 5 $\text {Ge} \text {V}$ to suppress the $J/\psi$ background. All leptons in the quadruplet are required to have an angular separation of $\Delta R>$ 0.1.

Each electron (muon) track is required to have a transverse impact parameter significance $|d_0/\sigma (d_0)|<5$ (3), to suppress the background from heavy-flavour hadrons. Reducible background from the Z+jets and $tt$ processes is further suppressed by imposing track-based and calorimeter-based isolation criteria on each lepton [131, 145]. A scalar $p_{\text {T}}$ sum (track isolation) is made from the tracks with $p_{\text {T}} > 500$ $\text {Me} \text {V}$ which either originate from the primary vertex or have $|z_0\sin \theta | < 3$ mm if not associated with any vertex and lie within a cone of $\Delta R =$ 0.3 around the muon or electron. Above a lepton $p_{\text {T}}$ of 33 $\text {Ge} \text {V}$, this cone size falls linearly with $p_{\text {T}}$ to a minimum cone size of 0.2 at 50 $\text {Ge} \text {V}$. Similarly, the scalar $E_{\text {T}}$ sum (calorimeter isolation) is calculated from the positive-energy topological clusters that are not associated with a lepton track in a cone of $\Delta R =$ 0.2 around the muon or electron. The sum of the track isolation and 40% of the calorimeter isolation is required to be less than 16% of the lepton $p_{\text {T}}$. The calorimeter isolation is corrected for electron shower leakage, pile-up and underlying-event contributions. Both isolations are corrected for track and topological cluster contributions from the remaining three leptons. The pile-up dependence of this isolation selection is improved compared with that of the previous measurements [17, 128, 146] by optimising the criteria used for exclusion of tracks associated with a vertex other than the primary vertex and by the removal of topological clusters associated with tracks. The signal efficiency of the isolation criteria is greater than 80%, improving the efficiency by about 5% compared with the previous analysis for the same background rejection.

The four quadruplet leptons are required to originate from a common vertex point. A requirement corresponding to a signal efficiency of better than 99.5% is imposed on the $\chi ^2$ value from the fit of the four lepton tracks to their common vertex.

If there is more than one decay final state per event with the priority quadruplet ($m_{12}$ closest to $m_Z$) satisfying the selection criteria, the quadruplet from the final state with highest selection efficiency, i.e. ordered 4$\mu$, 2e2$\mu$, 2$\mu$2e and 4e, is chosen as the Higgs boson candidate.

In the case of VH or $ttH$ production, there may be additional prompt leptons present in the event, together with the selected quadruplet. Therefore, there is a possibility that one or more of the leptons selected in the quadruplet do not originate from a Higgs boson decay, but rather from the V boson leptonic decay or the top quark semileptonic decay. To improve the lepton pairing in such cases, a matrix-element-based pairing method assuming the SM tensor structure is used for all events containing at least one additional lepton with $p_{\text {T}} >12$ $\text {Ge} \text {V}$ and satisfying the same identification, isolation and angular separation criteria as the four quadruplet leptons [17, 128]. For all possible quadruplet combinations that satisfy the selection, a matrix element for the Higgs boson decay is computed at LO using the MadGraph5_aMC@NLO [55] generator, with the reconstructed lepton momentum vectors as inputs to the calculation. The quadruplet with the largest matrix-element value is selected as the Higgs boson candidate. This method leads to a 50% improvement in correctly identifying the leptons in the quadruplet as those originating from a Higgs boson decay if an extra lepton is identified. The impact of the matrix element on the expected invariant mass distribution is shown in Fig. 2a.

To improve the four-lepton invariant mass reconstruction, the reconstructed final-state radiation (FSR) photons in Z boson decays are accounted for using the same strategy as the previous publications [17, 128]. Collinear FSR candidates are defined as candidates with $\Delta R <0.15$ to the nearest lepton in the quadruplet. Collinear FSR candidates are considered only for muons from the leading lepton pair, while non-collinear FSR candidates are considered for both muons and electrons from leading and subleading Z bosons.

Collinear FSR candidates are selected from reconstructed photon candidates and from electron candidates that share an ID track with the muon. Further criteria are applied to each candidate, based on the following discriminants: the fraction, $f_1$, of cluster energy in the front segment of the EM calorimeter divided by the total cluster energy to reduce backgrounds from muon ionisation; the angular distance, $\Delta R_{\text {cluster},\mu }$, between the candidate EM cluster and the muon; and the candidate $p_{\text {T}}$, which must be at least 1 $\text {Ge} \text {V}$. For all selected electron candidates and for photon candidates with $p_{\text {T}} \,{<}\, 3.5~\text {Ge} \text {V}$, a requirement of $f_1\,{>}\,0.2$ and $\Delta R_{\text {cluster},\mu }\,{<}\,0.08$ is imposed. The collinear photon candidates with $p_{\text {T}} \,{>}\, 3.5~\text {Ge} \text {V}$ are selected if $f_1\,{>}\,0.1$ and $\Delta R_{\text {cluster},\mu }\,{<}\,0.15$. Non-collinear FSR candidates are selected only from reconstructed isolated photons meeting the ‘tight’ criteria [129, 147] and satisfying $p_{\text {T}} \,{>}\, 10~\text {Ge} \text {V}$ and $\Delta R_{\text {cluster},\ell }\,{>}\,0.15$.

Only one FSR candidate is included in the quadruplet, with preference given to collinear FSR and to the candidate with the highest $p_{\text {T}}$. An FSR candidate is added to the lepton pair if the invariant mass of the lepton pair is between 66 and 89 $\text {Ge} \text {V}$ and if the invariant mass of the lepton pair and the photon is below 100 $\text {Ge} \text {V}$. Approximately 3% of reconstructed Higgs boson candidates have an FSR candidate and its impact on the expected invariant mass distribution is shown in Fig. 2b.

The Higgs boson candidates within a mass window of 115 $\text {Ge} \text {V}$ $<m_{4\ell }<$ 130 $\text {Ge} \text {V}$ are selected as the signal region. Events failing this requirement but that are within a mass window of 105 $\text {Ge} \text {V}$ $<m_{4\ell }{<}$ 115 $\text {Ge} \text {V}$ or 130 $\text {Ge} \text {V}$ $<m_{4\ell }{<}$ 160 (350) $\text {Ge} \text {V}$ are assigned to the sideband regions used to estimate the leading backgrounds as described in Sect. 6.

The selection efficiencies of the simulated signal in the fiducial region $|y_H|<2.5$, where $y_H$ is the Higgs boson rapidity, are about 33%, 25%, 19% and 16%, in the 4$\mu$, 2e2$\mu$, 2$\mu$2e and 4e final states, respectively.

Fig. 2

Impact on the expected invariant mass distribution of the selected Higgs boson candidates due to (a) matrix-element-based pairing for candidates with at least one extra lepton and (b) accounting for final-state radiation for candidates with an FSR candidate. For (a), the overflow events are included in the last bin

5 Event categorisation and production mode discrimination

In order to be sensitive to different production bins in the framework of simplified template cross-sections, the selected Higgs boson candidates in the mass window $115~\text {Ge} \text {V}< m_{4\ell } < 130~\text {Ge} \text {V}$ are classified into several dedicated reconstructed event categories. In addition, the events in the mass sidebands are also categorised for purposes of background estimation described in Sect. 6. In general, more than one production mode contributes to each reconstructed event category, as well as various background processes. For this reason, multivariate discriminants are introduced in most of the mutually exclusive reconstructed event categories to distinguish between these contributions.

5.1 Event categorisation

For signal events, the classification is performed in the order shown in the middle-right panel of Fig. 1 (from bottom to top) and as described below. First, those events classified as enriched in the $ttH$ process are split according to the decay mode of the two W bosons from the top quark decays. For semileptonic and dileptonic decays (ttH-Lep-enriched), at least one additional lepton with $p_{\text {T}}>$ 12 $\text {Ge} \text {V}$³ together with at least two b-tagged jets (with 85% b-tagging efficiency), or at least five jets among which at least one b-tagged jet (with 85% b-tagging efficiency) or at least two jets among which at least one b-tagged jet (with 60% b-tagging efficiency) is required. For the fully hadronic decay (ttH-Had-enriched), there must be either at least five jets among which at least two b-tagged jets (with 85% b-tagging efficiency) or at least four jets among which at least one b-tagged jet (with 60% b-tagging efficiency). Events with additional leptons but not satisfying the jet requirements define the next category enriched in VH production events with leptonic vector-boson decay (VH-Lep-enriched).

The remaining events are classified according to their reconstructed jet multiplicity into events with no jets, exactly one jet or at least two jets. Events with at least two reconstructed jets are divided into two categories: one is a ‘BSM-like’ category ($2j$-BSM-like) and the other ($2j$) contains the bulk of events with significant contributions from the $\mathrm {VBF}$ and VH production modes in addition to $\mathrm {ggF}$. The $2j$-BSM-like category requires the invariant mass $m_{jj}$ of the two leading jets to be larger than 120 $\text {Ge} \text {V}$ and the four-lepton transverse momentum, $p_{\mathrm {T}}^{4\ell }$, to be larger than 200 $\text {Ge} \text {V}$; the remaining events are placed in the $2j$ category.

Events with zero or one jet in the final state are expected to be mostly from the $\mathrm {ggF}$ process. Following the particle-level definition of production bins in Sect. 1.1, the 1-jet category is further split into four categories with $p_{\mathrm {T}}^{4\ell }$ smaller than 60 $\text {Ge} \text {V}$ ($1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Low}$), between 60 and 120 $\text {Ge} \text {V}$ ($1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Med}$), between 120 and 200 $\text {Ge} \text {V}$ ($1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {High}$), and larger than 200 $\text {Ge} \text {V}$ ($1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {BSM}$-$\mathrm {like}$).

The largest number of $\mathrm {ggF}$ events and the highest $\mathrm {ggF}$ purity are expected in the zero-jet category. The zero-jet category is split into three categories with $p_{\mathrm {T}}^{4\ell }$ smaller than 10 $\text {Ge} \text {V}$ ($0j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Low}$), between 10 and 100 $\text {Ge} \text {V}$ ($0j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Med}$) and above 100 $\text {Ge} \text {V}$ ($0j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {High}$). The first two categories follow the production bin splitting, and the last category improves the discrimination between VH ($V \rightarrow \ell \nu / \nu \nu$) and $\mathrm {ggF}$.

As illustrated in Fig. 1, there is a dedicated reconstructed event category for each production bin except for $\mathrm {gg2H}$-$2j$, $\mathrm {qq2Hqq}$-$\textit{VH}$ and $\mathrm {qq2Hqq}$-$\mathrm {VBF}$. These production bins are largely measured from the 2-jet reconstruction category, and to a lesser extent from the 1-jet categories, using multivariate discriminants (see Sect. 5.2). The $\mathrm {gg2H}$-$p_{\mathrm {T}}^H$-$\mathrm {High}$ production bin is measured simultaneously in all reconstructed event categories with high transverse momentum of the four-lepton system, independent of the reconstructed jet multiplicity.

The rightmost panel of Fig. 1 shows the background event classification. For estimating the tXX process from the mass sideband, a tXX-enriched sideband category (SB-$tXX$-enriched) is defined, which includes events with at least two jets including at least one tagged as a b-jet with 60% efficiency and $E_{\text {T}}^{\text {miss}}$ > 100 $\text {Ge} \text {V}$ in the $m_{4\ell }$ mass range 105–115 $\text {Ge} \text {V}$ or 130–350 $\text {Ge} \text {V}$. This region is dominated by $ttZ$ (87%) and has small contributions from $tt$, $tttt$, tWZ, $ttW$, $ttWW$, $ttWZ$, $ttZ\gamma$, $ttZZ$ and tZ. The tXX process is expected to give the largest contribution in ‘$ttH$-like’ categories. The large mass range for this category, larger than for the non-resonant ZZ as discussed next, allows better statistical precision for the estimate of this background.

For the estimation of non-resonant $ZZ^{*}$ production, events not meeting the criteria for the SB-$tXX$-enriched category and in the $m_{4\ell }$ mass range 105–115 $\text {Ge} \text {V}$ or 130–160 $\text {Ge} \text {V}$ are split according to the number of reconstructed jets: exactly zero jets (SB-$0j$), exactly one jet (SB-$1j$) or at least two jets (SB-$2j$). This mass range limits the contribution from the single-resonance process, $Z \rightarrow 4\ell$, and from the on-shell ZZ process. Similarly, events in the same mass range with an extra reconstructed lepton separately form the SB-$\textit{VH}$-Lep-enriched category, which is enriched with signal events containing leptons from the associated V leptonic decay or the top quark semileptonic decay. This category is mainly designed to improve the expected sensitivity for VH-Lep by about 5%, having a $\textit{VH}$ purity of about 19%.

The expected number of signal events is shown in Table 4 for each reconstructed event category separately for each production mode. The $\mathrm {ggF}$ and $bbH$ contributions are shown separately to compare their relative contributions, but both belong in the same ($\mathrm {ggF}$) production bin. The highest $bbH$ event yield is expected in the $0j$ categories since the jets tend to be more forward than in the $ttH$ process, thus escaping the acceptance of the $ttH$ selection criteria. The sources of uncertainty in these expectations are detailed in Sect. 7. The signal composition in terms of the Reduced Stage-1.1 production bins is shown in Fig. 3.

The separation of the contributions from different production bins, such as the $\mathrm {gg2H}$-$2j$, $\mathrm {qq2Hqq}$-$\textit{VH}$ and $\mathrm {qq2Hqq}$-$\mathrm {VBF}$ components contributing in categories with two or more jets, is improved by means of discriminants obtained using multivariate data analysis, as described in the following section.

Table 4 The expected number of SM Higgs boson events with $m_H=$ 125 $\text {Ge} \text {V}$ for an integrated luminosity of 139 fb$^{-1}$ at $\sqrt{s}=$ 13 $\text {Te} \text {V}$ in each reconstructed event signal ($115<m_{4\ell }<130$ $\text {Ge} \text {V}$) and sideband ($m_{4\ell }$ in 105–115 $\text {Ge} \text {V}$ or 130–160 $\text {Ge} \text {V}$ for $ZZ^{*}$, 130–350 $\text {Ge} \text {V}$ for tXX) category, shown separately for each production bin of the Production Mode Stage. The $\mathrm {ggF}$ and $bbH$ yields are shown separately but both contribute to the same ($\mathrm {ggF}$) production bin, and ZH and WH are reported separately but are merged together for the final result. Statistical and systematic uncertainties, including those for total SM cross-section predictions, are added in quadrature. Contributions that are below 0.2% of the total signal in each reconstructed event category are not shown and are replaced by ‘−’

Fig. 3

Standard Model signal composition in terms of the Reduced Stage-1.1 production bins in each reconstructed event category. The $bbH$ contributions are included in the $\mathrm {ggF}$ production bins

5.2 Multivariate production mode discriminants

To further increase the sensitivity of the cross-section measurements in the production bins (Sect. 1.1), multivariate discriminants using neural networks (NNs) [148] are introduced in many of the reconstructed signal event categories as observables used in the statistical fit, described in Sect. 8.2. The NN architecture and training procedure are defined using Keras with TensorFlow [149, 150]. These networks are trained using several discriminating observables, as defined in Table 5, on simulated SM Higgs boson signals with $m_{H} = 125~\text {Ge} \text {V}$ or non-Higgs-boson background. Due to the low number of signal events expected in the $0j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {High}$, $1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {BSM}$-$\mathrm {like}$ and ttH-Lep-enriched categories, only the observed yield is used as the discriminant in these categories.

Two types of NNs are used: feed-forward multilayer perceptron (MLP) and recurrent (RNN) [148,149,150,151,152]. Each NN discriminant combines two RNNs, one for the $p_{\text {T}}$-ordered variables related to the four leptons in the quadruplet and one for variables related to jets, and an MLP with additional variables related to the full event. The jet RNN accepts inputs from up to three jets. The outputs of the MLP and the two RNNs are chained into another MLP to complete an NN discriminant, which is trained to approximate the posterior probability for an event to originate from a given process. This is used in each reconstructed event category to discriminate between two or three processes, e.g. $\mathrm {ggF}$, $\mathrm {VBF}$ and ZZ background in the $1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Low}$ category. The variables used to train the MLP and RNNs for each category along with the processes being separated are summarised in Table 5.

Table 5 The input variables used to train the MLP, and the two RNNs for the four leptons and the jets (up to three). For each category, the processes which are classified by an NN, their corresponding input variables and the observable used are shown. For example, there are eight input variables for the Lepton RNN being trained if $p_{\mathrm {T}}^{\ell }$ and $\eta _{\ell }$ are listed. Leptons and jets are denoted by ‘$\ell$’ and ‘j’. See the text for the definitions of the variables

The NN training variables not previously defined are listed as follows. The kinematic discriminant $D_{ZZ^*}$[153], defined as the difference between the logarithms of the squared matrix elements for the signal decay (same as in Sect. 4) and squared matrix elements for the background process, is used to distinguish $\mathrm {ggF}$ from the non-resonant ZZ background. Three angles [7] are used to further distinguish these processes: the cosine of the leading Z boson’s production angle $\theta ^*$ in the four-lepton rest frame; the cosine of $\theta _1$ defined as the angle between the negatively charged lepton of the leading Z in the leading Z rest frame and the direction of flight of the leading Z in the four-lepton rest frame; and the angle $\phi _{ZZ}$, between the two Z decay planes in the four-lepton rest frame. The angular separation of the leading jet from the $4\ell$ system, $\Delta R_{4\ell j}$, is used to distinguish $\mathrm {VBF}$ or $ttH$ from $\mathrm {ggF}$. For categories with two or more jets, kinematic variables that also include the information from the two leading jets are used: the invariant mass, $m_{jj}$; the transverse momentum of the $4\ell$ and the 2-jet system, $p_{\text {T}}^{4\ell jj}$; and the Zeppenfeld variable, $\eta _{ZZ}^{\text {Zepp}}= \left| \eta _{4\ell } - \frac{\eta _{j1}+\eta _{j2}}{2}\right|$ [154]. The number of reconstructed jets, $N_\mathrm {jets}$, the number of b-tagged jets at 70% tagging efficiency, $N_{b\text {-jets,70\%}}$, and the scalar sum of the $p_{\text {T}}$ of all reconstructed jets, $H_{\text {T}}$, are used to identify the $ttH$ process.

Depending on the category and the number of processes being targeted, the NN has two or three output nodes. The value computed at each node represents the probability, with an integral of one, for the event to originate from the given process. For example, for the 0-jet category, two probabilities are evaluated, $\mathrm {NN_{ggF}}$ and $\mathrm {NN}_{ZZ}$. As these two values are a linear transformation of each other, only one output, $\mathrm {NN_{ggF}}$, is used as a discriminant in the fit model. In categories with three targeted processes, only two of the three corresponding output probabilities are independent. In a given category, a selection is applied on one of the three output probabilities to split the events in two subcategories. This output probability is then used as the discriminant for the subcategory of events passing the selection, while for the other subcategory one of the two remaining output probabilities is used. The selection criterion is chosen so as to provide the largest purity of the targeted process for events passing the selection. For example, in the 1-jet category, $\mathrm {NN_{VBF}}$ and $\mathrm {NN}_{ZZ}$ are used. The subcategory of events with $\mathrm {NN}_{ZZ}$ larger than 0.25 uses $\mathrm {NN}_{ZZ}$ as the discriminant in the fit model, while $\mathrm {NN_{VBF}}$ is used in the remaining subcategory. The subcategory definitions and observables used in all reconstructed event categories are summarised in Table 5.

6 Background contributions

6.1 Background processes with prompt leptons

Non-resonant SM $ZZ^{*}$ production via qq annihilation, gluon–gluon fusion and vector-boson scattering can result in four prompt leptons in the final state and constitutes the largest background for the analysis. While for the previous analyses [17, 128], simulation was exclusively used to estimate both the shape and normalisation, in this analysis the normalisation is constrained by a data-driven technique. This allows the systematic uncertainty to be reduced by removing both the theoretical and luminosity uncertainties contributing to the normalisation uncertainty.

As outlined in Sect. 5.1, to estimate the normalisation, sideband categories in the $m_{4\ell }$ mass region 105–115 $\text {Ge} \text {V}$ and 130–160 $\text {Ge} \text {V}$ are defined according to the jet multiplicity (SB-$0j$, SB-$1j$, SB-$2j$). The normalisation of the $ZZ^{*}$ background is simultaneously fitted with a common normalisation factor for signal region and sideband categories with the same jet multiplicity. For example, the $ZZ^{*}$ background is scaled by a common factor for $2j$, $2j$-BSM-like and SB-$2j$ categories. The background shape templates for $\mathrm {NN}$ discriminants and the expected fraction of events in relevant reconstructed signal-region event categories are obtained from simulation. As shown in Fig. 4a, good agreement is found between data and simulation for the shape of the NN observable. All expected distributions are shown after the final fit to the data for the Production Mode measurement (see Sect. 8) and are referred to as post-fit distributions in the following. The simulated distributions of the observables $p_{\mathrm {T}}^{4\ell }$ and $m_{jj}$ employed for the prediction of event fractions in each event category also agree with data, as seen in Figs. 4b, c respectively. The estimation of the $ZZ^{*}$ process in the jet multiplicity bins removes one of the leading theoretical uncertainties [155]. Due to the limited sensitivity and the low expected yield, the normalisation of $ZZ^{*}$ in $ttH$-like categories is estimated from simulation.

Fig. 4

The observed and expected (post-fit) distributions for an integrated luminosity of 139 fb$^{-1}$ at $\sqrt{\mathrm {s}}=$ 13 $\text {Te} \text {V}$ in the different background enriched regions: (a) $\mathrm {NN_{ggF}}$ in the SB-$0j$ sideband region, (b) $p_{\mathrm {T}}^{4\ell }$ in the sideband region combining the SB-$0j$, SB-$1j$ and SB-$2j$ categories, (c) $m_{jj}$ in the SB-$2j$ category, and (d) $N_\mathrm {jets}$ in the SB-$tXX$-enriched region. The SM Higgs boson signal is assumed to have a mass of $m_{H}$ = 125 $\text {Ge} \text {V}$. The uncertainty in the prediction is shown by the hatched band, calculated as described in Sect. 7. For comparison only, the hatched band includes the theoretical uncertainties in the SM cross-section for the signal and the background processes

Similarly, backgrounds affecting the $ttH$-like categories are estimated simultaneously from an enriched sample selected in a dedicated sideband region (SB-$tXX$-enriched), with the mass cut extended up to 350 $\text {Ge} \text {V}$ to improve the statistical precision of the estimate. The normalisation of the tXX process is simultaneously fitted across the ttH-Lep-enriched, ttH-Had-enriched and SB-$tXX$-enriched categories. The $N_\mathrm {jets}$ observable distribution, which is used to predict the event fractions in each category, is shown in Fig. 4d and agrees with data. In all other categories, the sensitivity of the tXX measurement is limited due to a small number of expected tXX events and its normalisation is estimated from simulation.

The contribution from VVV processes is estimated for all categories using the simulated samples presented in Sect. 3.

6.2 Background processes with non-prompt leptons

Other processes, such as $Z+\mathrm {jets}$, $tt$, and WZ, containing at least one jet, photon or lepton from a hadron decay that is misidentified as a prompt lepton, also contribute to the background. These ‘reducible’ backgrounds are significantly smaller than the non-resonant $ZZ^*$ background and are estimated from data using different approaches for the $\ell \ell +\mu \mu$ and $\ell \ell +ee$ final states [17, 128].

In the $\ell \ell +\mu \mu$ final states, the normalisation of the $Z+\mathrm {jets}$ and $tt$ backgrounds are determined by performing fits to the invariant mass of the leading lepton pair in dedicated independent control regions. The shape of the invariant mass distribution for each region is parameterised using simulated samples. In contrast to the previous analyses [17, 128], this fit is performed independently for each reconstructed event category, which removes the use of simulation to estimate the event fractions in these categories.

The control regions used to estimate this background are defined by closely following the requirements outlined in Sect. 4.2. The definition and modified requirements for each of the four control regions are:

• 1. an enhanced heavy-flavour control region with inverted impact-parameter and relaxed isolation requirements on the subleading lepton pair and relaxed vertex $\chi ^2$ requirements,

• 2. an enhanced $tt$ $e\mu +\mu \mu$ control region with an opposite-flavour leading lepton pair $e\mu$ and relaxed impact-parameter, isolation, and opposite-sign charge requirements on the subleading lepton pair $\mu \mu$, as well as relaxed vertex $\chi ^2$ requirements,

• 3. an enhanced light-flavour control region with inverted isolation requirements for at least one lepton in the subleading lepton pair, and

• 4. a same-sign $\ell \ell +\mu ^{\pm }\mu ^{\pm }$ control region with relaxed impact-parameter and isolation requirements.

The first two are the primary control regions used to estimate $Z+\mathrm {jets}$ and $tt$, and the latter two improve the estimate by reducing the statistical error of the fitted normalisation.

The background normalisations are obtained separately for the $Z+\mathrm {jets}$ and $tt$ background processes using the simultaneous fit in the four control regions. The normalisation $n_i^{CR}$ in each control region CR for the background process i is expressed as a fraction, $n_i^{CR} = t_i^{CR} \times N_i^{VR}$, of the normalisation $N_i^{VR}$ in a dedicated relaxed validation region (VR). $N_i^{VR}$ is used as the common parameter when fitting the normalisations in the different CRs. The transfer factor $t_i^{CR}$ is the ratio of the background contribution in the relaxed validation region and the given control regions. The relaxed validation region is defined by following the requirements outlined in Sect. 4.2 but by relaxing the impact-parameter and isolation requirements on the subleading lepton pair. This region contains a substantially larger number of events compared with the other four control regions, allowing a more reliable prediction of the shapes of the NN distributions. The shapes of the background NN distribution are then extrapolated together with the corresponding background normalisation from the relaxed validation to the signal region by means of additional transfer factors $T_i$. Transfer factors $t_i^{CR}$ and $T_i$ to extrapolate the background contributions from the control regions to the relaxed validation region and from there to the signal region are estimated from simulation and validated in several additional data control regions

The $\ell \ell +ee$ control-region selection requires the electrons in the subleading lepton pair to have the same charge, and relaxes the identification, impact parameter and isolation requirements on the electron candidate with the lowest transverse energy. This fake electron candidate, denoted by X, can be a light-flavour jet, an electron from photon conversion or an electron from heavy-flavour hadron decay. The heavy-flavour background is determined from simulation. Good agreement is observed between simulation and data in a heavy-flavour enriched control region.

The remaining background is separated into light-flavour and photon conversion background components using the sPlot method [156] which is performed on electron candidates X, separately for each reconstructed category in bins of the jet multiplicity and the transverse momentum of the electron candidate. The size of the two background components is obtained from a fit to the number of hits from the electron candidate X in the innermost ID layer in the $\ell \ell +ee$ data control region, where a hit indicates either a hadron track or an early conversion. A hit in the next-to-innermost pixel layer is used when the electron falls in a region that was either not instrumented with an innermost pixel layer module or where the module was not operating. The templates of the final discriminants for the mentioned fit of the light-flavour and photon conversion background components are obtained from simulated $Z+X$ events with an on-shell Z boson decay candidate accompanied by an electron X selected using the same criteria as in the $\ell \ell +ee$ control region. The simulated $Z+X$ events are also used to obtain the transfer factor for the X candidate for the extrapolation of the light-flavour and photon conversion background contributions from the $\ell \ell +ee$ control region to the signal region, after correcting the simulation to match the data in dedicated control samples of $Z+X$ events. The extrapolation to the signal region is also performed in bins of the electron transverse momentum and the jet multiplicity, separately for each reconstructed event category. A method similar to that for the $\ell \ell +\mu \mu$ final state is used to extract the NN shape, where the fractions of events from light-flavour jets and photon conversions are estimated from simulation and corrected transfer factors are used.

7 Systematic uncertainties

The systematic uncertainties are categorised into experimental and theoretical uncertainties. The first category includes uncertainties in lepton and jet reconstruction, identification, isolation and trigger efficiencies, energy resolution and scale, and uncertainty in the total integrated luminosity. Uncertainties from the procedure used to derive the data-driven background estimates are also included in this category. The second category includes uncertainties in theoretical modelling of the signal and background processes.

The uncertainties can affect the signal acceptance, selection efficiency and discriminant distributions as well as the background estimates. The dominant sources of uncertainty and their effect are described in the following subsections. The impact of these uncertainties on the measurements is summarised in Table 6.

Table 6 The impact of the dominant systematic uncertainties (in percent) on the cross-sections in production bins of the Production Mode Stage and the Reduced Stage 1.1. Similar sources of systematic uncertainties are grouped together: luminosity (Lumi.), electron/muon reconstruction and identification efficiencies and pile-up modelling (e, $\mu$, pile-up), jet energy scale/resolution and b-tagging efficiencies (Jets, flav. tag), uncertainties in reducible background (reducible bkg), theoretical uncertainties in $ZZ^{*}$ background and tXX background, and theoretical uncertainties in the signal due to parton distribution function (PDF), QCD scale (QCD) and parton showering algorithm (Shower). The uncertainties are rounded to the nearest 0.5%, except for the luminosity uncertainty, which is measured to be 1.7% and increases for the VH signal processes due to the simulation-based normalisation of the VVV background

7.1 Experimental uncertainties

The uncertainty in the combined 2015–2018 integrated luminosity is 1.7% [157], obtained using the LUCID-2 detector [158] for the primary luminosity measurements. This uncertainty affects the signal and the normalisation of the simulated background estimates when not constrained by the data sidebands.

The uncertainty in the predicted yields due to pile-up modelling ranges between 1% and 2% and is derived by varying the average number of pile-up events in the simulation to cover the uncertainty in the ratio of the predicted to measured inelastic cross-sections [159].

The electron (muon) reconstruction, isolation and identification efficiencies, and the energy (momentum) scale and resolution are derived from data using large samples of $J/\psi \rightarrow \ell \ell$ and $Z\rightarrow \ell \ell$ decays [129, 131]. Typical uncertainties in the predicted yields for the relevant decay channels due to the identification and reconstruction efficiency uncertainties are below $1\%$ for muons and 1%–2% for electrons. The uncertainty in the expected yields due to the muon and electron isolation efficiency is also taken into account, with the typical size being 1%. The uncertainties in the trigger efficiencies have a negligible impact. The uncertainties in the electron and muon energy and momentum scale and resolution are small and also have a negligible impact on the measurements.

The uncertainties in the jet energy scale and resolution are in the range 1%–3% [133]. The impact of these uncertainties is more relevant for the VH, $\mathrm {VBF}$ and $ttH$ production mode cross-sections (3%–5%) and for all the Reduced Stage-1.1 cross-section measurements, including the $\mathrm {ggF}$ process split into the different $N_\mathrm {jets}$ exclusive production bins (5%–20%).

The uncertainty in the calibration of the b-tagging algorithm, which is derived from dileptonic $tt$ events, amounts to a few percent over most of the jet $p_{\text {T}}$ range [138]. This uncertainty is only relevant in the ttH category, with its expected impact being approximately 1% in the $ttH$ cross-section measurement. The uncertainties associated with the $E_{\text {T}}^{\text {miss}}$ reconstruction have a negligible impact.

A shift in the simulated Higgs boson mass corresponding to the precision of the Higgs boson measurement, $m_{H} = 125.09 \pm 0.24$ $\text {Ge} \text {V}$ [160], is shown to have a negligible impact on the signal acceptance. A small dependency of the $\mathrm {NN_{ggF}}$ discriminant shape in the $0j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Low}$ and $0j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Med}$ categories on $m_{H}$ is observed for the signal (below 2% in the highest NN score bins) and is included in the signal model. This uncertainty affects the measurement of $\mathrm {ggF}$ production, as well as the measurements in other production bins with large $\mathrm {ggF}$ contamination.

For the data-driven measurement of the reducible background, three sources of uncertainty are considered: statistical uncertainty, overall systematic uncertainty for each of $\ell \ell +\mu \mu$ and $\ell \ell +ee$, and a shape systematic uncertainty that varies with the reconstructed event category. Since the yields are estimated by using a statistical fit to a control data region with large statistics, the inclusive background estimate has a relatively small (3%) statistical uncertainty. The systematic uncertainty for $\ell \ell +\mu \mu$ and the heavy-flavour component of $\ell \ell +ee$ is estimated by comparing the lepton identification, isolation and impact parameter significance efficiency between data and simulated events in a separate region, enriched with on-shell Z boson decays accompanied by an electron or a muon. For both the $\ell \ell +\mu \mu$ and $\ell \ell +ee$ estimates, the difference in efficiency is assigned as the uncertainty in the extrapolation of the yield estimate from the control region to the signal region. For the $\ell \ell +ee$ light-flavour component, the efficiency is derived from an enriched control region with a systematic uncertainty estimated by varying the assumed light- and heavy-flavour components. These inclusive uncertainties (6%) are treated as correlated across the reconstructed event categories. Finally, there are additional uncorrelated uncertainties (8%–70%) in the fraction of the reducible background in each event category due to the statistical precision of the simulated samples.

7.2 Theoretical uncertainties

The theoretical modelling of the signal and background processes is affected by uncertainties due to missing higher-order corrections, modelling of parton showers and the underlying event, and PDF$+\alpha _{\text {S}}$ uncertainties.

The impact of the theory systematic uncertainties on the signal depends on the kind of measurement that is performed. For signal-strength measurements, defined as the measured cross-section divided by the SM prediction, or interpretation of cross-section using the EFT approach, each source of theory uncertainty affects both the acceptance and the predicted SM cross-section. For the cross-section measurements, only effects on the acceptance need to be considered.

The impact of the theory systematic uncertainties on the background depends on the method of estimating the normalisation. If simulation is used, the uncertainties in the acceptance and the predicted SM cross-section are included. If the normalisation is estimated from a data-driven method, only the impact on the relative event fractions between categories is considered.

One of the dominant sources of theoretical uncertainty is the prediction of the $\mathrm {ggF}$ process in the different $N_\mathrm {jets}$  categories. The $\mathrm {ggF}$ process gives a large contribution in categories with at least two jets. To estimate the variations due to the impact of higher-order contributions not included in the calculations and migration effects on the $N_\mathrm {jets}$ $\mathrm {ggF}$ cross-sections, the approach described in Refs. [24, 161] is used, which exploits the latest predictions for the inclusive jet cross-sections. In particular, the uncertainty from the choice of factorisation and renormalisation scales, the choice of resummation scales, and the migrations between the 0-jet and 1-jet phase-space bins or between the 1-jet and $\ge 2$-jet bins are considered [24, 162,163,164]. The impact of QCD scale variations on the Higgs boson $p_{\text {T}}$ distribution is taken into account as an additional uncertainty. The uncertainty in higher-order corrections to the Higgs boson $p_{\text {T}}$ originating from the assumption of infinite top quark mass in the heavy-quark loop is also taken into account by comparing the $p_{\text {T}}$ distribution predictions to finite-mass calculations. An additional uncertainty in the acceptance of the $\mathrm {ggF}$ process in $\mathrm {VBF}$ topologies [165] due to missing higher orders in QCD in the calculation is estimated by variations of the renormalisation and factorisation scales using fixed-order calculations with MCFM [166]. An additional uncertainty in the Higgs boson $p_{\text {T}}$ distribution, derived by varying the renormalisation, factorisation and NNLOPS scale in the simulation, in the 0-jet topology is considered. This is particularly relevant when measuring the inclusive $\mathrm {ggF}$ cross-section using the $p_{\mathrm {T}}^{4\ell }$ categories for events with no jet activity. To account for higher-order corrections to $p_{\text {T}}^{Hjj}$, which is used as an NN input variable, the uncertainty is derived by comparing the predicted distribution obtained using Powheg NNLOPS and MadGraph5_aMC@NLO with the FxFx merging scheme.

For the $\mathrm {VBF}$ production mode, the uncertainty due to missing higher orders in QCD is parameterised using the scheme outlined in Ref. [23]. The migration effects due to the selection criteria imposed on the number of jets, transverse momentum of the Higgs boson, transverse momentum of the Higgs boson and the leading dijet system and the invariant mass of the two leading jets, used to define the full Stage 1.1 STXS production bins, are computed by varying the renormalisation and factorisation scales by a factor of two. The uncertainties are cross-checked with fixed-order calculations. Similarly, for the VH production mode with the associated V decaying leptonically, the scale variations are parameterised as migration effects due to the selection criteria imposed on the number of jets and the transverse momentum of the associated boson [167].

For the VH production mode with the associated V decaying hadronically and the $ttH$ production mode, the uncertainty due to missing higher orders in QCD is obtained by varying the renormalisation and factorisation scales by a factor of two. The configuration with the largest impact, as quantified by the relative difference between the varied and the nominal configuration, is chosen to define the uncertainty in each experimental category. These uncertainties are treated as uncorrelated among the different production modes. Due to the limited accuracy of the simulated samples, the uncertainties evaluated using this method for the total cross-sections are larger than those described in Ref. [24].

The uncertainties in the acceptance due to the modelling of parton showers and the underlying event are estimated with AZNLO tune eigenvector variations and by comparing the acceptance using the parton showering algorithm from Pythia 8 with that from Herwig 7 [168] for all signal processes. The uncertainty due to each AZNLO tune variation is taken as correlated among the different production modes while the difference between the parton showering algorithms is treated as an uncorrelated uncertainty. The uncertainties due to higher-order corrections to the Higgs boson decay are modelled using the PROPHECY4F [102, 105] and Hto4L [104, 169] generators. These corrections are below 2% and have a negligible impact on the results. A 100% uncertainty is assigned to heavy-flavour quark production modelling for the $\mathrm {ggF}$ contribution entering in the ttH category. This has a negligible impact on the results.

The impact of the PDF uncertainty is estimated with the thirty eigenvector variations of the PDF4LHC_nlo_30 Hessian PDF set following the PDF4LHC recommendations [40]. The modification of the predictions originating from each eigenvector variation is added as a separate source of uncertainty in the model. The same procedure is applied for the $\mathrm {ggF}$, $\mathrm {VBF}$, VH and $ttH$ processes, enabling correlations to be taken into account in the fit model.

The impacts of the theoretical uncertainties, as described above, on the shape of NN discriminants are also considered. For $\mathrm {ggF}$ production, a further cross-check is performed by comparing the NN shapes in the corresponding categories as predicted by Powheg NNLOPS and MadGraph5_aMC@NLO with the FxFx merging scheme. All the NN shapes from the two generators agree within the scale variations and, therefore, no additional shape uncertainty is included.

For signal-strength measurements, an additional uncertainty related to the $H\rightarrow ZZ^*$ branching ratio prediction [102, 105] is included in the measurement.

Since the normalisation of the $ZZ^{*}$ process in most reconstructed event categories is constrained by performing a simultaneous fit to sideband regions enriched in this contribution together with the signal regions, most of the theoretical uncertainty in the normalisation for this background vanishes. Nevertheless, uncertainties in the shapes of the discriminants for the $ZZ^{*}$ background and in the relative contribution of this background between the sidebands and the signal regions are taken into account. The uncertainties due to missing higher-order effects in QCD are estimated by varying the factorisation and renormalisation QCD scales by a factor of two; the impact of the PDF uncertainty is estimated by using the MC replicas of the NNPDF3.0 PDF set. Uncertainties due to parton shower modelling for the $ZZ^{*}$ process are considered as well. The impact of these uncertainties is below 2% for all production mode cross-sections measured. In addition, a comparison between Sherpa and Powheg is also taken as an additional source of systematic uncertainty. This model uncertainty is treated as uncorrelated among the different sideband-to-signal region extrapolations (in 0-jet, 1-jet and 2-jet categories).

The uncertainty in the gluon-initiated and the vector-boson-initiated $ZZ^{*}$ process is taken into account by changing the relative composition of the quark-initiated, the gluon-initiated and the vector-boson-scattered $ZZ^{*}$ components according to the theoretical uncertainty in the predicted cross-sections and the respective K-factors. In addition, the event yield and NN discriminant shapes in each event category are compared with the data in an $m_{4\ell }$ sideband around the signal region (105 $\text {Ge} \text {V}$ $<m_{4\ell }<$ 115 $\text {Ge} \text {V}$ or 130 $\text {Ge} \text {V}$ $<m_{4\ell }<$ 160 $\text {Ge} \text {V}$). Good agreement between the Sherpa predictions and the data is found.

For the tXX process, uncertainties due to PDF and QCD scale variations are considered in the relative fraction of events present in the $ttH$-like categories, in the SB-$tXX$-enriched control region and in the NN discriminant shape. Differences between MadGraph5_aMC@NLO and Sherpa are considered as an additional systematic uncertainty. For all other categories where this process is estimated from simulation, the impact of these uncertainties on the SM cross-section and acceptance are also considered.

Uncertainties in the PDF and in missing higher-order corrections in QCD are applied to the VVV background estimate, which is fully taken from MC simulation.

To probe the tensor structure of the Higgs boson coupling in the EFT approach, theoretical uncertainties due to PDF and QCD scale variations are assigned to the signal predictions based on the simulated highest-order SM signal samples. The same uncertainties are assigned to all corresponding BSM signal predictions, since it is shown using the MC signal samples simulated at LO accuracy that the uncertainties change negligibly as a function of the Wilson coefficients.

Fig. 5

The observed and expected (post-fit) four-lepton invariant mass distributions for the selected Higgs boson candidates, shown for an integrated luminosity of $139$ fb$^{-1}$ at $\sqrt{s}=13$ $\text {Te} \text {V}$. The SM Higgs boson signal is assumed to have a mass $m_{H}$ = 125 $\text {Ge} \text {V}$. The uncertainty in the prediction is shown by the hatched band, calculated as described in Sect. 7. For comparison only, the hatched band includes the theoretical uncertainties in the SM cross-section for the signal and the background processes

8 Measurement of the Higgs boson production mode cross-sections

8.1 Observed data

The expected and observed four-lepton invariant mass (post-fit) distributions of the selected Higgs boson candidates after the event selection are shown in Fig. 5.

The observed and expected (post-fit) distributions of the jet multiplicity, the dijet invariant mass, and the four-lepton transverse momenta in different $N_{\text {jets}}$ bins, which are used for the categorisation of reconstructed events, are shown in Fig. 6 for different steps of the event categorisation.

The expected numbers of signal and background events in each reconstructed event category are shown in Table 7 together with the corresponding observed number of events. The expected event yields are in good agreement with the observed ones. The observed and expected (post-fit) distributions of the NN discriminants are shown in Fig. 7 and in Figure 8. In addition, Fig. 8g, h show the observed and expected yields in the categories where no NN discriminant is used and in the mass sidebands used to constrain the $ZZ^{*}$ and tXX background, respectively. All distributions are in good agreement with the data.

The statistical interpretation of the results and compatibility with the SM are discussed in the following.

Fig. 6

The observed and expected distributions (post-fit) of (a) the jet multiplicity $N_{\text {jets}}$ after the inclusive event selection, the four-lepton transverse momenta $p_{\mathrm {T}}^{4\ell }$ for events with (b) exactly zero jets, (c) with exactly one jet and (d) with at least two jets and (e) the dijet invariant mass $m_{jj}$ for events with at least two jets. The SM Higgs boson signal is assumed to have a mass $m_{H}$ = 125 $\text {Ge} \text {V}$. The uncertainty in the prediction is shown by the hatched band, calculated as described in Sect. 7. For comparison only, the hatched band includes the theoretical uncertainties in the SM cross-section for the signal and the background processes

Table 7 The expected (pre-fit) and observed numbers of events for an integrated luminosity of 139 fb$^{-1}$ at $\sqrt{\mathrm {s}}=$ 13 $\text {Te} \text {V}$ in the signal region $115<m_{4\ell }<130$ $\text {Ge} \text {V}$ and sideband region $105<m_{4\ell }<115$ $\text {Ge} \text {V}$ or $130<m_{4\ell }<160$ $\text {Ge} \text {V}$ (350 $\text {Ge} \text {V}$ for tXX-enriched) in each reconstructed event category assuming the SM Higgs boson signal with a mass $m_{H} = 125$ $\text {Ge} \text {V}$. The sum of the expected number of SM Higgs boson events and the estimated background yields is compared with the data. Combined statistical and systematic uncertainties are included for the predictions. Expected contributions that are below 0.2% of the total yield in each reconstructed event category are not shown and replaced by ‘-’

8.2 Measurement of simplified template cross-sections

To measure the product $\vec {\sigma }\cdot$ $\mathcal {B}$ of the Higgs boson production cross-section and the branching ratio for $H \rightarrow ZZ^*$ decay for the production bins of the Production Mode Stage or the Reduced Stage 1.1, a fit to the discriminant observables introduced in Sect. 5.2 is performed using the likelihood function ${{{\mathcal {L}}}}( \vec {\sigma }, \vec {\theta })$ that depends on the Higgs boson production cross-section $\vec {\sigma } = \{ \sigma _{1}, \sigma _{2}, \ldots , \sigma _{N} \}$ where $\sigma _p$ is the cross-section in each production bin p and the nuisance parameters $\vec {\theta }$ accounting for the systematic uncertainties. The likelihood function is defined as a product of conditional probabilities over binned distributions of the discriminating observables in each reconstructed signal and sideband event category j,

$$\begin{aligned} \begin{aligned} {\mathcal {L}}(\vec {\sigma }, \vec {\theta })&= \displaystyle \prod _{j}^{N_{\mathrm {categories}}} \prod _{i}^{N_{\mathrm {bins}}} P\left( N_{i,j} ~| ~L\cdot \vec {{\sigma }} \cdot \mathcal {B} \cdot \vec {A}_{i, j}(\vec {\theta }) + B_{i,j}(\vec {\theta }) \right) \\&\quad \times \prod _{m}^{N_{\mathrm {nuisance}}} {\mathcal {C}}_{m}(\vec {\theta }) \;, \end{aligned} \end{aligned}$$

(3)

with Poisson distributions P corresponding to the observation of $N_{i,j}$ events in each histogram bin i of the discriminating observable given the expectations for each background process, $B_{i,j}(\vec {\theta })$, and for the signal, $S_{i,j}(\vec {\theta }) = L \cdot \vec {{\sigma }} \cdot \mathcal {B} \cdot \vec {A}_{i,j}(\vec {\theta })$, where L is the integrated luminosity and $\vec {A}_{i,j} = \{ A_{i,j}^{1}, A_{i,j}^{2}, \ldots , A_{i,j}^{N} \}$ is the set of signal acceptances from each production bin. The signal acceptance $A_{i,j}^{p}$ is defined as the fraction of generated signal events in the production bin p that satisfy the event reconstruction and selection criteria in the histogram bin i of the reconstructed event category j. For a given production bin p, the acceptance consists of $A_{i,j}^{p} = a^{p} \cdot \epsilon _{i,j}^{p}$, where $a^{p}$ is the particle-level acceptance in the fiducial region defined from requirements listed in Sects. 4 and 5 and $\epsilon _{i,j}^{p}$ is the reconstruction efficiency of these particle-level events. Constraints on the nuisance parameters corresponding to systematic uncertainties described in Sect. 7 are represented by the functions ${\mathcal {C}}_{m}(\vec {\theta })$. The cross-sections are treated as independent parameters for each production bin and correlated among the different reconstructed event categories. The test statistic used to perform the measurements is the ratio of profile likelihoods [170],

$$\begin{aligned} q(\vec {\sigma }) = -2 \ln \frac{{\mathcal {L}}(\vec {\sigma }, \hat{\hat{\vec {\theta }}}(\vec {\sigma })) }{{\mathcal {L}}(\hat{{\vec {\sigma }}}, \hat{{\vec {\theta }}} )} = -2 \ln \lambda (\vec {\sigma })\;, \end{aligned}$$

where $\vec {\sigma }$ represents only the cross-section(s) considered as parameter(s) of interest in a given fit. The likelihood in the numerator is the estimator of a conditional fit, i.e. with parameter(s) of interest ${\sigma _i}$ fixed to a given value, while the remaining cross-sections and nuisance parameters are free-floating parameters in the fit. The values of the nuisance parameters $\hat{\hat{\vec {\theta }}}(\vec {\sigma }))$ maximise the likelihood on the condition that the parameters of interest are held fixed to a given value. The likelihood in the denominator is the estimator of an unconditional fit in which all $\vec {\sigma }$ and $\vec {\theta }$ parameters are free parameters of the fit. The parameter of interest $\sigma$ in each production bin is alternatively replaced by $\mu \cdot \sigma _{\mathrm {SM}}(\vec {\theta })$, allowing an interpretation in terms of the signal strength $\mu$ relative to the SM prediction $\sigma _{\mathrm {SM}}(\vec {\theta })$.

Assuming that the relative signal fractions in each production bin are given by the predictions for the SM Higgs boson, the inclusive $H \rightarrow ZZ^{*}$ production cross-section for $|y_H|<2.5$ is measured to be:

$$\begin{aligned} \sigma \cdot \mathcal {B}\equiv & {} \sigma \cdot \mathcal {B}(H\rightarrow ZZ^*) \\= & {} 1.34 \pm 0.11(\mathrm {stat.}) \pm 0.04(\mathrm {exp.})\pm 0.04(\mathrm {th.})~\mathrm {pb} \\= & {} 1.34 \pm 0.12~ \mathrm {pb}, \end{aligned}$$

where the uncertainties are either statistical (stat.) or of experimental (exp.) or theoretical (th.) systematic nature.

The SM prediction is $(\sigma \cdot \mathcal {B})_{\mathrm {SM}}$  $\equiv$ $(\sigma \cdot \mathcal {B}(H\rightarrow ZZ^*))_{\mathrm {SM}}$ $= 1.33 \pm 0.08$ pb. The data are also interpreted in terms of the global signal strength, yielding

$$\begin{aligned} \mu= & {} 1.01\pm 0.08(\mathrm {stat.}) \pm 0.04(\mathrm {exp.})\pm 0.05(\mathrm {th.})\\= & {} 1.01\pm 0.11. \end{aligned}$$

The measured cross-section and signal strength are in an excellent agreement with the SM prediction, with a p-value of 98.6% for both compatibility tests.

The corresponding likelihood functions are shown in Fig. 9. The dominant systematic uncertainty in the cross-section measurement is the experimental uncertainty in the lepton efficiency and integrated luminosity measurements and theoretical uncertainties related to parton shower modelling affecting the acceptance. The signal-strength measurement is also affected by the theoretical uncertainty in the $\mathrm {ggF}$ cross-section due to missing higher-order corrections in QCD.

The expected SM cross-section, the observed values of $\sigma \cdot \mathcal {B}(H\rightarrow ZZ^*)$ and their ratio for the inclusive production and in each production bin of the Production Mode Stage and the Reduced Stage 1.1 are shown in Table 8.

The corresponding values are summarised in Fig. 10. In the ratio calculation, uncertainties in the SM expectation are not taken into account. The Production Mode Stage and Reduced Stage-1.1 measurements agree with the predictions for the SM Higgs boson. The p-values of the corresponding compatibility tests are 91% and 77%, respectively.

Fig. 7

The observed and expected NN output (post-fit) distributions for an integrated luminosity of 139 fb$^{-1}$ at $\sqrt{\mathrm {s}}=$ 13 $\text {Te} \text {V}$ in the different zero- and one-jet categories: (a) $\mathrm {NN_{ggF}}$ in $0j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Low}$, (b) $\mathrm {NN_{ggF}}$ in $0j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Med}$, (c) $\mathrm {NN_{VBF}}$ $1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Low}$ with $\mathrm {NN}_{ZZ}<0.25$, (d) $\mathrm {NN}_{ZZ}$ in $1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Low}$ with $\mathrm {NN}_{ZZ}>0.25$, (e) $\mathrm {NN_{VBF}}$ in $1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Med}$ with $\mathrm {NN}_{ZZ}<0.25$, (f) $\mathrm {NN}_{ZZ}$ in $1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {Med}$ with $\mathrm {NN}_{ZZ}>0.25$ and (g) $\mathrm {NN_{VBF}}$ in $1j$-$p_{\mathrm {T}}^{4\ell }$-$\mathrm {High}$. The SM Higgs boson signal is assumed to have a mass $m_{H}$ = 125 $\text {Ge} \text {V}$. The uncertainty in the prediction is shown by the hatched band, calculated as described in Sect. 7. For comparison only, the hatched band includes the theoretical uncertainties in the SM cross-section for the signal and the background processes. The bin boundaries are chosen to maximise the significance of the targeted signal in each category

Fig. 8

The observed and expected NN output (post-fit) distributions for an integrated luminosity of 139 fb$^{-1}$ at $\sqrt{\mathrm {s}}=$ 13 $\text {Te} \text {V}$ in the different categories: (a) $\mathrm {NN_{VBF}}$ in $2j$ with $\mathrm {NN}_{VH}<0.2$, (b) $\mathrm {NN}_{VH}$ in $2j$ with $\mathrm {NN}_{VH}>0.2$, (c) $\mathrm {NN_{VBF}}$ in $2j$-BSM-like, (d) $\mathrm {NN}_{ttH}$ in ttH-Had-enriched with $\mathrm {NN}_{tXX}<0.4$, (e) $\mathrm {NN}_{tXX}$ in ttH-Had-enriched with $\mathrm {NN}_{tXX}>0.4$ and (f) $\mathrm {NN}_{ttH}$ in VH-Lep-enriched. g Shows the categories where no NN discriminant is used while (h) shows the sidebands used to constrain the $ZZ^{*}$ and tXX backgrounds. The SM Higgs boson signal is assumed to have a mass $m_{H}$ = 125 $\text {Ge} \text {V}$. The uncertainty in the prediction is shown by the hatched band, calculated as described in Sect. 7. For comparison only, the hatched band includes the theoretical uncertainties in the SM cross-section for the signal and the background processes. The bin boundaries are chosen to maximise the significance of the targeted signal in each category

Fig. 9

Observed profile likelihood as a function of (a) $\sigma \cdot \mathcal {B}(H\rightarrow ZZ^*)$ normalised by the SM expectation and (b) the inclusive signal strength $\mu$; the scans are shown both with (solid line) and without (dashed line) systematic uncertainties

Table 8 The expected SM cross-section $(\sigma \cdot \mathcal {B})_{\mathrm {SM}}$, the observed value of $\sigma \cdot \mathcal {B}$, and their ratio $(\sigma \cdot \mathcal {B}) / (\sigma \cdot \mathcal {B})_{\mathrm {SM}}$ for the inclusive production and for each Production Mode Stage and Reduced Stage-1.1 production bin for the $H\rightarrow ZZ^*$ decay for an integrated luminosity of 139 fb$^{-1}$ at $\sqrt{\mathrm {s}}=13$ $\text {Te} \text {V}$. The $bbH$ (tH) contribution is included in the $\mathrm {ggF}$ ($ttH$) production bins. The uncertainties are given as (stat.)+(exp.)+(th.) for the inclusive cross-section and the Production Mode Stage, and as (stat.)+(syst.) for the Reduced Stage 1.1. The Reduced Stage-1.1 results are dominated by the statistical uncertainty and the impact of theory uncertainties is smaller than for the Production Mode Stage. The impact of the theory uncertainties for the Reduced Stage 1.1 is smaller than the least significant digit

Fig. 10

The observed and expected SM values of the cross-sections $\sigma \cdot \mathcal {B}$ normalised by the SM expectation $(\sigma \cdot \mathcal {B})_{\mathrm {SM}}$ for (a) the inclusive production and in the Production Mode Stage and (c) the Reduced Stage-1.1 production bins for an integrated luminosity of 139 fb$^{-1}$ at $\sqrt{\mathrm {s}}=13$ $\text {Te} \text {V}$. The fitted normalisation factors for the ZZ and tXX background are shown in the inserts. Different colours indicate different Higgs boson production modes (or background sources). The vertical band represents the theory uncertainty in the signal prediction. The correlation matrices between the measured cross-sections and the ZZ and tXX normalisation factors are shown for (b) the Production Mode Stage and (d) the Reduced Stage 1.1

For the $\mathrm {qq2Hqq}$-$\mathrm {VBF}$ bin, most of the sensitivity to the $\mathrm {VBF}$ production mode comes from the phase space with $m_{jj} > 350$ $\text {Ge} \text {V}$ and $p_{\text {T}} ^H < 200$ $\text {Ge} \text {V}$. To probe the $\mathrm {VBF}$ contribution more directly, the cross-sections in this and in the remaining phase space region of the $\mathrm {qq2Hqq}$-$\mathrm {VBF}$ bin are fitted separately to the data, simultaneously with the other Reduced Stage 1.1 bins, using the reconstruction categories described in Sect. 5. The cross-section in the $m_{jj} > 350$ $\text {Ge} \text {V}$ and $p_{\text {T}} ^H < 200$ $\text {Ge} \text {V}$  phase space is measured to be $0.060^{+0.025}_{-0.020}\ \text {pb}$ compared with the predicted cross-section of $0.0335 ^{+0.0007} _{-0.0011}\ \text {pb}$. This measurement has a correlation of 20% with the measurement in the $\mathrm {gg2H}$-$2j$ bin, while correlations with other bins are up to 50%.

The dominant contribution to the measurement uncertainty in the $\mathrm {ggF}$ Production Mode Stage bin originates from the same sources as in the inclusive measurement. For the $\mathrm {VBF}$ production bin, the dominant systematic uncertainties are related to parton showering modelling and jet energy scale and resolution uncertainties. The $\mathrm {VBF}$, VH and $ttH$ production bins are also affected by the theoretical uncertainties related to the modelling of the $\mathrm {ggF}$ process. For the Reduced Stage-1.1 bins, the dominant cross-section uncertainties are the jet energy scale and resolution, and parton shower uncertainties.

Figure 11 shows the likelihood contours in the ($\mathrm {ggF}$, $\mathrm {VBF}$), ($\mathrm {ggF}$, VH), ($\mathrm {VBF}$, VH) and ($\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {Low}$, $\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {High}$) planes. The other cross-section parameters are left free in the fit, i.e. they are not treated as parameters of interest. The compatibility with the SM expectation is at the level of 0.22, 0.25, 0.19 and 0.33 standard deviations, respectively.

Fig. 11

Likelihood contours at 68% CL (dashed line) and 95% CL (solid line) in the (a) ($\mathrm {ggF}$, $\mathrm {VBF}$), (b) ($\mathrm {ggF}$, VH), (c) ($\mathrm {VBF}$, VH) and (d) ($\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {Low}$, $\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {High}$) plane. The SM prediction is shown together with its theory uncertainty (filled ellipse). The VH parameter of interest is constrained to positive values

9 Constraints on the Higgs boson couplings in the $\kappa$-framework

The cross-sections measured at the Production Mode Stage are interpreted in the $\kappa$-framework described in Sect. 1.2. The relevant cross-sections and the branching ratio of Eq. (3) are parameterised in terms of the coupling-strength modifiers $\vec {\kappa }$. One interesting benchmark allows two different Higgs boson coupling-strength modifiers to fermions and bosons, reflecting the different structure of the interactions of the SM Higgs sector with gauge bosons and fermions. The universal coupling-strength modifiers $\kappa _{F}$ for fermions and $\kappa _{V}$ for vector bosons are defined as $\kappa _{V}= \kappa _{W}=\kappa _{Z}$ and $\kappa _{F}=\kappa _{t}=\kappa _{b}=\kappa _{c}=\kappa _{\tau }=\kappa _{\mu }$. It is assumed that there are no undetected or invisible Higgs boson decays. The observed likelihood contours in the $\kappa _{V}$–$\kappa _{F}$ plane are shown in Fig. 12 (only the quadrant $\kappa _{F}>0$ and $\kappa _{V}>0$ is shown since this channel is not sensitive to the relative sign of the two coupling modifiers). The best-fit value is $\hat{\kappa }_V=1.02\pm 0.06$ and $\hat{\kappa }_F=0.88\pm 0.16$, with the correlation of −0.17. The probability of compatibility with the Standard Model expectation is at the level of 75%.

10 Constraints on the tensor coupling structure in the EFT approach

To interpret the observed data in the framework of an effective field theory, an EFT signal model is built by parameterising the production cross-sections in each production bin of the Reduced Stage 1.1, as well as the branching ratio and the signal acceptances, as a function of the SMEFT Wilson coefficients introduced in Sect. 1.3. The constraints on the Wilson coefficients are then obtained from the simultaneous fit to the data in all reconstructed signal and sideband event categories. Due to the statistical precision of the data sample, the constraints are always set on one or at most two of the Wilson coefficients at a time, while the values of the remaining coefficients are assumed to be equal to zero.

10.1 EFT signal model

The EFT parameterisation of the production cross-sections in each production bin of the Reduced Stage 1.1 is obtained from Eq. (1) using simulated BSM samples introduced in Sect. 3. The contribution from the $gg\rightarrow Z(\rightarrow \ell \ell )H$ process is taken from the SM simulation and assumed to scale with BSM parameters in the same way as the $qq\rightarrow Z(\rightarrow \ell \ell )H$ processes. As in the case of simplified template cross-section measurements, $ttH$ and tH processes are combined into a single $ttH$ production bin. The cut-off scale is set to $\Lambda = 1$ $\text {Te} \text {V}$. Only LO computation of QCD and SM electroweak processes is provided, with LO effective couplings for the SM Higgs boson to gluon and to photon vertices. An assumption is made that higher-order corrections, applied in a multiplicative way, are the same for both the SM and the BSM LO predictions and therefore no changes in the parameterisation are expected due to higher-order effects [171]. With the current amount of data, the constraints from the $\mathrm {VBF}$, VH and $ttH$ production modes on the relevant Wilson coefficients still allow a rather large range of parameter values in which the quadratic term (the last term in Eq. (1)) cannot be neglected even though its contribution is suppressed by $\Lambda ^4$. Such dimension-six quadratic terms are therefore included in the EFT parameterisation. Since the linear terms from dimension-eight operators are suppressed by the same factor, they could in general also give similar non-negligible contributions. Dimension-eight terms are currently not available in the SMEFT model and are thus not taken into account.

Fig. 12

Likelihood contours at 68% CL (dashed line) and 95% CL (solid line) in the $\kappa _{V}$–$\kappa _{F}$ plane. The best fit to the data (solid cross) and the SM prediction (star) are also indicated

The branching ratio for the $H{\rightarrow }ZZ^* {\rightarrow }4\ell$ decay is parameterised in terms of Wilson coefficients following Eq. (2). The partial and total decay widths are calculated in MadGraph5_aMC@NLO. The total decay width is calculated by taking into account the dominant Higgs boson decay modes: $\gamma \gamma$, $Z\gamma$, $bb$, gg, WW and ZZ. Other decay modes are not affected by the probed Wilson coefficients. Their contribution to the total decay width is therefore given by the corresponding SM predictions.

The selection criteria for the four-lepton Higgs boson candidates, in particular the requirements on the minimum invariant mass $m_{34}$ of the subleading lepton pair, introduce an additional dependence of the signal acceptance on the BSM coupling parameters. The particle-level signal acceptance A, defined as the fraction of signal events satisfying the Higgs boson candidate selection criteria applied at particle-level, has therefore been simultaneously parameterised in terms of the three Wilson coefficients $c_{HW}$, $c_{HB}$ and $c_{HWB}$ ($c_{H{\widetilde{W}}}$, $c_{H{\widetilde{B}}}$ and $c_{H{\widetilde{W}}B}$) assuming that the values of CP-odd (CP-even) parameters vanish. The dependence of the acceptance on other EFT coupling parameters is shown to be negligible as these parameters have negligible or no impact on the $H\rightarrow ZZ^*$ decay. The acceptance correction relative to the SM prediction is described by a three-dimensional Lorentzian function with free acceptance parameters $\alpha _0$, $\alpha _1$, $\alpha _2$, $\beta _i$, $\delta _i$, $\delta _{(i,j)}$ and $\delta _{(i,j,k)}$,

$$\begin{aligned} \frac{A(\vec {c})}{A_{\text {SM}}}= & {} \alpha _0 + (\alpha _1)^2 \cdot \left[ \alpha _2 + \sum \limits _{i}\delta _i\cdot (c_i+\beta _i)^2 \right. \nonumber \\&\left. + \sum \limits _{\begin{array}{c} ij \\ i\ne j \end{array}} \delta _{(i,j)} \cdot c_i c_j + \underset{i\ne j\ne k}{\delta _{(i,j,k)}} \cdot c_{i} c_{j} c_{k}\right] ^{-1}, \end{aligned}$$

(4)

where indices i, j and k run over $({\textit{HW}}, {\textit{HB}}, {\textit{HWB}})$ in case of the acceptance correction for the set of CP-even parameters and over $({H{\widetilde{W}}}, {H{\widetilde{B}}}, {H{\widetilde{W}}}B)$ in case of the CP-odd parameters. A common parameterisation is used for all production bins since the differences between production bins are shown to be negligible. In addition, the reconstructed event categorisation criteria imposed on the selected Higgs boson candidates and the classification in bins of multivariate NN discriminant values do not impact the acceptance parameterisation. The impact of reconstruction efficiencies on the parameterisation is also negligible, such that Eq. (4) also holds for the ratio $A(\vec {c})/A_\text {SM}$ of reconstruction-level acceptances defined in Sect. 8. The resulting acceptance parameterisation curves are shown in Fig. 13 for the cases in which all but one of the Wilson coefficients are set to zero. For all cases, the acceptance correction is equal to one at the SM point. In the case of the $c_{HW}$ and $c_{HWB}$ Wilson coefficients, the acceptance corrections reach a maximum value slightly larger than one, leading to the shift of the maximum position from the SM point. This shift is compatible with the statistical accuracy of the fit and the impact of linear EFT terms which are not symmetric around the SM point.

The final parameterisation of signal yields relative to the SM prediction in each production bin of the Reduced Stage 1.1 is obtained as the product of the corresponding cross-section, branching ratio and acceptance parameterisations. The expected event yields normalised to the SM prediction are shown in Fig. 14 for each of the CP-even Wilson coefficients after setting all other coefficients to zero. Only production bins with the highest sensitivity to a given Wilson coefficient are shown. The impact of the quadratic terms in the EFT parameterisation can clearly be seen as a non-linear dependence on all but the $c_{HG}$ Wilson coefficient. For comparison, the predictions without the acceptance corrections ($\sigma \cdot \mathcal {B}$), and without both the acceptance and branching ratio corrections ($\sigma$) are also shown. Both the acceptance and the branching ratio parameterisations have a strong impact on the sensitivity to different Wilson coefficients, especially for the $c_{HW}$, $c_{HB}$ and $c_{HWB}$ parameterisations in gg2H production bins (Fig. 14a, b, c). Since these coefficients do not enter the $\mathrm {ggF}$ production vertex, the corresponding sensitivity is entirely driven by their impact on the decay and the acceptance of selected signal events. The acceptance corrections significantly degrade the sensitivity to the $c_{HW}$ coefficient (see Fig. 14a). Additional sensitivity to this coefficient can be gained from the qq2Hqq production bins as shown in Fig. 14d. The Wilson coefficients $c_{HG}$ and $c_{uH}$, on the other hand, do not affect the acceptance since they are not present in the decay vertex (Fig. 14e, f). The coefficient $c_{HG}$ still has a non-vanishing impact on the branching ratio through its contributions to the total decay width. Similar effects are also seen for the Wilson coefficients of CP-odd operators.

10.2 EFT interpretation results

The ratios of the expected signal yield for a chosen EFT parameter value to its SM prediction are shown in Fig. 15 in each production bin of the Reduced Stage 1.1, together with the corresponding measurement.

Fig. 13

The dependence of the signal acceptance normalised to the SM acceptance on the Wilson coefficients (a) $c_{HW}$ and $c_{H\tilde{W}}$, (b) $c_{HB}$ and $c_{H\tilde{B}}$, (c) $c_{HWB}$ and $c_{H\tilde{W}B}$ after setting all other coefficients to zero

Fig. 14

The expected event yields ($\sigma \cdot B \cdot A$) relative to the SM prediction as a function of the Wilson coefficient (a) $c_{HW}$, (b) $c_{HB}$ and (c) $c_{HWB}$ in the $\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {High}$ production bin, (d) $c_{HW}$ in the $\mathrm {qq2Hqq}$-$\mathrm {VBF}$ production bin, (e) $c_{HG}$ in the $\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {High}$ production bin and (f) $c_{uH}$ in the $ttH$ production bin. The dependence on only one Wilson coefficient is shown on each plot while setting all others to zero. For comparison, the predictions are also shown for the parameterisation without the acceptance corrections ($\sigma \cdot B$) and for the production cross-section only ($\sigma$) without the acceptance and the branching ratio corrections. The $\sigma$ parameterisations in (a), (b) and (c) coincide with the SM expectation at 1 as the coefficients $c_{HW}$, $c_{HB}$ and $c_{HWB}$ are not present in the $\mathrm {ggF}$ production vertex. Since the acceptance does not depend on the $c_{HG}$ and $c_{uH}$ parameters, no corresponding ($\sigma \cdot B \cdot A$) expectation is shown in (e) and (f). Similarly, no ($\sigma \cdot B$) expectation is shown in (f), since the $c_{uH}$ parameter has a negligible impact on the branching ratio. The bands indicate the expected precision of the cross-section measurement in a given production bin at the one standard deviation level

Fig. 15

The expected signal yield ratio for chosen (a) CP-even and (b) CP-odd EFT parameter values together with the corresponding cross-section measurement in each production bin of Reduced Stage 1.1. The parameter values correspond approximately to the expected confidence intervals at the 68% CL obtained from the statistical interpretation of data

Fig. 16

The observed and expected values of SMEFT Wilson coefficients from (a) CP-even and (b) CP-odd operators obtained for an integrated luminosity of $139$ fb$^{-1}$ at $\sqrt{s}=13$ $\text {Te} \text {V}$. Only one Wilson coefficient is fitted at a time while all others are set to zero. The values for the $c_{HG}$ and $c_{H\tilde{G}}$ coefficients are scaled by a factor of 100, and for the $c_{uH}$ and $c_{\tilde{u} H}$ coefficients by a factor of 0.05. The horizontal bands represent the expected measurement uncertainty

Table 9 The expected and observed confidence intervals at 68% and 95% CL on the SMEFT Wilson coefficients for an integrated luminosity of 139 $\hbox {fb}^{-1}$ at $\sqrt{s}=13 \ \text {Te} \text {V}$. Only one Wilson coefficient is fitted at a time while all others are set to zero

The EFT parameterisation of signal yields is implemented in the likelihood function of Eq. (3) using the BSM-dependent signal-strength parameters $\mu ^p(\vec {c})$ for each given production bin p,

$$\begin{aligned} \mu ^p(\vec {c}) = \frac{\sigma ^p(\vec {c})}{\sigma _{\text {SM}}} \cdot \frac{\mathcal {B} ^{4\ell }(\vec {c})}{\mathcal {B} ^{4\ell }_{\text {SM}}} \cdot \frac{A(\vec {c})}{A_{\text {SM}}}.&\end{aligned}$$

This is then fitted to the observed event yields. Default SM predictions at the highest available order are employed for the cross-sections and branching ratios multiplying the signal strengths in the likelihood function. Modifications of background contributions due to EFT effects are not taken into account.

The fit results with only one Wilson coefficient fitted at a time are summarised in Fig. 16 and in Table 9. The results are in good agreement with the SM predictions. The measurements are dominated by the statistical uncertainty. In the case of the CP-odd coupling parameters, each fit gives two degenerate minima since the corresponding EFT parameterisation contains only quadratic terms which are not sensitive to the sign of the fitted parameter. The fit of the CP-even coupling parameter $c_{uH}$ also results in two minima since the corresponding EFT parameterisation curve in the only sensitive $ttH$ production bin crosses the expected SM cross-section value at two different values of the $c_{uH}$ parameter (see Fig. 14f). The same is true also for the observed $ttH$ cross-section. The small degeneracies for other CP-even coupling parameters are removed by the combination of several sensitive production bins.

Fig. 17

Expected (dashed line) and observed (full line) 2D-fit likelihood curves at the 95% CL for the SMEFT Wilson coefficients of CP-even operators at an integrated luminosity of $139$ fb$^{-1}$ and $\sqrt{s}=13$ $\text {Te} \text {V}$. The best fit to the data (solid cross) and the SM prediction (star) are also indicated. Except for the two fitted Wilson coefficients, all others are set to zero

Fig. 18

Expected (dashed line) and observed (full line) 2D-fit likelihood curves at the 95% CL for the SMEFT Wilson coefficients of CP-odd operators at an integrated luminosity of $139$ fb$^{-1}$ and $\sqrt{s}=13$ $\text {Te} \text {V}$. The best fit to the data (solid cross) and the SM prediction (star) are also indicated. Except for the two fitted Wilson coefficients, all others are set to zero

The strongest constraint, driven mostly by the $\mathrm {ggF}$ reconstructed event categories, is obtained on the $c_{HG}$ coefficient related to the CP-even Higgs boson interactions with gluons. The highest sensitivity to this parameter is reached by the measurements in the $\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {Low}$ and $\mathrm {gg2H}$-$0j$-$p_{\mathrm {T}}^H$-$\mathrm {High}$ production bins due to the highest statistical precision. The sensitivity in the $\mathrm {gg2H}$-$p_{\mathrm {T}}^H$-$\mathrm {High}$ production bin, which is designed to target the BSM physics effects, is limited due to the small number of events observed in the corresponding reconstructed event category. Additional sensitivity in this bin may be provided by the two-loop interactions which are not implemented in the current simulation of the ggH vertex. The constrained range is stringent enough for the linear approximation to hold, i.e. the quadratic terms in the signal parameterisation are small compared with the linear ones (see Fig. 14e). The constraint on the $c_{H{\widetilde{G}}}$ parameter of the related CP-odd operator is worse by about a factor of three since the linear terms from CP-odd operators do not contribute to the total production cross-section. The constraints on the remaining EFT parameters are weaker, such that both the CP-even and CP-odd signals become dominated by the quadratic terms and are therefore comparable in size. The next-strongest constraints are obtained on the $c_{HB}$, $c_{HWB}$, $c_{HW}$, $c_{H{\widetilde{B}}}$, $c_{H{\widetilde{W}}B}$ and $c_{H{\widetilde{W}}}$ coefficients that mostly affect the $H\rightarrow ZZ^*$ decays. Due to the larger number of events in the 0-jet reconstructed event categories, the corresponding gg2H production bins provide the highest sensitivity to these decays. Additional smaller sensitivity is obtained from the production vertex of the $\mathrm {VBF}$ and VH production modes, with the dominant contribution from $\mathrm {qq2Hqq}$-$\mathrm {VBF}$ and $\mathrm {qq2Hqq}$-$\mathrm {BSM}$ bins. The latter one is designed to enhance the sensitivity to BSM physics. The qq2Hqq production bins improve in particular the sensitivity to the $c_{HW}$ and $c_{H{\widetilde{W}}}$ parameters that is otherwise significantly degraded by the acceptance corrections. Finally, looser constraints are set on the top-Yukawa coupling parameters $c_{uH}$ and $c_{{\widetilde{u}}H}$, driven by the measurements in the $ttH$ production bin.

To explore possible correlations between different Wilson coefficients, the simultaneous fits are also performed on two Wilson coefficients at a time. The corresponding results are shown in Fig. 17 for several combinations of two CP-even EFT parameters and in Fig. 18 for the corresponding CP-odd operators. The best-fit values as well as the deviation from the SM prediction are shown in Table 10. Good agreement with the SM predictions is observed for all such possible combinations.

Table 10 The best-fit values and the corresponding deviation from the SM prediction obtained from the two-dimensional likelihood scans of the CP-odd BSM coupling parameters performed with 139 $\hbox {fb}^{-1}$ of data at a centre-of-mass energy of $\sqrt{s}=13$ $\text {Te} \text {V}$. The limits are computed using the confidence-level interval method. Except for the two fitted BSM coupling parameters, all others are set to zero

The anti-correlation between the $c_{HW}$ and $c_{HB}$ coefficients, as well as between $c_{H{\widetilde{W}}}$ and $c_{H{\widetilde{B}}}$, is driven by their impact on the signal acceptance. The non-ellipsoidal shape is caused by the acceptance correction, which degrades the original branching ratio-driven sensitivity for increasing parameter values, in particular in the case of the $c_{HW}$ ($c_{H{\widetilde{W}}}$) coefficient. The sensitivity is, however, partially recouped by the $\mathrm {VBF}$ production vertex.

The ‘V’-shaped correlation between the $c_{HG}$ and $c_{HB}$ parameters is due to the interplay between the EFT parameterisation in the $\mathrm {ggF}$ production vertex and the parameterisation of the branching ratios and acceptances. The $\mathrm {ggF}$ production vertex provides the constraint on the $c_{HG}$ parameter alone, independently of $c_{HB}$. Due to the decay vertex with its acceptance corrections, this constrained range is shifted upward with increasing values of $c_{HB}$. Close to the SM point, the constrained $c_{HG}$ range remains approximately the same as without the decay constraints. An additional constraint on $c_{HB}$ is provided by the $\mathrm {VBF}$ production mode. Around the SM point, the $c_{HB}$ constraints correspond approximately to those from the one-dimensional parameter fit. Additional sensitivity to intermediate values of the $c_{HB}$ parameter is provided by the acceptance corrections, resulting in two additional allowed parameter regions that are disjoint from the region around the SM point. Similar arguments hold also for the CP-odd case with the $c_{H{\widetilde{G}}}$ and $c_{H{\widetilde{B}}}$ parameters. As opposed to the CP-even case, however, the likelihood contours are symmetric around the $c_{H{\widetilde{G}}}=0$ axis, since there are no linear terms contributing to the $\mathrm {ggF}$ production cross-section.

The correlation between the $c_{HG}$ and $c_{uH}$ ($c_{H{\widetilde{G}}}$ and $c_{{\widetilde{u}}H}$) parameters is introduced through the interference term in the $ttH$ vertex. However, the impact of this term on the final result is negligible since the $c_{HG}$ ($c_{H{\widetilde{G}}}$) parameter is already constrained to very small values compared with $c_{uH}$ ($c_{{\widetilde{u}}H}$). Therefore, the $ttH$ production vertex mainly constrains the $c_{uH}$ and $c_{{\widetilde{u}}H}$ parameters, while the $\mathrm {ggF}$ vertex constrains only the other two. The acceptance correction has no impact on these results. The CP-odd parameter range is less constrained than the CP-even one due to the missing linear $c_{{\widetilde{u}}H}$ terms in the cross-section parameterisation.

11 Conclusion

Higgs boson properties are studied in the four-lepton decay channel using $139$ fb$^{-1}$ of LHC proton–proton collision data at $\sqrt{s}=13$ $\text {Te} \text {V}$ collected by the ATLAS experiment. The Higgs boson candidate events are categorised into several topologies, providing sensitivity to different production modes in various regions of phase space. Additional multivariate discriminants are used to further improve the sensitivity in reconstructed event categories with a sufficiently large number of events. The cross-section times branching ratio for $H\rightarrow ZZ^*$ decay measured in dedicated production bins are in good agreement with the SM predictions. The inclusive cross-section times branching ratio for $H\rightarrow ZZ^*$ decay in the Higgs boson rapidity range of $|y_H|<2.5$ is measured to be $1.34 \pm 0.12$ pb compared with the SM prediction of $1.33\pm 0.08$ pb. Results are also interpreted within the $\kappa$-framework with coupling-strength modifiers $\kappa _V$ and $\kappa _F$, showing compatibility with the SM. Based on the product of cross-section, branching ratio and acceptance measured in Reduced Stage-1.1 production bins of simplified template cross-sections, constraints are placed on possible CP-even and CP-odd BSM interactions of the Higgs boson to vector bosons, gluons and top quarks within an effective field theory framework in the $H\rightarrow ZZ^*$ decay. The data are found to be consistent with the SM hypothesis.