Cross-section measurements for the production of a $Z$ boson in association with high-transverse-momentum jets in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

Cross-section measurements for a $Z$ boson produced in association with high-transverse-momentum jets ($p_{\mathrm{T}} \geq 100$ GeV) and decaying into a charged-lepton pair ($e^+e^-,\mu^+\mu^-$) are presented. The measurements are performed using proton-proton collisions at $\sqrt{s}=13$ TeV corresponding to an integrated luminosity of $139$ fb$^{-1}$ collected by the ATLAS experiment at the LHC. Measurements of angular correlations between the $Z$ boson and the closest jet are performed in events with at least one jet with $p_{\mathrm{T}} \geq 500$ GeV. Event topologies of particular interest are the collinear emission of a $Z$ boson in dijet events and a boosted $Z$ boson recoiling against a jet. Fiducial cross sections are compared with state-of-the-art theoretical predictions. The data are found to agree with next-to-next-to-leading-order predictions by NNLOjet and with the next-to-leading-order multi-leg generators MadGraph5_aMC@NLO and Sherpa.


Introduction
The measurement of -boson production 1 in association with jets, + jets, constitutes a powerful test of perturbative quantum chromodynamics (QCD) [1,2] and, in the case of high-energy jets, it provides a way to probe the interplay between QCD and higher-order electroweak (EW) processes [3][4][5][6]. The large + jets production cross section and the easily identifiable decays of the boson to charged-lepton final states offer a clean experimental signature which can be measured precisely. Such processes also constitute non-negligible backgrounds in measurements of the Higgs boson [7,8] and in searches for new phenomena [9][10][11], which often exploit the presence of high-T jets to enrich a data sample with potential signal. In those studies, predictions are used to extrapolate + jets backgrounds from control regions to the signal regions and to model the distributions of the final discriminants.
In the calculations of + 1-jet production at leading order (LO), the boson recoils against a quark or a gluon. At next-to-leading order (NLO), real and virtual QCD and EW effects play a role in + jets production, such as in topologies corresponding to dĳet events where a real boson is emitted from an incoming or outgoing quark leg [3][4][5][6]. Example Feynman diagrams for LO and NLO + jets production processes are shown in Figure 1. The latter case can lead to production rate enhancements proportional to s ln 2 ( T, 1 / ), where s is the strong coupling constant, T, 1 the transverse momentum of the leading jet, and the mass of the boson, and thus the effect can become very large for events with high-T jets. These events exhibit a collinear enhancement in the distribution of the angular distance between the boson and the closest jet. Although the enhancement can be probed in the region of small angular separation, this region also contains contributions where the boson is produced in association with larger numbers of jets, which must be included in the predictions. The measurements presented in this paper target QCD-only + jets production, treating EW + 2-jets (EW ) production [12] as a background. Measurements where the EW contribution is treated as signal and not subtracted as background are also performed and published in the HEPData entry [13] of this measurement.
The ATLAS Collaboration [14] at the Large Hadron Collider (LHC) [15] first measured angular distributions in high-T boson production with jets ( + jets) in the 8 TeV -collision data set [16]. The first similar measurement in + jets events was published by the CMS Collaboration and used a partial 13 TeV data set corresponding to 35.9 fb −1 [17]. Both measurements highlight the fact that the collinear region, where the angular separation between the / boson and the closest jet is small, represents a major challenge for contemporary Monte Carlo (MC) generators. The measurements presented in this paper include a wide range of new observables sensitive to the presence of high-T jets and to the collinear emission of a boson in dĳet events. The statistical power of the full LHC Run-2 data set makes it possible to tighten the collinear selection, and to measure key observables separately for collinear and non-collinear topologies.
This publication focuses on events that contain a -boson candidate reconstructed from either an + − or + − pair in association with hadronic jets defined as jets having transverse momentum greater than or equal to 100 GeV. The phase-space region with at least one associated jet is labelled as the inclusive region. In this region, the measured quantities are the transverse momentum of the leading jet ( T, 1 ), the transverse momentum of the boson ( T,ℓℓ ), the scalar sum of the transverse momentum of all selected jets and leptons ( T ), and the jet multiplicity. A high-T region is selected by requiring the presence of a jet with T ≥ 500 GeV. To test the prediction that this region is composed of two characteristic topologies, the soft radiation of a boson from a jet (collinear topology) and the hard scatter of a boson against a jet (back-to-back topology), the high-T region is split to cover different ranges of the angle between g q Z q g gq Z q Figure 1: Representative Feynman diagrams for the production of a boson in association with high-T jets. The + 1-jet events (left) are expected to populate the back-to-back region where the boson is balanced against a single high-T jet. In dĳet events (right), the boson is expected to be radiated from the quark leg, with kinematics leading to small values of the angular distance between the boson and the closest jet, Δ min , , and therefore populating the collinear region. the boson and the closest jet, and selected key observables are measured separately for each region. In the high-T region, the following observables sensitive to the presence of collinear -boson emission are studied: • Δ min , , the angular distance 2 between the boson and the closest jet. Real -boson radiation is expected to be enhanced at low values of Δ min , . At large values of Δ min , , the boson is balanced by a recoiling jet and large virtual EW corrections are expected. To enrich these two topologies, collinear and back-to-back regions are constructed by requiring Δ min , ≤ 1.4 and Δ min , ≥ 2.0, respectively.
• , , the ratio of the -boson T to the closest-jet T , defined as , ≡ T,ℓℓ T (closest jet) .
Collinear -boson radiation is expected to be dominated by soft bosons, resulting in very small values for this ratio.
• jets , the jet multiplicity. The back-to-back region is expected to be dominated by + 1-jet events, whereas the collinear region would be dominated by + 2-jets events.
The measurements of jet multiplicity and , are performed both in the full high-T region and separately in the collinear and back-to-back regions.
Measurements of jet multiplicity and Δ min , are also performed in an alternative high-energy region, constructed by requiring the scalar sum of the transverse momentum of all selected jets, T , to be at least 600 GeV. This alternative region probes high-energy events but does not depend on the presence of a single very energetic object. In this region, called the high-T region, a large fraction of the events have higher jet multiplicity. 2 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe. The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). The angular distance is defined as Δ ≡ √︁ (Δ ) 2 + (Δ ) 2 . When dealing with massive jets and particles, the rapidity = (1/2) ln[( + )/( − )] is used, where is the jet/particle energy and is the -component of the jet/particle momentum.
Predictions from the most recent generators combine NLO multi-leg matrix elements (ME) with a parton shower (PS) and hadronisation model [18][19][20][21]. Fixed-order parton-level theoretical predictions for + jets production at next-to-next-to-leading order (NNLO) are available for up to one associated jet [22][23][24][25]. In this paper, the cross-section measurements are compared with state-of-the-art multi-leg ME+PS generators and NNLO fixed-order + jets predictions from NNLOjet [24,25]. Virtual EW corrections were made available recently [26,27] and are included in one of the Sherpa [19] predictions studied in this paper.
The paper is organised as follows. Section 2 contains a brief overview of the ATLAS detector. The data and simulated samples, as well as additional predictions used in the analysis, are described in Section 3. The object definition and the event reconstruction at detector level are presented in Section 4, while Section 5 describes the background modelling and presents a comparison of measured and predicted yields at detector level. After background subtraction, the data are unfolded to particle level in a fiducial phase space with a procedure described in Section 6. The experimental and theoretical systematic uncertainties are estimated in Section 7. Section 8 presents the unfolded cross-section results and the comparisons with predictions. Conclusions are provided in Section 9.

