Study of Z bosons plus jets events using variables sensitive to double-parton scattering in pp collisions at 13 TeV

Double-parton scattering is investigated using events with a Z boson and jets. The Z boson is reconstructed using only the dimuon channel. The measurements are performed with proton-proton collision data recorded by the CMS experiment at the LHC at $\sqrt{s} =$ 13 TeV, corresponding to an integrated luminosity of 35.9 fb$^{-1}$ collected in the year 2016. Differential cross sections of Z + $\geq$ 1 jet and Z + $\geq$ 2 jets are measured with transverse momentum of the jets above 20 GeV and pseudorapidity $|\eta|$ $\lt$ 2.4. Several distributions with sensitivity to double-parton scattering effects are measured as functions of the angle and the transverse momentum imbalance between the Z boson and the jets. The measured distributions are compared with predictions from several event generators with different hadronization models and different parameter settings for multiparton interactions. The measured distributions show a dependence on the hadronization and multiparton interaction simulation parameters, and are important input for future improvements of the simulations.


Introduction
Proton-proton (pp) collisions at high energies result in many events with a Z boson produced in association with jets at large transverse momentum p T . The Z+jets process is an important background for many standard model (SM) measurements and searches for physics beyond the SM. Measurements of differential cross sections for Z boson production in association with jets can be used to test models of initial-state radiation, and multiparton interactions (MPI). Monte Carlo (MC) models are used to simulate collision processes and play a crucial role in understanding background for the new physics searches. The description of MPI is an important part of these simulations. However, MPI can not be completely described by perturbative quantum chromodynamics (QCD), so this requires a phenomenological description involving parameters that must be tuned with the help of data.
The ATLAS, CMS, and LHCb experiments have reported various measurements of Z+jets processes at centre-of-mass energies of 7 TeV [1-6], 8 TeV [7-9], and 13 TeV [10,11]. Z boson production in association with jets is a process with a clean experimental signature and is well understood theoretically. In this paper, measurements of Z+jets events are performed to explore observables sensitive to the presence of MPI. A possible contribution could come from the simultaneous occurrence of two parton-parton interactions and is usually, when the involved scale is large, called double-parton scattering (DPS). This analysis explores new signatures sensitive to DPS in events with a Z boson and one or more jets, and compares the data with MC simulations, using a variety of event generators simulating Z+jets processes with different MPI and hadronization models.
The ATLAS, ALICE, CMS, and LHCb Collaborations have previously reported measurements of DPS in various topologies, such as WW [12], W+jets [13,14], 4 jets [15], J/ψ production [16][17][18][19], double charm production [20, 21], but not in the Z+jets topology. Figure 1 shows typical diagrams of single parton scattering (SPS) and DPS production of Z+2 jet events. In SPS, the Z boson decaying into two muons and the two jets come from the same parton-parton interaction, whereas in the case of DPS, the Z boson and the two jets originate from two independent interactions. Events are categorized as Z+ ≥1 jet and Z+ ≥2 jets topologies, and the corresponding integrated cross sections are measured in a fiducial region. For consistency with previous DPS measurements at 7 TeV [13,14], jets are required to have a p T threshold of 20 GeV. Decreasing the jet p T threshold would increase the DPS contribution, but will also lead to a significant increase in the experimental uncertainties related to jet identification. With this threshold of 20 GeV, the DPS contribution is measured to be around 5.5% at 7 TeV [14] and, as estimated from simulation, at 13 TeV it is expected to increase by 20-30% because of the increased parton density. The differential cross sections are measured as functions of various observables based on the azimuthal (φ, in radians) separation and the p T balance between the Z boson and jets, as well as the p T balance between two jets in case of Z+ ≥2 jets events [13,14,22]. Since the two interactions in DPS are largely uncorrelated in the phase space studied, the shape of the distribution of these observables is expected to differ from that seen in SPS. Distributions of these observables normalized to the integrated cross section are measured because they have a lower systematic uncertainty due to the cancellation of the uncertainties correlated among bins. Previously measured distributions were used to extract parameters for DPS modelling in MC event generators [23,24].
The measurement is performed using events where the Z boson decays into two oppositely charged muons. The pp collision data used in the analysis were collected in 2016 at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 . Tabulated results are provided in the HEPDATA [25].

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 covering a pseudorapidity region of |η| < 2.5, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters, made of steel and quartz fibers, extend the η coverage provided by the barrel and endcap detectors to |η| < 5.0. Muons are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive plate chambers [26].
Events of interest are selected using a two-tiered trigger system. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of about 100 kHz within a latency of less than 4 µs [27]. 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 [28].
In addition to the Z+jets MG5 aMC samples with the CUETP8M1 tune described above, the measurements are compared with simulations using different PYTHIA8 tunes, such as CP5 [23], CP5 without MPI, and CDPSTP8S1-Wj [24]. The CP5 tune uses NNPDF3.1 PDF set at NNLO, with α S values of 0.118, and running according to NLO evolution. The CP5 tune is chosen since it is the standard tune obtained by fitting a large number of 1.96, 7, and 13 TeV measurements sensitive to soft and semihard multipartonic interactions [23]. The W+dijet DPS tune, CDPSTP8S1-Wj is derived from the parameters of PYTHIA8 tune 4C, with a variation of the impact parameter dependence, i.e. matter overlap function, which is the convolution of the matter distributions of the two incoming hadrons. The DPS simulation in MC is usually quantified in terms of a parameter known as effective cross section σ eff . In PYTHIA8, the value of σ eff is calculated by dividing the nondiffractive (ND) cross section by the so-called "enhancement factor", which depends on the parameters of the overlap matter distribution function and the limiting value of p T at a reference energy [52]. For central collisions, the enhancement factor tends to be large, translating to a lower value of σ eff and larger DPS contribution. For peripheral interactions, enhancement factors are small, giving large values of σ eff and a small DPS contribution.
The Z+jets events calculated at NLO with MG5 aMC are also interfaced with HERWIG7, apply- ing the CH3 tune for the underlying event description, hadronization, and showering [53][54][55][56]. The tune is derived by fitting measurements from pp collision data collected by the CMS experiment at √ s = 0.9, 7, and 13 TeV. The CH3 tune makes use of the NNPDF3.1 LO PDF set with α S = 0.130 for the simulation of MPI, but the NNLO PDF set with α S = 0.118 for all other components. The cross sections of Z+jets events simulated with different MC event generators and configurations are normalized to NNLO calculations with FEWZ v3.1 [57]. Table 1 summarizes all the event generators, PDF sets, and tunes used to produce both signals (including alternative samples), and the background processes.

