Measurement of $ZZ$ production in the $\ell\ell\nu\nu$ final state with the ATLAS detector in $pp$ collisions at $\sqrt{s} = 13$ TeV

This paper presents a measurement of $ZZ$ production with the ATLAS detector at the Large Hadron Collider. The measurement is carried out in the final state with two charged leptons and two neutrinos, using data collected during 2015 and 2016 in $pp$ collisions at $\sqrt{s} = 13$ TeV, corresponding to an integrated luminosity of 36.1 fb$^{-1}$. The integrated cross-sections in the total and fiducial phase spaces are measured with an uncertainty of 7% and compared with Standard Model predictions, and differential measurements in the fiducial phase space are reported. No significant deviations from the Standard Model predictions are observed, and stringent constraints are placed on anomalous couplings corresponding to neutral triple gauge-boson interactions.


Introduction
In the Standard Model (SM), the production of gauge boson pairs has a profound connection with the non-Abelian nature of the electroweak (EW) theory and with the spontaneous breaking of the EW gauge symmetry. In addition, a broad range of new phenomena beyond the SM (BSM) are predicted to reveal themselves through diboson production. The study of diboson production probes a cornerstone of the EW theory and possible BSM physics scenarios, and it constitutes a salient component of the physics programme at the Large Hadron Collider (LHC). Among all the diboson processes, the production of two on-shell Z bosons has the smallest cross-section, but is nevertheless quite competitive for measurements and searches, because of its generally good signal-to-background ratio for the fully leptonic decay channels. For instance, the ZZ process is a leading channel to search for anomalous neutral triple-gauge-boson couplings (aTGCs) [1] and to study the off-shell production of the Higgs boson [2, 3]. Figure 1 shows representative Feynman diagrams for ZZ production at the LHC. The dominant process is t-channel production with a quark and anti-quark initial state,  (b) q g Z Z q (c) Figure 1. Representative Feynman diagrams for ZZ production at the LHC: (a) lowest-order t-channel qqZZ production; (b) production of ZZ plus one parton through the qq initial state; (c) production of ZZ plus one parton through the qg initial state; (d) ggZZ production with a fermion loop; (e) ggZZ production involving an exchange of a virtual Higgs boson; (f) s-channel production with aTGCs.
hereafter denoted by the qqZZ process. Higher-order QCD corrections to the qqZZ process are found to be sizeable [4], and two tree-level diagrams concerning production of two Z bosons and one outgoing parton are shown. The gluon fusion process (ggZZ) includes two sub-processes, one with a fermion loop and the other involving a virtual Higgs boson. Although the ggZZ process only appears at O(α 2 S ), it nevertheless has a non-negligible contribution of O(10%) to the total ZZ production rate due to the large gluon flux at the LHC. The s-channel production is forbidden at the lowest order; however, the neutral TGCs can still acquire small values of O(10 −4 ) in the SM, due to the correction with a fermion loop [5]. The observation of aTGCs with larger values would hint at the existence of new physics.
Measurements of ZZ production at the LHC have been carried out in two decay final states, one with four charged leptons (4 ) and the other with two charged leptons and two neutrinos ( νν). Using LHC Run-1 and Run-2 data, multiple results [6-12] have been reported by the ATLAS and CMS experiments. The most precise results to date have been obtained from the 4 channel using 13 TeV data [8,12], where the integrated production cross-section has been measured to a precision of 5% and the upper bound on neutral aTGC parameters has been reduced to 10 −3 . The improved experimental precision -2 -
This paper presents a measurement of ZZ production using 36.1 fb −1 of data collected with the ATLAS detector in pp collisions at √ s = 13 TeV. This analysis is performed in the νν ( = e or µ) final state, which has a larger branching fraction but suffers from higher background contamination in comparison with the 4 channel. To ensure a good signal-tobackground ratio, the experimental selection requires one Z boson boosted against the other in the transverse plane, which results in a pair of high-p T isolated leptons and significant missing transverse momentum (E miss T ). The νν channel thus offers higher data statistics than the 4 channel for events with high-p T Z bosons, and offers competitive precision for integrated and differential measurements, as well as good sensitivity to aTGCs.
The dominant background arises from W Z production where the Z boson decays into a pair of charged leptons. About 60% of the W Z events which contribute to the νν final state have the W boson decaying leptonically (W → ν or W → τ ν → + 3ν, = e or µ), where the final-state lepton escapes detection. The remaining 40% W Z contribution is related to the W → τ ν decay with subsequent hadronic decays of the τ -lepton. Another important background comes from the processes that genuinely produce the νν final state but contain a lepton pair not originating directly from a Z-boson decay. This background, referred to as the non-resonant-background, includes W W , top-quark (tt and W t), and Z → τ τ production. The production of a Z boson in association with jet(s) (Z + jets) also constitutes a potentially large background source. The Z + jets events with large "fake" E miss T arise from heavy-flavour hadron decays in the accompanying jet(s), from jet mismeasurements in certain regions of the detector, and from the measurement resolution itself, owing to the additional pp collisions in the same or neighbouring proton bunch crossings (pile-up). The ZZ → 4 process yields a small contribution when one lepton pair misses detection, and it is considered as a background in this measurement. Finally, minor background contributions are expected from three-boson production (V V V with V = W or Z) and production of tt accompanied by one or two vector bosons (ttV ).
The integrated cross-section of ZZ production is measured in a fiducial phase space and then extrapolated to a total phase space. The determination of the fiducial (σ fid ZZ→ νν ) and total (σ tot ZZ ) cross-sections is obtained as shown in eq. (1.1): where C ZZ stands for an overall efficiency correction factor, A ZZ is the fiducial acceptance, and B is the branching fraction of the ZZ → νν ( = e, µ) decay. The signal yield N obs ZZ is determined through a fit to the observed E miss T spectrum, which leads to improved sensitivity compared with a simple event-counting method. The A ZZ (C ZZ ) factor is calculated as N exp,fid JHEP10(2019)127 events arising from the ZZ → τ τ νν decays with the subsequent τ → νν decays of both τ -leptons are considered as signal events at detector level but excluded in the fiducial measurements. Throughout this paper, "Z → " denotes the decays of a Z boson or a virtual photon into a charged-lepton pair. Furthermore, differential cross-sections are reported in the fiducial region for eight kinematic variables, which are sensitive to effects from higher-order corrections and possible BSM physics. These variables include the transverse momentum of the leading lepton (p 1 T ), the leading jet (p jet1 T ), the dilepton system (p T ), and the ZZ system (p ZZ T ), the transverse mass of the ZZ system (m ZZ T ), 1 the absolute rapidity of the dilepton system (|y |), the azimuthal angle difference between the two leptons (∆φ ), and the number of jets (N jets ). Since no significant deviations from the SM are observed, upper limits are placed on the aTGC parameters [1], which typically manifest themselves as a signal excess growing rapidly as the partonic centre-of-mass energy ( √ŝ ) increases. In this analysis, aTGCs are searched for using the p T spectrum in the fiducial region, motivated by the fact that p T is correlated with √ŝ and has a good experimental resolution.

