Azimuthal correlations in Z +jets events in proton–proton collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s} = 13\,\text {Te}\hspace{-.08em}\text {V} $$\end{document}s=13TeV

The production of Z bosons associated with jets is measured in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {p}\text {p}$$\end{document}pp collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13\,\text {Te}\hspace{-.08em}\text {V} $$\end{document}s=13TeV with data recorded with the CMS experiment at the LHC corresponding to an integrated luminosity of 36.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {fb}^{-1}$$\end{document}fb-1. The multiplicity of jets with transverse momentum \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\textrm{T}} > 30\,\text {Ge}\hspace{-.08em}\text {V} $$\end{document}pT>30GeV is measured for different regions of the Z boson’s \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\textrm{T}} (\text {Z })$$\end{document}pT(Z), from lower than 10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {Ge}\hspace{-.08em}\text {V}$$\end{document}GeV to higher than 100\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {Ge}\hspace{-.08em}\text {V}$$\end{document}GeV. The azimuthal correlation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta \phi $$\end{document}Δϕ between the Z boson and the leading jet, as well as the correlations between the two leading jets are measured in three regions of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\textrm{T}} (\text {Z })$$\end{document}pT(Z). The measurements are compared with several predictions at leading and next-to-leading orders, interfaced with parton showers. Predictions based on transverse-momentum dependent parton distributions and corresponding parton showers give a good description of the measurement in the regions where multiple parton interactions and higher jet multiplicities are not important. The effects of multiple parton interactions are shown to be important to correctly describe the measured spectra in the low \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\textrm{T}} (\text {Z })$$\end{document}pT(Z) regions.