ATLAS detector
The ATLAS experiment at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4 coverage in solid angle. It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer. The inner tracking detector covers the pseudorapidity range | | < 2.5. It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors. Lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity. A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (| | < 1.7). The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to | | = 4.9. The muon spectrometer surrounds the calorimeters and is based on three large superconducting air-core toroidal magnets with eight coils each. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. The muon spectrometer includes a system of precision tracking chambers and fast detectors for triggering. A two-level trigger system is used to select events. The first-level trigger is implemented in hardware and uses a subset of the detector information to accept events at a rate below 100 kHz. This is followed by a software-based trigger that reduces the accepted event rate to 1 kHz on average depending on the data-taking conditions. An extensive software suite [28] is used in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

Data set and simulated event samples
The data used in this analysis were recorded with the ATLAS detector from 2015 to 2018 in collisions at √ = 13 TeV (full Run-2 data set) and correspond to a total integrated luminosity of 139 fb −1 [29]. The mean number of interactions per bunch crossing, including the hard scattering and other interactions in the same and neighbouring bunch crossings (pile-up), was ⟨ ⟩ = 34.
MC simulation samples are used to estimate most of the contributions from background events, to unfold the data to particle level, and in comparisons with the unfolded data distributions. The generated samples were  [41], the Hessian NNPDF3.0nnlo PDF set [42] is used, and an analytic enhancement technique has been introduced [19]. The cross section in the high-T region has been considerably reduced relative to the prediction from the previous Sherpa version by switching to an improved matching scheme with a different treatment of unordered histories [40]. Sherpa 2.2.11 is also the only sample used in this paper which includes NLO virtual EW corrections [34,35]. The samples were produced using three options for the combination of NLO EW and QCD corrections: an additive, a multiplicative, and an exponential scheme. The nominal prediction is derived via the additive scheme; a systematic uncertainty band is derived from the envelope of all schemes. In contrast to virtual EW corrections, EW parton showers are not included in any of the generators used in this paper. The Sherpa 2.2.11 + jets samples are used for the nominal unfolding of the data distributions, to estimate the systematic uncertainties and in comparisons with the cross-section measurements.
A second + jets sample, referred to as MG5_aMC+Py8 FxFx [19], was produced by using the MadGraph5_aMC@NLO 2.6.5 [43] program to generate matrix elements at NLO accuracy in QCD for up to three additional partons in the final state. The NNPDF3.1nnlo set [42] was used in the generation. The parton showering and subsequent hadronisation was performed using Pythia 8.240 [21] with the A14 tune [44] and the NNPDF2.3lo PDF set [45]. The jet multiplicities were merged using the FxFx prescription [20]. This prediction is compared with the unfolded cross-section measurements.
A third sample of + jets events and an event sample from + jets processes were produced with the Sherpa 2. A fourth + jets sample, referred to as MG5_aMC+Py8 CKKWL, was generated using LO-accurate matrix elements with up to four final-state partons calculated by MadGraph5_aMC@NLO 2.2.2 [43]. The ME calculation employed the NNPDF3.0nlo PDF set and was interfaced to Pythia 8.186 [46] for the modelling of the parton shower, hadronisation, and underlying event. The overlap between matrix element and parton shower emissions was removed using the CKKW-L merging procedure [47,48]. The A14 tune of Pythia was used with the NNPDF2.3lo PDF set. This sample is used to validate the unfolding method and in comparisons with the unfolded cross-section measurements.
Two additional + jets samples were generated with the NNLOjet program [24,25], which computes fixed-order parton-level predictions for inclusive jet processes at higher orders in QCD. The NLO and NNLO predictions, referred to as NNLOjet@NLO and NNLOjet@NNLO, respectively, were calculated as higher-order corrections to the parton-level LO process of + 1-jet production. The NNPDF3.1nnlo set was used with a central scale choice of 0 = 1 2 ( T, + Σ ∈partons T, ) with T, = √︃ 2 ℓℓ + 2 T, . These samples are pure QCD predictions at parton level. To match the fiducial selection of the measurement (see Section 6), scale factors to correct from the Born level to the dressed-lepton level are computed and applied to these predictions. The slightly different overlap-removal procedure for jets and leptons used in these samples, due to the NNLOjet program design, is addressed by overlap-removal correction scale factors. Both sets of scale factors, deviating from unity at the percent level and computed separately for each bin of the measured observables, are published in the HEPData entry [13] of this measurement. Non-pertubative corrections are found to be consistent with zero when T, 1 exceeds 100 GeV and are not needed to match the fiducial selection of these measurements. These samples are used in comparisons with the unfolded cross-section measurements.
The EW process is defined by the -channel exchange of a weak boson and at tree level is calculated at O( 4 ew ) when including the decay of the boson [12]. In contrast, the strong process, which is covered by the + jets samples, has no weak boson exchanged in the -channel and at tree level is calculated at O( 2 ew 2 s ) when including the decay of the boson. The EW samples were produced in the vector-boson fusion (VBF) approximation with Herwig 7.1.5 [49,50] at NLO accuracy in the strong coupling, using VBFNLO 3.0.0 [51] to provide the loop amplitude. The MMHT2014lo PDF set [56] was used along with the default set of tuned parameters for parton showering, hadronisation and the underlying event. To account for the interference between strong and EW processes, a uniform modelling uncertainty of 25% in the EW cross section (40% in the collinear region), determined from simulation with MadGraph5_aMC@NLO 2.9.5, is applied [12].
The¯background in this measurement is derived with a data-driven method as described in Section 5. The MC¯events used for intermediate steps of the method were modelled using the Powheg Box v2 [52][53][54][55] generator at NLO with the NNPDF3.0nlo PDF set and the ℎ damp parameter 3 set to 1.5 top [57]. The events were interfaced to Pythia 8. 230 [21] to model the parton shower, hadronisation, and underlying event, with parameters set according to the A14 tune and using the NNPDF2.3lo set of PDFs. Thes ample is normalised to the cross-section prediction at NNLO accuracy, including the resummation of next-to-next-to-leading logarithmic (NNLL) soft-gluon terms calculated with Top++ 2.0 [58][59][60][61][62][63][64].
Single top quark production in the -channel, in the -channel, and in association with a boson ( ) was modelled using the Powheg Box v2 generator at NLO in QCD with the five-flavour scheme and the NNPDF3.0nlo set of PDFs. The diagram-removal scheme [65] was used to remove interference and overlap with¯production. The cross section is corrected to the theory prediction at approximate NNLO accuracy [66,67], while the -and -channel cross sections are corrected to the prediction at NLO accuracy [68,69].
Samples of diboson final states ( ) were produced with the Sherpa 2.2.1 or Sherpa 2.2.2 generator depending on the process, including off-shell effects and Higgs boson contributions where appropriate. Fully leptonic final states and semileptonic final states, where one boson decays leptonically and the other hadronically, were generated using matrix elements at NLO accuracy in QCD for up to one additional parton and at LO accuracy for up to three additional parton emissions. The matrix element calculations were matched and merged with the Sherpa parton shower as detailed above for Sherpa 2.2.1, and the NNPDF3.0nnlo set of PDFs was used.
The production of + final states was simulated with the Sherpa 2.2.8 [18] generator. Matrix elements at NLO QCD accuracy for up to one additional parton and LO accuracy for up to three additional parton emissions were matched and merged with the Sherpa parton shower as detailed above for Sherpa 2.2.1, and the NNPDF3.0nnlo set of PDFs was used.
Background events involving semileptonic decays of heavy quarks, hadrons misidentified as leptons, and, in the case of the electron channel, electrons from photon conversions are referred to collectively as 'multĳet events'. The multĳet background is estimated using data-driven techniques, as described in Section 5.
For bottom and charm hadron decays, the EvtGen 1.7.0 program [70] was used for MG5_aMC+Py8 FxFx samples, and EvtGen 1.2.0 was used for all other MadGraph and Powheg samples. The effect of multiple interactions in the same and neighbouring bunch crossings (pile-up) was modelled by overlaying the simulated hard-scattering event with inelastic events generated with Pythia 8.186 [46] using the NNPDF2.3lo PDF and the A3 tune [71]. The small differences in lepton reconstruction, isolation, and trigger efficiencies between simulation and data are corrected in the simulation on an event-by-event basis by applying efficiency scale factors for each lepton [72][73][74].