ATLAS detector
The ATLAS detector [17][18][19] is a large multi-purpose detector with a forward-backward symmetric cylindrical geometry and nearly 4π coverage in solid angle. 2 It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large superconducting toroidal magnets each having eight coils assembled radially and symmetrically around the beam axis. The inner-detector system (ID) is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range |η| < 2.5. A high-granularity silicon pixel detector covers the vertex region and usually provides four measurements per track. The pixel detector is followed by a silicon microstrip tracker which usually provides four measurement points per track. These silicon detectors are complemented by a transition radiation tracker, which enables radially extended track reconstruction and improved momentum measurements up to |η| = 2.0. The transition radiation tracker also provides electron identification information based on the fraction of hits (typically 30 hits in total) above a high-energy threshold designed for optimal electron-pion separation.
The calorimeter system covers the pseudorapidity range |η| < 4.9. Within the region |η| < 3.2, electromagnetic calorimetry is provided by barrel and endcap lead/liquid-argon (LAr) sampling calorimeters, with an additional thin LAr presampler covering |η| < 1.8, 2 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the centre 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). Table 1. Summary of Monte Carlo event simulation tools with their theoretical accuracy for each process, where "p" stands for parton(s). For the first two signal processes and the W Z process, the first (second) row describes the baseline (alternative) simulation. The theoretical accuracy of the normalisation used for the total production cross-section of each process is shown in the last column.
is about −4% at low m ZZ but has a larger impact at high m ZZ , which cancels out the positive QCD correction for m ZZ around 500 GeV. An alternative sample was generated with Sherpa2.2.2 [29] using the NNLO NNPDF3.0 PDF [30] and the same choice of QCD scales. The Sherpa generator and its associated UEPS algorithm has NLO QCD accuracy for inclusive observables and extended QCD precision for events with one or more outgoing partons (NLO for up to one parton, LO for two and three partons). The ggZZ events were simulated with the LO gg2vv3.1.6 [31,32] generator using the NNLO CT10 PDF, and then interfaced to Pythia8.186 using the A14 tune [33]. The production cross-section was corrected to NLO QCD accuracy using a K-factor of 1.7 reported in ref. [34]. An alternative modelling was provided by Sherpa2.1.1 [35] with the NLO CT10 PDF, which extended the LO QCD calculation to events with one parton. Both generators used m ZZ /2 for the QCD scales, and they incorporated both the fermion-loop and the Higgs processes, together with the interference between the two.
To study the effects of aTGCs, an additional sample for the SM qqZZ process was generated at NLO in QCD using Sherpa2.1.1 with the NLO CT10 PDF. The simulated sample was interfaced to a parton-level program [1] following the procedures detailed in ref. [36], and then event-by-event weights reflecting the relative change in the cross-sections due to any aTGCs were computed. A parameterisation of aTGC contributions as a function of any kinematic variable can be derived with this information. This procedure was adopted in the previous ZZ measurements [6-8].
Production of ZZ → 4 events was modelled in the same way as the signal events. The diboson background processes W Z and W W were generated with Powheg-Box v2 using the NLO ME calculation and the NLO CT10 PDF, and then interfaced to Pythia8.186 with the AZNLO tune. An alternative W Z sample was produced with Powheg-Box v2 interfaced to Herwig++ [37], for the study of UEPS uncertainties. The interference between the W W and ZZ processes in the νν final state was found to be negligible [16] and was therefore not considered in this analysis. Both the tt and W t events were simulated at NLO in QCD with Powheg-Box v2 [38,39] and interfaced to Pythia6.428 [40], and the production cross-sections were corrected to NNLO QCD precision [41,42]. Sherpa2.2.1 with the NNLO NNPDF3.0 PDF was used to model the Z + jets process. The production cross-section for the Z + jets process was calculated with NNLO QCD precision, while the simulation has NLO QCD precision for events with zero, one and two partons, and provided a LO QCD description for events with three to five partons. The rare V V V background, consisting of W W W , W W Z, W ZZ and ZZZ processes, was modelled with Sherpa2.1.1 with NLO QCD precision. MadGraph5 aMC@NLO [43] interfaced to Pythia8.186 was used to generate the ttV background events that account for ttW , ttZ and ttW W production processes. The ttV process was calculated at LO QCD accuracy, and its production cross-section was corrected to NLO QCD precision [43].
Generated events were then processed through the ATLAS detector simulation [44] based on GEANT4 [45] to emulate the response of the detector to the final-state particles. Pile-up was simulated with Pythia8.186 using the A2 tune [46] and overlaid on simulated events to mimic the real collision environment. The distribution of the average number of interactions per bunch crossing in the simulation was weighted to reflect that in data. Simulated events were processed with the same reconstruction algorithms as for the data. Furthermore, the lepton momentum scale and resolution, and the lepton reconstruction, identification, isolation and trigger efficiencies in the simulation were corrected to match those measured in data.