Introduction
In high-energy proton-proton (pp) collisions at the CERN LHC, the production of Z bosons is regarded as a standard measurement tool, because their properties can be measured very precisely in their leptonic decay channel, and the production cross section can be calculated with high precision.Although the production of Z bosons is a purely electroweak (EW) process, corrections from quantum chromodynamics (QCD) play an increasingly important role as the Z boson transverse momentum p T (Z) increases.At small p T (Z), where soft-gluon radiation is important, a resummation to all orders must be performed in order to obtain stable theoretical predictions [1][2][3][4] and to describe the measurements [5].When p T (Z) increases, hard partonic radiation becomes important and associated jets can be measured, allowing the study of QCD contributions to Z production.
This article describes a study by the CMS Collaboration of the production of Z bosons with associated jets at a center-of-mass energy of 13 TeV.We measure the multiplicity of jets with p T > 30 GeV in a pseudorapidity range of |η| < 2.4.In the region of low p T (Z), additional jets must balance the leading jet of p T > 30 GeV, whereas at large p T (Z) the Z boson is expected to balance the p T of the leading jet.We measure distributions in three (representative) p T (Z) regions: at low transverse momentum p T (Z) < 10 GeV; in the intermediate range of 30 < p T (Z) < 50 GeV; and in the large range of p T (Z) > 100 GeV.
The jet multiplicity, the azimuthal correlation ∆ϕ(Zj 1 ) between the Z boson and the leading jet, as well as the correlation ∆ϕ(j 1 j 2 ) between the two leading jets, is measured in these three ranges of p T (Z).At small p T (Z), a weak correlation between the Z boson and the leading jet is expected, whereas at large p T (Z) the azimuthal correlation is expected to be strong, since then the Z boson and the leading jet are most likely the highest p T objects in the event.The situation is opposite for ∆ϕ(j 1 j 2 ), where at small p T (Z) a strong correlation is expected, whereas at large p T (Z) the correlation will be weak.
The measurement of jet multiplicity as well as the measurements of the azimuthal correlations ∆ϕ(Zj 1 ), and ∆ϕ(j 1 j 2 ) in various ranges of p T (Z) provide an opportunity to make detailed comparisons with theoretical predictions.In particular, calculations of next-to-leading order (NLO) Z +jet production supplemented with parton shower (PS) and hadronization, as well as merged calculations with higher partonic jet multiplicity, can be studied.Of particular interest are the comparisons with predictions based on the parton branching (PB) method with transversemomentum dependent (PB-TMD) parton distribution functions (PDFs) [18][19][20] together with a TMD-based PS [21].A comparison with resummed calculations using the GENEVA [22][23][24][25] framework is also shown.

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 η coverage provided by the barrel and endcap detectors.Muons are de-tected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
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 around 100 kHz within a fixed latency of about 4 µs [26].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 [27].
The particle-flow algorithm (PF) [28] reconstructs and identifies each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector.The primary vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Section 9.4.1 of Ref. [29].
The energy of photons is obtained from the ECAL measurement.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.
The momentum resolution for electrons with p T ≈ 45 GeV from Z → e + e − decays ranges from 1.7-4.5%.It is generally better in the barrel region than in the endcaps, and also depends on the bremsstrahlung energy emitted by the electron as it traverses the material in front of the ECAL [30,31].The overall reconstruction efficiency is around 93% for electrons from Z decay.
The energy of muons is obtained from the curvature of the corresponding track.Muons are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers.The single-muon trigger efficiency exceeds 90% over the full η range, and the efficiency to reconstruct and identify muons is greater than 96%.Matching muons to tracks measured in the silicon tracker results in a relative p T resolution of 1% in the barrel and 3% in the endcaps for muons with p T up to 100 GeV.The p T resolution in the barrel is better than 7% for muons with p T up to 1 TeV [32].
The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response function of the calorimeters to hadronic showers.Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.
For each event, hadronic jets are clustered from these reconstructed particles using the infrared and collinear safe anti-k T algorithm [33,34] with a distance parameter of 0.4.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 whole p T spectrum and detector acceptance.Additional pp interactions within the same or nearby bunch crossings (pileup) contribute additional tracks and calorimetric energy depositions, increasing the apparent jet momentum.To mitigate this effect, tracks identified as originating from pileup vertices are discarded and a correction is applied to correct for any remaining contributions.
Jet energy corrections are derived from simulation studies so the average measured energy of jets becomes identical to that of particle level jets.In situ measurements of the momentum balance in dijet, photon+jet, Z+jet, and multijet events are used to determine any residual differences between the jet energy scale in data and in simulation, and appropriate corrections are made [35].Additional selection criteria are applied to each jet to remove jets potentially dominated by instrumental effects or reconstruction failures.The jet energy resolution amounts typically to 15-20% at 30 GeV, 10% at 100 GeV, and 5% at 1 TeV [35].
During the 2016 data-taking, a gradual shift in the timing of the inputs of the ECAL first-level trigger in the region at |η| > 2.0, referred to as prefiring, caused a specific trigger inefficiency.
For events containing an electron (a jet) with p T larger than 50 (100) GeV, in the region 2.5 < |η| < 3.0 the efficiency loss is about 10-20%, depending on p T , η, and timing.Correction factors were computed from data and applied to the acceptance evaluated by simulation.
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. [36].

