1 Introduction

Measurements of diboson production in proton–proton (\({\mathrm{p}} {\mathrm{p}} \)) collisions, such as Z boson pair (\({\mathrm{Z}} {\mathrm{Z}} \)) production, at the CERN LHC allow precision tests of the standard model (SM). In the SM, \({{\mathrm{Z}} {\mathrm{Z}} }\) production proceeds mainly through quark-antiquark t- and u-channel scattering diagrams. In calculations at higher orders in quantum chromodynamics (QCD), gluon-gluon fusion also contributes via box diagrams with quark loops. There are no tree-level contributions to \({{\mathrm{Z}} {\mathrm{Z}} }\) production from triple gauge boson vertices in the SM. Anomalous triple gauge couplings (aTGC) \({{\mathrm{Z}} {\mathrm{Z}} {\mathrm{Z}} }\) and \({{\mathrm{Z}} {\mathrm{Z}} \gamma }\) are introduced using an effective Lagrangian following Ref. [1]. In this parametrization, two \({{\mathrm{Z}} {\mathrm{Z}} {\mathrm{Z}} }\) and two \({{\mathrm{Z}} {\mathrm{Z}} \gamma }\) couplings are allowed by the electromagnetic gauge invariance and Lorentz invariance for on-shell \({{\mathrm{Z}}}\) bosons and are parametrized by two CP-violating (\({f_4^{{\mathrm{V}}}}\)) and two CP-conserving (\({f_5^{{\mathrm{V}}}}\)) parameters, where \({{{\mathrm{V}}} = ({\mathrm{Z}}, \gamma )}\). Nonzero aTGC values could be induced by new physics models such as supersymmetry [2]. The results can be also expressed in terms of parameters calculated within the effective field theory (EFT) framework, per convention used in Ref. [3] and references therein. In contrast to the anomalous couplings of electroweak (EW) vector bosons, the EFT framework allows an unambiguous calculation of loop effects and provides a simpler interpretation of the results than the aTGC framework.

