$ZZ \to \ell^{+}\ell^{-}\ell^{\prime +}\ell^{\prime -}$ cross-section measurements and search for anomalous triple gauge couplings in 13 TeV $pp$ collisions with the ATLAS detector

Measurements of $ZZ$ production in the $\ell^{+}\ell^{-}\ell^{\prime +}\ell^{\prime -}$ channel in proton-proton collisions at 13 TeV center-of-mass energy at the Large Hadron Collider are presented. The data correspond to 36.1 $\mathrm{fb}^{-1}$ of collisions collected by the ATLAS experiment in 2015 and 2016. Here $\ell$ and $\ell'$ stand for electrons or muons. Integrated and differential $ZZ \to \ell^{+}\ell^{-}\ell^{\prime +}\ell^{\prime -}$ cross sections with $Z \to \ell^+\ell^-$ candidate masses in the range of 66 GeV to 116 GeV are measured in a fiducial phase space corresponding to the detector acceptance and corrected for detector effects. The differential cross sections are presented in bins of twenty observables, including several that describe the jet activity. The integrated cross section is also extrapolated to a total phase space and to all Standard-Model decays of $Z$ bosons with mass between 66 GeV and 116 GeV, resulting in a value of $17.3 \pm 0.9$ [$\pm 0.6$ (stat.) $\pm 0.5$ (syst.) $\pm 0.6$ (lumi.)] pb. The measurements are found to be in good agreement with the Standard-Model predictions. A search for neutral triple gauge couplings is performed using the transverse momentum distribution of the leading $Z$-boson candidate. No evidence for such couplings is found and exclusion limits are set on their parameters.


Introduction
The study of the production of Z boson pairs in proton-proton (pp) interactions at the Large Hadron Collider (LHC) [1] tests the electroweak sector of the Standard Model (SM) at the highest available energies. Example Feynman diagrams of ZZ production at the LHC are shown in Figure 1. In pp collisions at a center-of-mass energy of √ s = 13 TeV, ZZ production is dominated by quark-antiquark (qq) interactions, with an O(10%) contribution from loop-induced gluon-gluon (gg) interactions [2,3]. The production of ZZ in association with two electroweakly produced jets, denoted EW-ZZ j j, includes the rare ZZ weak-boson scattering process. Study of ZZ production in association with jets is an important step in searching for ZZ weak-boson scattering, which has so far not been experimentally observed by itself.  Figure 1: Examples of leading-order SM Feynman diagrams for ZZ production in proton-proton collisions: (a) qq-initiated, (b) gg-initiated, (c) electroweak ZZ j j production, (d) electroweak ZZ j j production via weak-boson scattering.
The SM ZZ production can also proceed via a Higgs boson propagator, although this contribution is expected to be suppressed in the region where both Z bosons are produced nearly on-shell, as is the case in this analysis. Non-Higgs-mediated ZZ production is an important background in studies of the Higgs boson properties [4][5][6][7]. It is also a major background in searches for new physics processes producing pairs of Z bosons at high invariant mass [8][9][10][11] and it is sensitive to anomalous triple gauge couplings (aTGCs) of neutral gauge bosons, which are not allowed in the SM [12]. The SM does not have tree-level vertices coupling three neutral gauge bosons (ZZZ, ZZγ), because these would violate the underlying SU(2) L × U(1) Y symmetry. However, these couplings exist in some extensions of the SM, enhancing the ZZ production cross section in regions where the energy scale of the interaction is high.
An example Feynman diagram of ZZ production via aTGC is shown in Figure 2. Integrated and differential ZZ production cross sections were previously measured at √ s = 7 and 8 TeV by the ATLAS and CMS collaborations [13][14][15][16] and found to be consistent with SM predictions. The integrated pp → ZZ → + − + − cross section at √ s = 13 TeV was recently measured by the ATLAS Z/γ * q q Z Z Figure 2: Example Feynman diagram of ZZ production containing an aTGC vertex, here indicated by a red dot, which is forbidden in the SM.
[17] and CMS [18] collaborations, each analyzing data corresponding to an integrated luminosity of about 3 fb −1 . Searches for aTGCs were previously performed at lower center-of-mass energies by ATLAS [15], CMS [14,19], D0 [20], and by the LEP experiments [21]. This paper represents an extension of the ATLAS measurement, using a total of 36.1 ± 1.1 fb −1 of data collected with the ATLAS detector in the years 2015 and 2016.
In this analysis, candidate events are reconstructed in the fully leptonic ZZ → + − + − decay channel where and can be an electron or a muon. Throughout this analysis, "Z boson" refers to the superposition of a Z boson and a virtual photon in the mass range from 66 GeV to 116 GeV, as these are not strictly distinguishable when decaying to charged leptons. A fiducial phase space is defined, reflecting both the acceptance of the ATLAS detector [22,23] and the selections imposed on the reconstructed leptons in this analysis. Both the integrated and differential cross sections are measured, the latter with respect to 20 different observables. Ten of these directly measure the jet activity in the events. The observed event yields are unfolded to the fiducial phase space using simulated samples to model the detector effects. The integrated cross sections are inclusive with respect to jet production. For easier comparison to other measurements, the integrated fiducial cross sections determined in different leptonic channels are combined and extrapolated to a total phase space and to all Z boson decay modes. A search for aTGCs is performed by looking for deviations of the data from the SM predictions at high values of the transverse momentum of the leading-p T Z boson, which is one of the observables most sensitive to the energy scale of the interaction. 1 Differential fiducial cross sections are measured with respect to the following observables: • Transverse momentum of the four-lepton system, p T, 4 ; • Absolute rapidity of the four-lepton system, |y 4 |; • Separation in azimuthal angle between the two Z boson candidates, δφ Z 1 , Z 2 , defined such that it lies in the interval [0, π]; • Absolute difference in rapidity between the two Z boson candidates, |δy Z 1 , Z 2 |; • Transverse momentum of the leading-p T and the subleading-p T Z boson candidates, p T, Z 1 and p T, Z 2 ; • Transverse momentum of each of the four leptons; 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point in the center of the detector and the z-axis along the beam pipe. The x-axis points to the center of the LHC ring, and the y-axis points upward. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln[tan(θ/2)]. Transverse momentum p T is the projection of momentum onto the transverse plane.
• Number of jets with p T > 30 GeV and |η| < 4.5; • Number of jets with p T > 30 GeV and |η| < 2.4; • Number of jets with p T > 60 GeV and |η| < 4.5; • Scalar sum of the transverse momenta of all jets in the event with p T > 30 GeV and |η| < 4.5; • Absolute pseudorapidity of the leading-p T and the subleading-p T jets; • Transverse momentum of the leading-p T and the subleading-p T jets; • Absolute difference in rapidity between the two leading-p T jets, |δy jet 1 , jet 2 |; • Invariant mass of the two leading-p T jets, m jet 1 , jet 2 .
These measurements provide a detailed description of the kinematics in ZZ events and allow comparisons and validations of current and future predictions. Some of the differential measurements are particularly motivated: the transverse momentum of the four-lepton system directly measures the recoil against all other particles produced in the collision and therefore provides information about quantum chromodynamics (QCD) and electroweak radiation across the entire range of scales. The rapidity of the four-lepton system is sensitive to the z-component of the total momentum of the initial-state partons involved in the ZZ production. It may therefore be sensitive to the parton distribution functions (PDFs). The azimuthalangle separation and rapidity difference between the Z boson candidates probe their angular correlations and may help extract the contribution of double-parton-scattering ZZ production. The azimuthal-angle separation is also sensitive to radiation of partons and photons produced in association with the ZZ pair. The scalar sum of the transverse momenta of all jets provides a measure of the overall jet activity that is independent of their azimuthal configuration. The measurements of |δy jet 1 , jet 2 | and m jet 1 , jet 2 are particularly sensitive to the EW-ZZ j j process. They both tend to have larger values in weak-boson scattering than in other ZZ production channels, providing an important step towards the study of ZZ production via weak-boson scattering.