Theoretical predictions
The measured differential cross sections are compared with a variety of predictions.One of the NLO calculations uses MADGRAPH5 aMC@NLO [37] (version 2.2.2) event generator interfaced with PYTHIA8 [38] for PS and hadronization.The matrix element calculations include Z/γ * +0,1,2 jets at NLO.It is labeled as MG5 AMC+PY8 (≤ 2j NLO) in the following.
The measurements are also compared with predictions obtained from MADGRAPH5 aMC@NLO (version 2.6.9) with PB-TMD PDFs and the corresponding PS as implemented in CAS-CADE3 [21] (labeled as MG5 AMC+CA3).The matrix elements (MEs) are calculated at NLO for Z+1 and Z+2 partons separately.The parton density PB-NLO-HERAI+II-2018-SET2 [20], as well as the PB initial-state PS, follow angular ordering conditions [39][40][41][42].The advantage of the MG5 AMC+CA3 calculation is that the parameters of the PB-TMD initial-state PS are fixed by the PB-TMD PDFs.
In all calculations using MADGRAPH5 aMC@NLO the renormalization and factorization scales are set to µ r = µ f = 1/2 ∑ i H T i , where H T i is the scalar sum of the p T with i running over all final particles and partons in the ME calculation.The corresponding uncertainties are estimated as the envelope of the set of variations of µ r and µ f by factors of 2 and 1/2, in all possible combinations except the extreme cases (µ f , µ r ) = (2, 0.5), (0.5, 2).The PDF uncertainties are estimated as the standard deviation of observables when using weights from the replicas provided in the NNPDF 3.0 NLO [43] PDF set.
The corresponding versions of these generators at leading-order (LO) are also compared with the measurement.
The following calculations are used for comparison with the measurements (a summary is given in Table 1): • MG5 AMC+PY8 (≤ 2j NLO) is a fixed-order perturbative QCD calculation at NLO of up to 2 noncollinear high-p T partons for pp → Z+N, N = 0, 1, 2, supplemented with PS and multiparton interactions (MPIs) from PYTHIA8 (version 8.212).The parameters of the underlying event tune CUETP8M1 [44] are applied.The merging of PS and MEs is performed with the FxFx scheme [45] with the merging scale of 30 GeV and a minimal partonic p T for jets of p part T = 15 GeV.The NNPDF 3.0 NLO PDFs are used and α S (m Z ) = 0.118 is chosen, where m Z is the Z boson mass.The predictions from MG5 AMC+PY8 (≤ 2j NLO) are used to investigate the effect of MPI.
• MG5 AMC+PY8 (≤ 4j LO) includes MEs computed at LO for pp → Z+N partons, N = 0, 1 . . .4, using the k T -MLM [46] procedure to match the different parton multiplicities of the MEs to the PS, with the matching scale set to 19 GeV.The PYTHIA8 generator (version 8.212) is interfaced with MG5 AMC to include initial-and finalstate PS and hadronization, with settings defined by the CUETP8M1 tune [44].The Table 1: Description of the simulated samples used in the analysis.
• MG5 AMC+CA3 (Z+2 NLO) is a fixed-order perturbative QCD calculation at NLO of two noncollinear high-p T partons for pp → Z+2 with p part T > 15 GeV, supplemented with PB-TMD PDFs and parton showering and hadronization.The same PB-TMD distribution and PS as in MG5 AMC+CA3 (Z+1 NLO) is applied.
• MG5 AMC+CA3 (Z≤ 3j LO) uses MG5 AMC to generate Z+0,1,2,3 jet samples at LO with a partonic generation cut p part T > 15 GeV.The TMD merging [51] procedure for combining the TMD PS with the ME calculations is used.The same PB-TMD distributions and PS as in MG5 AMC+CA3 (Z+1 NLO) are applied.A merging scale value of 23 GeV is used, since it provides a smooth transition between ME and PS computations.An overall K-factor of 1.27 is applied to the prediction.MPI effects are not simulated.
The calculation uses the PDF4LHC15 NNLO PDF set [52] with α S (m Z ) = 0.118.The simulation of PS, hadronization and MPI is performed by PYTHIA8 (version 8.212) with the CUETP8M1 tune.

Simulated samples
Events generated by MG5 AMC+PY8 (≤ 2j NLO) are passed through a full detector simulation based on GEANT4 [55].The simulated events are reconstructed using standard CMS reconstruction packages.This sample is used for the simulation of the signal process to estimate efficiencies, systematic uncertainties and for the correction of the data for detector spreading effects and inefficiencies, the so-called unfolding procedure.
Other processes that can give a final state with two oppositely charged same-flavor leptons and jets are tt, single top, vector boson pair (VV) and W+jets.The tt and single top backgrounds are generated using POWHEG 2.0 [56][57][58][59][60][61] interfaced with PYTHIA8.The total cross section of tt production is normalized to the prediction with NNLO accuracy in QCD and next-to-nextto-leading logarithmic (NNLL) gluon radiation resummation calculated with TOP++ 2.0 [62].
The double vector boson productions are generated with MG5 AMC (WZ), POWHEG (WW), both interfaced to PYTHIA8, or with PYTHIA8 for ZZ.The total cross sections for the WZ and ZZ diboson samples are normalized to the NLO prediction calculated with MCFM 6.6 [63].The W+jets sample is generated by MG5 AMC at NLO accuracy, interfaced with PYTHIA8.The Z boson decay into τ + τ − is included in the signal simulation and considered as a background.

