Measurement of differential cross sections for Z boson pair production in association with jets at $\sqrt{s} =$ 8 and 13 TeV

This Letter reports measurements of differential cross sections for the production of two Z bosons in association with jets in proton-proton collisions at $\sqrt{s} =$ 8 and 13 TeV. The analysis is based on data samples collected at the LHC with the CMS detector, corresponding to integrated luminosities of 19.7 and 35.9 fb$^{-1}$ at 8 and 13 TeV, respectively. The measurements are performed in the leptonic decay modes ZZ $\to\ell^+ \ell^- \ell'^+ \ell'^-$, where $\ell,\ell' =$ e, $\mu$. The differential cross sections as a function of the jet multiplicity, the transverse momentum $p_\mathrm{T}$, and pseudorapidity of the $p_\mathrm{T}$-leading and subleading jets are presented. In addition, the differential cross sections as a function of variables sensitive to the vector boson scattering, such as the invariant mass of the two $p_\mathrm{T}$-leading jets and their pseudorapidity separation, are reported. The results are compared to theoretical predictions and found in good agreement within the theoretical and experimental uncertainties.


Introduction
The production of massive vector boson pairs is a key process for the understanding of both the non-Abelian gauge structure of the standard model (SM) and of the electroweak symmetry breaking mechanism. Thus, relevant information can be gathered measuring vector boson scattering [1] and triboson production processes that occur through the electroweak (EW) production of jets in association with bosons. Because of the very low cross sections for these processes compared to others leading to the same final state, a detailed understanding of the quantum chromodynamics (QCD) corrections to the associated production of vector boson pairs and jets is of paramount importance. The analysis presented in this Letter has been designed to provide such detailed understanding.
Both the ATLAS and CMS Collaborations have measured the inclusive production cross section of Z boson pairs and the differential cross sections as a function of Z boson pair observables [2][3][4][5][6][7][8]. In this Letter we present new measurements of differential cross sections for the production of two Z bosons in association with jets in proton-proton (pp) collisions at √ s = 8 and 13 TeV that extend the analyses of Refs. [6,8] to jet variables. The most recent publication from the AT-LAS Collaboration [4] includes jet variables as well. The decay modes of the Z boson to electron and muon ( = e, µ) pairs have been exploited. Reconstructed distributions are corrected for event selection efficiency and detector resolution effects by means of an iterative unfolding technique, which makes use of a response matrix to map physics variables at generator level onto their reconstructed values.
This Letter presents the dependence of the cross section on the jet multiplicity and the kinematic properties of the two p T -leading jets (where p T is the transverse momentum). Comparison with theoretical predictions provides an important test of the QCD corrections to ZZ production. Normalized differential cross sections as a function of the p T and pseudorapidity η of the two p T -leading jets, as well as their invariant mass (m jj ) and pseudorapidity separation (∆η jj ), are presented. The study of m jj establishes the basis for future multiboson final-state searches and for the investigation of phenomena involving interactions with four bosons at a single vertex, while the measurement of the ∆η jj distribution is instrumental in the study of vector boson scattering with the same final state as that of this analysis. The results presented in this paper form the ground on which the vector boson scattering analysis [9] has been built, and pave the road towards its extension to the full 13 TeV data set. All measurements are compared to predictions from recent Monte Carlo (MC) event generators. The data sets correspond to integrated luminosities of 19.7 and 35.9 fb −1 , collected by the CMS Collaboration at 8 and 13 TeV, respectively.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are silicon pixel and strip tracking detectors, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the η coverage provided by the barrel and endcap detectors up to |η| = 5. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid, using three different technologies: drift tubes for |η| < 1.2, cathode strip chambers for 0.9 < |η| < 2.4, and resistive plate chambers for |η| < 1.6. The silicon tracker measures charged particles within the range |η| < 2.5. For nonisolated particles in the range 1 < p T < 10 GeV and |η| < 1.4, the track resolutions are typically 1.5% in p T and 25-90  µm in the transverse (longitudinal) impact parameter [10].
The first level of the CMS trigger system [11], composed of custom hardware processors, uses information from the calorimeters and muon detectors to select the most interesting events within a time interval of less than 4 µs. The high-level trigger processor farm further decreases the event rate from around 100 kHz to less than 1 kHz, before data storage.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [12].