ATLAS detector
The ATLAS detector [22, 23] is a multipurpose particle detector with a cylindrical geometry. It consists of layers of inner tracking detectors, calorimeters, and muon chambers. The inner detector (ID) is immersed in a 2 T axial magnetic field generated by a thin superconducting solenoid and provides charged-particle tracking and momentum measurement in the pseudorapidity range |η| < 2.5. The calorimeter system covers the pseudorapidity range |η| < 4.9. Electromagnetic calorimetry is provided by high-granularity lead/liquid-argon calorimeters in the region |η| < 3.2 and by copper/liquid-argon calorimeters in the region 3.2 < |η| < 4.9. Within |η| < 2.47 the finely segmented electromagnetic calorimeter, together with the ID information, allows electron identification. Hadronic calorimetry is provided by the steel/scintillator-tile calorimeter within |η| < 1.7 and two copper (or tungsten)/liquid-argon calorimeters within 1.7 < |η| < 4.9. The muon spectrometer (MS) comprises separate trigger and high-precision tracking chambers. The precision chamber system covers the region |η| < 2.7 with three layers of monitored drift tubes, complemented by cathode strip chambers in the forward region, where the hit rate is highest. The muon trigger system covers the range |η| < 2.4 with resistive plate chambers in the central, and thin gap chambers in the forward regions. A two-level trigger system is used to select events of interest in real time [24]. The Level-1 trigger is implemented in hardware and uses a subset of detector information to reduce the event rate to a value of around 100 kHz. This is followed by a software-based high-level trigger system that reduces the event rate to about 1 kHz.  [33] (with α S = 0.118 at the Z pole mass), with the qq-initiated process simulated at next-toleading order (NLO) for ZZ plus zero or one additional jet and at leading order (LO) for two or three additional jets generated at the matrix-element level. A Sherpa 2.1.1 ZZ sample is generated with the loop-induced gg-initiated process simulated at LO using NLO PDFs, including subprocesses involving a Higgs boson propagator, with zero or one additional jet. The gg-initiated process first enters at nextto-next-to-leading order (NNLO) and is therefore not included in the NLO sample for the qq-initiated process. Due to different initial states, the gg-initiated process does not interfere with the qq-initiated process at NLO. The loop-induced gg-initiated process calculated at its LO (α 2 S ) receives large corrections at NLO (α 3 S ) [3]. The cross section of the sample is therefore multiplied by an NLO/LO K-factor of 1.67 ± 0.25 [3]. The EW-ZZ j j process is simulated using Sherpa 2.1.1 at its lowest contributing order in the electroweak coupling, α 6 (including the decays of the Z bosons). It includes the triboson subprocess ZZV → + − + − j j, where the third boson V decays hadronically. Sherpa also simulates parton showering, electromagnetic radiation, the underlying event, and hadronization in the above samples, using its default set of tuned parameters (tune). Throughout this paper, the prediction obtained by summing the above samples is referred to as the nominal Sherpa setup.

Simulated samples and theoretical predictions
An alternative prediction for the qq-initiated process is obtained using the Powheg method and framework [34,35] as implemented in Powheg-Box 2 [36], with a diboson event generator [37,38] used to simulate the ZZ production process at NLO. The simulation of parton showering, electromagnetic radiation, the underlying event, and hadronization is performed with Pythia 8.186 [39,40] using the AZNLO parameter tune [41]. This sample is used to estimate the systematic uncertainty due to modeling differences between the event generators.
Additional samples are generated to estimate the contribution from background events.  [43] and the A14 tune [44].
In all MC samples, additional pp interactions occurring in the same bunch crossing as the process of interest or in nearby ones (pileup) are simulated with Pythia using MSTW 2008 PDFs [45] and the A2 tune [46]. The samples are then passed through a simulation of the ATLAS detector [47] based on Geant 4 [48]. Weights are applied to the simulated events to correct for the small differences from data in the reconstruction, identification, isolation, and impact parameter efficiencies for electrons and muons [49,50]. Furthermore, the lepton momentum or energy scales and resolutions are adjusted such that data and simulation match [50,51]. NNLO cross sections for pp → ZZ → + − + − in the fiducial and total phase space are provided by Matrix [2], also in bins of the jet-inclusive unfolded distributions. They include the gg-initiated process at its lowest contributing order, which accounts for about 60% of the cross-section increase with respect to NLO [52]. The calculation uses a dynamic QCD scale of m 4 /2 and the NNPDF 3.0 PDFs, with NNLO PDFs being used also for the gg-initiated process. It uses the G µ electroweak scheme, in which the Fermi constant G µ as well as the pole masses of the weak bosons are taken as independent input parameters [53].
The NNLO calculation is also used for extrapolation of the integrated cross section from the fiducial to a total phase space. The estimation of PDF uncertainties with Matrix is currently unfeasible, because it would require repeating the entire calculation for each PDF variation, which is too computationally expensive. Therefore, these are estimated using an NLO (LO) calculation for the qq-initiated (gg-initiated) process from MCFM [54], taking the mass of the four-lepton system, m 4 , as the dynamic QCD scale. NLO PDFs are used for the gg-initiated process and its contribution is multiplied by the NLO/LO K-factor of 1.67 ± 0.25.
Electroweak corrections at next-to-leading order (NLO EW) [55,56] are calculated in the fiducial phase space, also in bins of the jet-inclusive unfolded distributions. The G µ scheme is used. The NLO/LO EW K-factor integrated across the entire fiducial phase space is about 0.95. The NLO EW corrections are calculated with respect to the qq-initiated process at LO in α S , meaning that they cannot be obtained differentially in observables that are trivial at LO in α S , e.g. the transverse momentum of the four-lepton system. Where a differential calculation is not possible, the integrated value in the fiducial phase space is used. The higher-order NNLO QCD and NLO EW corrections are applied to the predictions only where explicitly stated.
The NNLO calculations serve as the basis of a SM prediction incorporating the formally most accurate available predictions. The contribution of the gg-initiated process is multiplied by the NLO/LO K-factor of 1.67 ± 0.25. The NLO EW corrections are applied as multiplicative K-factors, differentially in the observable of interest if available, otherwise integrated over the fiducial phase space. They are never applied to the gg-initiated loop-induced process, as its topology is considered too different from the LO QCD predictions of the qq-initiated process for which the NLO EW corrections are calculated. The cross section of the EW-ZZ j j process calculated with Sherpa is added to the signal prediction.