Event selection
The particle flow (PF) algorithm [58] is used to reconstruct and identify individual particle candidates in an event, with an optimized combination of information from the various elements of the CMS detector. Energy deposits are measured in the calorimeters and charged particles are identified in the central tracking and muon systems.
Events are selected with a single-muon trigger requiring at least one isolated muon candidate with p T > 24 GeV. The primary vertex is the vertex candidate with the largest value of summed physics-object p 2 T . The physics objects are the track-only jets, clustered using the jet finding anti-k T algorithm [59,60] with the tracks assigned to candidate vertices as inputs, and the associated missing p T , which is the negative vector p T sum of all physics objects. Events are required to have at least two oppositely charged muons with p T > 27 GeV and |η| < 2.4. These muons are reconstructed by combining information from the inner tracker and the muon detector subsystems. The muon candidates are required to satisfy identification criteria based on the number of hits in each detector, the quality of the muon-track fit, and the consistency with the primary vertex, which is imposed by requiring the longitudinal and transverse impact parameters are less than 0.5 and 0.2 cm, respectively. The efficiency to reconstruct and identify the muons is greater than 96%. Matching muon candidates to tracks measured in the silicon tracker results in a relative p T resolution for muons with 20 < p T < 100 GeV of 1% in the barrel and 3% in the endcaps.
To suppress the multijet background, muons are required to be isolated. The relative isolation variable, I rel , for muons is defined as: (1) Here ∑ E neutral T and ∑ E γ T are the transverse energy sums of neutral hadrons and photons, respectively, within a cone of radius ∆R = √ (∆η) 2 + (∆φ) 2 = 0.4 around the muon track. The quantity ∑ p charged T represents the p T sum of the charged hadrons in the same cone around the muon associated with the selected vertex. Finally, ∑ p PU T is the p T sum of the charged hadrons in the same cone around the muon not associated with the selected vertex. A muon is considered isolated if I rel < 0.2. There are small residual differences in the trigger, identification and isolation efficiencies between data and simulation, which are measured using the "tag-and-probe" method [61], and included by applying scale factors to simulated events [26].
To select Z boson candidate events, the invariant mass of the two oppositely charged muons with highest p T is required to be close to the Z boson mass, 70 < m µ + µ − < 110 GeV. The Z boson candidate is required to be accompanied by at least one jet with p T > 20 GeV and |η| < 2.4. The overlap between the muons from the Z boson decay and the jets is removed by requiring a minimum ∆R distance of 0.4 between them.
Jets are clustered from PF candidates using the infrared-and collinear-safe anti-k T algorithm [59] with a distance parameter of 0.4, as implemented in the FASTJET package [60]. The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be, on average, within 5 to 10% of the true momentum over the entire p T spectrum and detector acceptance. Pileup can contribute additional tracks and calorimetric energy depositions, increasing the apparent jet momentum. To mitigate the effects of pileup, tracks identified as originating from pileup vertices are discarded, and a factor [62] is applied to correct for the remaining contributions. Jet energy corrections are derived from simulation studies so that the average measured response of jets becomes identical to that of particle-level jets. In situ measurements of the momentum balance in dijet, multijet, photon+jet, and Z+jets events with the Z boson decaying leptonically are used to correct residual differences between the jet energy scales in data and simulation [63,64]. The jet energy in simulation is spread to match the resolution observed in data. The jet energy resolution (JER) in data amounts typically to 15-20% at jet energy of 30 GeV, 10% at 100 GeV, and 5% at 1 TeV. Additional selection criteria are applied to remove jets potentially dominated by anomalous contributions from various subdetector components or reconstruction failures [65]. Jets identified as being likely to originate from pileup [66,67] are also removed by using a pileup jet ID discriminator for jets with 20 < p T ≤ 50 GeV. No pileup jet ID is applied to jets with p T ≥ 50 GeV, since here the pileup event contribution is negligible [67].
The selected events correspond to a sample of Z+jets events with a background of 2-5%, which is subtracted from the data before the unfolding. The dominant contribution is the production of tt pairs. The simulation of tt production is validated with data in a tt → eµ + X control region. The tt control region is constructed using data events by requiring an oppositely charged electron-muon pair with at least one jet. 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 originating from the electron track [68]. The longitudinal and transverse impact parameters for barrel (endcap) are required to be less than 0.10 (0.20) and 0.05 (0.10) cm, respectively. Electrons are required to have p T > 27 GeV and |η| < 2.4. The selection criteria for muons and jets are the same as discussed above. The invariant mass of the oppositely charged eµ pair is required to lie within the range of 70 < m eµ < 110 GeV. The difference in the data to simulation comparison in this control region is included as a part of the systematic uncertainty.

