Study of Z boson 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 s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV, corresponding to an integrated luminosity of 35.9 fb−1 collected in the year 2016. Differential cross sections of Z+ ≥1 jet and Z+ ≥2 jets are measured with transverse momentum of the jets above 20 GeV and pseudorapidity |η| < 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,   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].
All samples, except those based on sherpa, use pythia8 to model the initial-and finalstate parton showers and hadronization with the CUETP8M1 [24] or CUETP8M2T4 [46] tune. The CUETP8M1 tune includes the NNPDF 2.3 [47] LO PDF set with the strong coupling α S (m Z ) set to 0.1365 for space-and time-like shower simulation. The CUETP8M2T4 tune is based on the CUETP8M1 tune, which includes the NNPDF30_lo_as_0130 PDF set, but uses a lower value of α S = 0.1108 for the initial-state radiation component of the parton shower. Matching between the ME generators and the parton shower is done using the k T -MLM scheme [48,49] with the matching scale set at 19 GeV for the LO MG5_aMC sample, and the FxFx [50] scheme with the matching scale set to 30 GeV for the NLO MG5_aMC events. Generated events are processed through a full Geant4-based [51] CMS detector simulation and trigger emulation. The simulated samples include the effects of multiple interactions in each bunch crossing, referred to as pileup. The simulated events are reconstructed with the same algorithm used for the data.
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 her-wig7, applying 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: (4.1) -5 -

JHEP10(2021)176
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 -6 -JHEP10(2021)176 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 Object Selections 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.
• 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 production 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 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.

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 -10 -JHEP10(2021)176  Table 3. Uncertainty sources and their effect on the differential cross section distributions.  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, figure 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 figure 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 figure 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 50%.
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  Figure 2. Differential cross sections (left) 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. 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 area-normalized 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 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.
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 gen-  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. erators 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.  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.  Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.