Selection of νν events
This analysis selects a detector signature with a pair of high-p T isolated electrons (ee) or muons (µµ) and large E miss T . The ee and µµ channels are combined to obtain the final results. The event selection strategy was optimised to cope with the large background contamination. The selection requirements lead to a highly boosted Z boson back-to-back with the missing transverse momentum vector ( E miss T ). Backgrounds are further reduced by removing events with extra leptons or any jets containing b-hadrons ("b-jets"). Therefore, a precise understanding of the overall reconstruction and selection of leptons, jets, and E miss T is required in this measurement. Events are first required to have a collision vertex associated with at least two tracks each with p T > 0.4 GeV. The vertex with the highest scalar p 2 T sum of the associated tracks is referred to as the primary vertex.
Electrons are reconstructed from energy deposits in the EM calorimeter matched to a track reconstructed in the ID. The electron identification imposes selections on the number of hits in the ID and requirements on a likelihood discriminant, built from variables related to EM calorimeter shower shapes, track-cluster matching, track quality, and transition radiation [47]. Electrons must satisfy the "medium" identification criterion [47], which is about -7 -JHEP10(2019)127 90% efficient for electrons with p T ≈ 40 GeV. Candidate electrons must have p T > 7 GeV and pseudorapidity |η| < 2.47. Muons are reconstructed by combining all the hits associated with a pair of matched tracks reconstructed in the ID and MS, taking into account the energy loss in the calorimeter. Muons are identified by requiring a sufficient number of ID and MS hits, and good consistency between the ID and MS track measurements as well as good combined fit quality [48], and they must satisfy the "medium" identification criterion [48], which has an overall efficiency of 96%. Candidate muons are required to have p T > 7 GeV and |η| < 2.5. To further suppress misidentified lepton and cosmic-ray background contributions, the absolute value of the longitudinal impact parameter of leptons with respect to the primary vertex must be smaller than 0.5 mm, and the transverse impact parameter divided by its error must be less than 5 (3) for electrons (muons). In addition, the "loose" isolation criteria defined in refs. [47,48] are applied. The isolation selection imposes requirements on both the track-based and calorimeter-based isolation variables, and varies as a function of p T to maintain a uniform efficiency above 98% for prompt leptons.
Jets are reconstructed with the anti-k t algorithm [49] with radius parameter R = 0.4, using as input positive-energy topological clusters in the calorimeters [50-52]. The jet energy scale is calibrated using simulation and further corrected with in situ methods [51]. Candidate jets must have p T > 20 GeV and |η| < 4.5. Additional requirements using the track and vertex information inside a jet [53] are applied for jets with p T < 60 GeV and |η| < 2.5 to suppress pile-up contributions. Candidate b-jets (p T > 20 GeV and |η| < 2.5) are identified with an algorithm providing 85% signal efficiency and a rejection factor of 33 against light-flavour jets [54].
Leptons and jets may be close to each other or overlapping, even after implementing the full set of object selections. The appearance of such overlapping objects may lead to ambiguities in the event selection and in the energy measurement of the physics objects. A common procedure in the ATLAS experiment, as detailed in ref.
[55], is applied to resolve the ambiguities. This requirement helps to suppress the occurrence of two "problematic" scenarios, one with energy measurement of electrons biased due to nearby jets, and the other with a jet producing non-prompt muons through meson decays in flight.
The E miss T vector is computed as the negative of the vector sum of transverse momenta of all the leptons and jets, as well as the tracks originating from the primary vertex but not associated with any of the leptons or jets ("soft-term") [56]. The soft-term is computed in a way minimising the impact of pile-up in the E miss T reconstruction. Candidate events are preselected by requiring exactly two selected electrons or muons with opposite charges and p T > 20 GeV. The leading lepton is further required to have p T > 30 GeV, well above the threshold of the single-lepton triggers. To suppress the W Z background, events containing any additional lepton satisfying the "loose" rather than "medium" identification requirement, in addition to the other requirements, are rejected. The dilepton invariant mass (m ) is required to be in the range between 76 and 106 GeV, which largely reduces the contamination from the non-resonant-background. Figure 2 shows the observed and expected E miss T spectra after imposing the above requirements ("preselection"). The fractional experimental uncertainties in the expected spectra increase Data / Pred.  Figure 2. The E miss T distributions after the preselection for the ee (left) and µµ (right) channels. The expectation is derived from simulation and the hashed band accounts for the statistical and experimental uncertainties. The experimental uncertainties are described in section 7. The last bin in the distributions contains the overflow events. The arrows indicate that the differences between the data and the expectation in some bins exceed the y-axis scope of the bottom plots. The "other" background corresponds to the V V V and ttV processes. as a function of E miss T in the region dominated by the Z + jets process, as a result of the asymmetric migration effects along the steeply falling E miss T distribution and the large jet-related uncertainty for Z + jets events at high E miss T . The top-quark processes with genuine E miss T dominate the high E miss T region. For E miss T around 200 GeV, top-quark events generally contain less jet activity than Z +jets events: this leads to correspondingly smaller experimental uncertainties. The experimental uncertainties are elaborated in section 7. The data sample after the preselection is dominated by the Z + jets and non-resonantprocesses. To suppress these backgrounds, a further selection based on E miss T and event topology is applied.
Candidate events are required to have E miss T > 110 GeV and V T /S T > 0.65, where V T is the magnitude of the vector sum of transverse momenta of selected leptons and jets, and S T is the scalar p T sum of the corresponding objects. The variable V T /S T was found to be less sensitive to jet experimental uncertainties than similar variables such as E miss T /S T . To further reduce the impact of jet energy scale uncertainties, the calculation of V T and S T uses "hard jets" which are required to have p T > 25 GeV for the central region (|η| < 2.4) and p T > 40 GeV for the forward region (2.4 < |η| < 4.5). The E miss T cut suppresses the Z + jets contamination by many orders of magnitude, and the residual Z + jets events, which have large fake E miss T , are further suppressed by the V T /S T requirement. As the consequence of the combined E miss T and V T /S T requirement, the Z + jets process only constitutes a small fraction of the total background after the full selection.
Additional selection criteria based on angular variables are imposed to ensure the desired detector signature, which helps to further reject the Z + jets and non-resonantbackground events. The azimuthal angle difference between the dilepton system and E miss T , ∆φ( p T , E miss T ), must be larger than 2.2 radians, and the selected leptons must be close to each other, with the distance ∆R = (∆φ ) 2 + (∆η ) 2 < 1.9. Finally, events -9 -
containing one or more b-jets are vetoed to further suppress the tt and W t backgrounds.
The full event selection is summarised in table 2. Figure 3 gives the observed and simulated spectra for V T /S T , ∆R , ∆φ( p T , E miss T ), and the number of b-jets, where each plot is made with the implementation of all the cuts prior to the cut on that variable, according to the cut sequence in table 2.

