Measurement of differential cross sections for the production of a Z boson in association with jets in proton-proton collisions at $\sqrt{s}$ = 13 TeV

A measurement is presented of the production of Z bosons that decay into two electrons or muons in association with jets, in proton-proton collisions at a centre-of-mass energy of 13 TeV. The data were recorded by the CMS Collaboration at the LHC with an integrated luminosity of 35.9 fb$^{-1}$. The differential cross sections are measured as a function of the transverse momentum ($p_\mathrm{T}$) of the Z boson and the transverse momentum and rapidities of the five jets with largest $p_\mathrm{T}$. The jet multiplicity distribution is measured for up to eight jets. The hadronic activity in the events is estimated using the scalar sum of the $p_\mathrm{T}$ of all the jets. All measurements are unfolded to the stable particle-level and compared with predictions from various Monte Carlo event generators, as well as with expectations at leading and next-to-leading orders in perturbative quantum chromodynamics.


Introduction
The production of Z bosons in proton-proton (pp) collisions is described by the Drell-Yan process [1], where a quark and an antiquark from the colliding protons annihilate into a Z boson. At the CERN LHC, this Z boson is commonly produced with accompanying parton radiation via quantum chromodynamics (QCD), which provides a superb opportunity to gain a better theoretical understanding of both strong and electroweak physics in a jet environment. Specifically, events containing Z boson decays into oppositely charged lepton pairs (electrons and muons, but not taus) allow a sensitive evaluation of the accuracy of perturbative QCD [2][3][4] at the highest accessible energies for a broad range of kinematic configurations.
A precise understanding of the pp → Z(→ ℓ + ℓ − ) process is also critical in other standard model (SM) measurements, where it is an important background in studies of the properties of Higgs boson, and in analyses focusing on physics beyond the SM such as searches for dark matter and supersymmetric particles. The clean and readily identifiable signature and relatively large production rate for this process provide an opportunity to accurately constrain the parton distribution functions (PDFs) and probe the strong coupling strength α S .
In addition to these motivations, Z → ℓ + ℓ − + jets production serves as an important experimental benchmark. It is a key ingredient in calibrating the specific parts of the detector and the properties of reconstructed objects, e.g. the jet energy scale. Comparisons of Z+jets events with the expectations from Monte Carlo (MC) event generators and with reliable higher-order calculations can improve confidence in their predictions.
Measurements of differential cross sections for the production of Z bosons in association with jets have previously been reported by the ATLAS, CMS, and LHCb Collaborations in pp collisions at center-of-mass energies of 7 TeV [5-9], 8 TeV [10][11][12], and 13 TeV [13,14], and by the CDF and D0 Collaborations at the Fermilab Tevatron in proton-antiproton collisions at 1.96 TeV [15,16].
In this paper we present measurements of the differential cross sections for the production of Z bosons in association with jets recorded by the CMS Collaboration in 2016 with an integrated luminosity of 35.9 fb −1 . This is an update and an expansion of a previous CMS paper [13] that used 2015 data with an integrated luminosity of 2.19 fb −1 . The events with both electron and muon final states are combined and reconstructed as a pair of oppositely charged leptons that are required to have an invariant mass between 71 and 111 GeV. This mass range optimizes the signal acceptance, rejection of background, and the relative fraction of Z boson and virtual photon contributions.
The new analysis provides measurements of events with up to eight jets inclusively and five jets differentially, compared with the earlier measurement [13] of events with up to six jets inclusively and three jets differentially. Additionally, the ranges for all the observables are extended to larger values of transverse momentum (p T ), and the double-differential cross sections with respect to the leading jet and the Z boson are measured.
The cross sections are measured as a function of jet multiplicity (N jets ) and the individual jet kinematic variables: rapidity (y) and p T where the jets are ordered in decreasing p T . Jet kinematic variables are presented for events with jet multiplicities up to five jets. The term "inclusive" refers to distributions for events with at least N jets jets and the term "exclusive" for distributions where the events contain exactly N jets jets. The cross sections are also measured as a function of the scalar p T sum of the jets (H T ) for events having up to five jets.

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 a silicon pixel and strip tracker, 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 pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
The silicon tracker measures charged particles within the pseudorapidity range |η| < 2.5. During the LHC running period when the data used in this article were recorded, the silicon tracker consisted of 1440 silicon pixel and 15 148 silicon strip detector modules. For nonisolated particles of 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 [17].
The ECAL consists of 75 848 lead tungstate crystals, which provide the coverages |η| < 1.48 in the barrel region and 1.48 < |η| < 3.00 in two endcap regions. Preshower detectors consisting of two planes of silicon sensors interleaved with a total of 3X 0 of lead are located in front of each endcap ECAL. In the barrel section of the ECAL, an energy resolution of about 1% is achieved for unconverted or late-converting photons that have energies in the range of tens of GeV. The remaining barrel photons have a resolution of about 1.3% up to |η| = 1, rising to about 2.5% at |η| = 1.4. In the endcaps, the resolution of unconverted or late-converting photons is about 2.5%, while the remaining endcap photons have a resolution between 3 and 4% [18].
In the region |η| < 1.74, the HCAL cells have widths of 0.087 in pseudorapidity and 0.087 radians in azimuth (ϕ). In the η-ϕ plane, and for |η| < 1.48, the HCAL cells map on to 5 × 5 arrays of ECAL crystals to form calorimeter towers projecting radially outwards from close to the nominal interaction point. For |η| > 1.74, the coverage of the towers increases progressively to a maximum of 0.174 in ∆η and ∆ϕ. Within each tower, the energy deposits in ECAL and HCAL cells are summed to define the calorimeter tower energies, subsequently used to provide the energies and directions of hadronic jets. When combining information from the entire detector, the jet energy resolution amounts typically to 15% at 10 GeV, 8% at 100 GeV, and 4% at 1 TeV, to be compared to about 40, 12, and 5% obtained when the ECAL and HCAL calorimeters alone are used.
Muons are measured in |η| < 2.4, with detection planes made using three technologies, drift tubes, cathode strip chambers, and resistive plate chambers, used in conjunction with the tracker [19].
Events of interest are selected using a two-tiered trigger system [20]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a latency of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 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 is reported in Ref. [21].