Previous measurements of the production cross section for pairs of on-shell \({{\mathrm{Z}}}\) bosons at the LHC were performed by the CMS Collaboration with data sets corresponding to integrated luminosities of 5.1\(\,\text {fb}^{-1}\) at \(\sqrt{s} = 7\,\text {TeV} \) [4] and 19.6\(\,\text {fb}^{-1}\) at \(\sqrt{s} = 8\,\text {TeV} \) [5, 6] in the \({{\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\ell ^{\prime \prime }}\) and \({{\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\nu }\) decay channels, where \({\ell = {\mathrm{e}}}\) or \({{\upmu }}\) and \({\ell ^{\prime \prime } = {\mathrm{e}}, {\upmu }}\), or \({{\uptau }}\); and with an integrated luminosity of 2.6\(\,\text {fb}^{-1}\)  [7] and 35.9\(\,\text {fb}^{-1}\)  [8] at \(\sqrt{s} = 13\,\text {TeV} \) in the \({{\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\ell ' }\) decay channel, where \({\ell ' = {\mathrm{e}}}\) or \({{\upmu }}\). All of the results agree with SM predictions. The ATLAS Collaboration reported similar results at \(\sqrt{s} = 7\), 8, and 13\(\,\text {TeV}\)  [9,10,11,12,13,14], which also agree with SM predictions. These measurements are important to test predictions that were recently made available at next-to-next-to-leading order (NNLO) in QCD [15,16,17]. A comparison of these predictions with data for a range of center-of-mass energies provides an insight into the structure of the EW gauge sector of the SM.

This paper reports a study of the \({\mathrm{Z}} {\mathrm{Z}} \) production in the four-lepton decay channel (\({{\mathrm{p}} {\mathrm{p}} \rightarrow 2\ell 2\ell ' }\), where \({2\ell }\) and \({2\ell '}\) indicate pairs of opposite-sign electrons or muons) at \(\sqrt{s} = 13\,\text {TeV} \), with a data set corresponding to an integrated luminosity of \(137{\,\text {fb}^{-1}} \) recorded in 2016–2018. Both \({{\mathrm{Z}}}\) bosons are selected to be on-shell, defined as the mass range 60–120\(\,\text {GeV}\). Fiducial and total cross sections are measured, differential cross sections are presented as a function of different kinematic variables. The invariant mass distribution of the four-lepton system is used to search for anomalous \({{\mathrm{Z}} {\mathrm{Z}} {\mathrm{Z}} }\) and \({{\mathrm{Z}} {\mathrm{Z}} \gamma }\) couplings.

2 The CMS detector

A detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [18].

The central feature of the CMS apparatus is a superconducting solenoid of 6\(\text { m}\) internal diameter, providing a magnetic field of 3.8\(\text { T}\). Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter, which provide coverage in pseudorapidity \({| \eta | < 1.479}\) in a cylindrical barrel and \(1.479< | \eta | < 3.0\) in two endcap regions. Forward calorimeters extend the coverage provided by the barrel and endcap detectors to \(|\eta | < 5.0\). Muons are detected in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid in the range \(|\eta | < 2.4\), with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers.

Electron momenta are estimated by combining energy measurements in the ECAL with momentum measurements in the tracker. The momentum resolution for electrons with transverse momentum \(p_{\mathrm {T}} \approx 45\,\text {GeV} \) from \({{\mathrm{Z}} \rightarrow {{\mathrm{e}}}{}^{+} {{\mathrm{e}}}{}^{-}}\) decays ranges from 1.7% for nonshowering electrons in the barrel region to 4.5% for showering electrons in the endcaps [19]. Matching muons to tracks identified in the silicon tracker results in a \(p_{\mathrm {T}} \) resolution for muons with \(20<p_{\mathrm {T}} < 100\,\text {GeV} \) of 1.3–2.0% in the barrel and better than 6% in the endcaps. The \(p_{\mathrm {T}}\) resolution in the barrel is better than 10% for muons with \(p_{\mathrm {T}}\) up to 1\(\,\text {TeV}\)  [20, 21].

Events of interest are selected using a two-tiered trigger system [22]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100\(\text { kHz}\) within a time interval of less than 4 \(\upmu \)s. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1\(\text { kHz}\) before data storage.

3 Signal and background simulation

Several Monte Carlo (MC) samples are used in the analysis to optimize the selection, calculate the signal efficiency, and estimate background contamination. The pythia  8.226 and 8.230 [23, 24] packages are used for parton showering, hadronization, and the underlying event simulation with the CUETP8M1 tune [25] and the parton distribution function (PDF) NNPDF23_lo_as_0130 [26] for the 2016 data-taking period, and the CP5 tune [27] and the NNPDF 31_nnlo_as_0118 PDF for the 2017 and 2018 data-taking periods.

Signal events are generated with powheg 2.0 [24, 28,29,30,31] at next-to-leading order (NLO) in QCD for quark-antiquark (\({{\mathrm{q}} {{\mathrm{q}}}}\)) processes and leading order (LO) for quark-gluon processes. This includes \({{\mathrm{Z}} {\mathrm{Z}} }\), \({{\mathrm{Z}} \gamma ^*}\), \({{\mathrm{Z}}}\), \(\gamma ^*\gamma ^*\) with a constraint of \({m_{\ell \ell '} > 4\,\text {GeV}}\) applied to all pairs of oppositely charged leptons at the generator level to avoid infrared divergences. The gluon-gluon loop-induced process, \({{\mathrm{g}} {\mathrm{g}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\), is simulated at LO with mcfm  v7.0 [32]. It also includes interference with the SM Higgs off-shell production. The SM Higgs decay is modeled with jhugen 3.1.8 [33,34,35] at LO. The cross sections are scaled to correspond to cross section values calculated at NNLO in QCD for \({{\mathrm{q}} {{\mathrm{q}}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\) [15] (with a K factor of 1.1) and at NLO in QCD for \({{\mathrm{g}} {\mathrm{g}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\) [36] (K factor of 1.7). Electroweak \({\mathrm{Z}} {\mathrm{Z}} \) production in association with two jets is generated with MadGraph  [37] at LO. It amounts to approximately 1% of the total number of \({\mathrm{Z}} {\mathrm{Z}} \) events.

Simulated events for the irreducible background processes containing four prompt leptons in the final state, such as \({{{\mathrm{t}} {}{{\mathrm{t}}}} {\mathrm{Z}}}\), \({{\mathrm{W}} {\mathrm{W}} {\mathrm{Z}}}\), \({{\mathrm{W}} {\mathrm{Z}} {\mathrm{Z}}}\), and \({{\mathrm{Z}} {\mathrm{Z}} {\mathrm{Z}}}\), where the last three are combined and denoted as VVV, are generated with MadGraph 5_amc@nlo v2.4.2 [37] at NLO with zero or one outgoing partons in the matrix element calculation and merged with the parton shower using the FxFx scheme [38]. The same MC is used for \({{\mathrm{W}} {\mathrm{Z}}}\) simulation.

Event samples with aTGC contributions included are generated at LO with sherpa  v2.1.1 [39]. The distributions from the sherpa samples are normalized such that the total yield of the SM sample is the same as that of the powheg+mcfm sample. More details are discussed in Sect. 10.

The detector response is simulated using a detailed description of the CMS detector implemented with the Geant4 package [40]. The reconstruction in simulation and data uses the same algorithms. The simulated samples include additional interactions per bunch crossing, referred to as pileup. The simulated events are weighted so that the pileup distribution matches the data.

Results are also compared to fixed-order predictions produced via the Matrix framework [41], a parton-level MC generator that uses tree and one-loop amplitudes from OpenLoops 2 [42] and two-loop amplitudes from Ref. [43], capable of producing differential predictions at up to NNLO in QCD and NLO in EW, as implemented in Matrix  v2.0.0_beta1 [44]. The calculation is performed with the NNPDF31_nnlo_as_0118_luxqed [45] PDF set with dynamic renormalization (\(\mu _{\mathrm {R}}\)) and factorization scales (\(\mu _{\mathrm {F}}\)) set to the four lepton mass for the differential and fiducial predictions, and with fixed scale set to the nominal Z boson mass for the total cross section. The quark-induced processes are calculated at NNLO in QCD and NLO in EW. The gluon-induced contribution is calculated at NLO in QCD [46]. Photon-induced contributions are also included at up to NLO EW. The calculation uses massless leptons, which leads to a divergence at low dilepton mass. To avoid this divergence, we impose the requirement \({p_{\mathrm {T}} ^{\ell } > 5\,\text {GeV}}\) on the photon-induced component for total cross section predictions. With this condition, the photon-induced contribution is less than 1% of the total production rate. The quark-induced NNLO QCD and NLO EW contributions are combined multiplicatively, and the gluon- and photon-induced contributions are combined additively following the procedure described in Ref. [44]. The predictions reported here are consistent with those published in Refs. [15,16,17].

4 Event reconstruction

Individual particles – electrons, muons, photons, charged and neutral hadrons – in each collision event are identified and reconstructed with the CMS particle-flow (PF) algorithm [47] from a combination of signals from all subdetectors. Reconstructed electrons [19] and muons [20] are considered as lepton candidates if they have \({p_{\mathrm {T}} ^{{\mathrm{e}}} > 7\,\text {GeV}}\) and \({|\eta ^{{\mathrm{e}}} | < 2.5}\) or \({p_{\mathrm {T}} ^{{\upmu }} > 5\,\text {GeV}}\) and \({|\eta ^{{\upmu }} | < 2.4}\).

Lepton candidates are also required to originate from the primary vertex, defined as the reconstructed \({\mathrm{p}} {\mathrm{p}} \) interaction vertex with the largest value of summed physics object \(p_{\mathrm {T}} ^2\). The physics objects used in the primary vertex definition are the objects returned by a jet-finding algorithm [48, 49] applied to all charged tracks associated with the vertex. The distance of closest approach between each lepton track and the primary vertex is required to be less than 0.5\(\text { cm}\) in the plane transverse to the beam axis, and less than 1\(\text { cm}\) in the direction along the beam axis. Furthermore, the significance of the three-dimensional impact parameter relative to the primary vertex, \(\mathrm {SIP_{3D}}\), is required to satisfy \(\mathrm {SIP_{3D}} \equiv | \mathrm {IP} / \sigma _\mathrm {IP} | < 4\) for each lepton, where \(\mathrm {IP}\) is the distance of closest approach of each lepton track to the primary vertex and \(\sigma _\mathrm {IP}\) is its associated uncertainty.

Lepton candidates are required to be isolated from other particles in the event. The relative isolation is defined as

$$\begin{aligned} R_\text {iso}= & {} \left[ \sum _{\begin{array}{c} \text {charged} \\ \text {hadrons} \end{array}} p_{\mathrm {T}} \, + \, \max \left( 0, \sum _{\begin{array}{c} \text {neutral} \nonumber \\ \text {hadrons} \end{array}} \right. \right. p_{\mathrm {T}} \\&\left. \left. + \, \sum _{\text {photons}} \!\! p_{\mathrm {T}} \, - \, p_{\mathrm {T}} ^{\mathrm {PU}} \right) \right] \bigg / p_{\mathrm {T}} ^{\ell }, \end{aligned}$$
(1)

where the sums run over the charged and neutral hadrons, and photons identified by the PF algorithm, in a cone defined by \(\varDelta R \equiv \sqrt{\smash [b]{\left( \varDelta \eta \right) ^2 + \left( \varDelta \varphi \right) ^2}} < 0.3\) around the lepton momentum direction, where \(\varphi \) is the azimuthal angle in radians. To minimize the contribution of charged particles from pileup to the isolation calculation, charged hadrons are included only if they originate from the primary vertex. The contribution of neutral particles from pileup is \(p_{\mathrm {T}} ^{\mathrm {PU}}\). For electrons, \(p_{\mathrm {T}} ^{\mathrm {PU}}\) is evaluated with the “jet area” method described in Ref. [50]; for muons, it is half of the summed \(p_{\mathrm {T}} \) of all charged particles in the cone originating from pileup vertices. The average factor of one half accounts for the expected ratio of neutral to charged particle production in hadronic interactions. A lepton is considered isolated if \(R_\text {iso} < 0.35\).

The lepton reconstruction, identification, and isolation efficiencies are measured with a “tag-and-probe” technique [51] applied to a sample of \({{\mathrm{Z}} \rightarrow \ell ^+\ell ^-}\) data events. The measurements are performed in several bins of \({p_{\mathrm {T}} ^{\ell }}\) and \({|\eta ^\ell |}\). The electron reconstruction and selection efficiency in the ECAL barrel (endcaps) varies from about 85 (77)% at \({p_{\mathrm {T}} ^{{\mathrm{e}}} \approx 10\,\text {GeV}}\) to about 95 (89)% for \({p_{\mathrm {T}} ^{{\mathrm{e}}} \ge 20\,\text {GeV}}\), whereas in the barrel-endcap transition region this efficiency is about 85% averaged over electrons with \({p_{\mathrm {T}} ^{{\mathrm{e}}} > 7\,\text {GeV}}\). The muons are reconstructed and identified with efficiencies above \({\sim }98\%\) within \({|\eta ^{{\upmu }} | < 2.4}\).

5 Event selection

The primary triggers for this analysis require the presence of a pair of loosely isolated leptons of the same or different flavors [22]. The highest \(p_{\mathrm {T}}\) lepton must have \({p_{\mathrm {T}} ^{\ell } > 17\,\text {GeV}}\), and the subleading lepton must have \({p_{\mathrm {T}} ^{{\mathrm{e}}} > 12\,\text {GeV}}\) if it is an electron or \({p_{\mathrm {T}} ^{{\upmu }} > 8\,\text {GeV}}\) if it is a muon. The tracks of the triggering leptons are required to originate within 2 mm of each other in the plane transverse to the beam axis. Triggers requiring a triplet of lower-\(p_{\mathrm {T}}\) leptons with no isolation criteria, or a single high-\(p_{\mathrm {T}}\) electron or muon, are also used. An event is accepted if it passes any trigger regardless of the decay channel. The total trigger efficiency for events within the acceptance of this analysis is greater than 98%.

The four-lepton candidate selection is based on the one used in the recent CMS Higgs boson measurement  [52]. A signal event must contain at least two \({{\mathrm{Z}}/\gamma ^{*}}\) candidates, each formed from an oppositely charged pair of isolated electron or muon candidates. Among the four leptons, the highest \(p_{\mathrm {T}}\) lepton must have \(p_{\mathrm {T}} > 20\,\text {GeV} \), and the second-highest \(p_{\mathrm {T}}\) lepton must have \({p_{\mathrm {T}} ^{{\mathrm{e}}} > 12\,\text {GeV}}\) if it is an electron or \({p_{\mathrm {T}} ^{{\upmu }} > 10\,\text {GeV}}\) if it is a muon. All leptons are required to be separated from each other by \({\varDelta R \left( \ell _1, \ell _2 \right) > 0.02}\), and electrons are required to be separated from muons by \({\varDelta R \left( {\mathrm{e}}, \mu \right) > 0.05}\).

Within each event, all permutations of leptons giving a valid pair of \({{\mathrm{Z}}/\gamma ^{*}}\) candidates are considered separately. Within each four-lepton candidate, the dilepton candidate with an invariant mass closest to 91.2\(\,\text {GeV}\), taken as the nominal \({{\mathrm{Z}}}\) boson mass [53], is denoted \({{\mathrm{Z}} _1}\) and is required to have a mass greater than 40\(\,\text {GeV}\). The other dilepton candidate is denoted \({{\mathrm{Z}} _2}\) and is required to have a mass greater than 4\(\,\text {GeV}\). Both \(m_{Z_1}\) and \(m_{Z_2}\) are required to be less than 120\(\,\text {GeV}\). All pairs of oppositely charged leptons in the four-lepton candidate are required to have \({m_{\ell \ell '} > 4\,\text {GeV}}\) regardless of their flavor. In the rare case of further ambiguity, which occurs in less than 0.5% of events when five or more passing lepton candidates are found, the \({{\mathrm{Z}} _2}\) candidate that maximizes the scalar \(p_{\mathrm {T}} \) sum of the four leptons is chosen.

The \({{\mathrm{p}} {\mathrm{p}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\) cross section is measured using events where both \({m_{{\mathrm{Z}} _1}}\) and \({m_{{\mathrm{Z}} _2}}\) are greater than 60\(\,\text {GeV}\). Decays of the \({{\mathrm{Z}}}\) bosons to \({{\uptau }}\) leptons with subsequent decays to electrons and muons are heavily suppressed by the requirements on lepton \(p_{\mathrm {T}} \), and the contribution of such events is less than 0.5% of the total \({{\mathrm{Z}} {\mathrm{Z}} }\) yield. If these events pass the selection requirements of the analysis, they are considered signal, although they are not considered at generator level in the cross section measurement procedure. Thus, the correction for possible \(\tau \) decays is included in the efficiency calculation.

6 Background estimation

The requirement of four well-reconstructed and isolated lepton candidates strongly suppresses any background; therefore this analysis has very low background contributions, dominated by \({{\mathrm{Z}}}\) boson and \({{\mathrm{W}} {\mathrm{Z}}}\) diboson production in association with jets, and by t t production. In a small fraction of cases, particles from jet fragmentation satisfy both lepton identification and isolation criteria, and thus are misidentified as signal leptons. This background is estimated using control data samples, as decribed below.

The probability for such objects to be selected is measured from a sample of \({{\mathrm{Z}} +\ell _\text {candidate}}\) events, where \({{\mathrm{Z}}}\) denotes a pair of oppositely charged, same-flavor leptons that pass all analysis requirements and satisfy \({| m_{\ell ^+\ell ^-} - m_{{\mathrm{Z}}} | < 10\,\text {GeV}}\), where \({m_{{\mathrm{Z}}}}\) is the nominal \({{\mathrm{Z}}}\) boson mass. Each event in this sample must have exactly one additional object \({\ell _\text {candidate}}\) that passes relaxed identification requirements with no isolation requirements applied. The misidentification probability for each lepton flavor, measured in bins of lepton candidate \(p_{\mathrm {T}} \) and \(\eta \), is defined as the ratio between the number of candidates that pass the final isolation and identification requirements and the total number of candidates in the sample. The number of \({{\mathrm{Z}} +\ell _\text {candidate}}\) events is corrected for the contamination from \({{\mathrm{W}} {\mathrm{Z}}}\) production and for \({{\mathrm{Z}} {\mathrm{Z}} }\) events in which one lepton is not reconstructed. These events have a third genuine, isolated lepton that must be excluded from the misidentification probability calculation. The \({{\mathrm{W}} {\mathrm{Z}}}\) contamination is suppressed by requiring the missing transverse momentum \(p_{\mathrm {T}} ^\text {miss} \) to be below 25\(\,\text {GeV}\). The \(p_{\mathrm {T}} ^\text {miss} \) is defined as the magnitude of the missing transverse momentum vector \({\mathbf {p}}_{\mathrm {T}}^{\text {miss}} \), the projection onto the plane transverse to the beams of the negative vector sum of the momenta of all reconstructed PF candidates in the event, corrected for the jet energy scale. The transverse mass, calculated as \({m_{\mathrm {T}} \equiv \sqrt{\smash [b]{(p_{\mathrm {T}} ^\ell + p_{\mathrm {T}} ^\text {miss})^2 - ({\mathbf {p}}_{\mathrm {T}} ^{\ell } + {\mathbf {p}}_{\mathrm {T}}^{\text {miss}})^2}}}\), is required to be less than 30\(\,\text {GeV}\). The residual contribution of \({{\mathrm{W}} {\mathrm{Z}}}\) and \({{\mathrm{Z}} {\mathrm{Z}} }\) events, which can be up to a few percent of the events with \({\ell _\text {candidate}}\) passing all selection criteria, is estimated from simulation and subtracted.

To account for all sources of background events, two control samples are used to estimate the number of background events in the signal regions. Both are defined as samples that contain events with a dilepton candidate satisfying all requirements (\({{\mathrm{Z}} _1}\)) and two additional lepton candidates \({\ell ^{+}\ell ^{-}}\). In one control sample, enriched in \({{\mathrm{W}} {\mathrm{Z}}}\) events, one \({\ell }\) candidate is required to satisfy the full identification and isolation criteria and the other must fail the full criteria and instead satisfy only the relaxed ones; in the other, enriched in \({{\mathrm{Z}}}\)+jets events, both \({\ell }\) candidates must satisfy the relaxed criteria, but fail the full criteria. The additional leptons must have opposite charges and the same flavor (\({{\mathrm{e}} ^{\pm }{\mathrm{e}} ^{\mp }}\) and \({{\upmu } ^{\pm }{\upmu } ^{\mp }}\)). From this set of events, the expected number of background events in the signal region, denoted “\({{\mathrm{Z}}}\)+X” in the figures, is obtained by scaling the number of observed \({{\mathrm{Z}} _1+\ell ^{+}\ell ^{-}}\) events by the misidentification probability for each lepton failing the selection. The procedure is described in more detail in Ref. [54].

In addition to this reducible background, which contributes to approximately 1–2% of the \({\mathrm{Z}} {\mathrm{Z}} \) events, the \({{{\mathrm{t}} {}{{\mathrm{t}}}} {\mathrm{Z}}}\) and VVV processes with four prompt leptons are estimated from simulated samples to be around 1–1.5% of the expected \({{\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\ell '}\) yield. The total background contributions to the \({{\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\ell '}\) signal regions are summarized in Sect. 8.

7 Systematic uncertainties

The major sources of systematic uncertainty and their effect on the measured cross sections are summarized in Table 1. The lepton identification, isolation, and track reconstruction efficiencies in simulation are corrected with scaling factors derived with a tag-and-probe method and applied as a function of lepton \(p_{\mathrm {T}} \) and \(\eta \). To estimate the uncertainties associated with the tag-and-probe technique, the total yield is recomputed with the scaling factors varied up and down by the tag-and-probe fit uncertainties. The uncertainties associated with the lepton efficiency in the \({{\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\ell '}\) signal regions are 5% in the \({4{\mathrm{e}}}\), 3% in the \({2{\mathrm{e}} 2\mu }\), and 2% in the \(4\mu \) final states.

In both data and simulated event samples, trigger efficiencies are evaluated with a tag-and-probe technique. The ratio of the trigger efficiency estimated using data to the one estimated with simulation is applied to simulated events, and the size of the resulting change in the expected yield is taken as the uncertainty in the determination of the trigger efficiency. This uncertainty is around 1–2% of the final estimated yield.

Table 1 The contributions of each source of systematic uncertainty in the cross section measurements. The integrated luminosity uncertainty, and the PDF and scale uncertainties, are considered separately. All other uncertainties are added in quadrature into a single systematic uncertainty. Uncertainties that vary by decay channel are listed as ranges
Fig. 1
figure 1

Distributions of (left) transverse momentum and (right) pseudorapidity for individual leptons. Points represent the data with error bars showing the statistical uncertainties, histograms the expected SM predictions and reducible background estimated from data. The \(p_{\mathrm {T}}\) distributions includes overflow in the last bin

The largest uncertainty in the estimated background yield arises from differences in sample composition between the \({{\mathrm{Z}} +\ell _\text {candidate}}\) control sample used to calculate the lepton misidentification probability and the \({{\mathrm{Z}} +\ell ^+\ell ^-}\) control sample. An additional uncertainty arises from the limited number of events in the \({{\mathrm{Z}} +\ell _\text {candidate}}\) sample. A systematic uncertainty of 40% is applied to the lepton misidentification probability to cover both effects. Its impact varies by channel, but is of the order of 1% of the total expected yield.

The modeling of pileup relies on the total inelastic \({\mathrm{p}} {\mathrm{p}} \) cross Sect. [55]. The pileup uncertainty is evaluated by varying this cross section up and down by 5%.

Fig. 2
figure 2

Distributions of (upper left) \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\) for \({\mathrm{Z}} {\mathrm{Z}} \) events with \({60< m_{{\mathrm{Z}} _1, {\mathrm{Z}} _2} < 120 \,\text {GeV}}\); (upper right) mass of selected Z boson candidates; (lower left) transverse momentum of the \({{\mathrm{Z}} {\mathrm{Z}} }\) system; (lower right) transverse momentum of individual Z boson candidates. Points represent the data with error bars showing the statistical uncertainties, histograms the expected SM predictions and reducible background estimated from data. All \(p_{\mathrm {T}}\) and \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\) distributions include overflow in the last bin

Uncertainties because of factorization (\(\mu _\mathrm {F}\)) and renormalization (\(\mu _\mathrm {R}\)) scale choices on the \({{\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\ell '}\) acceptance are evaluated with powheg+mcfm by varying \(\mu _\mathrm {F}\) and \(\mu _\mathrm {R}\) up and down by a factor of two with respect to the default values \({\mu _\mathrm {F} = \mu _\mathrm {R} = m_{{\mathrm{Z}} {\mathrm{Z}} }}\), where \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\) is the invariant mass of the \({\mathrm{Z}} {\mathrm{Z}} \) system. All combinations are considered except those in which \(\mu _\mathrm {F}\) and \(\mu _\mathrm {R}\) differ by a factor of four. Parametric uncertainties (PDF+\(\alpha _\mathrm {S} \)) are evaluated according to the PDF4LHC prescription [56] in the acceptance calculation, and with NNPDF3.0 [57] in the cross section calculations. An additional theoretical uncertainty arises from scaling the powheg \({{\mathrm{q}} {{\mathrm{q}}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\) simulated sample from its NLO cross section to the NNLO prediction, and the mcfm \({{\mathrm{g}} {\mathrm{g}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\) samples from their LO cross sections to the NLO predictions. The change in the acceptance corresponding to this scaling procedure is about 1%.

The uncertainty in the integrated luminosity of the data samples is 2.5% (2016) [58], 2.3% (2017) [59], and 2.5% (2018) [60]. Since the luminosity uncertainty contains a significant uncorrelated portion, the relative luminosity uncertainty of the whole sample is smaller than for each individual year.

The same uncertainties are valid for both total and differential cross section measurements, but for the differential one there is also an additional uncertainty related to the unfolding procedure described in Sect. 9. It is estimated using MadGraph 5_amc@nlo instead of powheg+mcfm in unfolding. The unfolding uncertainty is included in the results and plots together with other uncertainites, but its effect is small compared to the statistical uncertainties of the measurement.

Table 2 Observed and expected prefit yields of \({{\mathrm{Z}} {\mathrm{Z}} }\) events, and estimated yields of background events, shown for each final state and combined. The statistical (first) and systematic (second) uncertainties are presented

8 Cross section measurement

The \(p_{\mathrm {T}}\) and \(\eta \) distributions for individual leptons are shown in Fig. 1. Both distributions contain four leptons per event. The invariant mass of the \({\mathrm{Z}} {\mathrm{Z}} \) system, the individual mass of reconstructed Z boson candidates in the \({\mathrm{Z}} {\mathrm{Z}} \) events, and their corresponding \(p_{\mathrm {T}}\) distributions are shown in Fig. 2. The last bins in \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\) and all \(p_{\mathrm {T}}\) distributions contain events from the overflow. The \({m_{\mathrm{Z}}}\) and Z \(p_{\mathrm {T}}\) distributions contain two Z candidates for each event. These distributions are shown for data and simulated events to demonstrate comparisons with SM expectations. The signal expectations include contributions from \({{\mathrm{Z}} {\mathrm{Z}} }\) production shown separately for \({{\mathrm{q}} {{\mathrm{q}}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\), \({{\mathrm{g}} {\mathrm{g}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\), and EW \({\mathrm{Z}} {\mathrm{Z}} \) processes in all figures and combined as “Signal” in Table 2. The EW \({\mathrm{Z}} {\mathrm{Z}} \) production contributes to approximately 1% of the total number of \({\mathrm{Z}} {\mathrm{Z}} \) events.

Table 3 Measured fiducial cross section for each data sample and combined. The first uncertainty is statistical, the second is experimental systematic, and the third is associated with the integrated luminosity
Table 4 Measured total \({\sigma ({\mathrm{p}} {\mathrm{p}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} )}\) cross section for each data sample and combined. The first uncertainty is statistical, the second is experimental systematic, the third is theoretical systematic. The fourth uncertainty is associated with the integrated luminosity

The irreducible background, which amounts to 1–1.5% of the total \({\mathrm{Z}} {\mathrm{Z}} \) yield, and reducible background are combined as “Background” in Table 2. The total background in this analysis is \(\approx 3\%\). The estimated yields agree well with the measured ones. The individual distributions are well described, except the \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\) distribution at high values of invariant masses and the \({p_{\mathrm {T}} ^{\mathrm{Z}}}\) distribution at high values of \(p_{\mathrm {T}} \). These are regions where the EW corrections may become important and will be discussed later in Sect. 10.

The measured yields are used to evaluate the \({{\mathrm{Z}} {\mathrm{Z}} }\) production cross section in the fiducial phase space. The signal acceptance is evaluated from simulation and corrected for each individual lepton flavor in bins of \(p_{\mathrm {T}} \) and \(\eta \) using factors obtained with the tag-and-probe technique. To include all final states in the cross section calculation, a simultaneous fit to the number of observed events in all decay channels is performed. The likelihood is composed as a combination of individual channel likelihoods for the signal and background hypotheses with the statistical and systematic uncertainties treated as scaling nuisance parameters. The combination of various data-taking periods is performed treating the theoretical uncertainties as fully correlated among various periods, whereas the experimental uncertainties are either correlated or uncorrelated, depending on their origin.

Fig. 3
figure 3

The total \({\mathrm{Z}} {\mathrm{Z}} \) cross section as a function of the proton–proton center-of-mass energy. Results from the CMS [4, 5] and ATLAS [9, 10, 14] experiments are compared to predictions from Matrix at NNLO in QCD and NLO in EW, and mcfm at NLO in QCD. The mcfm prediction also includes gluon-gluon initiated production at LO in QCD. The predictions use NNPDF31_nnlo_as_0118_luxqed and NNPDF3.0 PDF sets, respectively, and fixed factorization and renormalization scales \({\mu _\mathrm {F} = \mu _\mathrm {R} = m_{{\mathrm{Z}}}}\). Details of the calculations and uncertainties are given in the text. The ATLAS measurements were performed with a \({{\mathrm{Z}}}\) boson mass window of 66–116\(\,\text {GeV}\), instead of 60–120\(\,\text {GeV}\) used by CMS, and are corrected for the resulting 1.6% difference in acceptance. Measurements at the same center-of-mass energy are shifted slightly along the horizontal axis for clarity

The fiducial phase space for the \({{\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\ell '}\) cross section measurement is defined as: \({p_{\mathrm {T}} ^{\ell _1} > 20\,\text {GeV}}\), \({p_{\mathrm {T}} ^{\ell _2} > 10\,\text {GeV}}\), \({p_{\mathrm {T}} ^{\ell _{3,4}} > 5\,\text {GeV}}\), \({|\eta ^{\ell } | < 2.5}\), \({m_{2\ell } > 4\,\text {GeV}}\) (any opposite-sign same-flavor pair), \({60< m_{{\mathrm{Z}} _1}, m_{{\mathrm{Z}} _2} < 120\,\text {GeV}}\). The generator-level leptons used for the fiducial cross section calculation are “dressed” by adding the momenta of generator-level photons within \({\varDelta R\left( \ell ,\gamma \right) < 0.1}\) from the lepton momenta directions.

The measured \({{\mathrm{Z}} {\mathrm{Z}} }\) fiducial cross section presented in Table 3 can be compared to \(39.3^{+0.8}_{-0.7} \pm 0.6\text { fb} \) calculated with powheg+mcfm using the same settings as the simulated samples with K factors applied. The first uncertainty corresponds to the factorization and renormalization scales and the second to PDF, as described above. The powheg calculations used dynamic factorization and renormalization scales \({\mu _\mathrm {F} = \mu _\mathrm {R} = m_{2\ell 2\ell '}}\), whereas the contribution from mcfm is computed with dynamic scales \({\mu _\mathrm {F} = \mu _\mathrm {R} = 0.5 m_{2\ell 2\ell '}}\). It can also be compared to the prediction from Matrix  v2.0.0_beta1 of \(38.0^{+1.1}_{-1.0}\). The uncertainty in the Matrix prediction includes only the uncertainty due to the variation of \(\mu _\mathrm {F}\) and \(\mu _\mathrm {R}\).

The total \({{\mathrm{Z}} {\mathrm{Z}} }\) production cross section for both dileptons produced in the mass range 60–120\(\,\text {GeV}\) and \({m_{\ell ^+\ell ^{\prime -}} > 4\,\text {GeV}}\) is presented in Table 4. The nominal branching fraction \({\mathcal {B}({\mathrm{Z}} \rightarrow \ell ^+\ell ^-) = 0.03366}\) is used [53]. The measured total cross section can be compared to the theoretical value of \(16.9^{+0.6}_{-0.5} \pm 0.2\text { pb} \), calculated from powheg+mcfm with the same settings that is used for \({\sigma _{\mathrm {fid}} ({\mathrm{p}} {\mathrm{p}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} \rightarrow 2\ell 2\ell ')}\). It can also be compared to \(16.5^{+0.6}_{-0.5}\) \(\text { pb}\), calculated with Matrix  v2.0.0_beta1, or \(15.0^{+0.7}_{-0.6} \pm 0.2\) \(\text { pb}\), calculated with mcfm at NLO in QCD with additional contributions from LO \({{\mathrm{g}} {\mathrm{g}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\) diagrams and with the NLO NNPDF3.0 PDF set and fixed factorization and renormalization scales set to \({\mu _\mathrm {F} = \mu _\mathrm {R} = m_{{\mathrm{Z}}}}\).

Fig. 4
figure 4

Differential cross sections normalized to the fiducial cross section for the combined \({4{\mathrm{e}}}\), \({2{\mathrm{e}} 2{\upmu }}\), and \({4{\upmu }}\) decay channels as a function of \(p_{\mathrm {T}} \) for (left) all leptons, (right) all \({{\mathrm{Z}}}\) bosons in the event. The points represent the unfolded data with error bars showing the statistical uncertainties, the shaded histogram the powheg+mcfm \({\mathrm{Z}} {\mathrm{Z}} \) predictions, and the dashed curves correspond to the results of the Matrix and MadGraph 5_amc@nlo+mcfm calculations. The three lower panels represent the ratio of the measured cross section to the expected distributions from Matrix, powheg+mcfm and MadGraph 5_amc@nlo+mcfm. The shaded areas in all the panels represent the full uncertainties calculated as the quadratic sum of the statistical and systematic uncertainties, whereas the crosses represent only the statistical uncertainties

Fig. 5
figure 5

Differential cross sections normalized to the fiducial cross section for the combined \({4{\mathrm{e}}}\), \({2{\mathrm{e}} 2{\upmu }}\), and \({4{\upmu }}\) decay channels as a function of (left) \(p_{\mathrm {T}} \) of the \({{\mathrm{Z}} {\mathrm{Z}} }\) system, (right) the invariant mass of the \({\mathrm{Z}} {\mathrm{Z}} \) system. The points represent the unfolded data with error bars showing the statistical uncertainties, shaded histogram the powheg+mcfm \({\mathrm{Z}} {\mathrm{Z}} \) predictions, and the dashed curves correspond to the results of the Matrix and MadGraph 5_amc@nlo+mcfm calculations. The three lower panels represent the ratio of the measured cross section to the expected distributions from Matrix, powheg+mcfm and MadGraph 5_amc@nlo+mcfm. The shaded areas in all the panels represent the full uncertainties calculated as the quadratic sum of the statistical and systematic uncertainties, whereas the crosses represent only the statistical uncertainties

Fig. 6
figure 6

Differential cross sections normalized to the fiducial cross section for the combined \({4{\mathrm{e}}}\), \({2{\mathrm{e}} 2{\upmu }}\), and \({4{\upmu }}\) decay channels as a function of the azimuthal (left) and \(\varDelta \)R (right) separation of the two \({{\mathrm{Z}}}\) bosons. The points represent the unfolded data with error bars showing the statistical uncertainties, the shaded histogram the powheg+mcfm \({\mathrm{Z}} {\mathrm{Z}} \) predictions, and the dashed curves correspond to the results of the Matrix and MadGraph 5_amc@nlo+mcfm calculations. The three lower panels represent the ratio of the measured cross section to the expected distributions from Matrix, powheg+mcfm and MadGraph 5_amc@nlo+mcfm. The shaded areas in all the panels represent the full uncertainties calculated as the quadratic sum of the statistical and systematic uncertainties, whereas the crosses represent only the statistical uncertainties

Fig. 7
figure 7

Distribution of the reconstructed \({\mathrm{Z}} {\mathrm{Z}} \) mass for the combined \({4{\mathrm{e}}}\), \({2{\mathrm{e}} 2{\upmu }}\), and \({4{\upmu }}\) channels. Points represent the data with error bars showing the statistical uncertainties, the shaded histograms represent the SM prediction including signal and irreducible background from simulation, and the reducible background estimate from data. Dashed histogram represents an example of the aTGC signal. The last bin includes contribution from all events with mass above 1300\(\,\text {GeV}\)

The total \({{\mathrm{Z}} {\mathrm{Z}} }\) cross section is shown in Fig. 3 as a function of the \({\mathrm{p}} {\mathrm{p}} \) center-of-mass energy. Results from CMS [4, 5] and ATLAS [9, 10, 14] are compared to predictions from Matrix  v2.0.0_beta1 and mcfm. The uncertainties are statistical (inner bars) and statistical and systematic combined, as obtained from the fit (outer bars). The band around the Matrix predictions reflects scale uncertainties, while the band around the mcfm predictions reflects both scale and PDF uncertainties.

9 Differential cross sections

The differential distributions normalized to the fiducial cross sections are presented in Figs. 4, 5 and 6 for the combination of the \({4{\mathrm{e}}}\), \({2{\mathrm{e}} 2{\upmu }}\), and \({4{\upmu }}\) decay channels using the whole data sample. The fiducial cross section definition includes \({p_{\mathrm {T}} ^{\ell }}\) and \({|\eta ^{\ell } |}\) selections on each lepton, and the 60–120\(\,\text {GeV}\) mass requirement, as described in Sect. 4. Figure 4 shows the differential cross sections in bins of \(p_{\mathrm {T}} \) for: (left) all leptons in the event, (right) both \({{\mathrm{Z}}}\) bosons in the event, and in Fig. 5 (left) for the \(p_{\mathrm {T}} \) of the \({{\mathrm{Z}} {\mathrm{Z}} }\) system. Figure 5 (right) shows the normalized \({\mathrm{d}\sigma /\mathrm{d}m_{{{\mathrm{Z}} {\mathrm{Z}} }}}\) distribution. All \(p_{\mathrm {T}}\) and \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\) distributions include overflow in the last bin. Figure 6 shows the angular correlations between \({{\mathrm{Z}}}\) bosons. The data are corrected for background contributions and unfolded for detector effects using a matrix inversion method without regularization as described in Ref. [61], and compared with the theoretical predictions from powheg+mcfm, MadGraph 5_amc@nlo+mcfm and Matrix. The distributions include both Z boson candidates or all four leptons, where applicable, and are normalized to the numbers of objects in the event and to the fiducial cross section. The bottom part of each plot shows the ratio of the measured to the predicted values. The bin sizes are chosen according to the resolution of the relevant variables, trying also to keep the statistical uncertainties at a similar level for all the bins.

The distributions predicted by powheg+mcfm and MadGraph 5_amc@nlo+mcfm agree well with data, except for \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\). This distribution shows a small overestimate in the cross section at high invariant masses. The Matrix predictions describe this region better, which can be explained by the presence of the EW corrections in the Matrix calculations. The effect of EW corrections is in detail discussed in Ref. [44] and can reach 20–30% for \({m_{{\mathrm{Z}} {\mathrm{Z}} } = 1 \,\text {TeV}}\). On the other hand, the Matrix predictions show some deviation from the measurements as a function of \({p_{\mathrm {T}} ^{{\mathrm{Z}} {\mathrm{Z}} }}\) and for the azimuthal separation between the two Z bosons, which is not observed for powheg+mcfm and MadGraph 5_amc@nlo+mcfm predictions.

10 Limits on anomalous triple gauge couplings

The presence of aTGCs is expected to increase the event yield at high four-lepton masses. Figure 7 presents the distribution of the four-lepton reconstructed mass for the combined \({4{\mathrm{e}}}\), \({2{\mathrm{e}} 2{\upmu }}\), and \({4{\upmu }}\) channels, for the SM and an example of nonzero aTGC value with \(f_4^\gamma =0\), and \({f_4^{\mathrm{Z}} =0.0015}\). Limits on aTGCs are derived from fits to this distribution. The shaded histograms represent the SM predictions as described in the previous sections and the dashed curve shows the sherpa prediction. The sherpa SM predictions are normalized to the powheg+mcfm predictions including K factors and agree well with them in shape, as shown in Fig. 7. As a cross-check of the procedure, the sherpa SM distribution was also corrected bin-by-bin to the powheg+mcfm distribution, no difference was observed in the extracted limits. The presence of aTGC contribution increases the expected event yields at masses above 1300\(\,\text {GeV}\). In the fit, described below, this region is subdivided into two bins: 1300–2000\(\,\text {GeV}\) and above 2000\(\,\text {GeV}\). Typically 60–70% of the aTGC events have masses above 2000\(\,\text {GeV}\), whereas the expected SM contribution is approximately 2.4 and 0.2 events in the 1300–2000\(\,\text {GeV}\) and above 2000\(\,\text {GeV}\) bins, respectively.

The invariant mass distributions are interpolated from those obtained from the sherpa simulation for different values of the anomalous couplings in the range between 0 and 0.03. For each distribution, only one or two couplings are varied while all others are set to zero, thus creating a grid of points in the \({(f_{4}^{\mathrm{Z}}, f_{4}^\gamma )}\) and \({(f_{5}^{\mathrm{Z}}, f_{5}^\gamma )}\) parameter planes and the corresponding invariant mass distributions. In each \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\) bin, expected signal values are interpolated between the two-dimensional grid points using a second-order polynomial, since the cross section for the signal depends quadratically on the coupling parameters. A simultaneous fit to the values of aTGCs is performed for all lepton channels, see Ref. [62] for details. A profile likelihood method [53], Wald Gaussian approximation, and Wilks theorem [63] are used to derive one-(1D) and two-dimensional limits at 68 and 95% confidence levels (\(\text {CL}\)) on each of the aTGC parameters and combination of two of them, while all other parameters are set to their SM values. All systematic uncertainties are included by varying the number of expected signal and background events within their uncertainties. An additional 10% uncertainty is applied on the predictions of the SM and aTGC models to account for possible differences between model predictions and the interpolation used in the fit. No form factor [64] is used when deriving the limits; the results assume that the energy scale of new physics is very high. The constraints on anomalous couplings are displayed in Fig. 8. The curves indicate 68 and 95% \(\text {CL}\) contours; the dots indicate where the likelihoods reach their maximum. Coupling values outside the contours are excluded at the corresponding \(\text {CL}\). The crosses in the middle represent the observed 1D limits that are summarized in Table 5. The sensitivity is dominated by the statistical uncertainties.

Complete one-loop EW corrections to massive vector boson pair production [66, 67] were applied as a cross-check. The EW corrections to the \({\mathrm{Z}} {\mathrm{Z}} \) production cause the \({\mathrm{Z}} {\mathrm{Z}} \) mass spectrum to fall more rapidly at large masses. In addition, the overall cross section decreases by about 4%. The effect of NLO EW corrections is estimated by reweighting the SM sherpa sample as a function of \({m_{{\mathrm{Z}} {\mathrm{Z}} }}\) using weights derivedfrom the calculations described in Ref. [66]. This reweighting improves the expected limits by about 4–6%, whereas there is no effect on the observed limits. This is expected, since only the SM contribution is subject to the EW corrections; they are not applied on aTGCs. The limits are driven by the high mass tail above 1300 GeV. In this region the aTGC signal is much larger than the SM, and therefore the EW correction on the SM part has a very small effect on the predictions of the SM+aTGC model. This correction is much smaller than the uncertainty we apply in the fit procedure.

These results can be also expressed in terms of EFT parameters. The numerical relations between aTGCs and EFT parameters are given in Ref. [65]. The expected and measured limits in terms of EFT are presented in Table 5.

11 Summary

Fig. 8
figure 8

Two-dimensional observed (solid) and expected (dashed) contours exclusion limits at 95% \(\text {CL}\), and at 68 and 95% \(\text {CL}\), respectively, on the \({{\mathrm{Z}} {\mathrm{Z}} {\mathrm{Z}} }\) and \({{\mathrm{Z}} {\mathrm{Z}} \gamma }\) aTGCs. The plots show the exclusion contours in the \({(f_{4(5)}^{\mathrm{Z}}, f_{4(5)}^\gamma )}\) parameter planes. Dots show where the likelihoods reach their maximum. The coupling values outside the contours are excluded at the corresponding confidence level. The crosses in the middle represent the observed 1D limits. No form factor is used

Four-lepton final states have been studied in proton–proton collisions at \(\sqrt{s} = 13\,\text {TeV} \) with the CMS detector at the CERN LHC. The data sample corresponds to an integrated luminosity of 137\(\,\text {fb}^{-1}\), collected during 2016–2018. Themeasured \({{\mathrm{p}} {\mathrm{p}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} }\) total cross section is \(\sigma _{\text {tot}} ( {\mathrm{p}} {\mathrm{p}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} ) = 17.4 \pm 0.3 \,\text {(stat)} \pm 0.5 \,\text {(syst)} \pm 0.4 \,\text {(theo)} \pm 0.3 \,\text {(lumi)} \text { pb} \), where the \({{\mathrm{Z}}}\) boson masses are in the range \(60< m_{{\mathrm{Z}}} < 120\,\text {GeV} \). The results agree with the SM predictions, discussed in Sect. 8. The differential cross sections also agree well with the SM predictions. Improved limits on anomalous \({{\mathrm{Z}} {\mathrm{Z}} {\mathrm{Z}}}\) and \({{\mathrm{Z}} {\mathrm{Z}} \gamma }\) triple gauge couplings are established. These are the most stringent limits to date on anomalous \({{\mathrm{Z}} {\mathrm{Z}} {\mathrm{Z}}}\) and \({{\mathrm{Z}} {\mathrm{Z}} \gamma }\) triple gauge couplings and they improve the previous strictest results from CMS by \(\approx 30\)–40%.

Table 5 Expected and observed one-dimensional 95% \(\text {CL}\) limits on aTGC parameters. The corresponding constrains on EFT parameters are estimated using the transformation from Ref. [65]