Signal and background simulation
Several MC event generators are used to simulate the signal and background contributions. The MC simulation samples are employed to optimize the event selection, evaluate the signal efficiency and acceptance, estimate part of the background, and extract the unfolding response matrices used to correct for detector effects in the measured distributions.
For the 8 TeV data analysis, MADGRAPH5 1.3.3 [13,14] is used to simulate the production of the four-lepton final state at leading order (LO) in QCD with up to 2 jets included in the matrix element calculations. POWHEG 2.0 [15][16][17][18] is used for the simulation of the same process at next-to-leading-order (NLO). MADGRAPH5 aMC@NLO 2.3.3 (abbreviated as MG5 aMC@NLO in the following) [14,19], which simulates signal processes at NLO with 0 and 1 jet included in the matrix element calculations, is also used for comparison purposes at generator level. For the 13 TeV data analysis, the four-lepton processes are simulated at NLO in QCD with 0 or 1 jet included in the matrix element calculations with MG5 aMC@NLO and, for comparison, with POWHEG 2.0 at NLO. The latter is scaled by a factor of 1.1 to reproduce the total ZZ production cross section calculated at next-to-next-to-leading order (NNLO) [20] at 13 TeV. MG5 aMC@NLO, for both 8 and 13 TeV analyses, is used to evaluate the expected cross sections in bins of jet multiplicity at generator level. MG5 aMC@NLO and POWHEG 2.0, for both the 8 and 13 TeV analyses, include ZZ, Zγ * , Z, and γ * γ * processes, with the generator level constraint m + − > 4 GeV applied to all pairs of oppositely charged same-flavor leptons, to avoid infrared divergences.
The gg → ZZ processes, which occur via loop-induced diagrams, are generated at LO with MCFM 6.7 (7.0) [21] for the 8 (13) TeV analysis. The 13 TeV samples are scaled by a factor of 1.7 to match the cross section computed at NLO [22]. Electroweak production of four leptons and two jets is simulated at LO with PHANTOM [23]. This sample includes triboson processes, where the Z boson pair is accompanied by a third vector boson that decays into jets, as well as diagrams with quartic vertices.
Other diboson and triboson processes (WZ, Zγ, WWZ) as well as ttZ, tt, and Z+jets samples are generated at LO with MADGRAPH5 for the 8 TeV analysis, and at NLO with MG5 aMC@NLO, for the 13 TeV analysis.