Event samples
Candidate events are selected online using single-lepton triggers that require at least one isolated electron (muon) with p l T > 25(24) GeV and |η l | < 2.4. The total trigger efficiency for events within the acceptance of this analysis is greater than 90%. Simulated events for both signal and background processes are produced using various MC event generators, with the CMS detector response modelled using GEANT4 [22]. These events are then reconstructed with the same algorithms used to reconstruct collision data, and the simulated samples are normalized to the integrated luminosity of the data sample using their respective cross sections. For the simulation of the signal, we use a sample generated at next-to-leading order (NLO) with MAD-GRAPH5 aMC@NLO versions 2.2.2 and 2.3.2 (denoted MG5 aMC) [23] using the FxFx merging scheme [24,25]. Parton showering and hadronization are simulated with PYTHIA 8 (version 8.212) [26] using the CUETP8M1 tune [27]. The matrix element includes Z boson production with up to two additional jets generated at NLO with MG5 aMC, effectively yielding leading order (LO) accuracy for Z+3 jets.
The production of Z(→ ℓ + ℓ − ) + jets can be mimicked by various background sources: decays of W bosons resulting from top quark pair production (tt), diboson (WW, WZ, ZZ), triboson (ZZZ, WWZ, WZZ) production, and W bosons produced in association with jets, as well as Z + jets events in which the Z boson decays into (Z → τ + τ − ) + jets where the τ leptons decay leptonically. Background processes are split into two categories: resonant and nonresonant. Resonant background arises from events with genuine Z bosons (WZ, ZZ, tribosons, etc.) and is estimated using MC samples. The nonresonant background that comes from events that do not have a Z boson in the final state (such as tt) is estimated using data events with both an electron and a muon. Events with Z → τ + τ − are considered background and are estimated using the MG5 aMC signal sample at NLO.
Background samples corresponding to electroweak diboson and triboson production [28] are generated at NLO. The POWHEG BOX [29][30][31][32] is used for diboson samples with two leptonic decays (4ℓ, 3ℓν, and 2ℓ2ν) and MG5 aMC for all other diboson (2ℓ2q) and triboson samples. The MADSPIN [33] extension of MG5 aMC is used for diboson samples. For all samples, the NNPDF NLO 3.0 PDF [34] set is used and the generator is interfaced with PYTHIA 8 using the same CUETP8M1 tune as for the signal samples. The samples are normalized to the NLO cross sections calculated with MCFM 6.6 [35].
The simulated event samples include multiple pp collisions within a bunch crossing (pileup). Since the number of pileup interactions varies with the beam conditions, the samples are produced using an approximate pileup distribution. The actual distribution is measured in data and a weight is applied to each simulated event to correct for the difference. as inputs, and the associated missing transverse momentum, which is the negative vector sum of the p T of those jets.
The particle-level objects are defined with a lifetime of cτ > 1 cm (excluding neutrinos) and identified using the same algorithms as used for the data. Leptons are stable particles coming from Z boson decays, dressed by adding the momenta of all photons within R < 0.1 of their directions.
Electron candidates within the geometrical acceptance of |η| < 2.4, excluding the barrel-toendcap (1.444 < |η| < 1.566) transition regions of the ECAL, are reconstructed by combining the information from the ECAL and from the silicon tracker. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with the electron track. This "supercluster" [40] reconstruction efficiency for superclusters with an energy above 5 GeV is close to 100% [41]. To reduce the electron misidentification rate, electron candidates are subject to additional identification criteria that are based on the distribution of the electromagnetic shower in the ECAL, a matching of the trajectory of an electron track with the cluster in the ECAL, and compatibility of the track with the selected primary vertex.
Muon candidates within the geometrical acceptance of |η| < 2.4 are reconstructed with a global fit using both the inner tracking system and the muon spectrometer [19]. The momentum of the muons is obtained from the curvature of the corresponding track; for muons with 20 < p T < 100 GeV the resolution is 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 [19].
Jets are formed from the particles reconstructed by the PF algorithm using the FASTJET software package [39], and the anti-k T jet clustering algorithm [38] with a distance parameter R of 0.4. The jet four-momentum is defined as the sum of the four-momenta of its constituents. The technique of charged hadron subtraction [42] is used to reduce the pileup contribution by removing charged particles that originate from pileup vertices. The jet four-momentum is corrected for the difference observed in the simulation between jets built from PF candidates and generator-level particles. The jet mass and direction are kept constant when the corrections are applied. An offset correction is applied to jet energies to include the contribution from additional pp interactions within the same or previous bunch crossings. Further jet energy corrections are applied for differences between the observed and simulated number of pileup interactions, as obtained from zero-bias events and in the p T balance in dijet, Z + jet, and γ + jet events [43]. The jet energy corrections over the whole p T spectrum and detector acceptance are within 5 to 10% of the generator-level value. Tight identification quality criteria, based on the fraction of energy carried by charged and neutral hadrons, are applied to jets [44] to maximize the reconstruction efficiency while reducing the instrumental background. Jets are required to have |η| < 2.4, to be separated from all selected lepton candidates by at least a distance ∆R = √ (∆η) 2 + (∆ϕ) 2 = 0.4, and have a p T larger than 30 GeV for the single-differential cross sections and 20 GeV for the double-differential ones.
To compare the measured distributions with the theoretical predictions, various experimental corrections are applied after subtracting the total expected background from the observed number of events in each bin. The event acceptance and selection efficiency are estimated with simulation and are used to correct the data. Efficiency corrections are determined from the data using the "tag-and-probe" technique [45] to adjust for the efficiency differences between data and simulation for lepton reconstruction, identification, isolation, and trigger. A correction for the detector resolution effects is implemented using an unfolding technique (as discussed in Section 8).
After offline reconstruction, two leptons are required with the first having p T > 30 GeV and the second having p T > 20 GeV. We require that the two electrons (muons) with highest p T form a pair of oppositely charged leptons with an invariant mass in the range (91 ± 20 GeV). Electrons and muons are considered isolated based on the scalar p T sum of the nearby PF candidates with a distance ∆R < 0.4. For both electrons and muons, medium identification criteria are applied [46,47]. Corrections are applied to the muon momenta to adjust for a residual misalignment in the CMS detector between data and simulation [48].