Fiducial phase space
The fiducial phase space is defined using final-state particles, meaning particles whose average lifetime τ 0 satisfies cτ 0 > 10 mm [57]. A prompt lepton, photon, or neutrino refers to a final-state particle that does not originate from the decay of a hadron or τ lepton, or any material interaction (such as Bremsstrahlung or pair production) [57]. Particles other than leptons, photons, and neutrinos are never considered prompt in this analysis.
The requirements used to define the fiducial phase space mirror the selections applied to the reconstructed leptons. This is done to ensure that the extrapolation from the observed data to the fiducial phase space is as model-independent as possible, ideally depending only on detector effects. Events in the fiducial phase space contain at least four prompt electrons and/or prompt muons. The four-momenta of all prompt photons within ∆R = (∆η) 2 + (∆φ) 2 = 0.1 of a lepton are added to the four-momentum of the closest lepton. This dressing is done to emulate the effects of quasi-collinear electromagnetic radiation from the charged leptons on their experimental reconstruction in the detector [57]. Each dressed lepton is required to have transverse momentum p T > 5 GeV and absolute pseudorapidity |η| < 2.7.
All possible pairs of same-flavor opposite-charge dileptons are formed, referred to as quadruplets. In each quadruplet, the three highest-p T leptons must satisfy p T > 20 GeV, 15 GeV, and 10 GeV, respectively. If multiple selected quadruplets are present, the quadruplet minimizing |m a − m Z | + |m b − m Z | is selected, where m a, b is the mass of a given same-flavor opposite-charge dilepton and m Z = 91.1876 GeV is the Z boson pole mass [58]. All remaining requirements are applied to the leptons in the final selected quadruplet. Any two same-flavor (different-flavor) leptons i , ( ) j must be separated by ∆R( i , ( ) j ) > 0.1 (0.2). All possible same-flavor opposite-charge dileptons must have an invariant mass greater than 5 GeV, to match the same requirement in the selection of reconstructed events, which is introduced to reduce the background from leptonically decaying hadrons. If all leptons are of the same flavor, the dilepton pairing that minimizes |m a − m Z | + |m b − m Z | is chosen. The selected dileptons are defined as the Z boson candidates. Each is required to have an invariant mass between 66 GeV and 116 GeV. Based on the leptons in the chosen quadruplet, events are classified into three signal channels: 4e, 4µ, and 2e2µ.
Jets are used for several differential cross sections. They are clustered from all final-state particles except prompt leptons, prompt neutrinos, and prompt photons using the anti-k t algorithm [59] with radius parameter 0.4, implemented in FastJet [60]. Jets are required to have p T > 30 GeV and |η| < 4.5. Jets within ∆R = 0.4 of any selected fiducial lepton (as defined above) are rejected.
The fiducial selection is summarized in Table 1.

Signal-process definition
Some SM processes can pass the fiducial selection but are still excluded from the signal. They are considered irreducible backgrounds and are subtracted from the sample of selected candidate events. Any events containing four prompt leptons plus any additional leptons, neutrinos, or photons are considered irreducible backgrounds. An example is the triboson process ZZW + → + − + − + ν . In practice, predictions only exist for a subset of such processes. The irreducible backgrounds that are subtracted are discussed in Section 6. They are very small, approximately 1% of the predicted signal.
The fiducial phase space is inclusive with respect to jets, independently of their origin. Triboson (and higher boson-multiplicity) processes producing a ZZ pair decaying leptonically with any additional electroweak bosons decaying hadronically are included in the signal, as are any other SM processes of the pattern (ZZ → + − + − ) + (X → jets). In practice, only the process ZZV → + − + − j j is included in the theoretical predictions, in the EW-ZZ j j sample generated with Sherpa.
Production via double parton scattering in the same pp collision is included in the signal. Its contribution is not included in the theoretical predictions, but is expected to be smaller than 1% of the total signal yield. This estimate assumes incoherent double parton scattering and is based on a measurement of the effective area parameter at √ s = 7 TeV [61]. Various other measurements of the effective area parameter were made [62-69] and suggest no significant dependence on the center-of-mass energy nor the final state from which the area parameter was measured.  Pythia is used to calculate the fraction of produced events that fall in the fiducial phase space.

Event selection
The event selection begins with trigger and data-quality requirements. Candidate events are preselected by single-, di-, or trilepton triggers [24], with a combined efficiency very close to 100%. They must have at least one primary vertex [70] with two or more associated tracks with p T > 400 MeV. Events must satisfy cleaning criteria [71] designed to reject events with excessive noise in the calorimeters. The data are subjected to quality requirements to reject events in which detector components were not operating correctly.
Following this preselection, muons, electrons, and jets are selected in each event as described below. Based on these, the best lepton quadruplet is selected and required to satisfy further selection criteria.

Selection of muons, electrons, and jets
A muon is reconstructed by matching a track (or track segment) reconstructed in the MS to a track reconstructed in the ID [50]. Its four-momentum is calculated from the curvature of the track fitted to the combined detector hits in the two systems, correcting for energy deposited in the calorimeters. In regions with limited coverage from the MS (|η| < 0.1) or outside the ID acceptance (2.5 < |η| < 2.7), muons can also be reconstructed by matching calorimeter signals consistent with muons to ID tracks (calorimetertagged muons) or standalone in the MS, respectively. Quality requirements and the loose identification criteria are applied as described in Ref. [50]. Muons are required to have |η| < 2.7 and p T > 5 GeV. Calorimeter-tagged muons must have p T > 15 GeV.
An electron is reconstructed from an energy deposit (cluster) in the electromagnetic calorimeter matched to a high-quality track in the ID. Its momentum is computed from the cluster energy and the direction of the track and calibrated [51]. Electrons are required to have |η| < 2.47 and p T > 7 GeV.
Electrons can be distinguished from other particles using several identification criteria that rely on the shapes of electromagnetic showers as well as tracking and track-to-cluster matching quantities. Following the description in Ref.
[49], the output of a likelihood function taking these quantities as input is used to identify electrons, choosing the loose working point.
Leptons are required to originate from the hard-scattering vertex, defined as the reconstructed vertex [70] with the largest sum of the p 2 T of the associated tracks. The longitudinal impact parameter of each lepton track, calculated with respect to the hard-scattering vertex and multiplied by sin θ of the track, is required to be less than 0.5 mm. Furthermore, muons must have a transverse impact parameter calculated with respect to the beam line less than 1 mm in order to reject muons originating from cosmic rays. The significance of the transverse impact parameter 2 calculated with respect to the beam line is required to be less than three (five) for muons (electrons). Standalone muons are exempt from all three requirements, as they do not have an ID track.
Leptons are required to be isolated from other particles using both ID-track and calorimeter-cluster information. Muons Jets [73] are clustered from topological clusters in the calorimeters using the anti-k t algorithm [59] with radius parameter 0.4. Their energy is calibrated as described in Ref. [74]. They are required to have |η| < 4.5 and p T > 30 GeV, as in the fiducial definition. In order to reject jets originating from pileup interactions, they must either pass a jet vertex tagging selection [75,76] or have p T > 60 GeV.
In order to avoid the reconstruction of multiple electrons, muons, and/or jets from the same detector signature, all but one such overlapping objects are removed. Electron candidates sharing an ID track with a selected muon are rejected, except if the muon is only calorimeter-tagged, in which case the muon is rejected instead. Electron candidates sharing their track or calorimeter cluster with a selected higher-p T electron are rejected. Jets within ∆R = 0.4 of a selected lepton are rejected.