Event selection and particles reconstruction
The primary triggers for this analysis require the presence of two loosely isolated leptons of the same or of different flavor. The minimum p T for the first lepton is 17 GeV, while it is 8 (12) GeV for the second lepton in the 8 (13) TeV analysis. Triggers requiring a triplet of low-p T leptons with no isolation requirement and for, the 13 TeV analysis, isolated single-electron and singlemuon triggers, with minimal p T -thresholds of 27 and 22 GeV, respectively, help to increase the efficiency. The overall trigger efficiency for events that pass the ZZ selection is greater than 98%.
The offline event selection procedure is similar to that of the inclusive ZZ analyses [6-8] and is based on a global event description [32] that classifies particles into mutually exclusive categories: charged hadrons, neutral hadrons, photons, muons, and electrons. Events are required to have at least one vertex [10] within 24 cm of the geometric center of the detector along the beam direction, and within 2 cm in the transverse plane. Because of pileup the selected event can have several reconstructed vertices.
For the analysis at 8 TeV the vertex with the largest sum of the p 2 T of the tracks associated to it is chosen as the primary pp interaction vertex, while at 13 TeV the reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary vertex. The physics objects are the jets, clustered using the jet finding algorithm [33,34] with the tracks assigned to the vertex as inputs (including lepton tracks), and the associated missing p T , taken as the negative vector sum of the p T of those jets. Events with leptons are selected by requiring each lepton track to have a transverse impact parameter, with respect to the primary vertex, smaller than 0.5 cm and a longitudinal impact parameter smaller than 1.0 cm.
Electrons are measured in the range |η| < 2.5 by using both the tracking system and the ECAL. They are identified by means of a multivariate discriminant that includes observables sensitive to bremsstrahlung along the electron trajectory, the geometrical and momentum-energy agreement between the electron track and the associated energy cluster in the ECAL, the shape of the electromagnetic shower, and variables that discriminate against electrons originating from photon conversions [35]. The momentum resolution for electrons with p T ≈ 45 GeV from Z → e + e − decays ranges from 1.7% for nonshowering electrons in the barrel region to 4.5% for showering electrons in the endcaps [35].
Muons are reconstructed in the range |η| < 2.4 by combining information from the silicon tracker and the muon system [36]. The matching between the inner and outer tracks proceeds either outside-in, starting from a track in the muon system, or inside-out, starting from a track in the silicon tracker. The muons are selected among the reconstructed muon track candidates by applying minimal requirements on the track in both the muon system and the inner tracker system, and taking into account the compatibility with minimum-ionizing particle energy deposits in the calorimeters. In the intermediate range of 20 < p T < 100 GeV, matching muons to tracks measured in the silicon tracker results in a relative p T resolution of 1.3-2.0% in the barrel, and better than 6% in the endcaps. The p T resolution in the barrel is better than 10% for muons with p T up to 1 TeV [36].
Electrons (muons) are considered candidates for inclusion in the four-lepton final states if they have p T > 7 (5) GeV and |η | < 2.5 (2.4). In order to suppress electrons from photon conversions and muons originating from in-flight decays of hadrons, we place a requirement on the impact parameter computed in three dimensions. We require that the ratio of the impact parameter for the track and its uncertainty to be less than 4. To discriminate between prompt leptons from Z boson decay and those arising from electroweak decays of hadrons within jets, an isolation requirement for leptons is imposed. The relative isolation is defined as where the sums run over the charged and neutral hadrons, and photons, in a cone defined by ∆R ≡ √ (∆η) 2 + (∆φ) 2 around the lepton trajectory. The radius ∆R is set to be 0.4 and 0.3 in the 8 and 13 TeV data analyses, respectively. To minimize the contribution of charged particles from pileup to the isolation calculation, charged hadrons are included only if they originate from the primary vertex. The contributions of neutral particles from pileup to the activity inside the cone around a lepton is referred to as p PU T , and is obtained with different methods for electrons and muons. For electrons, p PU T is evaluated with the jet area method described in Ref. [37]. For muons, it is taken to be half the sum of the p T of all charged particles in the cone originating from pileup vertices. The factor of one-half accounts for the expected fraction of neutral to charged particles in hadronic interactions. A lepton is considered isolated if R iso < 0.4 (0.35) in the 8 (13) TeV data analysis.
The lepton momentum scales are calibrated in bins of p T and η using the decay products of known resonances decaying to lepton pairs. The measured lepton momentum scale is corrected with a Z → + − sample, by matching the peak of the reconstructed dilepton mass spectrum to the nominal value of m Z [38]. Muon momenta are calibrated by using J/ψ decays as well. We account for final-state radiation of leptons by correcting their momenta with photons of p T > 2 GeV and within a cone of ∆R = 0.5 around the lepton momentum direction [39,40]. The photons selected by this algorithm are excluded from the lepton isolation computation. The efficiency of the lepton reconstruction and selection is measured with the tag-and-probe technique [41] in bins of p T and η . This measurement is used to correct the simulation efficiency. Jets are reconstructed from particle candidates by means of the anti-k T clustering algorithm [33], as implemented in the FASTJET package [34], with a distance parameter of 0.5 (0.4) in the 8 (13) TeV data analysis. The jet energy resolution amounts typically to 15% at 10 GeV, 8% at 100 GeV, and 4% at 1 TeV.
Jet energy corrections are extracted from the data and the simulated events by combining several measurements and methods that account for the effects of pileup, non-uniform detector response, and residual data-simulation jet energy scale (JES) differences. The JES calibration [42, 43] relies on corrections parametrized in terms of the uncorrected p T and η of the jet, and are applied as multiplicative factors to the four-momentum vector of each jet.
In order to maximize the reconstruction efficiency while reducing the instrumental background and contamination from pileup jets, loose identification quality criteria [44] are imposed on jets, based on the energy fraction carried by charged and neutral hadrons, as well as charged leptons and photons. A minimum threshold of 30 GeV on the p T of jets is required to ensure that they are well measured and to reduce the pileup contamination. Jets are required to have |η| < 4.7 and to be separated from all selected lepton candidates by at least ∆R = 0.5 (0.4) in the 8 (13) TeV analysis.
A signal event must contain at least two Z/γ * candidates, each reconstructed from a pair of isolated electrons or muons of opposite charges. The highest-p T lepton must have p T > 20 GeV, and the second-highest lepton p e T > 10 (12) GeV if it is an electron, or p µ T > 10 GeV in case of a muon for the analysis at √ s = 8 (13) TeV. All leptons are required to be separated by ∆R ( , ) > 0.02, and electrons are required to be separated from muons by ∆R (e, µ) > 0.05.
Within each event, all permutations of oppositely charged leptons giving a valid pair of Z/γ * candidates are considered separately. For each 4 candidate, the lepton pair with the invariant mass closest to the nominal Z boson mass is denoted by Z 1 and the other dilepton candidate is denoted by Z 2 . Both Z 1 and Z 2 are required to have a mass between 60 and 120 GeV. All pairs of oppositely charged leptons in the 4 candidate are required to have m > 4 GeV regardless of their flavor to remove contributions from the decay of low-mass hadron resonances.
If multiple 4 candidates within an event pass this selection, the candidate with m Z 1 closest to the nominal Z boson mass is chosen. In the rare cases (0.3%) of further ambiguity, which may arise in events with more than 4 leptons, the Z 2 candidate that maximizes the scalar p T sum of the four leptons is chosen. The set of selection criteria just described is referred to as the ZZ selection, and gives a total of 288 (927) observed events at √ s = 8 (13) TeV. The corresponding number of expected signal events from MC prediction is about 271 (850).