Observables
In this paper, the cross sections are presented as a function of several kinematic and angular observables to characterize the production mechanisms of Z(→ ℓ + ℓ − ) + jets events. The cross sections are measured as a function of both the exclusive and inclusive jet multiplicities up to a total number of eight jets in the final state. In addition, they are measured as a function of the kinematic variables p T , y and H T for N JETS from one to five. Comparisons of jet multiplicity distributions with predictions from various MC generators show how accurately these generators describe different jet configurations.
The measurement of the distribution of Z boson transverse momentum p T (Z) for events with at least one jet is vital for understanding the balance of the p T between the jets and Z boson, and may be used for comparing theoretical predictions that produce multiple soft gluon emissions in different ways.
The rapidity of Z boson y(Z) is related to the momentum fraction x carried by the partons in the two colliding protons. Therefore, the y(Z) distribution directly reflects the PDFs of the interacting partons. At the LHC, the y(Z) distribution is symmetric around zero, therefore it is appropriate to measure the distribution of Z bosons as a function of |y|.
The distributions of jet H T and jet p T are important because they are sensitive to the effects of higher-order corrections, and provide an accurate estimation of the background from Z+jets process for SUSY searches.
The cross sections are measured as a function of the difference and sum of the rapidity of the leading two jets, |j i − j k |/2 and |j i + j k |/2, for events with at least two jets. For correlations between the Z boson and jets, the cross sections are measured as a function of the difference in rapidity, |y(Z) − y(j k )|/2, and the difference in the azimuthal angle ∆ϕ(Z, j k ). The rapidity sum, |y(Z) + y(j k )|/2, which is correlated with the event boost and is sensitive to the PDFs is also measured. The cross sections are also measured as a function of the difference in azimuthal angle ∆ϕ(j i , j k ) between the ith and kth jets from the p T -ordered list of jets in the event. Since the angular separation ∆ϕ between the Z boson and a jet is sensitive to soft gluon radiation, an advantage of studying the ∆ϕ distribution is that it depends on only the directions of the Z boson and a jet.
Lastly, double-differential cross sections are measured as functions of leading jet p T and y, leading jet and Z boson y, and Z boson p T and y. The measured cross sections are corrected for detector effects and compared with theoretical predictions at LO and NLO accuracy matched with the parton showering as implemented in MC generators.