Total and fiducial phase spaces
The definitions of the total and fiducial phase spaces are summarised in table 3. The total phase space is defined as in ref.
[8] for the ZZ → 4 measurement, requiring 66 < m , m νν < 116 GeV ( = e or µ), where the leptons and neutrinos originate from the Z-boson decays. The four-momenta of the leptons are defined at Born level, i.e. before any QED final-state radiation.
The fiducial phase space is defined with a set of criteria very close to that of the detector-level event selection (table 2). This strategy helps to reduce the amount of phase- Data / Pred.  Data / Pred.

ATLAS
Other bgds. Uncertainty space extrapolation in the fiducial measurements and therefore minimises the theoretical uncertainties of the results. The criteria are applied to "particle-level" physics objects, which are reconstructed from stable final-state particles, prior to their interactions with the detector. For electrons and muons, QED final-state radiation is partly recovered by adding to the lepton four-momentum the four-momenta of surrounding photons not originating from hadrons within an angular distance ∆R < 0.1 (dressed leptons). Particle-level jets are built with the anti-k t algorithm with radius parameter R = 0.4, using all finalstate particles as the input (excluding muons and neutrinos). As shown in table 3, the selection requirements for the fiducial phase space closely follow those in table 2. The E miss T vector is defined as the sum of transverse momenta of the two neutrinos from the Z-boson decays. This measurement requires large E miss T , which has a detector resolution of around 10 GeV [56] 15.7 ± 0.7 pb 13.9 ± 0.4 pb 1.8 ± 0.6 pb Table 4. Predictions for the signal yields at detector level, for the C ZZ and A ZZ coefficients defined in eq. (1.1), and for the cross-sections in the fiducial and total phase spaces. The first column gives the corresponding predictions for the total ZZ process, combined from those shown separately for the qqZZ and ggZZ sub-processes. The errors include both the statistical and systematic uncertainties (see section 7). The statistical, experimental, and theoretical uncertainties are also shown separately for the combined signal yields. the events selected at detector level, the E miss T threshold is therefore lowered to 90 GeV in the fiducial region. The efficiency of the b-jet veto is found to be 98% in the fiducial region and consistent between the Powheg+Pythia8 and Sherpa generators. No requirement is made on the number of b-jets in the fiducial selection. Table 4 gives the expected signal yields at detector level, the A ZZ and C ZZ factors, and the predicted cross-sections. The qqZZ and ggZZ processes have similar final-state kinematic distributions and their A ZZ and C ZZ factors are similar. The corresponding factors for the total ZZ process are averaged from that for the two sub-processes, weighted by the respective cross-sections. The cross-section predictions for the total phase space are corrected for the branching fraction of the ZZ → νν decays, 1.35% with a negligible uncertainty, obtained from refs. [57,58]. The expected fiducial and total cross-sections, σ exp,fid ZZ→ νν and σ exp,tot ZZ , are calculated from simulation, including the higher-order corrections detailed in section 3. The total uncertainties in these predictions are also provided in table 4, and the procedures used to derive these uncertainties are described in section 7.