Background estimation
The largest source of background arises from processes in which heavy-flavor jets produce secondary leptons, and from processes in which jets are misidentified as leptons. The main contributing processes are Z+jets, tt, and WZ+jets.
However, the lepton identification and isolation requirements reduce this background to a very small level compared to the signal. The residual contribution is estimated from data samples consisting of Z + events that are required to pass the ZZ selection described in Section 4, except that either one or both leptons belonging to the Z 2 candidate fail the isolation or identification requirements. Two control samples are selected, with one and two misidentified leptons, respectively. The background yield in the signal region is estimated by weighting the number of events in the control samples by the lepton misidentification rate measured in data in a dedicated control region. The procedure is identical to that of Refs. [7, 8] and is described in more detail in Ref. [39].
Another source of background arises from processes that produce four genuine high-p T isolated leptons, pp → ttZ and pp → WWZ. This contribution is small and is estimated by using the corresponding simulated samples.

Systematic uncertainties
The systematic uncertainties are estimated by varying the quantities that may affect the cross section and by propagating the changes to the analysis procedure. The systematic uncertainties from sources that may affect the differential cross section shapes have been estimated through the unfolding procedure by recomputing the response matrix, after varying each source of systematic uncertainty independently and in both directions, up and down. The systematic uncertainties are summarized in Table 1; uncertainties that are a function of the jet multiplicity are listed as a range.
The systematic uncertainty in the trigger efficiency is evaluated by taking the difference between the value obtained from the data and that from the simulated events, and it amounts to 1.5 (2.0)% of the differential cross sections measured with the 8 (13) TeV data. The uncertainties arising from lepton reconstruction and selection (identification, isolation, and impact parameter determination) depend on the jet multiplicity, are sensitive to statistical fluctuations, and range between 0.9 and 4.4%, in the 8 TeV analysis (3.7 and 4.5%, in the 13 TeV analysis). The JES uncertainty is the largest contribution affecting the differential cross section measurements. It increases with the jet multiplicity and reaches 9.2 (17.6)% when the number of jets exceeds two in the 8 (13) TeV analysis. Likewise, the uncertainty due to the jet energy resolution (JER) increases from 0.2 to 1.7% (2.1 to 8.4%) for the 8 (13) TeV samples. The larger JES and JER uncertainties for the 13 TeV sample reflect the increase in the number of soft jets (with p T close to the 30 GeV threshold) as a function of the center-of-mass energy.
The uncertainties in the Z+jets, WZ+jets, and tt background yields reflect their different relative fraction in the control samples used for this background estimation, and the statistical uncertainty of the control samples. The effect of these uncertainties increases with the jet multiplicity and amounts to 0.7-6.9% (0.5-2.5%) in the 8 (13) TeV measurement. The contribution to the uncertainty from the modeling of genuine four lepton background is smaller and varies between 0.1 and 2.0% (<0.1 and 1.2%) for the 8 (13) TeV data. The pileup uncertainty is evaluated by varying the pileup modeling in the MC samples within its uncertainty. The uncertainty in the integrated luminosity is 2.6 [45] and 2.5% [46] for the 8 and 13 TeV data, respectively.
The contribution of the unfolding procedure to the systematic uncertainty is obtained by comparing the results found with two different sets of MC samples: MADGRAPH5 + MCFM + PHANTOM (MG5 aMC@NLO + MCFM + PHANTOM) and POWHEG + MCFM + PHANTOM for the 8 (13) TeV measurement, and ranges from 0.2 to 3.7% (0.5 to 5.1%) at 8 (13) TeV. The impact of the relative contribution of the qq → ZZ and gg → ZZ processes in the response matrix definition is less than 1% and is evaluated by varying the corresponding cross section within their renormalization and factorization scale uncertainties. For 8 TeV, where no LO to NLO factor is applied to the MCFM cross section, the gg → ZZ cross section is varied by 100% of its value. The statistical uncertainties of the MC samples result in negligible contributions to the response matrix uncertainty. The systematic uncertainty arising from the choice of the PDF and the strong coupling strength α S has been evaluated using the PDF4LHC recommendations [47][48][49], using the CT10, MSTW08, and NNPDF2.3 [50] PDF sets, in the 8 TeV analysis, and the NNPDF3.0 set in the 13 TeV analysis.
The total systematic uncertainty is obtained by summing all the sources in quadrature, taking into account the correlations among the different channels.
For the normalized differential cross sections, only systematic uncertainties affecting the shape of the distributions are relevant. The uncertainties in the luminosity and trigger efficiency cancel out completely, as well as other contributions to the uncertainty in the total yield.

