Exploring high-purity multi-parton scattering at hadron colliders

Multi-parton interactions are a fascinating phenomenon that occur in almost every high-energy hadron--hadron collision, yet are remarkably difficult to study quantitatively. In this letter we present a strategy to optimally disentangle multi-parton interactions from the primary scattering in a collision. That strategy enables probes of multi-parton interactions that are significantly beyond the state of the art, including their characteristic momentum scale, the interconnection between primary and secondary scatters, and the pattern of three and potentially even more simultaneous hard scatterings. This opens a path to powerful new constraints on multi-parton interactions for LHC phenomenology and to the investigation of their rich field-theoretical structure.

At high-energy hadron colliders such as CERN's Large Hadron Collider (LHC), almost every event that gets studied is accompanied by multiple additional parton interactions (MPI) from the same proton-proton collision, cf.Fig. 1.For example in each proton-proton collision that produces a Z or Higgs boson (the "primary" process), models [30,31] suggest that there are about ten additional parton collisions that occur simultaneously, usually involving QCD scatterings of quarks and gluons.MPI are the subject of a rich array of studies, involving effects ranging from partonic correlations inside the proton to colour reconnections between final state quarks and gluons [1,2].Their modelling is an essential component of every major simulation tool [3][4][5].While often thought of as non perturbative, we shall see clearly below that MPI involve transverse momenta of up to 10 GeV and beyond, i.e. close to the scale of many of the primary processes regularly studied at the LHC.In a context where there is an ambitious worldwide effort to bring high precision in perturbative QCD calculations for those primary process [6][7][8], our current partial understanding of MPI scatters risks becoming a limiting factor across much of the LHC programme.
A major challenge in the experimental characterisation of MPI is the difficulty of unambiguously separating the MPI signal from the primary hard scattering.In this letter we propose an approach to investigating MPI scatters in Drell-Yan events that optimally suppresses the contamination from the primary hard scattering.This opens the path to a programme of experimental study of MPI that goes significantly beyond the current state of the art.Features that we will highlight include (a) the clarity of the MPI signal; (b) scope for quantitative investigations of the leading two hard scatters (2HS) that includes direct experimental sensitivity to perturbative interconnection [9][10][11][12][13] between the primary process and the second hard scatter, cf.Fig. 1b; and (c) the potential for observation of high-purity triple parton scattering (Fig. 1c) [14,15], as well as sensitivity to even more than three scatters.
The foundation of our approach is the well-known FIG. 1. Illustration of some MPI configurations that will be probed in this paper: (a) standard double hard scattering, producing a Z boson and a pair of jets; (b) perturbative interconnection between the partons involved in the two hard scatterings, where the q that produces the Z and the q that scatters to produce the dijet system both have a common origin in the perturbative splitting of a gluon; and (c) a process with three hard scatterings.
fact [16] that if one considers events where the Drell-Yan pair has a low transverse momentum, the amount of initial-state radiation (ISR) is strongly constrained.
To illustrate this quantitatively, we examine the average transverse momentum of the leading jet in Drell-Yan events as a function of the Z transverse momentum.Specifically, we consider Z → µ + µ − events and cluster all particles other than the muons with a jet algorithm (the anti-k t algorithm [17] with a jet radius of R = 0.7, as implemented in FastJet [18]).Experimentally, this observable could be studied using charged-track jets (see below), or possibly standard jets in a dedicated low-pileup run.For now, to help expose the basic physical dynamics and scales, we retain all particles in the jet clustering.one is from a resummed calculation (RadISH NNLL [19][20][21]), the other from a Monte Carlo simulation that uses a combination of MiNNLO [22,23] with POWHEG [24][25][26] and Pythia 8.3 [3] (with HepMC2 [27]), and the third is from Pythia alone.All Pythia results use the Monash tune [28].All three curves in Fig. 2 show the same features, namely that for almost the whole range of p tZ , the average leading jet p t is roughly proportional to p tZ (with a proportionality coefficient close to 1), a consequence of momentum conservation between the jet and the Z boson.For p tZ below about 2−3 GeV the average leading jet p t saturates.Events with very small p tZ mostly occur when the transverse recoil from one initial-state radiated gluon cancels with that from other initial-state radiation.In this region, the average leading jet p t has the parametric form ( [29], § A 1) where Λ is the scale of the Landau pole in QCD, M is the invariant mass of the Drell-Yan pair, β 0 = (11C A − 2n f )/(12π) and C F = 4/3, C A = 3, while n f is the number of light quark flavours.With n f = 5, this gives Λ 0.51 M 0.49 .In practice this simple scaling is accurate only for large values of M , and the result from a full NNLL resummation (green curve) can be read off as the intercept of the corresponding curve in Fig. 2, i.e. 2.5 GeV, which coincides well with the intercept of the simulations without MPI (blue curves).Next consider the red curves in Fig. 2, those with MPI.For high p tZ values, the leading jet p t again tracks p tZ .However for low p tZ values, the average leading jet p t saturates at a value of about 10 GeV, which is significantly above the MPI-off result.The interpretation is that in events with MPI, for low p tZ , the leading jet almost always comes from an MPI scatter, not from the hard scatter, and it has a characteristic scale of the order of 10 GeV.This may be surprising if one thinks of MPI as genuinely non-perturbative, but less so if one considers that Pythia simulates MPI as semi-hard scatterings [30,31].(We found similar results with Herwig 7.2's [4,32] implementation of a comparable model [33,34]).
Fig. 2 provides the foundation for the rest of this letter.Specifically, if we consider events with a stringent cut on p tZ , we ensure the near total absence of hadronic radiation from the primary scatter (defined as that producing the Z).Existing experimental work confirms that the relative MPI contribution is enhanced by choosing a low p tZ cut, for example using p tZ < 5 or 10 GeV [35][36][37][38][39]. From Fig. 2 we observe that if we choose a p tZ cut that corresponds to the onset of the low p tZ plateau of the MPI-off curves, i.e. p tZ < C Z = 2 GeV, we will obtain a near optimal selection for the study of MPI: if one takes C Z any higher, one increases contamination from hadronic activity due to the primary hard scatter; if one takes it any lower, there is no further advantage in terms of reducing primary hard-scatter contamination, but one loses cross section (and also reaches the limit of experimental lepton resolution).Our choice selects about 4−5% of the Drell-Yan events that pass the muon cuts, i.e. a cross section after the C Z cut of about 40 pb at √ s = 13.6 TeV.For an LHC Run 3 luminosity of 300 fb −1 , this would yield a sample of about 12 million events.
At first sight, Fig. 2 might suggest that MPI dynamics can be observed only at relatively low p ℓ tj ∼ 10 GeV.However after applying the p tZ cut, we can consider a much wider array of observables, some of which extend over a range of jet transverse momenta.The simplest is the cumulative inclusive jet spectrum, i.e. the average number of jets above some p tj,min , as a function of p tj,min , (2) To a good approximation this observable is given by a straight sum of the number of jets from the primary process and the number of jets from the MPI.The approximation is broken only by the potential overlap (in a cone of size R in rapidity and azimuth) of hadrons from the two scatters, and the approximation is exact in the limit of small R. Precisely for this reason, from here onwards we shall use R = 0.4 rather than the R = 0.7 of Fig. 2. All results (R = 0.4 and R = 0.7) use area subtraction [40,41] to further reduce the impact of such overlap, notably as concerns any underlying-event pedestal of transverse momentum from the softest part of the MPI.We include a jet rapidity cut, |y j | < 2, to mimic the central acceptance of the ATLAS and CMS detectors.
The cumulative inclusive jet spectrum is shown in Fig. 3.It is clear that the vast majority of jets come from  the MPI scatters rather than the primary scatter, even for the relatively large value of p tj,min = 50 GeV.Looking instead at moderately low p tj,min values, Fig. 3 indicates that on average there is one jet with p t ≳ 6 GeV, which is broadly consistent with the plateau at 10 GeV in Fig. 2, considering that Fig. 3 uses R = 0.4 instead of R = 0.7, and that it has a limited rapidity acceptance.Note that for large p tj,min , the sample without MPI is dominated by events where the Z is accompanied by two opposing jets.The Pythia8+MiNNLO sample includes the matrix element for that process at leading order (LO), while Pythia8 does not, thus explaining the observed difference between the two curves for p t ≳ 10 GeV.
It is useful to define the pure MPI contribution to the cumulative inclusive jet spectrum, In an actual experimental analysis, one might want to use a next-to-leading order (NLO) Z + 2jet sample to subtract the hard-event contribution.Let us now see how Eq. ( 3) connects with the widely used "pocket formula" for double-parton scattering.That formula states that the double-parton scattering cross section for two hard processes A and B is given by where σ eff for pp collisions is measured to be of the order of 15−20 mb [42][43][44][45][46][47][48] (for processes involving a vector boson) and is related to an effective area over which interacting partons are distributed in the proton.We take process A to be Z production with p tZ < C Z and process B to be inclusive jet production, and consider a p tj,min that is large enough for the pocket-formula to be valid, i.e. such that σ B /σ eff ≪ 1.This yields ( [29], § A 2) where dσjet dptj is the inclusive jet cross section for jet production, without any requirement that a Z be present in the event. 1 The right-hand-side of Eq. ( 5) does not involve C Z , and thus the pocket-formula prediction is that ⟨n(p tj,min )⟩ pure-MPI C Z should be independent of C Z .The pocket formula is, however, known to be an approximation.The difficulty of obtaining a pure MPI sample has so far limited the scope for investigating more sophisticated theoretical predictions.One particularly interesting effect not captured in the pocket formula relates to perturbative interconnection between the primary scattering and the secondary scattering, as in Fig. 1b, where at least some of the partons entering the two separate hard scattering processes (Z and dijet production) have a common origin, e.g. a perturbative g → q q splitting, with the q involved in Z production and the q involved in di-jet production.
Our procedure of constraining the Z transverse momentum means that the partons that annihilate to produce the Z will almost always have a low transverse momentum, which reduces the likelihood of their having been produced in a perturbative splitting.In contrast, if we relax the constraint on p tZ , we will allow for substantially more initial-state radiation from the partons that go on to produce the Z.The ISR partons can then take part in a separate hard scatter, i.e. increasing the interconnection contribution to 2HS, Fig. 1b.
To evaluate potential sensitivity to this effect, we examine the ratio between the 2HS rate with loose (C Z = 15 GeV) and tight (C Z = 2 GeV) constraints on p tZ , In each case the 2HS rate is normalised to the number of Z bosons that pass the selection cut.With the pocket formula the ratio should be 1, and so an experimental measurement of r 15/2 has the potential to provide powerful constraints on deviations from the pocket formula.6).Note the deviation from 1 when perturbative interconnection is turned on between the primary and secondary hard scatters, i.e. diagrams as in Fig. 1b.The dShower bands correspond to scale variation (see [29] §A 3 for further details).They include only the Zgg final state, which represents about 50% of independent 2HS, and so should be taken as qualitative.No jet rapidity cut is applied.
Fig. 4 shows the r 15/2 ratio evaluated in three ways.The Pythia8+MiNNLO curve corresponds to a full analysis, using Pythia8+MiNNLO curve itself (without MPI), to evaluate the no-MPI contribution for Eq. ( 3).Pythia8 does not include a perturbative interconnection mechanism (though it has correlations related to momentum conservation and colour reconnections [49]), and one sees a result consistent with r 15/2 = 1 to within statistical fluctuations.
Fig. 4 also includes curves from the dShower program [50,51].This is a state-of-the-art code that simulates a pure 2HS component, with the option of including interconnection effects according to Ref. [13].Rather than carrying out a full analysis (which would require a consistent merging with a 1HS component), we determine the r 15/2 ratio based on the truth Monte Carlo information about the transverse momentum of the hard outgoing partons in the 2 → 2 interaction, i.e. the second hard scattering.The pink curve is the result without interconnection (with MSTW2008 PDFs [52]), and is consistent with 1.The orange curve includes interconnection effects, and one clearly sees a 25-30% violation of the pocket formula.The scope for measuring this experimentally in a full analysis depends critically on the systematic errors associated with the subtraction of the no-MPI contribution in Eq. ( 3).The significance of such a signal is discussed in [29] §A 3, for various scenarios of uncertainties on the no-MPI term, and the conclusion is that reasonable assumptions lead to at least 2 standard deviations at low p tj,min , which would correspond to exclusion of the pocket formula.The significance can be raised by increasing the accuracy of the no-MPI predic- FIG. 5.The distribution of the absolute value of ∆ϕ12 between the two leading charged-track jets in events with ptZ < 2 GeV (cf.text for jet cuts).The plot shows a clear signal not just of 2HS (in the peak) but also of 3HS (plateau).
tions, e.g. with improved higher-order calculations.
The final question that we turn to is the sensitivity to more than two simultaneous perturbative scatterings.So far the only attempt to study this experimentally has been in triple charmonium production, where the measured cross section has a large uncertainty [15,53] and where generic difficulties in understanding charmonium production complicate the interpretation of the results.
Here we propose the study of charged-track jets, with moderately low p t cuts.To illustrate the study, we construct charged-track jets using charged particles with |η| < 2.4 and p t > 0.5 GeV.The use of charged particles enables the study of moderately low p t jets even in high-pileup runs, thus exploiting the full luminosity of the LHC.We order the jets in decreasing p t , and first study the two leading jets, with a "product" cut [54], √ p t1 p t2 > 9f chg GeV, and a ratio cut, p t2 > 0.6 p t1 .We quote the cuts in terms of a charged-to-neutral conversion ratio f chg = 0.65.The overall scale of the cuts ensures a non-negligible likelihood that each event contains at least one pair of jets.Fig. 5 shows results for the absolute difference in azimuthal angles between the two jets, ∆ϕ 12 .This observable is expected to peak around ∆ϕ 12 = π when the two jets come from the same hard partonic interaction, and to be uniformly distributed between 0 and π when the two jets come from distinct partonic interactions.The plot clearly shows both a peak and a continuum component.A parton-level based decomposition ([29] § A 5) of each histogram bin shows that the plateau is dominated by events with 3 hard scatterings (3HS), where each of the two leading jets comes from a different HS (each distinct from the one that produced the Z).The enhancement near ∆ϕ 12 = π originates mostly from 2HS where the two jets are from a single HS that is distinct from the one that produced the Z, cf.Fig. 1a.A measurement of ∆ϕ 12 would therefore provide clear and quantifiable indications not only of 2HS but also of 3HS.
With such an unambiguous signal of 3HS, one may wonder if it is possible to gain even further insight.One obvious question is whether one can identify a system with two back-to-back jets from one hard interaction and two further back-to-back jets from another hard interaction, all distinct from the Z hard interaction, cf.Fig. 1c.This appears to be on the edge of feasibility, but also brings sensitivity to 4HS ( [29], § A 4).
To conclude, the study of Drell-Yan events with a tight cut on p tZ opens the door to numerous new studies of multi-parton interactions, with high-purity 2HS samples, sensitivity to the perturbative quantum field theory effects that interconnect primary and secondary scatters, and the scope for extensive investigations into 3HS and perhaps even beyond.The studies outlined here are all possible with existing and Run 3 data.The subset of studies that extends to relatively low jet p t values should be feasible with charged-track jets.There is also ample scope for further exploration, for example in terms of the choices of jet cuts, or studies in other collision systems such as p Pb.We expect experimental results on these questions to have the potential for a significant impact not just on our intrinsic understanding of multi-parton interactions but also for the accurate modelling of hadron collisions that will be needed for the broad range of high precision physics that will be carried out at the highluminosity upgrade of the LHC and at potential future hadron colliders.
Appendix A: Supplementary material Here we provide additional material to help the reader reproduce the results described in the letter.
1.The average transverse momentum of the leading jet when ptZ → 0 In this appendix, we provide a derivation of Eq. (1).Our starting point is the following equation [20,21] describing the small-p t limit of the Drell-Yan spectrum.We proceed to evaluate Eq. (A1) at leading-logarithmic (LL) accuracy (radiation strongly ordered both in angle and transverse momentum).This is sufficient to extract the qualitative scaling given in Eq. (1).A more refined calculation with NNLL accuracy can be performed using the double-differential resummation of Ref. [21], as shown in Fig. 2. At LL, the variable k in Eq. (A1) represents the transverse momentum of the leading emission, which coincides with the leading jet at this logarithmic accuracy.The quantity R(k) is the Sudakov radiator (see e.g.Ref. [20]) and R ′ (k) ≡ dR(k)/d ln(M/k) with M being the invariant mass of the Drell-Yan pair.At LL it reads with being the one-loop coupling constant evaluated at M with β 0 = (11C A − 2 n f )/(12π).For the sake of simplicity we have neglected the effect of parton distribution functions, which are encoded in the Born cross section σ 0 .To proceed, we introduce the asymptotic moments of the leading jet p t as follows where To evaluate the quantity ⟨p ℓ tj ⟩ p tZ →0 we start by taking the limit p t → 0 of the Bessel function J 0 (p t b) → 1.We then focus on the integral in the exponent which reads In order to find an analytic estimate for ⟨p ℓ tj ⟩ p tZ →0 , we can now expand the hypergeometric function in Eq. (A6) to second order in powers of bk as shown in the r.h.s.This reflects the fact that the integral is dominated by finite and small bk ∼ 1, since the small p t limit is driven by radiation with non vanishing transverse momentum (hence finite b and small k).Defining y = bk we thus find where the approximation is obtained using the r.h.s. of Eq. (A6).We show the full numerical solution to the l.h.s. of Eq. (A7) and the r.h.s.approximation in Fig. 6 (left plot).While the r.h.s. of Eq. (A6) is formally accurate for large R ′ (k) we observe that it provides a reasonable approximation also at moderately low values of R ′ (k), even for R ′ (k) of order one, where one might have had concerns about the neglected 1/(R ′ ) 2 terms in Eq. (A7).
We then evaluate the remaining integrals over k, that is for n = 0, 1, using the saddle point method, which leads to Higher order corrections can be included in the above result, and modify the asymptotic scaling at most by a normalisation factor and subleading powers of 1/ ln M Λ .The asymptotic scaling in Eq. (A9) is compared to the full numerical calculation of ⟨p ℓ tj ⟩ p tZ →0 starting from Eq. (A1) as shown in Fig. 6 (right plot), which shows that the approximate solution captures the correct slope as a function of M in the asymptotically large M limit.