Quadruplet selection
As in the fiducial definition (Section 4.1), events must contain at least one quadruplet. All possible quadruplets in a given event are considered for further selection. At most one muon in each quadruplet may be a calorimeter-tagged or standalone muon. The three highest-p T leptons in each quadruplet must satisfy p T > 20 GeV, 15 GeV, 10 GeV, respectively. If multiple selected quadruplets are present, the best quadruplet is chosen as in the fiducial phase-space selection. Only the best quadruplet is considered further and the following requirements are imposed on the leptons in that quadruplet. Any two same-flavor (different-flavor) leptons i , ( ) j must be separated by ∆R( i , ( ) j ) > 0.1 (0.2). All possible same-flavor opposite-charge dileptons must have an invariant mass greater than 5 GeV, to reduce background from leptonic hadron decays. The two Z boson candidates, formed as in the fiducial definition, are required to have an invariant mass between 66 GeV and 116 GeV. Figure 3 shows the distribution of invariant masses of the Z boson candidates in the selected data events. Based on the leptons in the chosen quadruplet, events are classified into the 4e, 4µ, and 2e2µ signal channels.

Background estimation
The expected total background is very small, approximately 2% of the total predicted yield in each decay channel.
Irreducible backgrounds from processes with at least four prompt leptons in the final state are estimated with the simulated samples described in Section 3, including uncertainties from the cross-section predictions, luminosity measurement, and experimental effects, described in Section 7. Non-hadronic triboson processes (15% of the total background estimate) and tt Z processes with leptonic W/Z boson decays (19%) are considered. Simulated samples are also used to estimate the background from ZZ processes where at least one Z boson decays to τ leptons (8%), which is not an irreducible background as defined in Section 4.2.
Events from processes with two or three prompt leptons, e.g. Z, WW, WZ, tt, and ZZ events where one Z boson decays hadronically, can pass the event selection if associated jets, non-prompt leptons, or photons are misidentified as prompt leptons. This background is estimated using a data-driven technique as follows. A lepton selection that is orthogonal to the nominal selection in Section 5.1 is defined by reversing some of its requirements. Muons must fail the transverse impact parameter requirement or the isolation requirement, or both. Electrons must fail either the isolation requirement or the likelihood-based identification, but not both. A high-purity data sample of events containing a Z boson candidate decaying to a pair of electrons or muons is selected. The leptons forming the Z candidate must pass tight selection criteria, different from those used anywhere else in this analysis. Any additional reconstructed leptons in this sample are assumed to be misidentified, after the approximately 4% contamination from genuine third leptons from WZ and ZZ production has been subtracted using MC simulation. Using the observed rates of third leptons passing the nominal or the reversed selection, n l and n r , transfer factors f are defined as f = n l n r and measured in bins of p T and η of the third leptons. A background control sample of data events is then selected, satisfying all the ZZ selection criteria described in Section 5, except that one or two leptons in the final selected quadruplet are required to only satisfy the reversed criteria and not the nominal criteria. The number of observed events with one lepton (two leptons) satisfying only the reversed criteria is denoted N lllr (N llrr ). The events originate predominantly from processes with two or three prompt leptons. Using MC simulation, the contamination from genuine ZZ events is estimated to be approximately 36% of N lllr and approximately 1% of N llrr . The number of background events with one or two misidentified leptons can be calculated as where the superscript ZZ indicates the MC-simulated contributing events from ZZ production, w i indicates the simulated weight of the ith event, 3 and f i and f i are the transfer factors depending on p T and η of the leptons passing the reversed selection. In differential distributions, the yields in Eq. (1) are considered separately in each bin. Systematic uncertainties are applied to account for statistical fluctuations of the measured transfer factors, and for the simplification that the origins, fractions and selection efficiencies of misidentified leptons are assumed equal in the sample where the transfer factors are determined and in the background control sample. The latter uncertainties are estimated using transfer factors obtained from simulation of the different background processes and taking the difference between the result and the nominal method as the uncertainty. An additional uncertainty due to the modeling of the ZZ contamination in the background control sample is estimated by varying N ZZ lllr and N ZZ llrr up and down by 50%. The final total uncertainty is 100% (71%, 95%) in the 4e (2e2µ, 4µ) channel. The misidentified-lepton background amounts to 58% of the total background estimate. As a cross-check, the background is also estimated using an independent method in which ZZ events with one same-flavor same-charge lepton pair as one of the Z boson candidates are selected. The results are found to agree well with the nominal method, differing by less than one standard deviation in all channels.
Background from two single Z bosons produced in different pp collisions in the same bunch crossing is estimated by considering the Z boson production cross sections and the probability of the two primary vertices lying so close to each other that the detector cannot resolve them as separate vertices. It is found to be negligible (< 0.1% of the total signal prediction).
The observed and predicted event yields for signal and background are shown in Table 2. The prediction uncertainties are discussed in Section 7. Figure 4 shows the distributions of data and predictions for the mass and transverse momentum of the four-lepton system, the transverse momentum of the leading Z boson candidate, and the jet multiplicity. The data and the nominal Sherpa prediction agree well. The prediction using Powheg + Pythia to simulate the qq-initiated process tends to underestimate the normalization slightly, which can be understood from its lack of higher-order real-emission corrections that Sherpa implements. Powheg + Pythia also provides a worse description of high jet multiplicities, as it only describes one parton emission at matrix-element level.  Total prediction (Powheg + Pythia 193 ± 11 456 ± 24 286 ± 17 934 ± 50 with higher-order corrections, Sherpa) Table 2: Observed and predicted yields, using the nominal Sherpa setup for the signal predictions. All statistical and systematic uncertainties are included in the prediction uncertainties. An alternative total prediction is given, using Sherpa reweighted to the total NNLO prediction from Matrix with NLO EW corrections, adding the contribution of the EW-ZZ j j process generated with Sherpa, to predict the signal yield. A second alternative total prediction, identical to the nominal Sherpa setup, except using Powheg + Pythia with NNLO QCD and NLO EW corrections applied event by event to simulate the qq-initiated process, is shown at the bottom.   Events / GeV