Background estimation
After the event selection, the overall signal-to-background ratio is about 1.7. The W Z and non-resonant-backgrounds account for 72% and 21% of the total background contribution, respectively, and are estimated from control regions in data. The Z + jets background is largely suppressed, yielding a relative contribution of only 4%, and is estimated from data. Finally, the small remaining contributions from other processes, amounting in total to 3% of the total background, are estimated from simulation. The various background estimates and their uncertainties are described below.
To estimate the dominant resonant background from W Z production, a control region enriched in W Z events, with a purity of 90%, is defined using the preselection criteria, except that a third lepton with p T > 20 GeV and satisfying the medium identification criteria is required. Several further selections such as V T /S T > 0.3, b-jets veto, and m W T > 60 GeV, where m W T is constructed from the third lepton's transverse momentum and the E miss T vector, 3 are applied to suppress non-W Z contributions. A normalisation factor (f W Z ) is calculated in the control region as the number of observed events in data, subtracting the non-W Z contributions estimated from simulation, divided by the predicted W Z yield. The factor f W Z is found to be 1.26 ± 0.04 (stat), which is consistent with a recent W Z measurement [59], performed within a broader fiducial phase space and using a recent calculation of the W Z total cross-section at NNLO in QCD [60,76]. The statistical uncertainty of the data in the control region leads to a 3% uncertainty in the W Z estimate in the signal region. The systematic uncertainty is evaluated for the ratio of the W Z predictions in the signal and control regions. The experimental uncertainty in this ratio is 3.5%, and the theoretical uncertainty is 3.3%, calculated as the sum in quadrature of the PDF, scale, and UEPS uncertainties. The non-W Z contribution in the control region is less than 10%. The uncertainty related to the subtraction of the non-W Z contribution, estimated by imposing cross-section uncertainties for all the relevant processes, is found to be about 2%. The total uncertainty in the W Z estimate is about 6%. The kinematic distributions are estimated from simulation, with both the experimental and theoretical uncertainties considered.
To estimate the non-resonant-background, including W W , top-quark (tt and W t), and Z → τ τ production, a control region dominated by the non-resonant-processes (with a purity above 95%) is defined with all the event selection criteria in table 2, except that the final state is required to contain an opposite-sign eµ pair. The non-resonant-contribution in the ee (µµ) channel is calculated as one half of the observed data yield after subtracting the contribution from the other background processes in the control region, and then corrected for the difference in the lepton reconstruction and identification efficiencies between selecting an eµ pair and an ee (µµ) pair. The lepton efficiency correction is derived as the square root of the ratio of the numbers of µµ and ee events in data after the preselection. The choice of deriving the correction after preselection minimises the resulting statistical uncertainty. The total uncertainty in the non-resonant-estimate is about 16%, including the statistical uncertainty of the data in the control region (14%) and the method bias estimated from simulation (7%). The kinematic distributions for the nonresonant-background estimate in the signal region are predicted with simulation, and the assigned systematic uncertainty covers the experimental uncertainty in the simulated shape as well as the difference between data and simulation in the control region. Figure 4 gives two examples of comparing data and predictions in the W Z and nonresonant-background control regions. The left-hand figure is the m W T distribution in the W Z control region, where the normalisation factor f W Z is applied to the W Z simulation and good agreement between the observed and predicted shapes is found. The righthand figure is the E miss T distribution in the non-resonant-control region, where the W W and top-quark (tt and W t) production processes are dominant. Both the statistical and experimental uncertainties are included in the hashed bands in these figures.
A data-driven method is used to estimate the Z + jets background. This method defines three independent regions (labelled as B, C and D) which are enriched in Z + Regions B and D are dominated by the Z + jets process (with a purity greater than 95%), while its relative contribution in region C is only 70% because the tt contribution in this phase space region remains large. The derived Z + jets contribution is corrected for the closure factor (n A /n C × n D /n B ) estimated from simulation. This factor is found to be -14 -JHEP10(2019)127 0.9 and has a relative uncertainty of 48%, consisting of the statistical (40%), experimental (22%), and methodology uncertainties (15%). The experimental uncertainty in the closure factor is dominated by jet energy scale and resolution. The methodology uncertainty covers the variations obtained by changing the E miss T and V T /S T thresholds in the low-E miss T removal requirement by 40%. The Z +jets estimation is also subject to the statistical uncertainty of the data (5%) and the subtraction of non-Z contributions in the control regions (5%). The non-Z subtraction uncertainty is driven by the modelling uncertainty for the main non-Z process in region C (tt production), which is about 10-20% for E miss T above 100 GeV [61]. The total uncertainty on the Z + jets estimate is about 50%. The kinematic distributions for the Z + jets background in the signal region are derived from the data in region C, together with a systematic uncertainty assigned in a way similar to that described above for the non-resonant-background.
The ZZ → 4 , V V V and ttV (V ) backgrounds are estimated from simulation, and their contributions have a total uncertainty of 10-20%, including both the theoretical crosssection [8, 62, 63] and the experimental uncertainties.