Average number of jets with MPI and pocket formula
In this appendix we derive Eq. ( 5), which describes the average number of jets arising from 2HS in the pocketformula approximation.Our starting point is the pocket formula (4), which, for convenience, we recast as an explicit sum over the exclusive jet cross sections of a given multiplicity in the two processes A (Z production with p tZ < C Z ) and B (inclusive jet production).This reads where i, j ≥ 0 and σ (n) A,B denotes the cross section for producing exactly n jets in process A, B. The Drell-Yan cross section with a cut p tZ < C Z is then obtained as i σ (i) A few considerations about Eq. (A10) are in order.The pocket formula is a sensible approximation in situations where the total cross section for the production of one or more jets in process B is much smaller than the effective cross section σ eff .This is connected with the requirement that the total cross section for process A is preserved upon summing inclusively over the jet multiplicities.This translates to the following unitarity condition j=0 σ (j) where the value of σ (0) B term, which must be positive, is implicitly defined so as to ensure the equality.We choose to use the symbol ≃ in Eq. (A12) instead of an equal sign so as to reflect the fact that the expression is sensible only if the sum from j = 1 upwards is much smaller than σ eff .
To obtain ⟨n(p tj,min )⟩ pure-MPI C Z as defined in Eq. ( 3), we then calculate the differential inclusive jet spectrum from Eq. (A10), and get where dσ (n) X /dp tj is the differential inclusive jet spectrum in the subset of events of process X with n jets.We now make use of Eqs.(A11), (A12) and obtain where the use of ≃ in Eq. (A14) follows from Eq. (A12) and Eq. ( 5) then follows from Eq. (3).