The ZZ+jets differential cross section measurements
The distributions of the jet multiplicity combining the 4µ, 4e, and 2µ2e channels are shown in Fig. 1, together with the SM expectations, the estimated backgrounds, and the systematic uncertainty in the prediction.  The differential pp → ZZ → cross section is measured as a function of the jet multiplicity, the p T -leading jet transverse momentum (p j1 T ) and pseudorapidity (η j1 ) with the 8 and 13 TeV data. Because of the limited number of events with more than one jet at 8 TeV, the differential cross section as a function of the p T -subleading jet transverse momentum (p j2 T ) and pseudorapidity (η j2 ), as well as the invariant mass of the two p T -leading jets (m jj ) and their pseudorapidity separation (∆η jj ) are studied at 13 TeV only. For all measurements we consider jets with p j T > 30 GeV and |η j | < 4.7. For the jet multiplicity distribution we also present the measurements made with central jets (|η j | < 2.4) only. The measurements are performed for the two slightly different phase space regions adopted for the 8 [6] and 13 [8] TeV data, which are given in Table 2. The generator-level lepton momenta are corrected by adding the momenta of generator-level photons within ∆R ( , γ) < 0.1. The Z bosons are then selected with the same method adopted to extract the signal at the reconstruction level. TeV p e T > 7 GeV, |η e | < 2.5 p e T > 5 GeV, |η e | < 2.5 p µ T > 5 GeV, |η µ | < 2.4 p µ T > 5 GeV, |η µ | < 2.5 common definitions p 1 T > 20 GeV , p 2 T > 10 GeV m + − > 4 GeV (any opposite-sign same-flavor pair) 60 < (m Z 1 , m Z 2 ) < 120 GeV Each distribution is corrected for the event selection efficiency and the detector resolution effects by means of a response matrix that translates the physics variables at generator level into their reconstructed values. The correction procedure is based on the iterative D'Agostini unfolding method technique [51], as implemented in the ROOUNFOLD toolkit [52], and regularized by stopping after four iterations. The robustness of the result is tested against the singular value decomposition (SVD) [53] alternative unfolding method. For each measured distribution, a response matrix is evaluated using two different sets of generators: the first one includes MADGRAPH5 (qq → ZZ), MCFM (gg → ZZ) and PHANTOM (qq → ZZ + 2 jets) for the 8 TeV data set and MG5 aMC@NLO (qq → ZZ), MCFM (gg → ZZ) and PHANTOM (qq → ZZ + 2 jets) for the 13 TeV data set. In the second one, the POWHEG sample is instead used for the qq → ZZ process in both the 8 and 13 TeV data analyses. The former set is taken as the reference, while the latter is used for comparison and to estimate the systematic uncertainty due to the MC generator choice. After the unfolding, the cross sections for pp → ZZ + N jets → + N jets, for N = 0, 1, 2, and ≥3, are extracted.
The differential cross sections as a function of the jet multiplicity are shown in Fig. 2 for |η j | < 4.7 (upper) and for |η j | < 2.4 (lower). The ratios between the measured and expected distributions from the MADGRAPH5, MG5 aMC@NLO, and POWHEG set of samples for √ s = 8 TeV, and POWHEG and MG5 aMC@NLO for √ s = 13 TeV are also shown in the figures. Uncertainties in the MC predictions at the matrix element level are evaluated by varying the renormalization and factorization scales independently, up and down, by a factor of two with respect to the default values of µ R = µ F = m 4 for POWHEG and µ R = µ F = 1 2 ∑ p j T + ∑ p T for MG5 aMC@NLO. In the MCFM predictions, the uncertainty in the LO to NLO cross section scaling factor includes the renormalization and factorization scales uncertainty. The theoretical uncertainties also include the uncertainties in the PDF and α S . The measured and expected cross section values for |η j | < 4.7 are given in Tables 3 and 4.   the shape is included, which yields a smaller uncertainty compared to the unnormalized case. Figure 3 (top panels) shows the normalized differential cross section as a function of the jet multiplicity, with |η j | < 4.7. The observed fraction of events in the first bin with zero jets is larger than the predicted value, while for 1, 2, and ≥ 3 jets, the fraction is lower. Better agreement is observed for |η j | < 2.4 (Fig. 3, bottom panels). The measurements of the differential cross section as a function of the jet multiplicity are fairly well reproduced by the predictions both at 8 and 13 TeV when NLO matrix-element calculations are used in conjunction with PYTHIA 8 for parton showering, hadronization, and underlying event simulation. In the data, jets tend to have a lower p T value than in the simulations and therefore, on average, they are less likely to pass the 30 GeV threshold, thus increasing the number of events with no jets. The observation of fewer events than expected with at least one jet can be ascribed to a softer distribution of the transverse momentum of the hadronic particles recoiling against the diboson system. This explanation is supported by the measurement of a softer than expected p T distribution of the ZZ system [6, 7]. The observed discrepancy may be due to higher-order corrections to ZZ production, not included in MC samples used in this analysis, or to the parton shower modeling. Overall agreement is observed between data and theoretical predictions for all measurements related to the p T -leading and subleading jets. The ∆η jj distribution (Fig. 6, right) measured with 13 TeV data tends to be steeper than the MC predictions, but the differences are not statistically significant.