Systematic uncertainties
The measurement results and predictions are subject to theoretical and experimental uncertainties, as well as uncertainties related to the background estimation. The background uncertainties are explained in section 6. The statistical uncertainties of the simulated samples for both the signal and background processes are also taken into account wherever applicable. The systematic uncertainty sources for the signal process are detailed below.
The theoretical uncertainties for the dominant qqZZ signal sub-process are estimated with the Powheg+Pythia8 generator, since only the total cross-section has been calculated to NNLO QCD and NLO EW accuracies. The theoretical uncertainties originate from the PDF choice, the missing higher-order QCD calculation, and the UEPS modelling. The PDF uncertainty is calculated as the 68% confidence-level eigenvector uncertainty [26] of the nominal PDF used in the simulation. The uncertainty due to the QCD calculation, also referred to as the "scale" uncertainty, covers the variations of predictions from changing the QCD renormalisation and factorisation scales. The QCD scales are varied independently by factors ranging from one half to two, which in total yields seven different scale choices including the nominal one. The UEPS uncertainty is taken as the difference in the predictions between the Herwig++ and the default showering programs. The fractional theoretical uncertainty in A ZZ for the qqZZ process is about 1.8%, while the overall uncertainties in the cross-section predictions in the total and fiducial phase spaces are about 3% and 5%, respectively. The Sherpa generator is used to cross-check the nominal predictions, and the A ZZ factors from Powheg and Sherpa are consistent with each other within the uncertainty. The C ZZ predictions from the two generators are found to be consistent within the statistical uncertainty of 1%, and in this measurement, the theoretical uncertainty in C ZZ is neglected.
The understanding of the p T spectrum in the fiducial phase space is crucial for the study of aTGCs, and the predictions from the two generators differ by up to 10% for -15 -

JHEP10(2019)127
p T around 300 GeV, which is slightly above the theoretical uncertainty of the Powheg prediction. The Powheg prediction with the K-factors applied has better precision in terms of the EW calculation, while the Sherpa generator is expected to give a better description of ZZ production with extra QCD radiation. Finally, an uncertainty is applied to the p T prediction, as the sum in quadrature of the theoretical uncertainty estimated with Powheg and the difference between Powheg and Sherpa, which is about 5% for p T around 150 GeV and increases to about 11% for p T above 250 GeV.
A constant 30% uncertainty is assigned to the total ggZZ cross-section prediction, which covers the uncertainties concerning the NLO K-factor [34] and the potential missing higher-order contributions [64]. The A ZZ predictions for the ggZZ process from the gg2vv and Sherpa generators are found to be consistent, and the A ZZ uncertainty is estimated with Sherpa and found to be 4.6%. The theoretical uncertainty in C ZZ is neglected for the ggZZ process.
The major experimental uncertainties originate from the luminosity uncertainty, the momentum scale and resolution of leptons and jets, and the lepton reconstruction and selection efficiencies [47, 48, 51, 65]. Smaller experimental uncertainties are also considered, which include uncertainties due to the trigger selection efficiency, the b-jet identification efficiency, the calculation of the E miss T soft-term, and the variation of the average number of interactions per bunch crossing (hereafter referred to as pile-up uncertainty). The pileup uncertainty covers the uncertainty on the ratio between the predicted and measured inelastic cross-section in the fiducial volume defined by M X > 13 GeV where M X is the mass of the hadronic system [66]. Overall, the total experimental uncertainty on C ZZ is 3.1%, dominated by the jet and lepton components. The uncertainty in the combined 2015+2016 integrated luminosity is 2.1%. It is derived, following a methodology similar to that detailed in ref.
[67], and using the LUCID-2 detector for the baseline luminosity measurements [68], from calibration of the luminosity scale using x-y beam-separation scans.
The fractional uncertainties in A ZZ and C ZZ are summarised in table 5. In this analysis, the theoretical uncertainties are treated as uncorrelated between the qqZZ and ggZZ processes, while the experimental uncertainties are considered as fully correlated across the relevant processes and final-state channels. Table 6 lists separately for the ee and µµ channels the observed data yields and the expectations for the signal and background contributions after the final selection. Figure 5 shows for the combined ee and µµ channels the observed and expected E miss T distributions, which are in good agreement.