Event reconstruction
Events are used if they were recorded during stable beam conditions and if they satisfy detector and data-quality requirements [75]. They are required to have a primary vertex, defined as the vertex with the highest sum of track 2 T , with at least two associated tracks with T > 500 MeV [76]. Events are selected using triggers [77][78][79] that require at least two electrons or two muons, or the combination of at least one electron and one muon; the efficiencies for these triggers plateau in the region of T > 25 GeV.
Electron candidates are reconstructed from inner-detector tracks which come from the primary vertex and are matched to clusters of energy deposits in the EM calorimeter. To fulfil the primary-vertex condition, the electron track's transverse impact parameter significance must satisfy | 0 |/ ( 0 ) < 5.0, where 0 is the transverse impact parameter and ( 0 ) its uncertainty, and the longitudinal impact parameter 0 must satisfy | 0 sin( )| < 0.5 mm, where is the angle of the track to the beamline. Electron candidates must satisfy the 'Medium' likelihood-based identification requirements [72] based on EM shower shapes, track quality, and track-cluster matching. They must also satisfy the 'PflowLoose' [72] isolation requirement. Electron candidates are used in the analysis if they have T ≥ 25 GeV and | | < 1.37 or 1.52 < | | < 2.47.
Muon candidates are identified by matching inner-detector tracks from the primary vertex to either full tracks or track segments reconstructed in the muon spectrometer. The candidates must satisfy the following primary-vertex requirements: the transverse impact parameter significance must satisfy | 0 |/ ( 0 ) < 3.0 and the longitudinal impact parameter must satisfy | 0 sin( )| < 0.5 mm, where 0 , ( 0 ), 0 and are as defined above for the electrons. Muons are required to pass 'Medium' identification requirements [73,74] based on quality criteria applied to the inner-detector and muon-spectrometer tracks. Muon candidates with T ≥ 300 GeV must satisfy tighter identification requirements in the muon spectrometer in order to improve the muon-T resolution. Muons must also satisfy the 'PflowLoose' isolation requirement, built from tracking and calorimeter information, with a muon-T -dependent variable cone size Δ [74]. Muon candidates are used in the analysis if they have T ≥ 25 GeV and | | < 2.4.
Jets of hadrons are reconstructed using a particle-flow algorithm [80] based on noise-suppressed positiveenergy topological clusters in the calorimeter. Energy deposited in the calorimeter by charged particles is subtracted and replaced by the momenta of tracks which are matched to those topological clusters. The jets are clustered using the anti- [81] algorithm implemented in the FastJet package [82] with a radius parameter = 0.4. They are further calibrated according to in situ measurements of the jet energy scale [83]. Analysis jets are required to have a calibrated T ≥ 100 GeV and | | < 2.5.
Electrons, muons and jets are reconstructed and identified independently. An overlap-removal procedure is then applied to uniquely identify these objects in an event. For the lepton-jet overlap removal, softer jets with T ≥ 30 GeV and | | < 2.5 are considered. Preselected jets with a high probability to have been initiated by an electron or a radiated photon such that Δ between the jet and a lepton is smaller than 0.2 are removed. In a second step, leptons closer than Δ = 0.4 to any remaining jet are removed.
Events are selected if they contain a -boson candidate reconstructed from two same-flavour, oppositecharge leptons (ℓ = , ) and with dilepton invariant mass 71 GeV ≤ ℓℓ ≤ 111 GeV. The selected events are also required to contain at least one analysis jet. Events which satisfy the above selection requirements define the inclusive + jets region. A dedicated high-T region is created by requiring the leading jet to have T, 1 ≥ 500 GeV. This latter region is split into the collinear region, where the angular distance between the boson and the closest jet must be Δ min , ≤ 1.4, and the back-to-back region that requires Δ min , ≥ 2.0. An alternative phase space is also defined by requiring T ≥ 600 GeV, labelled as the high-T region, where T is defined as the scalar sum of the T of the analysis jets.