Notes on the extraction of r 15/2
Here we further comment on the elements that enter into Fig. 4. In the Pythia8+MiNNLO case, we calculate ⟨n(p t,min )⟩ pure-MPI C Z defined in Eq. ( 3) by taking the difference between the simulation with and without inclusion of the full tower of MPI.For p tj, min ≥ 20 GeV, as considered in Fig. 4, Pythia8 indicates that the dominant contribution stems from events with at most two hard scatters.This observation justifies the comparison between the Pythia8+MiNNLO result and a simulation of just double parton scattering, as is provided by the dShower code.We stress that dShower, unlike Pythia, simulates events with just a fixed number of scatterings, i.e. it provides only the pure 2HS component of the full MPI ladder.In particular, this means that no events with a single hard scattering are present in the simulation.The lack of this 1HS process strongly affects the relative fraction of jets from the Z(+jets) process as compared to the 2 → 2 process.Therefore, instead of performing a full jet analysis, as done in Pythia8+MiNNLO, we extract r 15/2 by running dShower with different generation cuts (p tj, min ) on the jets produced by the second hard scattering (pp → jj) accompanying the pp → Z process.We note that this procedure is not equivalent to the experimental definition of r 15/2 but allows us to estimate the magnitude of interconnection effects with a state-of-the-art calculation.
The dShower simulation relies on 3 light active flavours.Specifically, for the pocket formula simulation we adopt the n f = 3 MSTW2008 set [52] and an effective cross section σ eff = 18.5 mb, while we use the n f = 3 DGS set [13,50] to account also for the interconnection between the two scatterings.Another point to note is that the dShower simulation includes only a subset of partonic channels, specifically those contributing to a Zgg final state (thus the actual interconnected diagram differs from that illustrated in Fig. 1b).In independent 2HS, for the p tj,min range shown in Fig. 4, the gg channel represents 50−60% of the total second-hard scatter dijet rate.We expect that other partonic channels would see a comparable degree of interconnection, resulting in a similar value for r 15/2 across all channels, however this point clearly deserves to be verified with an explicit calculation.Thus the dShower code gives an indication of the order of magnitude of perturbative interconnection effects, but does not, as yet, predict the full quantitative picture.
We estimate the theoretical uncertainty on the dShower prediction by varying for each of the two processes: (i) in the case when interconnection effects are on, the cutoff on the impact parameter, i.e. the ν-scale in Ref. [50], and (ii) the renormalisation scale (µ R ), which also affects the definition of the shower starting scale.More precisely, we    4 and 7 but for a pt cut on the Z of ptZ < CZ = 10 GeV.This shows higher MPI fraction but lower interconnection effects across all values of pt,min.The Z selection results in a Pythia8+MiNNLO cross section of σp tZ <10 GeV ≃ 340 pb.set dShower's ParamNu parameter to be either 1 or 0.5, corresponding to an impact parameter cutoff of the order of the hard scale or half of it, respectively.Regarding µ R , we set the flag MuRisPT to be either True or False, corresponding to setting the renormalisation scale to be of the order of either the invariant mass of the two scatterings or the transverse momentum scale of the outgoing particles (the leptons in the DY case).In all cases, we run with UnequalScale=True so that the shower starting scales in the two hard scatterings are independent.The envelope of all these variations constitutes the uncertainty band displayed in Fig. 4.
One general concern in the extraction of the r 15/2 ratio, Eq. ( 6), is whether the relative MPI contribution to ⟨n(p t,min )⟩ 15 , shown in Fig. 7, is sufficiently large that one can reliably determine ⟨n(p t,min )⟩ pure-MPI 15 after subtracting the no-MPI contribution, including its uncertainties.In particular, one might wonder whether it would be beneficial to lower the upper p tZ cut.Fig. 8 demonstrates the higher purities that can be achieved by lowering the loose cut to 10 GeV, but at the cost of a reduced impact of interconnection effects.Alternatively, one can select a bin in p tZ , e.g. 10 < p tZ < 15 GeV as shown in Fig. 9.In this case, we find an enhanced signal of interconnection effects but low purities.At these low purities, even a small offset in jet energies between the no-MPI and the MPI samples (e.g.due to imperfections of the area subtraction) may result in enhanced systematic errors on the r x/2 determination.This, combined with an enhanced sensitivity to statistical fluctuations, may be the cause of the apparent deviation of the Pythia+MiNNLO r x/2 result from one at high p tj,min .
An analysis of the uncertainties can help us understand which choice of cuts might give the most significant determination of deviations from the pocket formula.In particular one should examine how r x/2 would be determined  experimentally, with j exp being the measured jet rate and h th the theoretically determined no-MPI rate.The latter is subject to a theoretical uncertainty, which we write as ∆h th = f h th , with f being a fractional error.The measured jet rate is affected both by statistical and systematic uncertainties.We assume that the experimental systematic uncertainties would largely be correlated and so cancel in r x/2 .We have checked that the statistical uncertainty is much smaller than the theory uncertainty for an integrated luminosity of 300 fb −1 . 3Therefore, in what follows, we only consider the propagation of ∆h th into the r x/2 uncertainty.We find where ρ quantifies the correlation between the theory error at C Z = x GeV and at C Z = 2 GeV.We remind the reader that the pocket-formula Eq. ( 5), corresponds to r x/2 = 1.In the following, we use as our estimate of r x/2 the dShower result.The theoretical uncertainty is expressed in terms of the MPI fraction or "purity": p = n pure-MPI /(n pure-MPI + h th ), as shown in the lower panels of Figs. 3 and 7. (The purity is defined as that obtained when r x/2 = 1, and is always as obtained with a |y j | < 2 cut).The quantity of interest is the significance (number of σ) for observing effects that go beyond the pocket-formula.It can be obtained as (r x/2 − 1)/∆r x/2 .We use a dShower-like signal as a baseline.We plot the significance in Fig. 10.The three bands show distinct choices of the upper p tZ cut, 15 GeV (as used in the main text), 10 GeV and 10 < p tZ < 15 GeV.The columns show three values of the ρ parameter when evaluating the theoretical systematic uncertainty via Eq.(A17): 0 (f x and f 2 are fully uncorrelated), 0.5 and 1 (full correlation).The different rows show a range of assumptions about the theoretical errors on the hard component: the choices f x = f 2 = 5 − 10% mimic the expected accuracy of the NNLO calculations of the Z + 2-jet rate that should hopefully become available in the next few years.The f x = 10% and f 2 = 20% reflects the uncertainty associated with scale variations in the Pythia8+MiNNLO samples.The observed significance of the deviation from the pocket formula depends strongly both on the assumptions for f 2 and f x and on the value of the correlation parameter ρ, as well as on p tj,min .In the optimistic f x = f 2 = 5% scenario, at the lower end of the p tj,min range, we see at least 4σ significance even with ρ = 0, and over 5σ if some correlation is assumed.Other scenarios still all give at least 2σ at low p tj,min , which would be sufficient to exclude the pocket-formula in the presence of an r x/2 effect of the size suggested by dShower.The much larger significances in the lower left-hand plot are an artefact of an almost exact cancellation of correlated uncertainties between different p tZ cuts.Overall, we see the choice of upper p tZ cut is not too critical.A final comment is that the generally improved significance with small f 2 and f x may provide an additional motivation for Z + 2-jet calculations at NNLO and beyond. 3The impact of the statistical uncertainty should be properly assessed in the case of a dedicated low-pileup run.