Integrated cross-section results
The integrated fiducial and total cross-sections (σ fid ZZ→ νν and σ tot ZZ ) are determined by binned maximum-likelihood fits to the E miss T distributions. As shown in figure 5, the signal-to-background ratio increases as E miss T becomes larger. The use of E miss T improves the precision of the measured fiducial cross-section relatively by 5% compared with the case where no kinematic information is used.  Total expected 347 ± 5 ± 15 372 ± 6 ± 16 Table 6. Observed data yields and expected signal and background contributions, shown separately for the ee and µµ channels. The errors shown for the expected yields correspond to the statistical and systematic contributions in that order. The expected background and signal+background yields are shown in the last two rows, where the uncertainties are computed as the sum in quadrature of those from the individual processes.
The expected yield in each channel i and in each E miss T bin j is given by: where L is the integrated luminosity, N ij bkg the expected background yield, B the branching fraction for the ZZ → νν decay ( = e or µ), and f ij ZZ is the fraction of signal events   in bin j with respect to the total distribution. The number of events follows a Poisson distribution in each bin, and the systematic uncertainties are treated as Gaussian nuisance parameters, θ k , in the fit. For each source of systematic uncertainty, k, a single nuisance parameter is used for all the processes and channels where this uncertainty matters. The statistical uncertainty due to the limited size of simulated samples is treated as uncorrelated among bins and channels. The binned likelihood function is built over all bins as follows: where N ij obs is the observed data yield in each bin. Table 7 summarises the main sources of uncertainty in the measured combined fiducial cross-section, where individual sources of a similar nature are grouped together. The statistical and total systematic uncertainties in the measurement are of similar sizes. Table 8 shows the measured fiducial cross-sections, separately for each channel and for their combination, together with the breakdown of their uncertainties. The ee and µµ channel cross-sections are compatible within their respective statistical uncertainties. The measured combined fiducial cross-section has a total uncertainty of 7%, which is significantly better than the previous measurement [7], and comparable in size to that obtained in the ZZ → 4 channel [8,12]. Table 8 also shows the combined measured total cross-section, as well as the predictions for the cross-sections, as taken from  Table 7. Relative contributions to the measured combined fiducial cross-section from the main sources of uncertainty after the fit procedure. The total uncertainty includes the data statistical and systematic components. For the systematic uncertainty, the individual sources of a similar nature are grouped together for simplicity. "Sim. stat." indicates the uncertainty source corresponding to the limited size of the simulation samples for the signal and background processes.  Table 8. Measured and predicted integrated cross-sections in the fiducial and total phase spaces, together with the breakdown of their uncertainties. The luminosity uncertainty is quoted separately from the other systematic uncertainties. The measurements are also shown separately for the ee and µµ channels in the case of the fiducial cross-section. measurement is about 13% higher than the prediction, which is not significant given the size of the measurement and prediction uncertainties.

Differential cross-section results
Differential cross-sections are measured in the fiducial phase space by counting data events observed in each bin of the observables of interest, after subtracting the expected background contribution, and correcting for the detector effects with the unfolding procedure, chosen here to be the iterative Bayesian unfolding method of ref. [69]. The unfolding process takes into account fiducial corrections (correcting for events outside the fiducial phase space but passing the detector-level selections), bin-to-bin migrations due to detector resolution, and detector inefficiencies. An optimal number of two iterations is used for this analysis, as a balance between the size of the statistical uncertainty in the measurement and residual biases from the method. The residual bias is in almost all bins below 1%, as estimated by comparing the results obtained using different prior distributions (constant, expected, observed) in the unfolding process. The experimental uncertainties for the measurement results are evaluated by varying the response matrices for the unfolding according to the ±1σ effects of each uncertainty source, and by comparing the resulting unfolded results with the nominal one. The background uncertainties are considered at the stage of the background subtraction. The statistical uncertainty of the data is estimated by repeating the unfolding procedure with 2000  sets of pseudo-data and then taking the root mean square of the deviations of the resulting spectra from the data spectrum. The response matrices are also subject to the statistical uncertainty of the simulated samples, which is estimated using a similar approach. Figures 6 and 7 present the measured differential cross-sections for the eight observables of interest defined in section 1. The binning for each variable is chosen to minimise the bin-to-bin migrations while preserving a sufficient number of events per bin. The N jets spectrum is measured only for hard jets as defined in section 4, and the p jet1 T distribution is obtained in the fiducial phase space of events containing at least one hard jet. The predictions from Powheg+gg2vv and Sherpa are also shown in figures 6 and 7, and are found to be in agreement with the measurements within uncertainties. The electroweak production of ZZ associated with two jets is not taken into account in the predictions due to its negligible contribution. The differential measurements are largely dominated by the statistical uncertainty on the data, but the systematic uncertainties contribute significantly in certain regions of phase space. As an example, the uncertainties from the various sources for the differential measurement of the p T distribution are listed in table 9.
Comparing with the Run-1 results of ref.
[7], this measurement is obtained from a larger dataset with highly improved accuracy and for a wider range of observables.