Background estimation
Backgrounds from single-boson, diboson and single-top-quark processes are estimated using the MC samples described in Section 3, while top-pair production and 'multĳet event' contributions from semileptonic decays of heavy quarks, hadrons misidentified as leptons, and electrons from photon conversions are estimated from data. A summary of the composition and relative importance of the backgrounds in the various signal regions is given in Table 2. The overall purity of the + jets selections (defined as the expected fraction of signal events after the final selection) ranges from 94% in the inclusive region to 87%, 86%, 87% and 88% in the high-T , collinear, back-to-back and high-T regions, respectively. Backgrounds are dominated by¯, diboson and EW processes, with fractions of 2%-5% , 2%-6% and 1%-5%, respectively. The fraction of events arising from → + − , + jets, + , and multĳet backgrounds is below the percent level in all signal regions. The¯background is evaluated with a data-driven methodology. A¯-enriched control region is constructed with the same event selection as the signal region, but with ± ∓ final states instead of the same-flavour + − or + − pairs. This control region contains only percent-level contributions from + jets, + jets, diboson, single-top and → + − events. The prediction for the¯distributions in the signal region is obtained by multiplying the corresponding measured distributions in the control region (after subtracting the non-¯contributions) by / and / scale factors [84]. These factors are computed bin by bin in the signal region for each distribution using simulation.
Diboson backgrounds are dominated by two contributions: semileptonic and final states, and fully leptonic diboson final states where the decay products of a gauge boson are reconstructed as jets. The measured kinematics of the boson decay products, as well as the production of one or two additional jets, agrees with predictions from Sherpa, as demonstrated in Refs. [85-89] within the modelling uncertainties described in Section 7. The simulation of the EW events done with Herwig+VBFNLO agrees with measurements performed in ( → ℓ + ℓ − )-enriched phase spaces [12]. Due to their good performance in previous measurements, simulations are used to describe the diboson and EW backgrounds in this analysis.
Multĳet events are assessed with a data-driven approach using a template fit of the ℓℓ distribution. The ℓℓ template for this background is derived from data in a multĳet-enriched control region, which is defined by either inverting or dropping the lepton selection requirements associated with isolation, identification and charge. The sub-percent contributions to the multĳet template that do not originate from the multĳet background are evaluated and subtracted using simulation. The fit is performed over the range 51 GeV ≤ ℓℓ ≤ 151 GeV. The contribution from multĳet events in the analysis is then estimated in the invariant-mass interval of the signal region (71 GeV ≤ ℓℓ ≤ 111 GeV). The resulting fraction of multĳet events is at the sub-percent level and so is neglected in this analysis. Figure 2 shows the data and predicted event yields as a function of T, 1 in the electron channel and as a function of Δ min , in the high-T region in the muon channel. The Sherpa 2.2.11 predictions agree with the data in general, but do not describe it precisely in the full range of the measurements. The distributions are discussed in more detail in Section 8.

