Search for jet quenching effects in high-multiplicity pp collisions at √ s = 13 TeV via di-jet acoplanarity

The ALICE Collaboration reports a search for jet quenching effects in high-multiplicity (HM) proton– proton collisions at √ s = 13 TeV, using the semi-inclusive azimuthal-difference distribution ∆ ϕ of charged-particle jets recoiling from a high transverse momentum (high-p T , trig ) trigger hadron. Jet quenching may broaden the ∆ ϕ distribution measured in HM events compared to that in minimum bias (MB) events. The measurement employs a p T , trig -differential observable for data-driven suppression of the contribution of multiple partonic interactions, which is the dominant background. While azimuthal broadening is indeed observed in HM compared to MB events, similar broadening for HM events is observed for simulations based on the PYTHIA 8 Monte Carlo generator, which does not incorporate jet quenching. Detailed analysis of these data and simulations show that the azimuthal broadening is due to bias of the HM selection towards events with multiple jets in the final state. The identification of this bias has implications for all jet quenching searches where selection is made on the event activity.


Introduction
Strongly-interacting matter at extreme temperature and density forms quark-gluon plasma (QGP), a state in which the dominant degrees of freedom are sub-hadronic [1,2].The QGP filled the early universe, and it is generated today in the collision of heavy atomic nuclei at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory and at the Large Hadron Collider (LHC) at CERN.Measurements at these facilities, and their comparison with theoretical calculations based on viscous relativistic hydrodynamics, show that the QGP exhibits complex collective behavior, flowing with very low specific shear viscosity [3].
A key question in the experimental study of the QGP is the limit of its formation in terms of the size of the colliding nuclei.To explore this question, measurements have been carried out for "small collision systems", in which one of the projectiles is a proton or light nucleus and the other a heavy nucleus, additionally with selection of high event activity (EA) such as produced particle multiplicity or forward neutron production (see e.g.Ref. [4]).Final-state hadronic distributions in small systems exhibit experimental signatures which are associated with production of the QGP in heavy-ion collisions [5], including collective flow [6][7][8][9] and enhancement in the production of strange hadrons [10].
Jets are the hadronic remnants of hard (high momentum transfer Q 2 ) interactions of quarks and gluons from the hadronic projectiles.Measurements of jet production and jet structure in elementary collisions are well described by theoretical calculations based on perturbative QCD (pQCD) [11][12][13].In heavy-ion collisions, jets interact with the QGP, resulting in modifications of jet production rates, structure, and angular distributions, which provide incisive probes of the QGP ("jet quenching") [14,15].
Jet quenching is a necessary consequence of QGP formation in small systems, though its effects are expected to be small [16][17][18][19][20][21][22].A common signature of jet quenching is the suppression of inclusive hadron or jet yield measured in heavy-ion collisions compared to that expected by scaling the corresponding yield measured in minimum-bias (MB) pp collisions, using a Glauber modeling of the collision geometry [14,15,23].Inclusive jet yield suppression has been reported using this approach in EA-selected d-Au collisions at RHIC (EA based on forward multiplicity) [24] and in p-Pb collisions at the LHC (EA based on forward transverse energy E T ) [25].However, Glauber modeling using these EA metrics in small systems is subject to significant non-geometric bias [4,22,[26][27][28][29][30][31][32].An alternative choice for the EA metric, based on zero-degree neutron measurements (nZDC), is found to be less biased, though the scaling still has model-dependent assumptions and uncertainties [4].Such biases and uncertainties limit the sensitivity of the measurement of jet quenching effects in small systems using Glauber-scaled inclusive yield observables.At present there is no significant evidence, beyond experimental uncertainties, of jet quenching in small systems using this approach [33][34][35].
The PHENIX Collaboration has recently searched for jet quenching in d-Au collisions at a center-ofmass energy per nucleon-nucleon collision √ s NN = 200 GeV by the measurement of both π 0 and direct photon (γ dir ) inclusive yields in collisions selected by EA [32].Jet quenching is studied using the γ dir inclusive yield to estimate empirically the rate of hard processes, which does not depend upon Glauber scaling.Strong suppression in the ratio of π 0 and γ dir inclusive yields is observed for high-EA relative to MB collisions, though with absolute value that is consistent with unity within the normalization uncertainty.
An alternative approach has been proposed to search for medium-induced inclusive yield suppression in small systems utilizing MB collisions of light nuclei [19,20], which likewise does not require Glauber modeling for the yield scaling.However, this approach cannot be applied to EA-selected event populations in small collision systems, where the strongest signals characteristic of collective flow have been observed [6][7][8][9][10].
The comprehensive search for jet quenching effects in small collision systems therefore also requires Search for jet quenching effects in HM pp collisions at √ s = 13 TeV ALICE Collaboration approaches based on coincidence observables, which are self-normalized and likewise do not require a Glauber model calculation of nuclear geometry for an EA-selected population.The measurement of jets recoiling from a high-p T hadron trigger in p-Pb collisions at √ s NN = 5.02 TeV has set a limit of 0.4 GeV/c (90% confidence) for medium-induced energy transport to angles greater than 0.4 radians relative to the jet axis, for high-EA collisions selected with criteria based both on forward multiplicity and on nZDC [36].Measurement of the distribution of hadrons recoiling from a high-p T jet trigger in EA-selected (nZDC) p-Pb collisions at √ s NN = 5.02 TeV likewise finds no significant jet quenching signal within uncertainties [37].The measurement of azimuthal anisotropy of high-p T hadrons finds small but non-zero second Fourier coefficient v 2 for events selected by EA based on forward E T [38], though such effects cannot be attributed solely to jet quenching.
It therefore remains an open question whether the collective effects observed in small systems are indeed due to QGP formation, or whether they arise from other phenomena [22,[39][40][41].New searches for jet quenching effects in small systems are required to resolve this issue.
In this article we present a novel search for jet quenching effects in EA-selected high-multiplicity (HM) pp collisions at √ s = 13 TeV.Since "collision geometry" is ill-defined for EA-selected pp collisions, inclusive observables are not appropriate for such a search.Rather, we utilize the semi-inclusive hadron+jet acoplanarity observable [36,[42][43][44][45], i.e. the distribution of the azimuthal angle ∆ϕ between a high-p T hadron trigger and correlated recoil jets, comparing ∆ϕ measurements in HM-selected and MB populations.Jets are reconstructed from charged particles using the anti-k T algorithm [46] with resolution parameter R = 0.4.Jet quenching in the QGP is expected to broaden the ∆ϕ distribution relative to that in vacuum, due to in-medium multiple scattering [47][48][49][50][51][52].Indeed, significant in-medium acoplanarity broadening has been observed in central Pb-Pb collisions for large-R recoil jets at low p T , though the physical mechanisms underlying this broaden remain an open question [44,45].However, at present there is no theoretical guidance for the magnitude of the jet transport parameter q [14] or alternative characterizations of jet quenching in EA-selected pp collisions, and this is therefore entirely an experiment-driven search.
The analysis is based on the ∆ recoil observable developed for semi-inclusive coincidence measurements of jets recoiling from a high-p T hadron trigger [42].Precise suppression of jet yield uncorrelated with the trigger particle (uncorrelated background yield) is crucial in this analysis, since the uncorrelated yield of jets generated by multiple partonic interactions (MPI) from independent high-Q 2 processes can mimic azimuthal broadening arising from jet quenching.The ∆ recoil observable provides data-driven suppression of uncorrelated background yield through the difference of trigger hadron-normalised recoil jet distributions in two exclusive trigger p T intervals (Sec.5).Selection of the HM population utilizes a large data sample recorded by ALICE with an online HM trigger during the 2016-2018 LHC pp runs at √ s = 13 TeV.The ∆ϕ distributions from the HM-selected and MB event populations are compared, revealing a striking acoplanarity broadening in HM-selected events.However, similar broadening is also observed in PYTHIA-based simulations.The physical origin of the broadening is elucidated through a detailed study of the rapidity dependence of jet production and the number of recoil jets in HM-selected and MB-selected events.
The paper is organized as follows: Sec. 2 presents the data set and offline analysis; Sec. 3 presents characterization of event activity using forward multiplicity; Sec. 4 presents jet reconstruction; Sec. 5 presents the coincidence observable ∆ recoil ; Sec. 6 presents data corrections; Sec.7 presents systematic uncertainties; and Sec. 8 presents the physics results and their interpretation.

Data set and offline analysis
The ALICE detector and its performance are described in Refs.[53,54].Data for this analysis were recorded during the 2016, 2017, and 2018 LHC runs with pp collisions at √ s = 13 TeV.Events were selected online using signals in the V0 detectors [55], which are plastic scintillator arrays covering the pseudorapidity ranges 2.8 < η < 5.1 (V0A) and −3.7 < η < −1.7 (V0C).The V0 signal is proportional to the total number of charged particles (multiplicity) in the detector acceptance.Two different V0 trigger configurations were employed, called minimum bias (labelled "MB") and high multiplicity ("HM").The MB trigger required the in-time coincidence of V0A and V0C signals, while the HM trigger required the sum of V0A and V0C signal amplitudes (denoted as V0M) to be at least five times larger than the mean signal amplitude in MB events (denoted as ⟨V0M⟩).The HM trigger selected 0.1% of MB events with the largest value of V0M.
The EA is characterized offline by the scaled V0 signal, V0M/⟨V0M⟩ = (V0A + V0C)/ ⟨V0A + V0C⟩, which is insensitive to changes in V0 gain in the different data-taking periods due to scintillator aging.It also provides ordering of events in terms of EA without the need for precise calculation of the absolute V0 signal in model calculations, for a well-defined comparison of such models with data.The value of ⟨V0M⟩ is calculated separately for each data-taking run lasting a few hours, as a function of the collision vertex position along the beam axis.The HM selection is further constrained in the offline analysis to the range 5 < V0M/⟨V0M⟩ < 9.The lower bound of 5 is determined by the online HM trigger threshold, while the upper bound of 9 is determined by the range over which the V0M/⟨V0M⟩ distributions for the three different measurement periods are consistent; higher values may be affected by residual, uncorrected pileup effects.
In the offline analysis, jets are measured at midrapidity using charged particles reconstructed with the ALICE central barrel detectors, covering the range |η| < 0.9.Track reconstruction is based on space points measured by the Inner Tracking System (ITS) and Time Projection Chamber (TPC) [54].Primary event vertices are reconstructed offline based on global tracks, which are required to have space points in the Silicon Pixel Detector (SPD) forming the two innermost layers of the ITS.Accepted events are required to have the primary vertex within |z vtx | < 10 cm, where z vtx is the location of the vertex along the beam axis relative to the nominal center of the ALICE detector.
For MB-triggered events, the pileup rate due to multiple hadronic pp collisions in the same LHC bunch crossing is less than 3.5%.The pile-up contribution is suppressed offline by rejecting events with multiple reconstructed event vertices.Monte Carlo studies estimate that the residual pileup contribution following this cut is negligible for the MB populations and about 1% for the HM population.As discussed in Sec. 5, the observable ∆ recoil used in the analysis provides data-driven suppression of uncorrelated background yield, which also includes residual pileup events that are not rejected by the multiple-vertex algorithm.
After event selection, the data sets have an integrated luminosity of 32 nb −1 for the MB trigger and 10 pb −1 for the HM trigger.
During the data taking, the ITS had non-uniform efficiency, and the analysis therefore utilizes hybrid tracks [13,56] to achieve azimuthally uniform tracking response.Hybrid tracks consist of good quality global tracks with at least one hit in the SPD, and complementary tracks without SPD signals.To ensure good momentum resolution, the momentum of these complementary tracks is determined using the primary vertex as a constraint.Reconstructed tracks with |η| < 0.9 and p T > 0.15 GeV/c are accepted for the analysis.Hybrid track reconstruction efficiency is 0.85 at p T = 1 GeV/c, 0.82 at p T = 10 GeV/c, and 0.74 at p T = 50 GeV/c.Tracking efficiencies for MB and HM events are similar.Primary-track momentum resolution is 0.7% at p T = 1 GeV/c, 1.3% at p T = 10 GeV/c, and 3.7% at p T = 50 GeV/c.
Simulations are utilized for data corrections and for comparison to theoretical calculations.The simulations are based on the PYTHIA 8 event generator [57] with Monash tune [58], and a detailed GEANT3 model [59] of the ALICE detector response, which includes production of secondary particles and realistic hit digitization.Events generated by PYTHIA 8 without detector effects are denoted "particle-level," and such events passed through GEANT3 are denoted "detector-level." For particle-level events, the V0A and V0C responses are determined by counting the number of charged Search for jet quenching effects in HM pp collisions at √ s = 13 TeV ALICE Collaboration particles in their acceptance.The coincidence requirement of the online trigger is modeled by requiring particle-level events to have particles in both V0A and V0C, while detector-level events are required to have GEANT-generated hits in both V0A and V0C.The MB events are used to calculate the V0M distributions at both detector and particle level.Figure 1 shows the V0M/⟨V0M⟩ probability distribution for MB pp collisions at √ s = 13 TeV.The lower limit for HM event selection, V0M/⟨V0M⟩ = 5, is indicated by the red dashed-dotted line.The figure also shows the V0M/⟨V0M⟩ probability distribution at the particle and detector level for PYTHIA 8generated events for MB pp collisions at √ s = 13 TeV.These distributions differ because of a large contribution at forward angles of secondary particles generated in detector material [60].The particlelevel distribution falls more rapidly than that observed in data in the range V0M/⟨V0M⟩ > 4.However, the detector-level distribution qualitatively reproduces that observed in data, with probability densities which lie within a factor ∼ 2 of each other over a range of nine orders of magnitude in V0M/⟨V0M⟩.

V0M/⟨V0M⟩ distributions
In order to compare PYTHIA 8 particle-level HM-selected distributions to data, we assume that the secondary particle yield due to interactions in detector material is on average proportional to the primary multiplicity in the V0 acceptance.The HM selection for the particle-level distribution is therefore chosen to select the same fraction of the MB cross section (0.1%) as the HM selection V0M/⟨V0M⟩ > 5.0 used for data.This selection corresponds to particle-level V0M/⟨V0M⟩ > 4.4, as indicated by the dark dashed line in Fig. 1.
We note, however, that such selections of the same cross section fraction at the particle and detector level only select the same fraction of the distributions.They cannot select the same population eventby-event, due to significant fluctuations in the correlation between particle-and detector-level events.For the detector model used in this study, about 35% of events passing the HM selection at the detector level also do so at the particle level.This is a generic issue, with similar features expected for any event selection in small systems based on forward EA.
Search for jet quenching effects in HM pp collisions at √ s = 13 TeV ALICE Collaboration

Jet reconstruction
Jet measurements in this coincidence analysis are corrected for two distinct background effects: uncorrelated jet yield, which is corrected using the ∆ recoil observable discussed in Sect.5; and p T -smearing of correlated jets due to overlap with uncorrelated background components, as detailed in this section.These corrections are carried out in distinct analysis steps.
Several types of reconstructed jets are consequently used in the analysis, which are distinguished using the notation defined in Refs.[36,42].For data, p raw,ch T,jet refers to the raw output of the jet reconstruction algorithm; p reco,ch T,jet refers to p raw,ch T,jet after subtraction of the event-wise estimate of the background contribution to p T,jet , ρ × A jet , which is the product of median jet p T density in the event and the jet area (Eq.1); and p ch T,jet refers to the fully corrected jet transverse momentum.For simulations, p part T,jet refers to charged-particle jets at the particle level, and p det T,jet charged-particle jets at the detector level; both quantities are corrected by ρ × A jet using the following procedure.
Jets are reconstructed from accepted charged-particle tracks.Particles are assumed to be massless and their four-momenta are combined with the boost-invariant p T recombination scheme [61].Jet reconstruction is carried out twice for each event.The first reconstruction pass uses the k T algorithm [61] with R = 0.4 and accepts jets with |η jet | < 0.9 − R. The first-pass jet population is used to determine ρ, the event-wise estimate of the background energy density [62], where p raw,ch,i T,jet and A i jet are the raw jet p T and the area [63] of the i th jet in the event.Jet area is calculated using the ghost area method of FastJet, with a ghost area of 0.005 [63].The two hardest jets in the event are excluded from the median calculation.The most probable value of ρ is zero in both the MB and HM populations, while the mean value of ρ is 0.09 GeV/c for MB and 1.24 GeV/c for HM-selected events.For events containing a charged track in |η| < 0.9 with 20 < p T < 30 GeV/c, the mean value of ρ is 0.39 GeV/c for MB and 1.62 GeV/c for HM events.
The second reconstruction pass uses the anti-k T algorithm [61] with R = 0.4.The acceptance for the second pass is likewise |η jet | < 0.9 − R over the full azimuth.
The jet p T obtained from the second pass is then adjusted for the median background p T density ρ according to [62] The jet p T scale and jet p T resolution are the same as in Ref. [64].

Observables and raw data
The analysis utilizes a differential observable based on the semi-inclusive distribution of charged-particle jets recoiling from a high-p T trigger ("h+jet") [42] (see also [36,43]).The key components of this approach are summarized in this section.
The goal of the analysis is the search for broadening of the ∆ϕ distribution in HM-selected events due to medium-induced jet scattering, by comparison to the MB population.A significant source of uncorrelated background yield to this process arises from MPIs, in which multiple uncorrelated high-Q 2 partonic interactions occur in the same pp collision, with one such interaction generating a trigger hadron and an-other generating a recoil jet in the acceptance.The ∆ϕ distribution of such MPI pairs is by definition uniform on average, thereby limiting the measurement sensitivity to broadening of the ∆ϕ distribution from medium-induced scattering of correlated recoil jets.
Precise background yield correction must be carried out in a fully data-driven way, without model dependence.We therefore employ the ∆ recoil observable [42], which is the difference between semi-inclusive recoil jet distributions for two ranges of p T,trig , both normalized to the number of trigger hadrons, where TT denotes "trigger track."In this analysis, TT Sig = TT{20, 30} specifies the range 20 < p T,trig < 30 GeV/c for the Signal trigger distribution, and TT Ref = TT{6, 7} specifies the range 6 < p T,trig < 7 GeV/c for the Reference trigger distribution.These intervals were chosen to optimize the opposing requirements of obtaining high statistical precision and limiting the kinematic range for more precise comparison of the same observables with different EA selections.The number of trigger hadrons measured in each TT class is denoted N trig .The azimuthal difference ∆ϕ between TT and recoil jet is defined to have the range [0, π] radians.
The ∆ recoil distribution is a two-dimensional function of p T,jet and ∆ϕ [36,42].We define its onedimensional projections, ∆ recoil (p T,jet ) and ∆ recoil (∆ϕ), onto the p reco,ch T,jet and ∆ϕ axes respectively, for restricted ranges in the other kinematic variable.These projections are shown in Figs. 2 and 3 both prior to and after corrections, as indicated by the functional argument (e.g.∆ recoil (p reco,ch T,jet ) or ∆ recoil (p ch T,jet )).The scaling factor c Ref , which is applied to the Reference distribution (second term in Eq. 3), accounts for the different phase space available to observe uncorrelated yield in the Signal and Reference distributions [42,43].In this analysis, the value of c Ref is determined from the ratio of trigger-normalized Signal and Reference recoil jet yields in the bin 0 < p reco,ch T,jet < 1 GeV/c, which is expected to be dominated by uncorrelated background yield.This gives values c Ref = 0.95 ± 0.03 (syst.) in MB events and c Ref = 0.94 ± 0.03 (syst.) in HM events.
Since the uncorrelated yield is by definition independent of TT, it therefore contributes with equal magnitude to the two terms in Eq. 3 and is therefore removed by the subtraction.While ∆ recoil is a differential observable and not an absolutely normalized yield, its two terms are nevertheless calculable perturbatively [65].Measurements of ∆ recoil in minimum-bias pp collisions are well described by PYTHIA 8 [36,42].
In the analysis of both MB and HM-selected events, the dataset is divided into two statistically independent subsets, with the Signal distribution determined using 95% of all events, and the Reference distribution determined using the remaining 5%.These fractions were chosen to provide an equal number of trigger hadrons in the two TT classes, in order to optimize the statistical precision of ∆ recoil .The statistical error due to the N trig normalization is negligible.
If two hadrons in an event satisfy the TT condition, one is chosen at random.This selection ensures that the p T -differential TT distribution has the same shape as the inclusive charged-hadron yield, which is an essential requirement for a semi-inclusive measurement [42].

Corrections
Particle-level jets are clustered from all final-state charged particles generated with PYTHIA 8 as described in Sec. 2, except for weak decay daughters [66,67].The measured (detector-level) ∆ recoil distribution is corrected to the particle level.
Subtraction of the scaled Reference distribution in Eq. 3 accurately removes the uncorrelated jet yield [42].However, the resulting ∆ recoil distribution is still subject to instrumental effects which smear p T,jet and ∆ϕ, whose correction is discussed in this section.These smearing effects are encoded in the fourdimensional response matrix R instr , which maps the true (particle-level) distribution of ∆ recoil onto the measured distribution, The p T,jet and angular smearing due to instrumental effects is very similar in the inclusive and recoil jet populations in these simulations.In addition, the inclusive population has a substantially larger sample in the simulated dataset.The matrix R instr is therefore constructed using the inclusive jet population, by matching particle-level and detector-level jets in phase space within ∆R = (∆η) 2 + (∆ϕ) 2 < 0.3.
The response matrix has p T bins of width 1 GeV/c, so that the difference in the spectrum shape of the inclusive and recoil jet populations has negligible effect on this procedure.
Track reconstruction efficiency, which is the dominant instrumental effect, is found to be the same for MB and HM events.The MB and HM analyses therefore use the same response matrix, obtained using MB events.Unmatched particle-level jets are tabulated and their rate is applied as an efficiency correction, following the unfolding correction discussed below.Jet matching efficiency for both MB and HM events is 0.98 at p ch T,jet = 10 GeV/c and consistent with unity at higher p ch T,jet .The corrected ∆ recoil distribution is calculated by regularized inversion of Eq. 4, using two-dimensional iterative Bayesian unfolding [68] implemented in the RooUnfold package [69].The input distribution ∆ Meas recoil is specified in the range 10 < p det T,jet < 100 GeV/c and π/2 < ∆ϕ det < π rad.The prior distribution for initiating the unfolding is the ∆ recoil particle-level spectrum calculated using PYTHIA 8 simulations.For unfolding the MB dataset, the prior is calculated using MB PYTHIA 8 events, whereas for unfolding the HM dataset the detector-level HM selection is applied to the simulated data, as shown in Fig. 1.
Regularization is optimized by requiring that the unfolded distributions from successive iterations exhibit a mean change of less than 2%, averaged over all p T bins.The optimum number of iterations using this criterion is found to be 5, for both the MB and HM analyses.
While correction of the ∆ recoil distribution for instrumental effects is carried out by regularized unfolding using the full instrumental response matrix R instr , insight into the unfolding procedure can also be gained by parametric characterization of the main detector-level effects.Detector-level effects in p T,jet are characterized by the jet energy scale shift (JES) between particle-level and detector-level, ⟨(p det T,jet − p Azimuthal angular smearing is likewise found to have the same magnitude in MB and HM events.Azimuthal smearing of tracks contributing to TT{20, 30} has RMS ≈ 2 × 10 −3 rad, while azimuthal smearing of the jet axis has RMS ≈ 26 × 10 −3 rad for 10 < p ch T,jet < 15 GeV/c and 11 × 10 −3 rad for 60 < p ch T,jet < 80 GeV/c.The corresponding distributions of ∆ϕ have RMS of 34 × 10 −3 rad for 10 < p ch T,jet < 15 GeV/c and 13 × 10 −3 rad for 60 < p ch T,jet < 80 GeV/c.Subtraction of the reference distribution in the ∆ recoil observable corrects the uncorrelated background yield.However, the resulting two-dimensional distribution as a function of p ch T,jet and ∆ϕ is still smeared by residual fluctuations of the uncorrelated background component, which causes the underlying event density to deviate locally from ρ. Correction for such fluctuations requires model assumptions (see e.g.Refs.[35,64]) and is therefore not included in the unfolding; rather, its magnitude is assessed for two common model choices and the variation in unfolded distributions is included in the systematic uncertainty, as discussed in Sec.7.4.
The unfolding procedure is validated by a closure test in which the input is PYTHIA 8 detector-level events, the full analysis chain, including unfolding, is carried out, and the output is compared to the particle-level PYTHIA 8 spectrum.Good agreement between both distributions is found within statistical uncertainties, thereby confirming the robustness of the applied corrections.

Systematic uncertainties
The main sources of systematic uncertainty in the measurement of ∆ recoil are related to track reconstruction efficiency; track p T resolution; the unfolding procedure; the determination of the scaling factor c Ref ; and residual background p T -density fluctuations relative to ρ.The systematic uncertainty due to each source is estimated by varying the appropriate parameters and rerunning the full analysis chain.
However, the finite statistical precision of the data imposes a limit on the precision with which the spectrum can be unfolded.This limit is taken into account in the determination of the systematic uncertainties by using the following procedure, which is described in detail in Ref. [36].For each parameter variation the spectrum is unfolded 10 times, with the central values of the spectrum varied randomly and independently using a Poisson distribution corresponding to the statistical error of each data point [36].The ratio of the unfolded spectrum from each iteration to that of the primary analysis is calculated, and at each point the median of this set of ratios is assigned as the systematic uncertainty from this source.Tables 1 and 2 show representative values of the systematic uncertainty in ∆ recoil (p ch T,jet ) and ∆ recoil (∆ϕ) measurements for MB and HM-selected event populations, respectively.

Tracking efficiency and track p T resolution
The tracking efficiency uncertainty is 3% [56].To assess the corresponding uncertainty in the ∆ recoil distribution, a variation of R instr is constructed in which 3% of all tracks are randomly discarded.While it is in practice not possible to generate R instr with 3% higher tracking efficiency, this uncertainty is expected to be symmetric.The resulting uncertainty ranges from < 1% at low p ch T,jet to 7% at high p ch T,jet , with only minor dependence on EA.
To assess the systematic uncertainty due to track p T resolution, two different instances of R instr are generated, with p T -resolution corresponding to that of either real data or detector-level MC data hybrid tracks.This source makes negligible contribution to the total systematic uncertainty.

Unfolding
The unfolding procedure has several parameters whose values influence the corrected ∆ recoil distribution: number of iterations; choice of prior spectrum; and range and binning of the raw input distribution.Each source was varied independently: -The regularization condition was varied by ±1 iterations with respect to the optimized value of 5.
The corresponding uncertainty is found to be small, since unfolding converges rapidly to a stable result.
-Variations in the prior spectrum were obtained using the particle-level ∆ recoil spectra generated by the POWHEG MC event generator [70,71] matched to PYTHIA 8 for parton shower and hadronization, and with different choices of regularization and factorization scale.
-The p T,jet binning was varied by shifting the bin boundaries by 1-2 GeV/c, and by changing the lower bound of the input spectrum from 10 GeV/c to 6 GeV/c.The binning in ∆ϕ was not varied, since ∆ϕ smearing effects are small.
The systematic uncertainty attributed to unfolding is the maximum deviation in the ∆ recoil spectrum from varying these parameters, relative to the ∆ recoil spectrum using the primary analysis parameters.
For ∆ recoil (p ch T,jet ), the resulting relative systematic uncertainty is about 0.6% for MB and 0.9% for HM events, with a weak dependence on p ch T,jet .For ∆ recoil (∆ϕ), the relative systematic uncertainty is smallest at ∆ϕ = π and increases monotonically towards ∆ϕ = π/2, for both MB and HM events.

Scaling factor c Ref
The value of the c Ref scaling factor in Eq. 3 was varied in the range [0.9, 1].This range brackets the c Ref values obtained by changing the p reco,ch T,jet bin in which it is evaluated from (0, 1) GeV/c to (−1, 0) GeV/c, and by its variation with ∆ϕ.Different choices of c Ref modify the ∆ recoil spectrum relative to that of the primary analysis result, with uncertainty decreasing as a function of p ch T,jet .Representative values are provided in Tables 1 and 2.

Background fluctuations
As discussed in Sec.6, no correction is applied directly for the effect of residual background fluctuations; rather, a model-dependent estimate of its magnitude contributes to the systematic uncertainty.For this estimate, a response matrix which encodes the effect of residual background fluctuations, R bkgd , is convoluted with the instrumental response matrix R instr in Eq. 4. The matrix R bkgd is determined for events selected with TT{20, 30}, using two methods: -Calculate the sum of track p T in a cone R = 0.4 placed randomly in the acceptance, excluding overlap with the leading and sub-leading jets, and the TT.This sum is corrected for the median background density, where πR 2 is the cone area and ρ is defined in Eq. 1.The matrix R bkgd is the distribution of δ p RC T .
-Embed a high-p T track perpendicular in azimuth to the TT and at the same value of η, with p T distributed uniformly in the range 0-20 GeV/c.Jet reconstruction is then carried out, the jet candidate containing the embedded track is identified, and the quantity δ p ET T is calculated as where p emb T is transverse momentum of the embedded track.The matrix R bkgd is the distribution of δ p ET T .
Unfolding is then carried out for both choices, and the assigned systematic uncertainty is the maximum difference of the two unfolded ∆ recoil distributions from that of the primary analysis.

Total systematic uncertainty
The total systematic uncertainty of the unfolded ∆ recoil distribution is the quadrature sum of the contribution from each source.The systematic uncertainties from all sources except unfolding are correlated between the HM and MB analyses.This correlation is accounted for in the systematic uncertainty of ratios of ∆ recoil distributions by estimating the systematic uncertainty directly from the spread of the ratios calculated for each variation of the analysis configuration.
8 Results Figure 3 shows fully-corrected distributions of ∆ recoil (p ch T,jet ) and ∆ recoil (∆ϕ) measured in MB and HMselected pp collisions at √ s = 13 TeV, together with calculations based on PYTHIA 8, Monash tune, and a pQCD calculation at LO with Sudakov broadening [51,72,73] (the latter only for ∆ recoil (∆ϕ) in MB collisions).Since the pQCD calculation is LO, there are large uncertainties in its normalization; the ∆ recoil (∆ϕ) distributions from the calculation are therefore scaled to the integrated yield of the data in the same p ch T,jet bin in order to compare their shapes.Both PYTHIA 8 and the pQCD calculation (∆ recoil (∆ϕ) shape only) are consistent with the distributions measured in MB events, within experimental uncertainties.

Fully corrected distributions
Comparison of the MB and HM ∆ recoil (p ch T,jet ) distributions reveals a yield suppression in HM collisions that is largely independent of p T,jet , though there is a hint of a harder recoil jet spectrum for HM events.The mean value of the yield ratio HM/MB in the left panel of Fig. 3 is 0.78.The ∆ recoil (∆ϕ) distributions show that the jet-yield suppression in HM events occurs predominantly in the back-to-back configuration, with the yield ratio HM/MB in the bin ∆ϕ ∼ π measured to be 0.69 for 20 < p ch T,jet < 40 GeV/c (Fig. 3, middle panel) and 0.78 for 40 < p ch T,jet < 60 GeV/c (Fig. 3, left panel).The total yield is suppressed, while the azimuthal distribution is broadened; such broadening could arise from jet quenching, i.e. mediuminduced jet scattering occurs preferentially in HM events.Notably, however, PYTHIA 8 particle-level distributions likewise exhibit jet yield suppression and azimuthal broadening for HM-selected events, accurately reproducing the measured distributions.Since PYTHIA 8 does not incorporate jet quenching, this disfavors jet quenching as the predominant effect generating the broadening seen in data.
8.2 Origin of HM induced TT-jet acoplanarity in PYTHIA 8 Figure 4: PYTHIA 8 particle-level simulation of the probability distribution of the yield of charged-particle jets (R = 0.4) recoiling from a high-p T hadron (TT{20, 30}) as a function of η jet for various p ch T,jet intervals, in pp collisions at √ s = 13 TeV.Left: MB events; right: HM events.V0A and V0C acceptances are also shown.
To clarify the origin of the broadening, a more detailed investigation is carried out using both data and particle-level distributions from PYTHIA 8 simulations.
Figure 4 shows the calculated pseudorapidity (η jet ) distribution of charged-particle jets recoiling from a high-p T trigger hadron with TT{20, 30} in MB and HM-selected event populations, for various lower thresholds in p ch T,jet .Since the per-trigger yield varies with p ch T,jet , these distributions are normalized to unit integral to enable a direct comparison of their shapes.The acceptances of V0A and V0C are also shown; note that V0C covers smaller values of |η| than V0A.While the η jet distribution for MB events is symmetric, the η jet distribution for HM events is highly asymmetric, with significant enhancement in the relative rate of recoil jets in the V0C acceptance.The enhancement is largest for the highest value of p ch T,jet .
Figure 5 shows similar distributions for recoil jets with p ch T,jet > 25 GeV/c, in various intervals of V0M/ ⟨V0M⟩.A striking enhancement is observed in the per-trigger jet yield within the V0C acceptance for the largest values of V0M/ ⟨V0M⟩.Since HM events are selected based on the value of V0M/ ⟨V0M⟩, it is evident that the HM selection induces a bias which enhances the rate of hard recoil jets in the V0C acceptance.In both data and simulations, the probability to observe a single jet is suppressed in HM relative to MB events by 2-3%, while the probability to observe multiple jets is significantly enhanced.
The effect of the HM selection bias shown in Figs. 4, 5 and 6 on the acoplanarity distributions in Fig. 3 can be understood as follows.The HM requirement preferentially selects events with a jet in the V0C acceptance.For the leading-order (LO) di-jet channel, in which the TT and recoil jet are azimuthally back-to-back, this depletes both the per-trigger rate for recoil jets observed in the central barrel (|η jet | < 0.9 − R) at ∆ϕ ∼ π and the relative rate of single recoil jets in the acceptance.The recoil jets observed in the ALICE central barrel in HM events therefore arise at enhanced rate from higher-order production processes, with multiple recoil jets in the final state.These jets from higher-order processes have broader distribution in ∆ϕ than jets from LO production.This picture is confirmed by carrying out the PYTHIA 8 analysis for much larger recoil jet acceptance, |η jet | < 5.6.In that case, no back-to-back yield depletion is observed.
The yield suppression and azimuthal broadening seen in Fig. 3 may therefore arise from the effect of different phase space being used for EA characterization (V0A and V0C) and recoil jet measurement (central barrel), combined with the HM selection bias towards events with a jet in the V0C acceptance.It therefore cannot be attributed uniquely to jet quenching in HM-selected events.This conclusion about interpretability in terms of jet quenching of the effects seen in Fig. 3 is more general than this specific analysis, however.The bias of the HM event selection identified here is generic, and must also be taken into account for the interpretation in terms of QGP formation of other phenomena observed in small systems.

Summary
This article reports a search for jet quenching effects in high multiplicity pp collisions at √ s = 13 TeV, based on the semi-inclusive azimuthal distribution of charged-particle jets recoiling from a high-p T hadron trigger in the ALICE central barrel acceptance.Significant azimuthal broadening is observed in events selected to have high multiplicity in forward detectors, which may arise from jet quenching.However, similar broadening is also observed in simulations with the PYTHIA 8 event generator, which does not incoporate jet quenching.
Detailed analysis of the data and simulations reveals that the high-multiplicity event selection, which is based on the V0 detectors at forward and backward rapidities, preferentially selects events with an energetic jet in the forward detector, consequently biasing the acoplanarity distribution measured in the central region.We conclude that this observable is therefore not uniquely sensitive to jet quenching effects in small collision systems.This bias is generic, however, and it should likewise be taken into account in the interpretation of other measurements which probe the formation of a quasi-equilibrated quark-gluon plasma in small collision systems.

Figure 1 :
Figure 1: Probability distribution of V0M/⟨V0M⟩ in MB pp collisions measured at √ s = 13 TeV, and in simulated MB pp events generated by PYTHIA 8 at the particle and detector level.The vertical dashed lines indicate the lower bound for HM selection for data and particle-level simulations.

Figure 2 Figure 2 :
Figure 2 shows selected ∆ recoil (p reco,ch T,jet ) and ∆ recoil (∆ϕ) distributions, together with their corresponding Signal and Reference distributions.The Signal and Reference distributions have similar magnitude only in the region p reco,ch T,jet < 20 GeV/c, while at larger values of p reco,ch T,jet the Reference distribution falls below the Signal distribution and ∆ recoil is similar in magnitude to the Signal distribution.
⟩, while jet energy resolution (JER) is its width, σ (p det T,jet )/p part T,jet , where σ (p det T,jet ) is RMS of the p det T,jet distribution.For p part T,jet = 10, 40, and 70 GeV/c, JES is −10%, −14%, and − 18%, and JER is 25%, 22%, and 23%, respectively.Other systematic effects which are not encoded in the response matrix, such as the precision of the ALICE central barrel B-field scale (∼ 10 −3 ), make negligible contribution to JES and JER.

Figure 6
Figure 6 provides additional insight into the bias induced by the HM event selection.The figure shows the probability to observe a specific number of jets with R = 0.4 and p ch T,jet > 15 GeV/c recoiling against a high-p T hadron trigger (TT{20, 30}) in the ALICE central barrel acceptance (|η jet | < 0.5), for MB and HM pp collisions at √ s = 13 TeV.Distributions are shown for both data and PYTHIA 8 calculations.

Figure 5 :
Figure 5: Same as Fig. 4, but for recoil jets with p ch T,jet > 25 GeV/c for various intervals in V0M/ ⟨V0M⟩.Distributions are normalized per trigger.

Figure 6 :
Figure 6: Distribution of probability per trigger hadron of number of jets with R = 0.4 and p ch T,jet > 15 GeV/c recoiling from a high-p T hadron (TT{20, 30}) at midrapidity (|η jet | < 0.5), for MB and HM pp collisions at √ s = 13 TeV.Left panels: ALICE data; right panels: simulations at the particle level with PYTHIA 8 Monash tune.Lower panels show the ratio HM/MB.The insert in the lower panels has magnified vertical scale to show the ratio of probabilities to observe a single jet.

Table 1 :
Main sources of systematic uncertainty and total uncertainty in ∆ recoil (p ch T,jet ) and ∆ recoil (∆ϕ) in representative bins, for MB events.