Observables and the unfolding
Events are categorized in the Z+ ≥1 jet and Z+ ≥2 jets subsets, and the differential cross sections are measured as function of the following observables: 1. Z+ ≥1 jet events • Azimuthal angle between the Z boson and the highest p T (leading) jet: ∆φ(Z, j 1 ).
For SPS events, the Z boson and the leading jet will balance each other, hence this observable will peak around π, whereas for the DPS process the distribution is expected to be flat because of the absence of a correlation between the Z boson and the leading jet produced from the two different scatterings. • Relative p T imbalance between the Z boson and leading jet is This observable is expected to be close to zero for SPS events, whereas in the case of DPS this observable will have higher values.

Z+ ≥2 jets events
• Azimuthal angle between the Z boson and dijet system: ∆φ(Z, dijet). Here, dijet means the resulting three-momentum of the leading and subleading jets. For SPS events, the dijet system p T will balance the Z boson p T , therefore, this observable will peak around π. In the case of DPS production, the distribution is expected to be flat since the Z boson and the two jets are originating from two independent scatterings. • Relative p T imbalance between the Z boson and dijet system is For SPS events this observable is expected to be 0, whereas for DPS events this observable will have higher values. • Relative p T imbalance between the leading (j 1 ) and subleading (j 2 ) jets is For DPS events, the two jets are expected to balance each other, therefore this observable will be around 0. For SPS events, the two jets are correlated with the Z boson and not expected to balance each other.
The reconstructed distributions are corrected for the event selection efficiency and detector resolution using an unfolding technique that employs a response matrix to map the reconstructed observables onto the particle-level values. The unfolding is done with 20 detector-level and 10 particle-level ∆φ or ∆ rel p T bins. The unfolding is performed using the TUNFOLD package [69], which is based on a least squares fit with a possible Tikhonov regularization term [70]. Since the effect of regularization is minimal on the reported observables, the unfolding is performed without regularization.
The Z+jets events simulated with MG5 aMC + PYTHIA8 with tune CUETP8M1 are used to construct the response matrix. At particle level, events are required to have at least two oppositely charged muons with p T > 27 GeV and |η| < 2.4. The particle-level definition of a muon corresponds to a generator-level muon coming from the Z boson decay, "dressed" by adding the momenta of all photons within ∆R < 0.1 around both muon directions to account for the FSR effects [71]. The particle-level jets with p T > 20 GeV and |η| < 2.4 are formed from stable particles (cτ > 1 cm), except neutrinos, using the same anti-k T jet algorithm as for reconstructed jets. A possible overlap between particle-level jets and a pair of muons from the Z boson decay is removed by requiring a minimum distance of 0.4 between them. The distributions are unfolded to the particle level in the fiducial region defined in Table 2.
In simulation, the reconstructed jets and a pair of muons are spatially matched to the corresponding particle-level objects by requiring that they are within ∆R of 0.4 from one another. Events that have reconstructed objects without matched particle-level objects are included in the background category and are excluded from the sample. This contribution includes the events where the selected jets originate from pileup. The contribution of background with no jet at particle level, but at least one jet at the reconstructed level is about 4.1%. The contribution of background with no jet (1 jet) at particle level and at least 2 jets at the reconstruction level is around 1.1 (4.8)%. The simulation of pileup jets is validated in a control region enriched in pileup jets, obtained by inverting the criteria used to reject the pileup jets. The simulation describes the data well, within the uncertainties, in the pileup-enriched control region, validating the simulation of background events from pileup.
Events that have particle-level objects in the fiducial volume, but no matching reconstructed objects, are accounted for with acceptance and efficiency corrections. In addition, there may be events in which the particle-level jet passing the fiducial selection does not lead to a reconstructed level jet that passes the fiducial selection. Events of this type are considered as background at the reconstruction level and are not considered. However, at the generator level, these are genuine signal events missed because of the detector and reconstruction inefficiencies. These missed events are accounted for via a signal acceptance correction.
Finally, the unfolded distributions are scaled with the inverse of the integrated luminosity to obtain the differential cross section.