Data analysis
The differential cross section of Z bosons with associated jets is measured in bins of p T (Z), as functions of the jet multiplicity, the azimuthal angles ∆ϕ(Zj 1 ) and ∆ϕ(j 1 j 2 ), where j 1 is the leading jet and j 2 is the second-leading jet.

Event selection
The data samples recorded in 2016 correspond to an integrated luminosity of 36.3 fb −1 .Events with a pair of leptons (µ + µ − or e + e − ) consistent with the decay of a Z boson and with jets reconstructed from PF candidates are selected as Z+jet events.Those events are required to pass a series of selection criteria to reduce the background contributions.An event is selected if the double muon (electron) trigger with 18 and 7 (23 and 12) GeV thresholds in p T or a single muon trigger with a threshold of 24 GeV is satisfied.In the offline selection, the leading (subleading) electron and muon candidates must have transverse momenta of p T > 25 (20) GeV in a range of |η| < 2.4.Only events with pairs of oppositely charged muons (electrons) with an invariant mass in the range 91 ± 15 GeV are accepted.
Muon candidates are required to be isolated from other particles, as specified by an isolation criteria, I ISO : where the sums run over the corresponding particles inside a cone of radius ∆R = √ (∆η) 2 + (∆ϕ) 2 = 0.4 around the muon candidate considering separately charged hadrons (charged), neutral hadrons (neutral), photons (EM), and charged particles from pileup (PU).
Electrons are required to be isolated from other particles, as specified by an isolation criteria, I ISO : where the sums run over the corresponding particles inside a cone of radius ∆R = 0.3.The term ρA eff represents a correction for pileup effects, where ρ corresponds to the amount of p T added to the event per unit area and A eff is the area of the isolation region weighted by a factor that accounts for the dependence of the pileup transverse energy density on the electron η [31,64].
Jets are required to have a minimum p T of 30 GeV to ensure that they are well measured and to reduce the pileup contamination.Jets are limited to a rapidity range of |y| < 2.4, and are required to be isolated from the lepton candidates by ∆R ℓ,j > 0.4.To keep only charged particles originating from the Z boson vertex, charged particles identified as originating from pileup vertices are discarded.As discussed in Section 2, jet energy corrections are applied to data and simulation.The jet energy resolution (JER) in simulation is further spread to match that in data.
The simulated events are reweighted such that their pileup distribution matches the measured one in each data-taking period.
Several corrections for leptons are applied to the simulation yields to compensate for the measured differences between the efficiencies in data and simulation.These corrections are applied as trigger, lepton identification, and lepton isolation scale factors.The values of the scale factors are close to one.An additional trigger inefficiency correction due to the prefiring effect is included.The exclusive jet multiplicity in different regions of p T (Z) for muon and electron channels is shown in Fig. 1.

Correction for the detector effects
Detector effects, like inefficiencies and the spreading of the particle momentum, energy and angle, are corrected using the an unfolding procedure, which is applied after background subtraction.The iterative D'Agostini method as implemented in RooUnfold [65,66] is used.The iteration is affected by fluctuations that increase with the number of iterations.The fluctuations are studied for each distribution and the procedure of unfolding is stopped before the fluctuations become significant with respect to the statistical uncertainty, following the method used in [16].Through the unfolding procedure the cross section at the stable-particle level is obtained.Particles are considered stable if their proper lifetime is above 10 mm/c.Neutrinos are not included.The momentum of the leptons is calculated including photons in a cone of a radius of ∆R = 0.1 ("dressed" leptons).The phase space definition for the final cross sections is given in Table 2.