Z plus four-jet study
Figure 11 shows a Z plus four-jet study that is intended to help examine the structure of 3HS, in particular the situation where the Z and each of the two pairs of jets arises from distinct hard scatterings.We apply the usual p tZ < C Z = 2 GeV requirements, and the same cuts for the two highest p t jets as in Fig. 5, but with an additional constraint of ∆ϕ 12 > 3π/4, so as to enhance the contribution from the situation where the two leading jets are from the same hard interaction.We then apply product and ratio cuts to a second pair of jets, jets 3 and 4, √ p t3 p t4 > 9f chg GeV, p t4 > 0.6 p t3 .We also apply a rapidity cut |∆y i,j | > 1, with i = {1, 2} and j = {3, 4}, to reduce the likelihood that a jet in the first pair and a jet in the second pair originate from the fragmentation of a single hard parton.Finally, we plot the distribution of ∆ϕ 34 in Fig. 11.We see some degree of peak around ∆ϕ 34 = π, for the most part a consequence of the 3HS that we were trying to isolate.Meanwhile the plateau region receives contributions from a mix of 3HS, 4HS and even some 5HS, illustrating the considerable potential of such a Z + 4-jet analysis.Clearly it would be interesting, in both the dijet and 4-jet studies, to further investigate the structure of different numbers of interactions, for example by varying the jet p t cuts so as to modify the relative contributions from different numbers of hard interactions.