Systematic uncertainties
The measurements have various sources of systematic uncertainties.
• Jet energy resolution and scale: The effect of the uncertainty in jet energy scale (JES) or JER [63,64] is evaluated by varying the JES or JER within the associated uncertainty and performing the unfolding procedure with the modified distribution. The variations of JES corrections within their uncertainties change the differential cross sections by 2-8%, whereas the area-normalized distributions are affected by up to 4%. The variations of JER corrections within their uncertainty change the differential cross section distributions by 1-7% and area-normalized distributions up to 5%.
• Pileup jet identification: Simulated events are corrected for the differences in the jet identification efficiency between data and simulated events. The uncertainties in these corrections affect the measurements by up to 1-2% for differential cross section distributions and less than 0.5% for area-normalized distributions.
• Closure uncertainty: The effect of model dependence is evaluated by comparing unfolded results obtained using response matrices constructed with the MG5 aMC + PYTHIA8 with tune CUETP8M1 and SHERPA generators having different ME and parton showering, as discussed in Section 3. The calculated uncertainty is then symmetrized.
The effect of scale uncertainties is estimated using a set of generator weights that correspond to variations of renormalization (µ F ) and factorization (µ R ) scales up and down by factors of 2 from their nominal values, excluding the pair of extreme variations. The unfolded distributions are obtained for all such combinations and their envelope is quoted as the uncertainty.
The uncertainty in the PDFs is estimated using the 100 replicas of the NNPDF 3.0 PDF set [72]. The unfolded distributions are reproduced using the weights of the replicas and a standard deviation is computed on a bin-by-bin basis [72].
These sources, added in quadrature, affect the differential cross section by 1-5% and area-normalized distributions up to 1-4% in the case of Z+ ≥1 jet events, whereas the effect is within 9% for the differential cross section and up to 7% for the areanormalized distribution in the case of Z+ ≥2 jet events. • Integrated luminosity: It is determined with a 2.5% [73] uncertainty for differential cross section distributions, but completely cancels in the area normalized distributions.
• Pileup weighting: The distribution of the mean number of interactions per bunch crossing of the simulated samples is weighted to match that of the data. The uncertainty related to pileup weighting is estimated by varying the total inelastic cross section by ±4.6% [74]. The effect is negligible for both the differential and areanormalized distributions. After the pileup weighting, the vertex multiplicity in simulation shows overall good agreement with data, but there is an overestimation of around 40% for high vertex multiplicities. To investigate the effect of this residual discrepancy, the simulated events are additionally corrected to reproduce the vertex multiplicity distributions observed in data. The data are unfolded with the weighted response matrix, and the results are compared with the unfolded results without weighting. These sources, added in quadrature, affect the differential cross section by 1.0-1.5% and area-normalized distributions up to 1%.
• Muon selection: The systematic uncertainty related to various muon selection criteria such as muon identification, isolation and trigger scale factors is less than 1% for both differential cross section and area-normalized distributions.
The effect of the muon momentum corrections on the measurement is very small (<0.1%), therefore no additional systematic uncertainty is assigned.
• Background modelling: There is a small contribution from tt events, whereas the contribution from other processes such as dibosons, W+jets, and QCD multijet pro-duction is smaller than 1%. To calculate the systematic uncertainty related to the simulated background contribution, the cross sections of the background samples are varied by their uncertainties. The systematic uncertainty related to the tt process is obtained from the differences between data and simulation in a tt → eµ + X control region. The variation of the background contribution within the uncertainties affects the differential cross section less than than 0.2% for Z+ ≥1 jet and less than 0.6% for Z+ ≥2 jets events. The effect of this variation on the area-normalized distributions is less than 0.2%.
Tables 3 and 4 summarize the effect of various systematic uncertainties for the differential cross section and the normalized distributions. These systematic uncertainties are considered uncorrelated and are added in quadrature. Table 3: Uncertainty sources and their effect on the differential cross section distributions.