Background estimation
The contributions from background processes are estimated using MC-based simulations, described in Section 3.1, and are subtracted from the measured distributions.The dominant background, tt, is verified with data control samples, using the same criteria as for the measurement, but requiring the two leptons to have different flavours (eµ instead of µ + µ − or e + e − ).
The effect of mismodeling of top quark distributions is covered by the MC uncertainties.Therefore, no additional correction or uncertainty is applied [14].The Z → τ + τ − decays are considered as a background, and their contribution is estimated from simulation and subtracted during the unfolding procedure.

Uncertainties
The statistical uncertainties from the measured observables are propagated to the final results via the unfolding procedure.The systematic uncertainties originate from the following sources: • Jet energy scale: Variations of the jet energy scale corrections [35,67] are applied as functions of p T and η for individual run periods; this affects the differential cross sections by 3-7%.
• Jet energy resolution: The JER [35,67] uncertainty is obtained by varying the spreading factor to match the simulated jet energy resolution to data by one standard deviation around its central value, resulting in an uncertainty of up to 1-2%.
• Efficiency correction: The uncertainty coming from the measurements of trigger efficiency, lepton reconstruction, and lepton identification is estimated by varying the scale factors by their uncertainties, as described in Ref. [5].The resulting uncertainty in the differential cross section measurement is less than 1%.
• Luminosity: The uncertainty in the integrated luminosity is 1.2% [68].It is applied as a global scale factor to the cross section as well as to the normalization of the background samples.
• Pileup: The determination of the simulated pileup profile is based on a total inelastic pp cross section of 69.2 mb [69].Alternative pileup profiles are generated by varying this cross section by 5% affecting the measurement by 1-2%.
• Prefiring: The prefiring uncertainty is estimated by up and down variations of the prefiring weight.The uncertainty is less than 0.5%.
• Background: The theoretical uncertainty in the calculation of background processes is used to estimate the uncertainty in the background modeling.The main source of background is tt production on this process.An uncertainty of around 6% is estimated using the TOP++2.0program, which includes scale and PDF variations.The resulting uncertainty is less than 0.2%.Systematic uncertainties stemming from other background processes are negligible.
• Unfolding and model: The uncertainty of unfolding and modeling is estimated by reweighting the simulated signal event sample to match the data and using this as an alternative model for unfolding.This gives an uncertainty of about 2%.The uncertainty coming from the finite size of the simulation sample that is used to correct the data for detector effects results in an uncertainty of 2-8%.
A summary table of these uncertainties is given in Table 3.All the systematic uncertainties are quadratically summed assuming independent uncertainty sources.