Determination of hard-scattering jet permutations
Several steps are required in order to obtain the breakdown into numbers of hard scatterings shown in Figs. 5 and  11.
In a parton-level Monte Carlo simulation with Pythia8, it is possible to associate each parton with a specific underlying hard scattering.To do so, we use the event record as represented through the HepMC2 package [27].In identifying the hard scattering association of each parton, some care is required, for example, to make sure that initial-state radiation (and its subsequent showering) is correctly treated.This is illustrated in Fig. 12, which shows an event with three hard scatterings (each represented in red).All partons in the event (both intermediate and final) are colour-coded according to their associated hard scatter.One sees that final partons (those with no further vertices emanating from them) may have their origin both before and after the hard scattering.
Given a hard scattering association for each parton, the next step is to obtain a hard-scattering association for a given jet.Our approach is, for each jet, to identify the fraction of the jet's momentum that comes from each of the hard scatters.We declare the hard scatter that contributes the most to the jet to be the main source of that jet.
The next issue is that of how to transfer the information to hadron-level analyses such as those in Figs. 5 and 11.Ideally, one would want to be able to identify, for each hadron, which MPI it came from.However hadrons may come from more than one MPI, for example due to colour reconnections [49].Therefore to obtain Figs. 5 and 11, we carry out two analyses: one at hadron level, which determines the normalisation of each bin, and one at parton level, which determines the relative contributions of different hard-scattering permutations to each bin.
If the hadron-level analysis uses all hadrons, we believe the above procedure to be adequate.However Figs. 5 and 11 use only charged hadrons, which introduces extra fluctuations (for example changing the p t ordering of the jets).To reflect this in the parton-level analysis, we adopt a heuristic approach that splits each parton collinearly into three or four pieces (with equal probability), distributes the parton's momentum randomly between the different pieces, and then assigns each piece a non-zero charge with a 61% probability (we do not impose charge conservation). 4 To test the ability of such a procedure to correctly simulate charged-to-full fluctuations, we take two samples of simulated jets, one at hadron level, the other at parton level.In each case we require the full jet to have a minimum p t of 7 GeV.
In the hadron-level sample, we examine the distribution of the ratio of the charged-hadrons' total p t in each jet to the full jet p t .In the parton-level sample, we examine the distribution of the ratio of the "charged" parton pieces' total p t in each jet to the full jet p t .The two distributions are shown in Fig. 13 and can be seen to be remarkably similar.In determining the relative fractions of different hard-scattering permutations for a given bin of Figs. 5 and 11, we use jets obtained from the clustering of just the charged parton pieces.We have verified that the histograms (summed over all numbers of hard scatters) in the 2-jet and 4-jet analyses have similar shapes in the charged-hadron hadron-level jets partonic scheme FIG. 13.Illustration that the procedure described in the text to split partons into collinear pieces, some of which are taken charged ("partonic scheme"), correctly reproduces the hadron-level distribution of charged-to-full transverse momentum fraction in jets.The pt,j > 7 GeV cut is applied on the full jets.See text for further details.
and charged-parton-piece analyses.We do, however, find a difference in overall normalisation, by a factor of about 1.5−2.5 which is expected, because full hadron-level jets tend to have less energy than the full parton level jets, and the splitting of partons into collinear pieces does not correct for that.We do not expect this to significantly modify the relative fractions of different hard-scattering permutations in Figs. 5 and 11.