Systematic uncertainties
The sources of systematic uncertainty are introduced below. Their effects on the predicted integrated signal yields after event selection are shown in Table 3.
For leptons and jets, uncertainties of the momentum or energy scale and resolution are considered [50, 51, 74]. Uncertainties of the lepton reconstruction and identification efficiencies [49, 50] as well as the efficiency of the jet vertex tagging requirements [75,76] in the simulation are taken into account. All of the above depend on the p T and η of the lepton or jet. The electron efficiency uncertainties contain contributions associated with the basic reconstruction, the identification, and the isolation. In addition to correlated components, each is split into O(10) uncorrelated components to take into account the partial decorrelation between individual electrons in different η-p T regions. For muons, the efficiency uncertainties associated with individual muons are treated as fully correlated. This leads to a larger uncertainty for muons than for electrons. As the selection is fully jet-inclusive, jet uncertainties do not affect the integrated yields and are therefore not shown in Table 3.
The pileup modeling uncertainty is assessed by varying the number of simulated pileup interactions. The variations are designed to cover the uncertainty of the ratio of the predicted and measured cross section of non-diffractive inelastic events producing a hadronic system of mass m X > 13 GeV [77]. The uncertainty of the integrated luminosity is 3.2%. It is derived from a preliminary calibration of the luminosity scale using a pair of x-y beam-separation scans performed in August 2015 and May 2016, following a methodology similar to that detailed in Ref.  Table 3. A predicted theoretical modeling uncertainty is taken into account in the unfolding of differential cross sections. It is estimated by using Powheg + Pythia instead of Sherpa to generate the qq-initiated subprocess, and taking the absolute deviation of the result obtained with this setup from the one obtained with the nominal Sherpa setup as an uncertainty, symmetrizing it with respect to the nominal value. This contribution is not shown in Table 3, because it is never applied to yields, where it would be dominated by cross-section normalization differences rather than differences in the reconstruction efficiencies. A further source of uncertainty are statistical fluctuations in the MC samples. The associated uncertainty of the measured differential cross sections is < 1% in most bins, reaching 3% in rare cases.
In the search for aTGCs, an additional uncertainty due to the factorization approximation of NNLO QCD and NLO EW corrections is applied as follows. Following a criterion motivated in Ref.
[81], events are classified as having high QCD activity if i p T, i > 0.3 i | p T, i |, where the sums are over fiducial leptons. In events with high QCD activity, the NLO EW K-factors are in turn not applied and applied with doubled deviation from unity, as 1 + 2(K-factor − 1). The deviations from the nominal result are taken as uncertainties.
The uncertainty of the misidentified-lepton background is described in Section 6. A 30% normalization uncertainty is applied for triboson and tt Z backgrounds with four genuine leptons to account for the crosssection uncertainty. The background uncertainties are considered uncorrelated with other sources.
The propagation of uncertainties in the unfolding as well as the estimation of unfolding-specific uncertainties is described in Section 9.  Table 3: Relative uncertainties in percent of the predicted integrated signal yields after event selection, derived using the nominal Sherpa setup. All uncertainties are rounded to one decimal place.

Integrated cross section
The integrated fiducial cross section σ fid is determined by a maximum-likelihood fit in each channel separately as well as for all channels combined. The expected yield in each channel i is given by where L is the integrated luminosity, and N bkg is the expected background yield. The factor C ZZ is applied to correct for detector inefficiencies and resolution effects. It relates the background-subtracted number of selected events to the number in the fiducial phase space. C ZZ is defined as the ratio of generated signal events satisfying the selection criteria using reconstructed objects to the number satisfying the fiducial criteria using the particle-level objects defined in Section 4.1. It is determined with the nominal Sherpa setup. The C ZZ value and its total uncertainty is determined to be 0.494 ± 0.015 (0.604 ± 0.017, 0.710 ± 0.027) in the 4e (2e2µ, 4µ) channel. The dominant systematic uncertainties come from the uncertainties of the lepton reconstruction and identification efficiencies in the simulation, the choice of MC event generator, QCD scales and PDFs, and the modeling of pileup effects. Other smaller uncertainties come from the scale and resolution of the lepton momenta as well as statistical fluctuations in the MC sample. Table 4 gives a breakdown of the systematic uncertainties.
The likelihood function to be minimized in the cross-section fit is defined as  where is the probability of observing N obs events given that the yield follows a Poisson distribution with mean N exp , and L corr and L uncorr are products of Gaussian nuisance parameters corresponding to the uncertainties of L, C ZZ , and N bkg . The term L corr contains the nuisance parameters that are fully correlated between channels, i.e. all except the statistical uncertainties, while L uncorr contains those that are uncorrelated, i.e. the statistical uncertainties of C ZZ and N bkg in each channel. Nuisance parameters corresponding to different sources of systematic uncertainty are considered uncorrelated. In the combined cross-section fit, the product over channels is taken in the likelihood function shown in Eq. (2), fixing the relative contributions of the signal channels to their theoretically predicted values. Table 5 shows the integrated fiducial cross sections for each channel as well as all channels combined, along with a theoretical prediction. Measurements and predictions agree within approximately one standard deviation, except for the 4e channel, where they agree within approximately 2.5 standard deviations. The sum of the 4e and 4µ cross sections is not equal to the 2e2µ cross section. This is because of interference in the 4e and 4µ channels, the bias caused by the pairing prescription in the fiducial definition, as well as other small differences in the fiducial selection (different ∆R( i ( ) j ) requirement, m > 5 GeV for any same-flavor opposite-charge pair). Figure 5 shows the ratio of measured over predicted cross sections. The goodness of the combined cross-section fit is assessed, taking as hypothesis that the relative contributions of the channels are as predicted. This assumes lepton universality in Z → + − , which is experimentally confirmed to high precision [82,83]. Using the maximum likelihood for the observed yields, L obs , and for the expected yields, L exp , the ratio −2 ln(L obs /L exp ) is found to be 8.7. The p-value is calculated as the fraction of 10 5 MC pseudoexperiments giving a larger ratio than the fit to data, and is found to be 2.3%. This relatively low p-value is driven by the compatibility of the 4e channel with the other two channels.   Figure 5: Comparison of measured integrated fiducial cross sections to a SM prediction based on an NNLO calculation from Matrix with the gg-initiated contribution multiplied by a global NLO correction factor of 1.67. A global NLO EW correction factor of 0.95 is applied, except to the gg-initiated loop-induced contribution, and the contribution of around 2.5% from EW-ZZ j j generated with Sherpa is added. For the prediction, the QCD scale uncertainty is shown as one-and two-standard-deviation bands.