Phenomenological models and theoretical calculations
We compare the measured Z + jets differential cross sections with three predictions: MG5 aMC at NLO, MG5 aMC at LO, and the GENEVA MC program. The two MG5 aMC calculations (version 2.2.2) [25] are interfaced with PYTHIA 8 (version 8.212) [26]. For the LO MG5 aMC, the generator calculates LO matrix elements (MEs) for five processes: pp → Z + N jets with N = 0, . . . , 4. The NNPDF 3.0 LO PDF is used and α S (m Z ) is set to 0.130. The NLO MG5 aMC prediction includes NLO ME calculations for pp → Z+N jets with N up to 2. The NNPDF 3.0 NLO PDF set is used and α S (m Z ) is set to 0.118. Both predictions use PYTHIA 8 to model the initial-and final-state radiation, parton showers and hadronization with the CUETP8M1 [27] tune that is done with the NNPDF 2.3 [49] LO PDF and α S (m Z ) = 0.130. The ME and parton shower matching is done using the k T -MLM [23,50] scheme with the matching scale set at 19 GeV for the LO MG5 aMC and the FxFx [24] scheme with the matching scale set to 30 GeV for the NLO MG5 aMC.
The third calculation uses the GENEVA 1.0-RC3 (GE) simulation program [51,52], where a nextto-leading-order (NNLO) calculation for Drell-Yan production is combined with higher-order resummation. Logarithms of the 0-jettiness resolution variable, τ, also known as beam thrust and defined in Ref. [53], are resummed at next-to-next-to-leading logarithmic (NNLL) accuracy, including part of the next-to-NNLL corrections. The accuracy refers to the τ dependence of the cross section and is denoted as NNLL ′ τ . The PDF set NNPDF3.1 NNLO [54] is used for this calculation and α S (m Z ) is set to 0.118. The resulting parton-level events are further combined with parton showering and hadronization provided by PYTHIA 8 using the same tune as for MG5 aMC.
In this analysis, uncertainties in the ME calculation for the MG5 aMC and GENEVA predictions are estimated using the procedure recommended by the authors of the respective generators. For the MG5 aMC prediction, the factorization (µ F ) and renormalization (µ R ) scales are varied by a factor of 0.5 and 2 to estimate the uncertainty coming from missing higher-order terms in the fixed-order calculation. An envelope of the six variations is used with the two extremes (one scale varied by a factor 0.5 and the other by a factor 2) excluded. For the GENEVA sample, µ F and µ R are simultaneously varied by 0.5 and 2, leading to two combinations, their uncertainties are symmetrized by using the maximum of the up and down uncertainties for both cases. The uncertainty from the resummation in GENEVA is estimated using six profile scales [55,56], as described in Ref. [51], and added in quadrature to the scale uncertainty. The PDF uncertainty in the MG5 aMC sample is estimated using the set of 100 replicas of the NNPDF 3.0 NLO PDF and the uncertainty in the α S value used in the ME calculation is estimated by varying it by ±0.001. The PDF and α S uncertainties are added in quadrature to the ME calculation uncertainties. For both MG5 aMC and GENEVA, all these uncertainties are obtained using the reweighting method [51,57] implemented in these event generators.