Results
The production cross section of Z+jets is measured in the phase space given in Table 2.The Z boson is identified via its leptonic decay channel.The results of the muon and electron decay In Fig. 2, the exclusive jet multiplicity is shown for three different ranges of p T (Z).At low p T (Z), the majority of events have no jet with p T > 30 GeV and only about 1% of the events have one or more jets.This suggests that the p T (Z) is mainly compensated by softer (p T < 30 GeV) radiation at low p T (Z).Events with higher jet multiplicity indicate that the dominant hard process is essentially a jet production process, and the Z boson is radiated from a quark, as an electroweak correction to a pure QCD process.At high p T (Z), the majority of events have at least one jet with a tail towards higher jet multiplicities, which indicates that the hardest process is indeed Z+jet, and additional jets originate from higher-order QCD corrections.
In Fig. 2, the measurement is compared with the generator MG5 AMC+PY8 (≤ 2j NLO) with and without multiparton interactions.The MPI contribution is important in the low p T (Z) region, but also at higher p T (Z) and higher jet multiplicities MPI plays a role.The prediction of MG5 AMC+PY8 (≤ 2j NLO) including MPI agrees with the measurement, even for high jet multiplicities.This behaviour is consistent with the prediction of the dependence of MPI effects on event kinematics from the p T (Z) reported in [72].
In Fig. 3, a comparison of the measurement with predictions from MG5 AMC+CA3 (Z+1 NLO), MG5 AMC+CA3 (Z+2 NLO) and GENEVA (Z+0 NNLO) is shown.Both MG5 AMC+CA3 (Z+1 NLO) and MG5 AMC+CA3 (Z+2 NLO) predictions are multiplied by a factor 1.2 to account for the normalization of PB TMD set 2 (as discussed in Section 3).For p T (Z) > 30 GeV the Z+1 (Z+2) predictions describe well the one (two) jet multiplicities, whereas at higher multiplicities a deviation from these measurement is observed, which can be attributed to the missing MPI contributions (as shown in Fig. 2).The GENEVA (Z+0 NNLO) predictions, which include MPI, are in agreement for low jet multiplicities for low p T (Z), whereas higher jet multiplicities are not well described because of missing higher order contributions in the ME calculations.
In Fig. 4, the measurement is compared with predictions from MG5 AMC+PY8 (≤ 4j LO) and MG5 AMC+CA3 (Z≤ 3j LO).The prediction from MG5 AMC+PY8 (≤ 4j LO) describes the measurements in all p T (Z) regions.The MG5 AMC+CA3 (Z≤ 3j LO) prediction agrees with the measurements in all p T (Z) ranges, except in the second bin at low p T (Z) values where MPI plays a significant role.
In Fig. 5, the azimuthal correlation, ∆ϕ(Zj 1 ), between the Z boson and the leading jet is shown for three different ranges of p T (Z).In the range p T (Z) < 10 GeV, the Z boson is only very weakly correlated with the leading jet, thus the distribution is almost uniform.In the region p T (Z) > 100 GeV, the Z boson is highly correlated with the leading jet and the cross section       falls more than two orders of magnitude from the back-to-back region to the small ∆ϕ(Zj 1 ) region.The systematic uncertainty in the low p T (Z) range is O(10%).In Fig. 5, the predictions from MG5 AMC+PY8 (≤ 2j NLO) with and without MPI are compared with the measurement.In the low p T (Z) range, MPI contributes about 40%, and even in the region of 30 < p T (Z) < 50 GeV, the contribution from MPI could be about 20% in the small-∆ϕ(Zj 1 ) region.The prediction of MG5 AMC+PY8 (≤ 2j NLO) including MPI describes the measurements.
In Fig. 6, the measurement is compared with MG5 AMC+CA3 (Z+1 NLO), MG5 AMC+CA3 (Z+2 NLO) and GENEVA (Z+0 NNLO).In low p T (Z) range, the MG5 AMC+CA3 (Z+1 NLO) and MG5 AMC+CA3 (Z+2 NLO) predictions differ from the measurements due to the missing contribution of MPI.In the high p T (Z) region the predictions agree better with the measurements (the region ∆ϕ(Zj 1 ) → π is not accessible in the Z+2 calculation).The GENEVA (Z+0 NNLO) prediction agrees with the measurement at low p T (Z), whereas at larger p T (Z) the prediction differs from the measurement because of missing higher order contributions.
In Fig. 7, the predictions from MG5 AMC+PY8 (≤ 4j LO) and MG5 AMC+CA3 (Z≤ 3j LO) are compared with the measurement.The MG5 AMC+PY8 (≤ 4j LO) is in agreement with the measurement.The MG5 AMC+CA3 (Z≤ 3j LO) prediction is too low in the low p T (Z) region, due to the missing MPI contribution, whereas other p T (Z) ranges are described.
In Fig. 8, the azimuthal correlation ∆ϕ(j 1 j 2 ) between the two leading jets is shown for the three different ranges of p T (Z).A strong correlation between the two leading jets is observed at small p T (Z), whereas only a weak correlation is seen at large p T (Z).This indicates that at low p T (Z) the process is dominated by a jet production process and that the Z boson is radiated from a quark (EW correction) and therefore the jets are correlated.On the contrary, at large p T (Z) the process is dominated by Z+jet production, with higher-order QCD corrections in form of additional jets, which are only weakly correlated.The measurement is compared with MG5 AMC+PY8 (≤ 2j NLO), with and without MPI.Except for the highest p T (Z) region, the contribution from MPI is significant, especially in the small ∆ϕ(j 1 j 2 ) range.The prediction obtained with MG5 AMC+PY8 (≤ 2j NLO) including MPI describes the measurement well over the whole range.
The contribution from MPI is significant in the low p T (Z) regions and becomes negligible when p T (Z) > 100 GeV.The calculation MG5 AMC+PY8 (≤ 2j NLO) describes the measurements within the scale uncertainties, if an appropriate tune for PS and underlying event parameters are applied (here, the CUETP8M1 tune).
The predictions of MG5 AMC+CA3 using PB-TMD PDFs and initial-state PB PS come remarkably close to the measurements in regions of phase space where they are applicable.The predic-                tion of GENEVA (Z+0 NNLO) for Z+jet observables describe the measurements in the regions of low jet multiplicity but show differences when two or more jets are selected since higher order contributions are not included.The prediction of the merged LO calculations MG5 AMC+PY8 (≤ 4j LO) and MG5 AMC+CA3 (Z≤ 3j LO) describe the measurements quite well keeping in mind the role of MPI at low p T (Z).
In Ref. [73], calculations using TMDs and the "winner-takes-all" jet recombination scheme for the azimuthal angular decorrelation in Z+jet are described.The calculations reported here do not change significantly if the "winner-takes-all" recombination scheme in the anti-k T algorithm is applied.
In this paper the differential cross section measurements are presented in three representative regions of the Z boson transverse momentum, the results for intervals 10 < p T (Z) < 30 GeV and 50 < p T (Z) < 100 GeV can be also found in HEPData [74].