Unfolding of detector effects
The cross-section measurements presented in this paper (see Section 1) are performed within the fiducial acceptance region defined by the following requirements: • Two same-flavour, opposite-charge leptons with T ≥ 25 GeV and | | < 2.5 • At least one jet, where jets must have T ≥ 100 GeV and | | < 2.5.
The high-T , collinear, back-to-back, and high-T signal regions are defined in analogy to the detector-level definitions.
The cross sections are defined at particle level, corresponding to 'dressed' electrons and muons. A dressed lepton is defined as the four-vector combination of a prompt lepton (that does not originate from the decay of a hadron or a -lepton, or from a photon conversion) and all prompt photons within a surrounding cone of size Δ = 0.1. The particle level also includes jets found by applying the anti-algorithm with radius   parameter = 0.4 to final-state particles with decay length > 10 mm, excluding dressed -boson decay products. Overlap removal is also applied at particle level: jets with T ≥ 30 GeV within Δ = 0.2 of a dressed lepton are removed, followed by the removal of leptons within Δ = 0.4 of the remaining jets. This overlap removal is applied at particle level in order to best match the detector response, especially in the collinear region where the detector is not able to discriminate easily between nearby objects.
The fiducial cross sections are evaluated from the reconstructed kinematic observables for events that pass the selection described in Section 4. The expected background components, as described in Section 5, are subtracted from the distributions in data.
An iterative unfolding technique [90] with two iterations, as implemented in the RooUnfold package [91], is used to unfold the background-subtracted data to the particle level, thereby accounting for the impact of detector inefficiencies and resolution [72-74, 80, 83]. Before entering the iterative unfolding, the background-subtracted data are corrected for the expected fraction of events passing the detector-level selection but not the particle-level selection. The unfolding is carried out with the response matrices constructed from the Sherpa 2.2.11 + jets samples. The unfolded event yields are divided by the integrated luminosity of the data sample to provide the final fiducial cross sections [92]. The electron and muon channels are unfolded separately and then combined to measure the production cross section for a boson decaying into a single charged-lepton flavour ( → ℓ + ℓ − ).
The binning of all observables is optimised to keep the statistical uncertainty below 10% and to maximize the purity, or the fraction of reconstructed events where the reconstructed and truth values fall in the same bin. The latter is kept above 60%, with typical values of 70 − 90%. In order to mitigate the modelling uncertainty due to migration effects across the lower edge of the T, 1 distribution, this observable is unfolded using two underflow bins within 60 GeV ≤ T, 1 ≤ 100 GeV. In a similar fashion, an underflow bin is added for the T,ℓℓ , T and jet multiplicity distributions for events where the leading jet does not pass the T ≥ 100 GeV selection but instead has 60 GeV ≤ T, 1 ≤ 100 GeV. The unfolding of , is performed in two dimensions using three bins for the T of the closest jet for each bin of , .