Extrapolation to total phase space and all Z boson decay modes
Extrapolation of the cross section is performed to a total phase space for Z bosons with masses in the range from 66 GeV to 116 GeV and any SM decay. The total phase space is the same as the fiducial phase space (Section 4.1), except that no p T , η, and ∆R requirements are applied to the leptons. The ratio of the fiducial to total phase-space cross section is determined using the Matrix setup described in Section 3 and found to be A ZZ = 0.58 ± 0.01, where the uncertainty includes the following contributions. A similar value is found when the calculation is repeated with the nominal Sherpa setup, and the difference between these (1.0% of the nominal value) is included in the uncertainty of A ZZ . Other included uncertainties are derived from PDF variations (0.4%, calculated with MCFM) and QCD scale variations (0.8%).
To calculate the extrapolated cross section, the combined fiducial cross section is divided by A ZZ and by the leptonic branching fraction 4 × (3.3658%) 2 [58], where the factor of four accounts for the different flavor combinations of the decays. The result is obtained using the same maximum-likelihood method as for the combined fiducial cross section, but now including the uncertainties of A ZZ as additional nuisance parameters. The used leptonic branching fraction value excludes virtual-photon contributions. Based on a calculation with Pythia, including these would increase the branching fraction ZZ → + − + − by about 1-2%.
The extrapolated cross section is found to be 17.3 ± 0.9 [±0.6 (stat.) ±0.5 (syst.) ±0.6 (lumi.)] pb. The NNLO prediction from Matrix, with the gg-initiated process multiplied by a global NLO correction factor of 1.67 [3] is 16.9 +0.6 −0.5 pb, where the uncertainty is estimated by performing QCD scale variations. A comparison of the extrapolated cross section to the NNLO prediction as well as to previous measurements is shown in Figure 6.  [84][85][86], and to pure NNLO predictions from Matrix (with no additional higher-order corrections applied). The total uncertainties of the measurements are shown as bars. Some data points are shifted horizontally to improve readability.