Summary
We have measured the Z+jet production cross section in proton-proton collisions at the LHC at a center-of-mass energy of 13 TeV.The associated jet multiplicity for various regions of the transverse momentum of the Z boson, p T (Z), was measured.At p T (Z) < 10 GeV only about 1% of the events have jets with p T > 30 GeV, with nonnegligible cross sections at high jet multiplicity.At 30 < p T (Z) < 50 GeV, most of the events have at least one jet, with a significant tail to higher jet multiplicities.The azimuthal angle ∆ϕ(Zj 1 ) between the Z boson and the leading jet, as well as the azimuthal angle ∆ϕ(j 1 j 2 ) between the two leading jets, was measured for the three p T (Z) regions.At low p T (Z), the Z boson is only loosely correlated with the jets, but the two leading jets are strongly correlated.At large p T (Z), the Z boson is highly correlated with the leading jet, but the two leading jets are only weakly correlated.
The measurement shows that at low p T (Z) the Z boson appears as an electroweak correction to high-p T jet production, whereas at large p T (Z) the dominant process is Z+jet production.
The predictions of MG5 AMC+CA3 (Z+1 NLO) and MG5 AMC+CA3 (Z+2 NLO) using the parton branching method with transverse-momentum dependent (PB-TMD) parton densities, which do not include MPI effects, and the corresponding PS agree with the measurements in the regions where MPI effects are negligible.The prediction from GENEVA (Z+0 NNLO) using matrix elements at next-to-next-to-leading order for Z production, supplemented with resummation, PS and MPI from PYTHIA8, agrees with the measurements in the low jet multiplicity region.
The leading order prediction of MG5 AMC+PY8 (≤ 4j LO), including merging of jet multiplicities, describes the measurements well.The prediction of MG5 AMC+CA3 (Z≤ 3j LO) using PB-TMD parton densities and PS with merging of jet multiplicities agrees well with the measurements in the regions where MPI is negligible.
In summary, Z+jet measurements challenge theoretical predictions; a good agreement can be achieved by including contributions of multiparton interactions, parton showering, parton densities, as well as multijet matrix element merging.The differential measurements provided here help to disentangle the various contributions and illustrate where each contribution be-perMicro Corporation; the Welch Foundation, contract C-1845; and the Weston Havens Foundation (USA).