Systematic uncertainties Theory modelling uncertainties
Theoretical modelling uncertainties from the MC predictions are considered when unfolding the data and in the comparisons with the cross-section measurements.
Modelling uncertainties are taken into account by varying the QCD scales, the PDFs and, in the case of Sherpa 2.2.11, the virtual EW corrections. The effect of QCD scale uncertainties is defined by the envelope of variations resulting from changing the renormalisation ( r ) and factorisation ( f ) scales by factors of two with an additional constraint of 0.5 ≤ r / f ≤ 2. Uncertainties due to the PDF parameterisation are evaluated using sets of PDF variations [93]. The PDF uncertainties also include a comparison with the nominal MMHT2014nnlo [56] PDF and the CT14nnlo [94] PDFs. For Sherpa + jets and diboson processes, the PDF uncertainties also include a consistent variation of s in the PDF and in the hard scatter based on NNPDF3.0nlo [42]. The prediction from Sherpa 2.2.11 also considers uncertainties related to the NLO virtual EW corrections, derived from the envelope of the additive, multiplicative and exponentiated EW correction schemes [19]. The uncertainties associated with the virtual EW corrections are maximal and amount to 5% where the EW corrections are largest: in the back-to-back region with large T,ℓℓ and Δ min , ≈ . In comparison, the effect of QCD scale uncertainties on Sherpa 2.2.11 predictions ranges between 10% and 60%, with average values near 25%. The corresponding range in MG5_aMC+Py8 FxFx is between 5% and 20%. The differences between these two generators and their uncertainties are further explored in Ref. [19]. In the NNLOjet predictions, the QCD scale uncertainties are typically in the range between 5% and 10% and constitute the dominant systematic component. Due to computational limitations in the NNLOjet program, the predictions do not include PDF uncertainties.
The diboson predictions used in the background subtraction from data include PDF and scale uncertainties. The EW prediction includes the effects of scale, PDF and interference uncertainties, which amount to normalisation uncertainties of 9%, 2% and 25%-40% respectively (see Section 3). The effects of scale and PDF uncertainties on the single-top predictions amount to a total normalisation uncertainty of about 4%, primarily from the normalisation to theory predictions at NNLO and NLO accuracies.

Systematic uncertainties in cross-section measurements
Systematic uncertainties in the measured cross sections stem from experimental, MC-modelling and unfolding uncertainties. The uncertainties are propagated to the data cross sections by varying the subtracted background and the MC inputs to the unfolding procedure (response matrix, fraction of unmatched events, reconstruction efficiency). They are treated as being correlated over kinematic regions, over distributions of observables and, where applicable, over channels and between signal and background processes.
Experimental uncertainties specific to each leptonic final state ( → + − and → + − ): Systematic uncertainties in the lepton-candidate selection are related to the reconstruction, identification, isolation, and trigger [72,73]. Uncertainties in the lepton calibrations can cause changes in the acceptance, owing to the migration of events across the T threshold and the ℓℓ boundaries. The uncertainties in the electron energy scale and resolution are taken into account [95], as are those related to the muon momentum scale, inner-detector and muon-spectrometer resolution, and sagitta-bias correction [73].
Experimental uncertainties common to the electron and muon final states: Systematic uncertainties associated with jet reconstruction are addressed via jet-energy-scale (JES) variations in a 29-nuisanceparameter scheme and jet-energy-resolution (JER) variations in a 13-nuisance-parameter scheme [96,97]. Imperfect modelling of the effects of pile-up leads to acceptance changes for different jet multiplicities. To assess this uncertainty, the average number of pile-up interactions is varied in simulation. The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [98], obtained using the LUCID-2 detector [99] for the primary luminosity measurements.
Modelling uncertainties: Distribution-shape variations from PDF, scale and EW uncertainties in the Sherpa 2.2.11 + jets simulation, computed as described above, are propagated to the unfolded cross sections via the response matrices and associated unmatched-events and efficiency corrections. Although the uncertainties in the simulation can be large, their effect on the cross-section measurement is minimised by the unfolding technique used. The systematic uncertainties in the modelling of background MC samples are propagated to the unfolded cross section via the background subtraction in the signal regions. The background-modelling uncertainty comes mainly from the diboson and EW backgrounds.
Systematic uncertainties associated with the unfolding procedure: Systematic uncertainties account for possible residual biases in the unfolding procedure, such as those due to the modelling of the signal events or the finite bin width used in each distribution. The limited size of a simulation sample can also create biases in the distribution used in the unfolding procedure. The following uncertainties from the unfolding procedure are considered: • The statistical uncertainties of the MC inputs to the unfolding procedure are propagated to the unfolded cross sections with a 'toy' simulation method based on 1000 ensembles (pseudo-experiments) of unfolding inputs.
• The effects of the mismodelling of the data by the MC simulation on the results of unfolding procedure are derived by reweighting the Sherpa 2.2.11 + jets MC simulation at particle level for each unfolded observable, such that the MC simulation distribution matches the backgroundsubtracted data at the reconstruction level. The reweighted MC simulation is unfolded with the non-reweighted response matrix and the uncertainty is obtained by comparing the unfolded result against the reweighted distribution at particle level (non-closure).
• An additional uncertainty is derived to account for more subtle differences between the Sherpa 2.2.11 and MG5_aMC+Py8 CKKWL generators (e.g. hadronisation models, additional soft objects, distributions in other kinematic dimensions) which are not accounted for by the previous method.
A non-closure test is performed where the MG5_aMC+Py8 CKKWL samples are first reweighed to match Sherpa 2.2.11 particle-level distributions for each observable in turn and subsequently unfolded with Sherpa 2.2.11. The uncertainty is evaluated by comparing the unfolded result and the reweighted distribution at particle level.
The total fractional uncertainties of the unfolded differential cross sections in T, 1 and Δ min , for the combined → ℓ + ℓ − measurement, performed as detailed in Section 8, are shown in Figure 3. Table 3 shows the breakdown of the total relative statistical and systematic uncertainties in the measured integrated cross sections for + jets production in the five kinematic search regions. In the inclusive and high-T regions, jet uncertainties dominate with a relative contribution of 2.6% and 2.8%, respectively. In the high-T region and its collinear and back-to-back subregions the jet and data statistical uncertainties are largest, with relative contributions of 3.2% and 2.1% in the high-T region, respectively, 2.8% and 2.9% in the collinear region, and 3.6% and 2.7% in the back-to-back region.  Figure 3: Fractional uncertainties in the measured differential cross-sections in T, 1 in the inclusive region (left) and in Δ min , in the high-T region (right) from the combined → ℓ + ℓ − measurement.