Search for aTGCs
The search for aTGCs is carried out using the unfolded p T distribution of figure 6 in the fiducial phase space. The contribution due to aTGCs is introduced using an effective vertex function approach [1]. It includes two coupling parameters 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 . The sensitivity range of this search is found to be within the unitarity bounds, so no form factors [1] are applied to the coupling parameters in this analysis. Furthermore, the coupling parameters are used to extract information about the dimension-eight operators of the effective field theory [70].
The relative change in the SM production cross-section in each p T bin is parameterised as a quadratic function depending on the coupling parameters. The parameterisation for the qqZZ process is derived from simulation as described in section 3, while the impact of the aTGCs is neglected for the ggZZ process. Since the sensitivity to possible aTGCs lies in the high-p T region, only the bins with p T > 150 GeV are considered in the search. Figure 8 compares the measured p T spectrum in this region of phase space to the SM prediction alone and to the SM prediction augmented with aTGCs corresponding to different values of the coupling parameters described above. The highest-p T bins are required to contain a minimum of ten observed events at the detector level, which arises from the fact that the limits on aTGCs derived from the unfolded distribution rely on a Gaussian approximation for the statistical fluctuations from the data. The total uncertainty in the measured cross-section in the last bin is about 40%, while the corresponding SM prediction has an uncertainty of 15%. The sensitivity to aTGCs is thus still limited by the statistical uncertainty of the data.
Since the observation is consistent with the SM prediction, 95% confidence intervals (CIs) for the coupling parameters are derived. The test statistic is based on a profile loglikelihood ratio [71], and the likelihood function is described by a multivariate Gaussian distribution with the following form: where x represents the measured p T spectrum, µ the expected one, and C is the expected covariance matrix in the differential measurement, with the term |C| standing for the determinant of the covariance matrix. Both x and µ have five elements, corresponding to the five p T bins used in this search. The expectation µ depends on the coupling parameters and a set of Gaussian nuisance parameters, θ s , which encapsulate the theoretical uncertainties in the prediction. The expected covariance matrix C is first derived by unfolding the predicted SM spectrum at detector level and then modified to account for the fact that the prediction µ, varied in the statistical test, differs from the SM prediction µ 0 in the fiducial region in the case of non-zero coupling parameters. The relative difference between µ and µ 0 can be expressed using a vector k, with k i = µ i /µ i 0 for bin-i. Each matrix element C ij is the sum of a statistical uncertainty component (C ij stat ) and systematic uncertainty components arising from SM ZZ process (C ij ZZ-syst ) and background contribution (C ij bkg-syst ); the dependence of C on µ is implemented such that C ij stat and C ij ZZ−syst scale with k i × k j and k i × k j , respectively.
The 95% CIs are derived by scanning the parameter space and using a frequentist method with the CL s formalism [72]. In the first step, the one-dimensional CI for each aTGC parameter is derived with all other parameters set to zero, and the results are shown in table 10. These intervals are found to be more stringent than those obtained  Table 10. One-dimensional 95% confidence intervals for the aTGC parameters described in the text. The limits on each parameter are derived with the other parameters set to zero. from the ZZ → 4 channel based on the same ATLAS data sample [8]. In the second step, two-dimensional CIs for each pair of coupling parameters are derived with the other two parameters set to zero. Figure 9 presents the two-dimensional CI contours for the six possible pairs of aTGC coupling parameters. Finally, in the context of effective field theories, a one-dimensional 95% CI can be placed on the CP-even dimension-eight operator CB W /Λ 4 [70]. This is obtained by extrapolating the results of table 10 following the conversion formula in ref. [70]. A 95% CI of [−4.0, 4.0] in units of TeV −4 is set for CB W /Λ 4 , which is improved by 30% with respect to the ATLAS results in the ZZ → 4 channel [8].

Conclusion
This paper presents a measurement of ZZ production using the νν final state, based on 36.1 fb −1 of data collected by the ATLAS detector at the LHC in pp collisions at √ s = 13 TeV. The fiducial cross-section of the combined ee and µµ channels is measured to be σ fid ZZ→ νν = 25.4 ± 1.4 (stat) ± 0.9 (syst) ± 0.5 (lumi) fb, in agreement with the SM prediction of 22.4 ± 1.3 fb. The integrated cross-sections in the total and fiducial phase spaces are measured with a total uncertainty of 7%, which is significantly better than the previous measurement using the 8 TeV data. The measured cross-sections are slightly larger than the SM predictions, but the difference is not significant considering the measurement and prediction uncertainties. Differential cross-sections are reported for eight different kinematic variables in the fiducial phase space, and no significant deviation from the expectations is found. The measured p T spectrum in the fiducial phase space is used to set limits on aTGCs, and the obtained 95% confidence intervals for the aTGC parameters are more stringent than those derived from the 4 channel.   Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.    [57] Particle Data Group, The review of particle physics, Phys. Rev. D 98 (2018) 010001.