Fig. 2 FIG. 2 .
FIG.2.The average leading jet transverse momentum ⟨p ℓ tj ⟩ as a function of the Z transverse momentum, in Z → µ + µ − events, with muon selection cuts as indicated in the plot.A radius of R = 0.7 is used here to reduce the loss of transverse momentum from the jet due to final-state radiation and hadronisation and so more accurately track the transverse momentum of the underlying parton.

FIG. 6 .
FIG.6.Left: comparison of the exact solution of the integral in Eq. (A7) to its analytic approximation in the r.h.s. of the equation.Right: comparison of the exact numerical scaling for ⟨p ℓ tj ⟩p tZ →0 to the analytic asymptotic estimate given in Eq. (A9).

FIG. 7 .
FIG.7.The analogue of Fig.3in the letter but for a pt cut on the Z of ptZ < CZ = 15 GeV.This shows that at moderately high pt, the MPI fraction is still reasonable, ∼ 25%.As discussed in the text, such an MPI fraction should still allow for a quantitatively reliable extraction of ⟨n(pt,min)⟩ pure-MPI 15 , as needed for the evaluation of the r 15/2 .The Z selection results in a Pythia8+MiNNLO cross section of σp tZ <15 GeV ≃ 450 pb.

FIG. 9 .
FIG. 9. Same as Figs.4 and 7 but for a pt cut on the Z of 10 < ptZ < 15 GeV.This shows lower MPI fraction but higher interconnection effects across all values of pt,min.The Z selection results in a Pythia8+MiNNLO cross section of σ10<p tZ <15 GeV ≃ 110 pb.