Results
The integrated and differential fiducial cross sections are measured in the electron and muon channels separately, and the compatibility of the results from the two channels is tested. These results are then combined using the Best Linear Unbiased Estimate (BLUE) method [100]. In the combination, all systematic uncertainties except the lepton-related experimental uncertainties are treated as correlated between the electron and muon channels. Data and MC statistical uncertainties are treated as uncorrelated between channels. The combined measurements are consistent with both individual decay channels for the full set of observables. In general, the uncertainties in the measured cross sections are dominated by the systematic uncertainties and are smaller than the uncertainties in the predictions, except for NNLOjet@NNLO, which matches or exceeds the precision of the measurements in some kinematic regions.
Integrated fiducial cross sections for + jets production are evaluated in the inclusive, high-T , collinear, back-to-back and high-T signal regions (see Section 6) by summing over the respective unfolded jet-multiplicity distributions. The -boson and jet transverse momenta (two correlated quantities) are fundamental observables of the + jets process and probe pertubative QCD over a wide range of scales. Moreover, understanding the kinematics of jets in events with vector bosons produced in association with several jets is essential for the modelling of backgrounds for other SM processes and searches beyond the SM. Figure 5 shows the differential cross section as a function of T,ℓℓ and T, 1 . The high-T,ℓℓ region is dominated by the back-to-back topology and receives significant negative corrections due to EW effects. In contrast, events with a high-T jet typically result in both back-to-back and collinear topologies. The Sherpa 2.2.1 and MG5_aMC+Py8 CKKWL generators predict a harder T, 1 distribution than seen in the data, resulting in an overestimation of the cross section for high T, 1 . In contrast, Sherpa 2.2.11 and MG5_aMC+Py8 FxFx show significantly better modelling of the T, 1 spectrum and are also in good agreement with the measured T,ℓℓ spectrum. The smaller cross section from Sherpa 2.2.11 relative to Sherpa 2.2.1 in the high-T region is attributed to the improved matching scheme with a different treatment of unordered histories [19]. The prediction from MG5_aMC+Py8 FxFx models the data more precisely, due to the inclusion of NLO matrix elements with three partons. The NNLOjet@NNLO predictions describe the data very precisely, except for very large values of T,ℓℓ (and T, 1 ), where negative NLO virtual EW corrections of 10%-20% are expected and the QCD-only calculation overestimates the data.
Jet-multiplicity distributions provide an excellent probe of QCD. Whereas the + 1-jet bin is most sensitive to PDF effects, those with higher jet multiplicities are more sensitive to perturbative QCD effects [101]. Jet-multiplicity distributions also probe the validity of predictions in the presence of jet vetoes, which are frequently used in searches that require a specific number of jets in the selection. These vetoes create additional logarithmic terms, which are not explicitly included in the theoretical predictions presented in this paper. Figure 6 shows the differential cross section as a function of the jet multiplicity in the inclusive and high-T regions. As expected [101], the jet multiplicity in the inclusive phase space follows a downward staircase pattern, whereas in the high-T phase space, the cross section increases between 1-jet and 2-jet events followed by a downward pattern for higher jet multiplicities. While Sherpa 2.2.1 and Pred. / data MG5_aMC+Py8 CKKWL tend to overestimate the cross section for higher jet multiplicities, Sherpa 2.2.11 and MG5_aMC+Py8 FxFx agree with the data for both the inclusive and high-T regions, the latter generator at a higher level of precision. The NNLOjet predictions agree with the data in the inclusive and high-T phase spaces at high precision when it is expected from the order of the calculation, and at lower precision when the order of the calculation is exceeded.
The angular distance between the boson and the closest jet, Δ min , , and the ratio of the -boson transverse momentum to the closest-jet transverse momentum, , , provide a way to distinguish between the presence of collinear -boson emission and back-to-back topologies. In the high-T selection, the collinear region is sensitive to logarithmic enhancements in production proportional to s ln 2 ( T, 1 / ), whereas the back-to-back region receives non-negligible virtual EW corrections. Figure 7    accumulation of events with low values of , . This hypothesis is validated in Figure 8, which shows the measurement of the , distribution for the subset of collinear events defined by Δ min , ≤ 1.4 where only the accumulation of low-, events is observed. In contrast, the measurement of the , distribution for the back-to-back selection defined by Δ min , ≥ 2.0 is populated by events with , ≈ 1. The jet multiplicity is also measured separately for the collinear and back-to-back topologies as shown in Figure 9. It is found that the collinear region is dominated by dĳet events whereas the back-to-back region is dominated by + 1-jet events. FxFx and Sherpa 2.2.11 with data even for very collinear events indicates that resummation of the additional logarithmic terms, e.g. via EW showers, is not needed at the present level of theoretical and experimental precision in the experimentally accessible kinematic regime. The NNLOjet@NNLO prediction models the data distribution in both the collinear and the back-to-back regions with high precision. The good performance in the collinear phase space is remarkable, as this region contains a large fraction of events with at least three jets, which are simulated only at LO in NNLOjet@NNLO and not at all by NNLOjet@NLO. The QCD-only calculation overestimates the cross section for exact back-to-back events in the high-T  region, dominated by high-T,ℓℓ events, consistent with the pattern observed in Figure 5.
An alternative way to select a high-energy phase space is by requiring a large value of T or T . The former is often used as a dynamical scale choice, whereas the latter is more suited to selecting a phase space similar to the high-T region. Figure 10 shows the differential cross section as a function of T .
The central values from Sherpa 2.2.1 increasingly overestimate the cross section with increasing T , while still marginally agreeing with the data within modelling uncertainties of up to 50%. The prediction from MG5_aMC+Py8 CKKWL shows a similar trend. In contrast, the state-of-the-art predictions from Sherpa 2.2.11, MG5_aMC+Py8 FxFx and NNLOjet@NNLO model the data well, the last with higher precision. Figure 11 shows the differential cross section as a function of Δ min , and jet multiplicity in the high-T region. The high-T region probes high-energy events where the energy is typically shared by several jets. Compared to the high-T region, the jet multiplicity distribution is shifted towards higher values, and the back-to-back peak in the Δ min , distribution is suppressed relative to events where a jet is in closer proximity to the boson. As in the high-T region, Sherpa 2.2.1 is marginally consistent with the data within a large theoretical uncertainty, with the overestimation of the central values most pronounced in the collinear region and for higher jet multiplicities, with MG5_aMC+Py8 CKKWL showing a similar trend. Both Sherpa 2.2.11 and MG5_aMC+Py8 FxFx are consistent with the data, the latter at higher precision. The NNLOjet@NNLO prediction models this region well and with high precision, even though the fraction of events with at least three jets is larger than in the high-T phase space.          Figure 11: Differential cross section as a function of the angular distance Δ min , between the boson and the closest jet (left) and of the jet multiplicity (right) in the high-T region. In the jet-multiplicity distribution, the last bin is inclusive in higher jet multiplicities. Further details are provided in the caption of Figure 5.