Summary
The differential cross sections for the production of Z pairs in the four-lepton final state in association with jets in proton-proton collisions at √ s = 8 and 13 TeV have been measured. The data correspond to an integrated luminosity of 19.7 (35.9) fb −1 for a center-of-mass energy of 8 (13) TeV. Cross sections are presented for the production of a pair of Z bosons as a function of the number of jets, the transverse momentum p T , and pseudorapidity of the p T -leading and subleading jets. Distributions of the invariant mass of the two p T -leading jets and their separation in pseudorapidity are also presented. Good agreement is observed between the measurements and the theoretical predictions when next-to-leading order matrix-element calculations are used together with the PYTHIA parton shower simulation. Cross sections for ZZ  Figure 2: Differential cross sections of pp → ZZ → 4 as a function of the multiplicity of jets with |η j | < 4.7 (top panels) and |η j | < 2.4 (bottom panels), for the 8 (left) and 13 (right) TeV data. The measurements are compared to the predictions of MG5 aMC@NLO, POWHEG, and MADGRAPH5 (8 TeV only) sets of samples. Each MC set, along with the main MC generator, includes the MCFM and PHANTOM generators. PYTHIA 6 and PYTHIA 8 are used for parton showering, hadronization, and underlying event simulation, for the 8 and 13 TeV analysis, respectively, with the sole exception of MG5 aMC@NLO, which is always interfaced to PYTHIA 8. The total experimental uncertainties are shown as hatched regions, while the colored bands display the theoretical uncertainties in the matrix element calculations.  [GeV]   production in association with jet have been measured with a precision ranging from 10 to 72% (8 to 38%) at 8 (13) TeV, for jet multiplicities ranging from 0 to ≥ 3. The systematic uncertainty is of the same size, or smaller, than the statistical one. Analyses using future, larger data sets, with smaller statistical uncertainties, will allow the theoretical prediction of ZZ+jets to undergo more stringent tests.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus      [8] CMS Collaboration, "Measurements of the pp → ZZ production cross section and the Z → 4 branching fraction, and constraints on anomalous triple gauge couplings at √ s = 13 TeV", Eur. Phys. J. C 78 (2018) 165, doi:10.1140/epjc/s10052-018-5567-9, arXiv:1709.08601.