Figure 1 :
Figure 1: The exclusive jet multiplicity distribution before unfolding in three different regions of p T (Z): p T (Z) < 10 GeV (upper), 30 < p T (Z) < 50 GeV (middle), p T (Z) > 100 GeV (lower) for the µ + µ − channel (left) and the e + e − channel (right).The error bars around the data points represent the statistical uncertainties.

Figure 2 :
Figure 2: Jet multiplicity in three different regions of p T (Z): p T (Z) < 10 GeV (upper left), 30 < p T (Z) < 50 GeV (upper right), p T (Z) > 100 GeV (lower).The error bars on the data points represent the statistical uncertainty of the measurement, and the hatched band shows the total statistical and systematic uncertainties added in quadrature.Predictions using MG5 AMC+PY8 (≤ 2j NLO) with and without MPI are shown.

Figure 4 :
Figure 4: Jet multiplicity in three different regions of p T (Z): p T (Z) < 10 GeV (upper left), 30 < p T (Z) < 50 GeV (upper right), p T (Z) > 100 GeV (lower).The error bars on the data points represent the statistical uncertainty of the measurement, and the hatched band shows the total statistical and systematic uncertainties added in quadrature.Predictions from MG5 AMC+PY8 (≤ 4j LO) and MG5 AMC+CA3 (Z≤ 3j LO) are shown.Different normalization factors are applied, as described in the text.

Figure 5 :
Figure 5: Cross section as a function of ∆ϕ(Zj 1 ) between the Z boson and the leading jet in the three p T (Z) bins: p T (Z) < 10 GeV (upper left), 30 < p T (Z) < 50 GeV (upper right), p T (Z) > 100 GeV (lower).The error bars on the data points represent the statistical uncertainty of the measurement, and the hatched band shows the total statistical and systematic uncertainties added in quadrature.Predictions using MG5 AMC+PY8 (≤ 2j NLO) with and without multiparton interactions are shown.

Figure 7 :
Figure 7: Cross section as a function of ∆ϕ(Zj 1 ) between the Z boson and the leading jet in three p T (Z) bins: p T (Z) < 10 GeV (upper left), 30 < p T (Z) < 50 GeV (upper right), p T (Z) > 100 GeV (lower).The error bars on the data points represent the statistical uncertainty of the measurement, and the hatched band shows the total statistical and systematic uncertainties added in quadrature.Predictions from MG5 AMC+PY8 (≤ 4j LO) and MG5 AMC+CA3 (Z≤ 3j LO) are shown.Different normalization factors are applied, as described in the text.

Figure 8 :
Figure8: Cross section as a function of ∆ϕ(j 1 j 2 ) between two leading jets in three p T (Z) regions: p T (Z) < 10 GeV (upper left), 30 < p T (Z) < 50 GeV (upper right), p T (Z) > 100 GeV (lower).The error bars on the data points represent the statistical uncertainty of the measurement, and the hatched band shows the total statistical and systematic uncertainties added in quadrature.Predictions using MG5 AMC+PY8 (≤ 2j NLO) with and without multiparton interactions are shown.

Figure 10 :
Figure 10: Cross section as a function of ∆ϕ(j 1 j 2 ) between two leading jets in three p T (Z) regions: p T (Z) < 10 GeV (upper left), 30 < p T (Z) < 50 GeV (upper right), p T (Z) > 100 GeV (lower).The error bars on the data points represent the statistical uncertainty of the measurement, and the hatched band shows the total statistical and systematic uncertainties added in quadrature.Predictions from MG5 AMC+PY8 (≤ 4j LO) and MG5 AMC+CA3 (Z≤ 3j LO) are shown.Different normalization factors are applied, as described in the text.

Table 3 :
Systematic uncertainties on the unfolded differential cross section