Differential cross sections
Differential cross sections are obtained by counting candidate events in each bin of the studied observable, subtracting the expected background, and unfolding to correct for detector effects. The unfolding takes into account events that pass the selection but are not in the fiducial phase space (which may occur due to detector resolution or misidentification), bin migrations due to limited detector resolution, as well as detector inefficiencies. To minimize the model dependence of the measurement, the unfolding corrects and extrapolates the measured distributions to the fiducial phase space, rather than extrapolating to nonfiducial regions. For each given observable distribution, all of the above detector effects are described by a response matrix R whose elements R i j are defined as the probability of an event in true bin j being observed with the detector in bin i. The response matrix therefore relates the true distribution t and the background-subtracted measured distribution m, Two examples of response matrices are shown in Figure 7. The purity, defined as the fraction of events that are reconstructed in their true bin, is greater than 70% for jet-inclusive observables, except in very few bins. In jet-exclusive observables, the purity is greater than 60% in most bins, but drops to as low as 35% in some bins. This is due to contribution from jets originating from or contaminated by pileup interactions, as well as worse jet energy resolution and poorer knowledge of the jet energy scale than is the case for leptons.
The unfolding is performed by computing the inverse of the response matrix, using regularization to numerically stabilize the solutions, decreasing their statistical uncertainty at the cost of a small regularization bias. An iterative unfolding method based on Bayes' theorem [87] is used, which combines the measured distribution with the response matrix to form a likelihood and takes the predicted true distribution as prior. It applies Bayes' theorem iteratively, using the posterior distribution as prior for the next iteration, each iteration decreasing the dependence on the initial prior. Depending on the observable, either two or three unfolding iterations are performed. In each case, the number is optimized to minimize the overall uncertainty. More iterations lead to higher statistical uncertainty and fewer iterations to higher unfolding method uncertainty due to stronger dependence on the theoretical prediction of the underlying distribution.
The nominal response matrices, corrections, and priors are obtained using the nominal Sherpa setup. Reconstructed objects in the MC simulation are not required to have a corresponding generated object.  The statistical uncertainty due to fluctuations in the data is estimated by generating 2000 sets of random pseudodata following a Poisson distribution in each bin whose expectation value is the number of observed data events in that bin. The unfolding is repeated with the sets of pseudodata, taking the root mean square of the deviation of the resulting unfolded spectrum from the actual unfolded data as the statistical uncertainty in each bin. Another uncertainty due to statistical fluctuations in the MC simulations used to obtain the response matrix is obtained the same way, repeating the unfolding using randomly fluctuated copies of the response matrix.
Experimental and theoretical-modeling uncertainties are estimated by repeating the unfolding with the varied response matrix and taking the deviation from the nominal of the resulting unfolded distribution as the uncertainty. Background uncertainties are estimated by subtracting the varied background predictions from the data before unfolding.
The uncertainty due to imperfect modeling of the observable by MC simulation as well as the inherent bias of the unfolding stemming from regularization is estimated using a data-driven method [88]. The initial priors are reweighted by a smooth polynomial function such that there is very good agreement between the prior folded with the response matrix and the observed data. The folded reweighted prior is unfolded using the nominal response matrix. The deviations of the obtained unfolded distribution from the reweighted prior are used as the unfolding bias uncertainty in each bin. This uncertainty is smaller than 1% of the cross section in almost all bins, but reaches 22% in individual bins (such as the first bin of the mass of the two leading-p T jets, where the modelling of the data is poor).
The unfolding is repeated using Powheg + Pythia instead of Sherpa to model the qq-initiated process and the difference between the unfolded distributions obtained in this way is assigned as an additional systematic event generator uncertainty.
The statistical uncertainty of the data is typically in the range of 5-40% of the cross section. It dominates the total uncertainty in most bins. In jet-inclusive observables, the largest systematic uncertainty comes from the modeling of the response matrix (up to approximately 25%). In jet-exclusive observables, the jet energy scale uncertainty is an additional large contribution (3-23%). Figure 8 shows detailed bin-by-bin uncertainties for selected observables.    Figure 9(a) shows the transverse momentum of the four-lepton system, p T, 4 . The cross section has a peak around 10 GeV and drops rapidly towards both lower and higher values. At low p T, 4 , the resummation of low-p T parton emissions is important and fixed-order descriptions are inadequate. For this reason, the fixed-order predictions are not shown in the first two bins, 0-5 GeV and 5-15 GeV. The region below p T, 4 = 60 GeV is modeled slightly better by predictions that include a parton shower, again suggesting the importance of resummation. Above 60 GeV, the fixed-order NNLO predictions describe the data slightly better. Figure 9(b) shows the absolute rapidity of the four-lepton system, which drops gradually towards high values. This distribution is potentially sensitive to a different choice of PDF, describing the momentum distribution of the incoming partons. Fixed-order calculations and predictions including a parton shower model this observable reasonably well, within the statistical and systematic uncertainties. The predictions tend to slightly underestimate the cross sections for small values of |y 4 |. : Measured and predicted differential cross sections for (a) the transverse momentum and (b) the absolute rapidity of the four-lepton system. The statistical uncertainty of the measurement is shown as error bars, and shaded bands indicate the systematic uncertainty and the total uncertainty obtained by summing the statistical and systematic components in quadrature. The ratio plots only show the total uncertainty. A pure NNLO calculation from Matrix is shown with no additional corrections applied. The best SM prediction is based on this NNLO calculation, with the gg-initiated contribution multiplied by a global NLO correction factor of 1.67. For the p T, 4 distribution in (a), the NLO EW correction is applied as a global factor of 0.95 as a differential calculation is not available. For the |y 4 | distribution in (b), an NLO EW correction factor is applied in each bin. The contribution from EW-ZZ j j generated with Sherpa is added. For the fixed-order predictions, the QCD scale uncertainty is shown as a shaded band. Parton-showered Sherpa and Powheg + Pythia predictions are also shown. For better visualization, the last bin is shown using a different x-axis scale where indicated by the dashed vertical line. Figure 10(a) presents the azimuthal angle separation between the two Z boson candidates. The fixed-order predictions only describe the shape of the gg-initiated process at LO and therefore predict a distribution that is more peaked at π than those from Sherpa and Powheg + Pythia, where the parton shower shifts some events towards lower values. Figure 10(b) shows the distribution of the absolute rapidity difference of the two Z boson candidates, which drops towards high values and is modeled by all calculations to within the uncertainties. A pure NNLO calculation from Matrix is shown with no additional corrections applied. The best SM prediction is based on this NNLO calculation, with the gg-initiated contribution multiplied by a global NLO correction factor of 1.67. For the δφ Z 1 , Z 2 distribution in (a), the NLO EW correction is applied as a global factor of 0.95 as a differential calculation is not available. For the |δy Z 1 , Z 2 | distribution in (b), an NLO EW correction factor is applied in each bin. The contribution from EW-ZZ j j generated with Sherpa is added. For the fixed-order predictions, the QCD scale uncertainty is shown as a shaded band. Parton-showered Sherpa and Powheg + Pythia predictions are also shown. For better visualization, the last bin is shown using a different x-axis scale where indicated by the dashed vertical line. Figure 11 shows the transverse momentum of the leading-p T and subleading-p T Z boson candidates, exhibiting a wide peak around 50 GeV and 30 GeV, respectively. Anomalous triple gauge couplings (as discussed in Section 10) would manifest as an excess in the cross section at large values of the transverse momentum of the Z bosons, which is not observed in these differential cross-section distributions (the last bin in each distribution is consistent with the SM predictions). The discrepancies at p T of about 50 GeV, 90 GeV in the leading Z boson candidate are related to the excesses seen in Figure 4(c). The local significance of these excesses with respect to the Sherpa prediction is estimated to be 2.3 and 2.0 standard deviations respectively. This estimate is based on the corresponding bins in the measured distribution before unfolding, as the statistical treatment is simpler due to the statistical uncertainties being Poissonian.
In the estimation, both the predicted and observed yields are taken to be Poissonian and the systematic uncertainties of the prediction are taken into account.  Figure 11: Measured and predicted differential cross sections for the transverse momentum of (a) the leading-p T and (b) the subleading-p T Z boson candidates. The statistical uncertainty of the measurement is shown as error bars, and shaded bands indicate the systematic uncertainty and the total uncertainty obtained by summing the statistical and systematic components in quadrature. The ratio plots only show the total uncertainty. A pure NNLO calculation from Matrix is shown with no additional corrections applied. The best SM prediction is based on this NNLO calculation, with the gg-initiated contribution multiplied by a global NLO correction factor of 1.67. An NLO EW correction factor is applied in each bin. The contribution from EW-ZZ j j generated with Sherpa is added. For the fixed-order predictions, the QCD scale uncertainty is shown as a shaded band. Parton-showered Sherpa and Powheg + Pythia predictions are also shown. For better visualization, the last bin is shown using a different x-axis scale. The scale change is indicated by the dashed vertical line. Figure 12 presents the transverse momenta of the leptons in the final selected quadruplet. From the highest-p T to the lowest-p T lepton, the distribution becomes less peaked and more symmetric about the peak, while the position of the peak shifts from ∼60 GeV to ∼50 GeV, then ∼35 GeV, and finally ∼25 GeV. All lepton p T distributions agree well with the predictions.  Figure 12: Measured and predicted differential cross sections with respect to the transverse momenta of the leptons in the final selected quadruplet, in descending order of transverse momentum. A pure NNLO calculation from Matrix is shown with no additional corrections applied. The best SM prediction is based on this NNLO calculation, with the gg-initiated contribution multiplied by a global NLO correction factor of 1.67. An NLO EW correction factor is applied in each bin. The contribution from EW-ZZ j j generated with Sherpa is added. For the fixed-order predictions, the QCD scale uncertainty is shown as a shaded band. Parton-showered Sherpa and Powheg + Pythia predictions are also shown. For better visualization, the last bin is shown using a different x-axis scale. The scale change is indicated by the dashed vertical line. Figure 13 shows the jet multiplicity distributions as well as the scalar sum of the transverse momenta of all selected jets. Powheg + Pythia shows a clear trend towards underestimating the cross section at jet multiplicities greater than one and large jet scalar p T sum, which is expected, because in Powheg + Pythia only the hardest parton emission is included at the matrix-element level. Sherpa, however, includes up to three parton emissions at the matrix-element level, and exhibits good agreement with the measurements for these higher jet multiplicities. The central-jet multiplicity in Figure 13(b) is an exception, as Powheg + Pythia describes it slightly better than Sherpa. The most significant observed disagreement is the deficit in the bin 60 GeV < p T < 90 GeV of the jet scalar p T sum. It has a local significance of 2.3 standard deviations with respect to the Sherpa prediction, estimated from the corresponding bins in the measured distribution before unfolding.  Figure 14 shows the transverse momentum and absolute pseudorapidity of the leading-p T and subleadingp T jets. Within the relatively large uncertainties, Sherpa provides a good description of the kinematics. Powheg + Pythia also describes the shapes of the |η| distributions well, while its normalization is too low for the subleading-p T jet. Powheg + Pythia does not describe the p T distribution of either the leading-p T or subleading-p T very well. A deficit of events is observed in the bin 2.0 < |η| < 2.5 of the subleading-p T jet. Its local significance with respect to the Sherpa prediction is estimated to be 3.2 standard deviations, based on the corresponding bins in the measured distribution before unfolding.  Figure 15 shows the rapidity difference and invariant mass of the two leading-p T jets. The EW-ZZ j j production process predicted by Sherpa is shown separately, in addition to the process-inclusive predictions from Sherpa and Powheg + Pythia. This contribution falls much less steeply towards higher values of the presented observables. The contribution from this process in the last bins in each distribution improves the agreement between prediction and measurement, demonstrating the importance of this process at these ends of the kinematic phase space.  Figure 15: Measured and predicted differential cross sections for (a) the absolute difference in rapidity between the two leading-p T jets and (b) the invariant mass of the two leading-p T jets. The statistical uncertainty of the measurement is shown as error bars, and shaded bands indicate the systematic uncertainty and the total uncertainty obtained by summing the statistical and systematic components in quadrature. The ratio plots only show the total uncertainty. Parton-showered Sherpa and Powheg + Pythia predictions are shown. In addition to the process-inclusive predictions, the EW-ZZ j j production process predicted by Sherpa is shown separately. It is also shown normalized to the process-inclusive Sherpa prediction to facilitate shape comparisons between the EW-ZZ j j subprocess and the process-inclusive ZZ production. For better visualization, the last bin is shown using a different x-axis scale. The scale change is indicated by the dashed vertical line.