Background estimation
Two categories of background events are considered: resonant and nonresonant. The resonant background, which consists mainly of multiboson events with at least one Z boson produced in the final state, is estimated using simulation. The background from nonresonant events containing two leptons primarily from W boson decays such as those appearing in tt is estimated from data events. The decay Z → τ + τ − is considered as a background and is estimated from the MG5 aMC signal MC sample. The backgrounds from events where one or two jets are misidentified as a lepton, such as W+jets or multijets, is negligible.
The method used for estimating the nonresonant background uses a control region in data containing events with one electron and one muon e ± µ ∓ passing all other signal region criteria. This control region is used to estimate the nonresonant background in the signal region by applying a transfer factor to account for cross section and lepton efficiency differences between channels. Assuming lepton flavor symmetry, the cross section for the e ± µ ∓ and either e + e − or µ + µ − channel differs only by a factor of 2. The difference in the efficiency between electrons and muons is estimated using the total yields of the two channels. Resonant signal and background are estimated in the control region with the same signal simulation and subtracted to avoid double counting.
The background-subtracted Z → e + e − and Z → µ + µ − event numbers are 10 and 18 million respectively. The kinematic properties of the Z boson and the leading jet, and measurement of jet multiplicity are shown in Figs. 1-3 together with the results of the simulation. Background samples corresponding to diboson electroweak production is denoted as "VV", and nonresonant background samples are denoted as "NRB" in the figure legends. The fraction of background events is small compared with the signal and amounts to approximately 1% for ≥0 jets increasing to 10% at 5 or more jets. For p T variables, the background increases from 1% below 100 GeV to 10% in the high-p T tails.