Results
The production cross sections in the fiducial region defined in Table 2 are measured to be 158.5 ± 0.3 (stat) ± 7.0 (syst) ± 1.2 (theo) ± 4.0 (lumi) pb for Z+ ≥1 jet events and 44.8 ± 0.4 (stat) ± 3.7 (syst) ± 0.5 (theo) ± 1.1 (lumi) pb for Z+ ≥2 jet events. The measured cross sections are described, within the uncertainties, by different simulations, except for the MG5 aMC + PYTHIA8 with CP5 tune MPIOFF and the DPS-specific CDPSTP8S1-WJ tune. The cross section of the DPS-specific tune is predicted to be 10% higher than the measured cross section. The cross section of CP5 tune without MPI underestimates the measured cross section up to 10% for Z+ ≥1 jet events and up to 16% for Z+ ≥2 jet events. Table 5 summarizes the measured and predicted cross sections for the Z+ ≥1 jet and Z+ ≥2 jet processes.
For Z+ ≥ 1 jet events, Fig. 2 (3) shows the differential cross section measurements (left) and the area-normalized distributions (right) as a function of ∆φ(Z, j 1 ) (∆ rel p T (Z, j 1 )), respectively. Different MC event generators (except for the MG5 aMC + PYTHIA8 with the DPS-specific tune CDPSTP8S1-WJ) describe, within the uncertainties, the overall differential cross section as a function of ∆φ and ∆ rel p T , apart from a few discrepancies in some specific regions of these observables. MG5 aMC + PYTHIA8 generator prediction with the DPS-specific tune CDPSTP8S1-WJ overestimates the cross section up to 10-20%, but correctly describe the shapes of the ∆φ and ∆ rel p T distributions. The MG5 aMC + PYTHIA8 prediction (with CP5 tune) overestimates (up to 20%) the measurement in the lower-∆ rel p T , where SPS is expected to be dominant. The prediction of MG5 aMC + HERWIG7 describes, within the uncertainties, the shape of the ∆ rel p T distribution, but deviates up to 15-20% in the lower-∆φ region where DPS is expected to contribute more.
The Z+jets calculation of MG5 aMC + PYTHIA8 without MPI does not describe the measurement and is lower than the measurement by up to 50% in both the lower ∆φ and higher ∆ rel p T regions, where the MPI contribution is expected to be the largest.
For Z+ ≥2 jet events Fig. 4 (5) shows the differential cross section measurements (left) and the area-normalized distributions (right) as a function of ∆φ(Z, dijet) (∆ rel p T (Z, dijet)), respectively. The differential cross section, as a function of ∆φ, is reasonably well described by the different predictions. The predictions of MG5 aMC + PYTHIA8 (with the CDPSTP8S1-WJ tune) and SHERPA best describe the shape of the measured distributions, whereas the MG5 aMC + PYTHIA8 (with CP5 tune) prediction deviates by up to 15% in the lower-∆φ region, where DPS production of Z+jets is expected to contribute more. The shape of the ∆ rel p T distribution is best described by the predictions of MG5 aMC + PYTHIA8 (with the CDPSTP8S1-WJ tune) and SHERPA. The MG5 aMC + PYTHIA8 (with CP5 tune) prediction overestimates, up to 20%, in the lower-∆ rel p T region.
The differential cross section as a function of ∆ rel p T between two jets, as shown in Fig. 6, is well described by the different predictions except for MG5 aMC + PYTHIA8 (with the CDPSTP8S1-WJ tune), which overestimates the differential cross section measurements up to 15%. The shape of the ∆ rel p T distribution is described well by the predictions presented, except some deviations shown by MG5 aMC + HERWIG7 mainly in the higher region of the distribution. If MPI is not included in the simulation, predictions underestimate the differential cross section and fail to describe the shape of all the observables for Z+ ≥2 jet events with deviations up to  and area-normalized distributions (right) as functions of ∆φ between the Z boson and the leading jet for Z+ ≥1 jet events. The uncertainties in the predictions are shown as coloured bands around the theoretical predictions including statistical, PDF, scale, and tune uncertainties for the NLO MG5 aMC + PYTHIA8 (with CP5 tune) and the statistical uncertainty only for the LO MG5 aMC + PYTHIA8 (with CP5 tune), NLO MG5 aMC + PYTHIA8 (with CDPSTP8S1-WJ tune, CP5 tune with MPI-OFF), NLO MG5 aMC + HERWIG7 (with tune CH3), and SHERPA predictions. In the top panel, the vertical bars on the data points represent statistical uncertainties, whereas in the bottom panels, the total uncertainty in data is indicated by the solid yellow band centred at 1. In the legend, the χ 2 per degree of freedom is given to quantify the goodness of fit of the model to the data.  Figure 3: Differential cross sections (left) and area-normalized distributions (right) as functions of the p T imbalance between the Z boson and the leading jet for Z+ ≥1 jet events. The uncertainties in the predictions are shown as coloured bands around the theoretical predictions including statistical, PDF, scale, and tune uncertainties for the NLO MG5 aMC + PYTHIA8 (with CP5 tune) and the statistical uncertainty only for the LO MG5 aMC + PYTHIA8 (with CP5 tune), NLO MG5 aMC + PYTHIA8 (with CDPSTP8S1-WJ tune, CP5 tune with MPI-OFF), NLO MG5 aMC + HERWIG7 (with tune CH3), and SHERPA predictions. In the top panel, the vertical bars on the data points represent statistical uncertainties, whereas in the bottom panels, the total uncertainty in data is indicated by the solid yellow band centred at 1. In the legend, the χ 2 per degree of freedom is given to quantify the goodness of fit of the model to the data.   Figure 4: Differential cross sections (left) and area-normalized distributions (right) as functions of ∆φ between the Z boson and the dijet for Z+ ≥2 jet events. The uncertainties in the predictions are shown as coloured bands around the theoretical predictions including statistical, PDF, scale and tune uncertainties for the NLO MG5 aMC + PYTHIA8 (with CP5 tune) and the statistical uncertainty only for the LO MG5 aMC + PYTHIA8 (with CP5 tune), NLO MG5 aMC + PYTHIA8 (with CDPSTP8S1-WJ tune, CP5 tune with MPI-OFF), NLO MG5 aMC + HERWIG7 (with tune CH3), and SHERPA predictions. In the top panel, the vertical bars on the data points represent statistical uncertainties, whereas in the bottom panels, the total uncertainty in data is indicated by the solid yellow band centred at 1. In the legend, the χ 2 per degree of freedom is given to quantify the goodness of fit of the model to the data.
The model parameters for different MPI and hadronization models are mostly derived by fitting minimum bias and low-p T ( 3 GeV) MPI measurements. From the above data, it is clear that our measurements can test the accuracy of various predictions and the simulation of MPI. Most of the predictions follow a similar trend describing the differential cross sections and area-normalized distributions with a few exceptions.
The PYTHIA8 CP5 tune (with MPI) describes the differential cross section measurements within the uncertainty, but deviate (up to 15-20%) from the measurements in the lower-region of ∆ rel p T (Z, j 1 ) and ∆ rel p T (Z, dijet), where SPS is expected to be dominant. In the case of areanormalized distributions, the prediction from the MG5 aMC + PYTHIA8 CP5 tune (with MPI) describes the shape of the ∆ rel p T (j 1 , j 2 ) distribution within the uncertainty, but overestimates in the lower-region of ∆ rel p T (Z, j 1 ) and ∆ rel p T (Z, dijet), but underestimate otherwise. The LO calculations with MG5 aMC + PYTHIA8 provide a similar agreement as obtained by the NLO calculation. The MG5 aMC + HERWIG7 prediction describes the measurements within the uncertainty except some deviations in describing the shapes of the ∆φ(Z, j 1 ) and ∆ rel p T (j 1 , j 2 ) distributions. The SHERPA prediction is in reasonable agreement with the measurements within the uncertainty.  Figure 5: Differential cross sections (left) and area-normalized distributions (right) as functions of the p T imbalance between the Z boson and the dijet for Z+ ≥2 jet events. The uncertainties in the predictions are shown as coloured bands around the theoretical predictions including statistical, PDF, scale, and tune uncertainties for the NLO MG5 aMC + PYTHIA8 (with CP5 tune) and the statistical uncertainty only for the LO MG5 aMC + PYTHIA8 (with CP5 tune), NLO MG5 aMC + PYTHIA8 (with CDPSTP8S1-WJ tune, CP5 tune with MPI-OFF), NLO MG5 aMC + HERWIG7 (with tune CH3), and SHERPA predictions. In the top panel, the vertical bars on the data points represent statistical uncertainties, whereas in the bottom panels, the total uncertainty in data is indicated by the solid yellow band centred at 1. In the legend, the χ 2 per degree of freedom is given to quantify the goodness of fit of the model to the data.  Figure 6: Differential cross sections (left) and area-normalized distributions (right) as functions of the p T imbalance between leading and subleading jets for Z+ ≥2 jet events. The uncertainties in the predictions are shown as coloured bands around the theoretical predictions including statistical, PDF, scale, and tune uncertainties for the NLO MG5 aMC + PYTHIA8 (with CP5 tune) and the statistical uncertainty only for the LO MG5 aMC + PYTHIA8 (with CP5 tune), NLO MG5 aMC + PYTHIA8 (with CDPSTP8S1-WJ tune, CP5 tune with MPI-OFF), NLO MG5 aMC + HERWIG7 (with tune CH3), and SHERPA predictions. In the top panel, the vertical bars on the data points represent statistical uncertainties, whereas in the bottom panels, the total uncertainty in data is indicated by the solid yellow band centred at 1. In the legend, the χ 2 per degree of freedom is given to quantify the goodness of fit of the model to the data.
The predictions based on the DPS-specific tune CDPSTP8S1-WJ describe the shape of the distributions within the uncertainty. Since the parameters were derived by fitting only 7 TeV measurements, this suggests that the collision energy dependence of the MPI parameters is well modelled in this tune.