1 x 1 x 1 x
FIG.10.Statistical significance of the detection of the breaking of the pocket-formula with the r x/2 observable for three different values of ρ in Eq. (A17) (one per column) and three different assumptions for the fractional uncertainties, fx and f2, on the no-MPI cross section (one per row).See main text for further details.

pp s = 13 . 6
FIG. 11.The distribution of ∆ϕ34 from the four-jet study in the text, illustrating the rich decomposition into different numbers of MPI.
FIG.12.A graph representation of the branching and scattering in a parton-level event, as simulated with Pythia8.Each red region corresponds to a hard partonic scattering.The showering associated with production of the Drell-Yan pair is shown in cyan, while the showering associated with each of the two other hard scatters (which are both 2 → 2 processes) is shown respectively in magenta and yellow.Typical Pythia8 events contain significantly more than two additional hard scatters, but the number has been restricted in this graph for simplicity, and simulation of final-state radiation has been turned off for the same reason.
Note that with C Z = 15 GeV, the pure-MPI jet fraction is predicted by Pythia8+MiNNLO to be about 25% at p tj,min = 40 GeV ([29], § A 3), which should be adequate for a quantitative extraction of r 15/2 .FIG. 4. Pythia8+MiNNLO and dShower results for the r 15/2 ratio of Eq. (