Unfolding procedure
In this analysis, unfolding is performed to remove detector effects and estimate the particle level distributions in data. The MG5 aMC MC sample is used to extract the detector transformation, called the response matrix, that feeds into the unfolding algorithm. The unfolding procedure consists of performing a least-squares fit with optional Tikhonov regularization [58], as implemented in the TUnfold software package [59]. In this analysis the best value for the regularization parameter is chosen using the L-curve method [60]. Closure tests are performed by checking the unfolded distributions with the original data.
The momenta of the leading leptons are summed to obtain the particle-level Z boson momentum. The particle-level objects are required to pass the same kinematic selections as at detector level.

Systematic uncertainties
The sources of experimental uncertainty are divided into the following categories: Jet energy scale (JES) and jet energy resolution (JER); lepton efficiencies (identification, isolation, and track reconstruction); lepton energy scale (LES) and resolution (LER); trigger efficiency; luminosity; pileup; background and unfolding uncertainties. The uncertainties listed above are assumed to be independent such that each is computed individually and added in quadrature to obtain a total uncertainty. To compute the systematic uncertainty from each source, the analysis is repeated using the source values increased and decreased by 1σ from the central value. This results in bin-by-bin uncertainty contributions from each source in the unfolded distributions.
The JES uncertainty originates mainly from the uncertainty in the single-particle response. It is the dominant source of systematic uncertainty. It affects the reconstruction of the transverse energy of the selected jets. In this analysis, jet energy corrections (JEC) are applied to include inefficiencies, nonlinearities and finite resolution in energy and position of the reconstructed jets. The effect of the JES uncertainty is studied by scaling the reconstructed jet energy up and down by p T and η-dependent scale factors. A similar procedure is followed for the JER. The Data  Figure 1: The Z boson candidate p T (upper) and |y| (lower) for events with at least one jet. The muon (left) and electron (right) channels are shown separately. The background is estimated from both simulation and data driven methods (such as nonresonant background, NRB) as described in Section 7. The error bars around the data points represent the statistical uncertainties. The distribution ratio of simulation to data is shown in the bottom frames, with error bars that represent the total statistical uncertainties from the data and simulation samples. Data ee Data ee →  Data uncertainties in the JES and JER vary from 1-11% as a function of jet multiplicity.
Scale factors for lepton efficiencies are applied on an object-by-object basis so that the simulation samples reflect the inefficiencies observed in data. The lepton identification, isolation, track reconstruction and trigger efficiencies in simulation are corrected with scaling factors derived with a tag-and-probe method and applied as a function of lepton p T and η. To estimate the uncertainties, the total yield is recomputed with the scaling factors varied up and down by the fit uncertainties. The uncertainty associated with lepton efficiency in the electron (muon) channel is 1% (0.5%).
The LES and LER uncertainties make a small contribution to the overall lepton uncertainties of ∼1% for each channel.
A normalization uncertainty of 2.5% is assigned to account for the imperfect knowledge of the integrated luminosity [61]. This uncertainty is propagated to the measured differential cross sections.
The uncertainty coming from the pileup model is estimated by varying the amount of pileup events by 4.6% up and down [62] when reconstructing the response matrices from the simulation. The difference in the unfolded data is the uncertainty.
The uncertainty in the unfolding procedure comes from: (1) the statistical uncertainty in the response matrix coming from the finite size of the MC sample used to compute it; and (2) the possible event generator dependence of the response matrix itself. Because of the finite binning, a different distribution could lead to a different response matrix. This uncertainty is estimated by weighting the MC to agree with the data in each distribution to be unfolded and building a new response matrix. The weights are extracted from the data-to-MC ratio of a fine-binned histogram at the reconstruction level. The fine binning allows us to account for the effect of the distribution of events within each measurement bin. The difference between the nominal results and the results unfolded using the alternative response matrix is assumed the systematic uncertainty. Statistical fluctuations in the response matrix are propagated analytically in the TUnfold package.
Lastly, the background samples are varied by their corresponding cross section uncertainty before being subtracted from data prior to unfolding.