Summary
The CMS Collaboration has measured the differential cross sections for Z+ ≥1 jet and Z+ ≥2 jet events using proton-proton collision data at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 collected in the year 2016. The Z boson is reconstructed using the dimuon channel. This is the first measurement performed to explore observables sensitive to the presence of multi-parton interaction (MPI) using the Z+jets process at 13 TeV. Within the fiducial region, the production cross sections of Z+ ≥1 jet and Z+ ≥2 jet events are measured to be 158.5 ± 0.3 (stat) ± 7.0 (syst) ± 1.2 (theo) ± 4.0 (lumi) pb and 44.8 ± 0.4 (stat) ± 3.7 (syst) ± 0.5 (theo) ± 1.1 (lumi) pb, respectively. The measured integrated cross section in the fiducial region with jets is well described by the event generators SHERPA, MG5 aMC + PYTHIA8, and MG5 aMC + HERWIG7 predictions. The prediction obtained with MG5 aMC + PYTHIA8 with the double-parton scattering (DPS) specific tune CDPSTP8S1-WJ overestimates the measurements by 10-15%, but correctly describes the shape of all the observables. The prediction from MG5 aMC + PYTHIA8 with the CP5 tune, derived by fitting soft quantum chromodynamics (QCD) measurements, describes the differential cross section and area-normalized distributions. However, there are parts of the distributions that are not well described, such as single parton scattering that dominates lower regions of transverse momentum imbalance ∆ rel p T distributions. Predictions with other MPI models describe the measurements well (SHERPA) or reasonably well (MG5 aMC + HERWIG7) except a few deviations in describing the shapes of the ∆φ(Z, j 1 ) and ∆ rel p T (j 1 , j 2 ) distributions. The measured distributions show a significant sensitivity to MPI. A proper simulation of MPI is essential to describe the shape of the measured distributions and hence these results are a useful input to further improve DPS-specific tunes and a global tune in combination with other soft QCD measurements.

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 centres and personnel of the Worldwide LHC Computing Grid and other centres 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   [4] CMS Collaboration, "Measurements of jet multiplicity and differential production cross sections of Z+jets events in proton-proton collisions at √ s = 7 TeV", Phys. Rev. D 91 (2015) 052008, doi:10.1103/PhysRevD.91.052008, arXiv:1408.3104. [6] LHCb Collaboration, "Study of forward Z + jet production in pp collisions at √ s = 7 TeV", JHEP 01 (2014) 033, doi:10.1007/JHEP01(2014)033, arXiv:1310.8197.