Conclusions
This paper presents measurements of cross sections for a boson produced in association with hightransverse-momentum jets and decaying into a charged-lepton pair. The data were collected from collisions at √ = 13 TeV with the ATLAS detector at the LHC during 2015-2018 and correspond to an integrated luminosity of 139 fb −1 . Measurements were performed on events that contain a -boson candidate reconstructed from an + − pair or a + − pair in association with hadronic jets defined as jets having transverse momentum greater than or equal to 100 GeV. Primarily, only the QCD component of + jets production is measured, treating the EW processes as a background. The paper focuses on selections with very high leading-jet T ( T, 1 ≥ 500 GeV) and very high T ( T ≥ 600 GeV), which are used to study two populations of events -collinear events and back-to-back events -with distinct patterns in distributions of Δ min , , , , and jet multiplicity. The data distributions were unfolded to the particle level and compared with state-of-the-art generator predictions and fixed-order calculations. Both Sherpa 2.2.1 and MG5_aMC+Py8 CKKWL overestimate the cross sections for large values of T, 1 , T and T . The predictions from MG5_aMC+Py8 FxFx, with matrix elements for up to three partons at NLO, offer a significant improvement over MG5_aMC+Py8 CKKWL (which models up to four partons at LO) and in general match the data with high precision. Similarly, Sherpa 2.2.11, with the addition of a fifth parton at LO in the matrix element, the addition of NLO virtual EW corrections, and a different treatment of unordered histories in the parton shower, shows a significant improvement over Sherpa 2.2.1 and agrees with the data. The NNLO calculations at fixed order from NNLOjet describe the data cross sections at a very high level of precision, including in high-T regions where a sizeable fraction of the events have more than two jets. The calculation exhibits a harder T,ℓℓ spectrum than the data in a region where larger negative EW corrections are expected. [7] ATLAS Collaboration,      The ATLAS Collaboration