Results
The measurements from the electron and muon channels are consistent with each other within the statistical and systematic uncertainties, and hence they are combined. To combine the two channels, a hybrid method based on the weighted mean and the best linear unbiased estimates method [63,64] is used to calculate the cross section values and uncertainties. The covariance matrix of the combination is calculated assuming that all uncertainty sources are correlated between channels except the statistical components and those associated with the lepton reconstruction and identification. Figure 4 shows the measured cross sections as a function of N jets for a total number of up to eight jets in the final state. The trend of the jet multiplicity represents the expectation of the perturbative QCD prediction for an exponentially falling spectrum with the number of jets. The agreement is good for the exclusive distributions for all the theoretical estimations, remaining within the uncertainties and going up to the maximum number of final-state partons included in the ME, namely four in the MC generators used here. The GENEVA generator predicts a steeper spectrum than observed due to the lack of hard jets at ME level beyond two.
The size of the 2016 data samples allows us to determine the differential cross sections for jet multiplicities up to eight jets, and to measure the cross sections as a function of several kinematic observables up to five jets. The combined single-differential cross sections are shown in Figs. 5-22, while double-differential cross sections are given in Figs. 23-25. All results are compared with theoretical predictions from MG5 aMC at LO and at NLO. Since the GENEVA predictions are effectively LO in QCD at two jets, only the results with at least one or two jets are compared with GENEVA.
The jet y and p T up to five leading jets can be seen in Figs. 5-9. For both quantities, the data distributions are well reproduced by the simulations. The MG5 aMC at LO and NLO, describe the data well in general. The GENEVA prediction shows good agreement for the measured p T and y of the first jet, although it underestimates the cross section at low p T in the second jet.
In addition, the inclusive jet differential cross sections as a function of H T for events with at least one, two and three jets, respectively, are presented in Fig. 10. The MG5 aMC predictions at both LO and NLO are compatible with the measurement. The cross section at higher values of H T is slightly overestimated, but the discrepancy is compatible with the theoretical and experimental uncertainties. The slopes of the distributions for the first two jet multiplicities predicted by GENEVA do not describe the data.
The measured cross section as a function of the dijet mass is also shown in Fig. 11. The three predictions considered here agree with the measurement within the experimental uncertainties, except for a dijet mass below ∼50 GeV, where the predictions made with GENEVA show a deficit with respect to the measurement. The MG5 aMC at NLO generator has better agreement with the measurement in this region. The MG5 aMC at LO generator predicts a distribution that falls more steeply for a dijet mass above ∼100 GeV.
The rapidity distributions of the Z boson and jets are reasonably well modelled by the predictions, but the correlations between the rapidities, which have been shown by measuring multidimensional differential cross sections and distributions of rapidity differences and sums (Figs. 13-17), are not well described by the multileg LO calculation. We have shown that the NLO multileg event generator reproduces the rapidity difference distributions well. The rapidity sum is also successfully described. For this variable the discrepancy with the LO calculation could be due to a different choice of the parton distribution functions. The azimuthal angles between the Z boson and the jet (Figs. [18][19][20] and between the jets (Figs. [21][22] are well described by the predictions including the LO one. The results for the double-differential cross sections are presented in Figs. 23-25 and are compared with the predictions described in Section 6. The double-differential cross sections are shown for at least one jet as a function of the leading jet p T and y (Fig. 23), leading jet and Z boson y (Fig. 24), Z boson p T and y (Fig. 25). In general, all the predictions are in agreement with the data, and the NLO MG5 aMC prediction provides a better description than the LO MG5 aMC and GENEVA predictions for double-differential cross sections in Fig. 25. In the low-p T region GENEVA gives good description as expected by the resummation.
Overall, the MG5 aMC at NLO predictions describe the data within theoretical uncertainties over a range of kinematic variables. In the regions of NLO accuracy, such as the first and second jet p T and y, the agreement is within 10% up to the TeV scale.
The differential cross section results with covariance matrices are presented in HEPData [65].      for events with at least three jets. The measurement statistical (resp. systematic) uncertainties are represented with vertical error bars (resp. hashed areas). The measurement is compared to the NLO and LO MG5 aMC predictions described in Section 6. The predictions uncertainties, estimated as explained in this section, are represented by colored areas in the ratio plots (light color for the statistical part and darker color for the total). Only the statistical uncertainty is displayed for the LO prediction.     10: The measured differential cross section as a function of H T for events with at least one (left), two (right), and three (bottom) jets. The measurement statistical (resp. systematic) uncertainties are represented with vertical error bars (resp. hashed areas). The measurement is compared to the NLO MG5 aMC, LO MG5 aMC, and GENEVA (for N jets ≥ 1 and N jets ≥ 2) predictions described in Section 6. The predictions uncertainties, estimated as explained in this section, are represented by colored areas in the ratio plots (light color for the statistical part and darker color for the total). Only the statistical uncertainty is displayed for the LO prediction.          Figure 19: The measured differential cross section as a function of the Z boson and subleading jet azimuthal difference for events with at least two (left) and three (right) jets. The measurement statistical (resp. systematic) uncertainties are represented with vertical error bars (resp. hashed areas). The measurement is compared to the NLO MG5 aMC, LO MG5 aMC, and GENEVA (for N jets ≥ 2) predictions described in Section 6. The predictions uncertainties, estimated as explained in this section, are represented by colored areas in the ratio plots (light color for the statistical part and darker color for the total). Only the statistical uncertainty is displayed for the LO prediction.   Figure 21: The measured differential cross section as a function of the leading and subleading jet azimuthal difference for events with at least two (left) and three (right) jets. Details on the presentation of the results are given in Fig. 19.    theo. unc.

Measurement Prediction
Figure 25: Double differential cross sections as a function of Z boson p T and |y| for events with at least one jet. Details on the presentation of the results are given in Fig. 4.

Summary
The production of Z bosons, decaying into a pair of electrons or muons, in association with jets has been studied in proton-proton collisions at a center-of-mass energy of 13 TeV at the LHC in 2016 by the CMS experiment using a data set with an integrated luminosity of 35.9 fb −1 . Differential cross sections have been measured for Z bosons decaying to electrons or muons with transverse momentum p T > 25 GeV and pseudorapidity |η| < 2.4 requiring at least one jet with p T > 30 GeV and |η| < 2.4. Differential cross sections have been measured as a function of the exclusive and inclusive jet multiplicities (N jets ), the p T of the Z boson, and kinematic variables that include jet transverse momenta, the scalar sum for up to five inclusive N jets , rapidity, dijet invariant mass (M jj ) and their sum values.
The results, corrected for detector effects through unfolding, are compared with three theoretical predictions: (1) the expectations are computed from particle-level simulations using merged leading-order (LO) calculations with the k T -MLM parton-showers and matrix-element matching scheme; (2) using next-to-leading-order (NLO) calculations and the FxFx merging scheme; and (3) the GENEVA MC program, where a next-to-NLO (NNLO) calculation for Drell-Yan production is combined with higher-order resummation.
High precision is achieved in measuring the cross sections using the latest experimental methods and larger sets of data than were available in previous CMS publications. The increased number of events allows to extend the kinematic range to higher values of p T and mass. The measurements presented in this paper provide a detailed description of the topological structure of Z → ℓ + ℓ − + jets events that is complementary to the existing measurements of rates and associated jet multiplicities.
The kinematics of Z+jets events is studied in detail. The measured differential cross sections and N jets distributions are within the experimental and theoretical uncertainties. Some deviations are observed for N jets > 3. Such discrepancies offer the possibility of using these data to further improve the modeling. The results also indicate that multiparton NLO calculations can be used to estimate the Z → ℓ + ℓ − + jets contributions to measurements and searches at the LHC.

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 and other centers 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, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies: BMBWF and FWF  [5] ATLAS Collaboration, "Measurement of the production cross section of jets in association with a Z boson in pp collisions at √ s = 7 TeV with the ATLAS detector", JHEP 07 (2013) 032, doi:10.1007/JHEP07(2013)032, arXiv:1304.7098.