Search for anomalous triple gauge couplings
The search for aTGCs uses the reconstructed transverse momentum of the leading-p T Z boson candidate (p T, Z 1 ) to look for deviations of the data from the SM, as this variable is found to provide the highest sensitivity to their predicted effects. The four-lepton mass provides similar sensitivity, but is not used, because no dedicated calculation of NLO EW corrections for pp → ZZ → + − + − production binned in the four-lepton mass is available for the analysis.
The considered aTGC signal model uses an effective vertex function approach [89]. It includes two coupling strengths that violate charge-parity (CP) symmetry, f γ 4 and f Z 4 , as well as two CP-conserving ones, f γ 5 and f Z 5 . No unitarizing form factor is used, as the sensitivity of the measurement is well within the unitarity bounds.
The expected number of events N in the aTGC search is parameterized in terms of the coupling strengths, on which it depends both linearly and quadratically, where N SM is the SM expectation and the N i j are yield coefficients that depend on the final-state particle momenta. To determine the coefficients N i j , 2 × 10 5 events with aTGC are generated at LO with one fixed reference set of coupling strengths using Sherpa and the CT10 PDF set. Based on the kinematic properties of each event, the coefficients N i j are extracted using a framework [90] based on the BHO program [91].
The yield for all other values of the coupling strengths can then be calculated using Eq. (3).
The SM prediction N SM is constructed separately using the highest-order calculations available. The nominal Sherpa setup is used, except that the qq-initiated process is generated with Powheg + Pythia and each event reweighted by NNLO QCD and NLO EW corrections binned in p T, Z 1 . The SM ZZ predictions, estimated backgrounds, as well as observed yields are shown in Table 6 as a function of p T, Z 1 . These contributions are also shown in Figure 16 together with two different aTGC predictions. The considered systematic uncertainties of the predictions are the same as in the integrated cross-section measurement. An additional uncertainty associated with the combination of NNLO QCD and NLO EW corrections for the SM ZZ → + − + − process is assigned as described in Section 7, ranging from ∼1% in the lowest to ∼10% in the highest p T, Z 1 bin. The p T, Z 1 bins are optimized using the predictions to maximize the expected sensitivity.
The data are found to be consistent with the SM predictions, and no indication of aTGCs is observed. Confidence intervals of aTGC parameters are determined using the expected and observed yields in bins of p T, Z 1 as reconstructed by the detector.
A frequentist method [92] is used to find the 95% confidence level (CL) intervals for the aTGC parameters. The predicted and observed yields are assumed to follow Poissonian probability density functions, while the systematic uncertainties are treated as nuisance parameters constrained by Gaussian functions. The 9.2 ± 2.8 0.43 ± 0.13 0.15 ± 0.05 0.078 ± 0.028 Misid. lepton background 12 ± 8 0.17 ± 0.11 < 0.1 < 0.1 Table 6: Observed and predicted yields in bins of the transverse momentum of the leading-p T Z boson candidate. All statistical and systematic uncertainties are included in the prediction uncertainties, including the uncertainty associated with the combination of NNLO QCD and NLO EW corrections for the SM ZZ → + − + − process.  expected confidence intervals and their one-and two-standard-deviation confidence bands are established using many independent sets of randomly generated pseudodata following a Poisson distribution whose expectation value is the SM prediction in each bin.
Confidence intervals are set for each coupling strength individually, setting all others to zero, using 2500 sets of pseudodata. The expected and observed 95% CL intervals are listed in Table 7. The onedimensional confidence intervals are more stringent than those derived in previous measurements by the ATLAS and CMS collaborations [14,15,19] and at the Tevatron and LEP colliders [20,21]. In addition, two-dimensional 95% CL intervals are obtained by allowing pairs of aTGC parameters to vary simultaneously, while setting the others to zero, using 26000 sets of pseudodata. They are shown in Figure 17.  Confidence intervals are also provided for parameters of the effective field theory (EFT) in Ref. [93], which includes four dimension-8 operators describing aTGC interactions of neutral gauge bosons. The coefficients of the operators are denoted CB W /Λ 4 , C BW /Λ 4 , C WW /Λ 4 , and C BB /Λ 4 , where Λ is the energy scale of the new physics described by the EFT. They can be linearly related to the parameters f γ 4 , f Z 4 , f γ 5 , and f Z 5 as described in Ref. [94]. Thus Eq.
(3) can be reformulated in terms of the EFT coefficients and confidence intervals set in the same way as for the coupling strengths. The resulting one-dimensional EFT confidence intervals are given in Table 8.
Two-dimensional EFT confidence intervals are shown in Figure 18.  Table 8: One-dimensional expected and observed 95% CL intervals of EFT parameters using the transformation from Ref. [94]. Each limit is obtained setting all other EFT parameters to zero.   Figure 18: Observed and expected two-dimensional 95% CL intervals in planes of different pairs of EFT parameters using the transformation from Ref. [94]. The EFT parameters other than those shown are set to zero. The black straight lines indicate the observed one-dimensional confidence intervals at 95% CL.

Conclusion
The production of pairs of Z bosons is studied in the ZZ → + − + − channel in 13 TeV proton-proton collisions produced at the LHC and recorded with the ATLAS detector, using data corresponding to an integrated luminosity of 36.1 ± 1.1 fb −1 . Integrated fiducial cross sections are measured separately in the three decay channels 4e, 2e2µ, and 4µ as well as in their combination. They are found to agree well with NNLO SM predictions with NLO QCD corrections for the gg-initiated production process as well as NLO EW corrections applied. The combined cross section is extrapolated to a total phase space and all SM Z boson decays. Differential cross sections are measured for 20 observables. They are compared to NLO predictions with parton shower, to fixed-order NNLO predictions, and to fixedorder predictions from calculations at the highest known orders for the different subprocesses (NNLO pp → ZZ, NLO gg → ZZ, NLO EW corrections, electroweak pp → ZZ j j). In general, the predictions describe the observables reasonably well. Using the transverse momentum of the leading-p T Z boson candidate, confidence intervals are obtained for parameters of aTGCs forbidden at tree level in the SM, both parameterized as aTGC coupling strengths and in an effective field theory approach. No significant deviations from the SM are observed.