Predicting parton energy loss in small collision systems

Medium induced parton energy loss is not conclusively established neither in very peripheral heavy-ion collisions nor in proton-ion collisions. However, the standard interpretation of azimuthal momentum anisotropies in theses systems implies some partonic rescattering. The upcoming lightion runs at the LHC provide a unique opportunity to search for parton energy loss in different systems of similar size. Here, we make predictions for the expected parton energy loss signal in the charged hadron spectra in a system size scan at LHC. We test a large set of model assumptions against the transverse momentum and centrality dependence of the charged hadron nuclear modification factor in lead-lead and xenon-xenon collisions at the LHC. We then attempt to make a model agnostic prediction for the charged hadron nuclear modification factor in oxygen-oxygen collisions.


I. INTRODUCTION
The observed factor 5 suppression of the charged hadron nuclear modification factor R h AA in central √ s N N = 130 GeV Au-Au collisions at RHIC marked the start of experimental energy loss studies two decades ago [1,2]. Pb+Pb collision data from the LHC showed that this quenching increases mildly with center-of-mass energy, and that nuclear modifications remain visible in hadron spectra up to the transverse momentum p ⊥ ∼ O (100 GeV) [3][4][5][6]. An important early finding at RHIC was that (within experimental uncertainties) quenching disappears in d+Au collisions where no dense medium was expected to interact with high-p ⊥ partons in the final state [7][8][9]. This finding was later corroborated at LHC where quenching is absent in TeV-scale pPb collisions [10][11][12][13].
A reassessment of the conclusions drawn from these data in small systems may be needed in the light of the recent LHC discovery of strong collectivity ("flow") in soft multi-hadron correlations [14][15][16][17], and its confirmation in the subsequent analysis of small collision systems at RHIC [18,19]. According to the standard phenomenological interpretation, v n measurements indicate significant final state interactions between colored degrees of freedom in small collision systems. This raises the questions of why high-p ⊥ energy loss effects have escaped so far experimental detection in small systems and how such effects could be revealed in future experiments. To address this, our work develops and documents parton energy loss models that extend to the smallest hadronic collision systems.
Parton energy loss in the QCD medium was predicted in the pioneering works of Bjorken [20] and of Gyulassy, Pluemer and Wang [21,22]. It was given a first QCDbased treatment by Baier, Dokshitzer, Mueller, Peigné and Schifff (BDMPS) [23,24], and by Zakharov (Z) [25,26], with later refinements by others [27][28][29]. These works calculate, for an arbitrary number of interactions with the medium, the non-abelian Landau-Pomeranchuk-Migdal (LPM) effect that underlies medium induced parton splitting. The same LPM effect was found independently by Arnold, Moore and Yaffe (AMY) when developing an effective kinetic transport formulation of hard degrees of freedom in QCD finite temperature field theory [30,31]. Spurred by the measurement of quenched jets (as opposed to quenched high-p ⊥ hadrons) at the LHC, much subsequent theoretical work aimed at extending the BDMPS-Z formalism to multi-parton final states, either by encoding jet quenching in Monte Carlo simulations [32][33][34][35][36][37][38] or by extending the BDMPS-Z formalism to higher order in α s and thus to higher number of medium induced gluons in the final state [39][40][41][42].
In the present work we focus on modeling the suppression of high momentum hadron spectra. Our starting point is a particularly clean and simple reformulation of the BDMPS-Z formalism due to Arnold [43] from which we determine the probability distribution of parton energy loss ("quenching weight") and the resulting hadron nuclear modification factor following Ref. [44]. There have been several model comparisons to quenched hadron spectra with the systematic study of the centrality dependence of the nuclear modification factor [45][46][47][48][49][50][51][52]. These works focus on the centrality in PbPb, XeXe (AuAu) collisions at the LHC (at RHIC). Our aim is to validate an energy loss model on this centrality dependence and to use it for predicting nuclear modification factors in the foreseen TeV-scale minimum bias collisions of lighter nuclei, i.e., in oxygen-oxygen (OO) and argon-argon (ArAr) collisions [53]. For the heavy quark nuclear modification in small systems, a similar approach has been followed in Ref. [54]. For the nuclear modification of jets, several studies of pPb at the LHC [55][56][57][58] (see also [59]) arrived at quenching effects that are larger than the current bounds set by experiments. A community-wide study of future physics opportunities for high-density QCD at the LHC [53] asked for further modeling efforts, noting that current Monte Carlo models of parton energy loss [34] may somewhat over predict medium effects in argon and xenon collisions.
Theoretical uncertainties in applying the BDMPS-Z formalism to quenched hadron spectra have been analyzed in a community-wide study [60], and they have been included in subsequent extractions of the jet transport coefficientq from data [47,61]. In addition, there are known event selection and geometry biases that in peripheral AA collisions complicate the model comparison of nuclear modification factors [62]. One qualitative conclusion of the present study will be that an energy loss model based on the BDMPS-Z formalism and consistent with experimental data in PbPb and XeXe collisions can result in sufficiently small nuclear modifications in OO collisions that a high accuracy baseline is needed to detect medium induced energy loss. In our companion paper [63] we show that this is indeed possible.
In Section II, we shall provide a description of different building blocks of a parton energy loss model. We also comment on the system size dependence of theoretical uncertainties. Section III presents our results on momentum and system size dependence of the charged hadron nuclear modification factor. Because our simplified model does not take into account all the details of modeling soft QCD medium evolution in heavy ion collisions, we vary various model assumptions to test the robustness of our predictions. Although it is not the main focus of our paper, we also checked the model predictions for high momentum hadron v 2 (p ⊥ ). Our conclusions are given in Section IV.

II. SIMPLE PARTON ENERGY LOSS MODEL
Most formulations of parton energy loss for single inclusive hadron production start from the framework of collinearly factorized perturbative QCD. In this framework, a generic hadronic cross section can be schematically written as where the perturbatively computable hard partonic (gluon (g) and quark (q)) cross sections σ vac g/q are convoluted with the universal process-independent parton distribution functions (PDFs) that describe the parton content of the hadrons and with the hadronic fragmentation functions (FFs). This starting point provides a systematically improvable baseline for calculating the spectra in the absence of medium effects.
Nuclear effects in Eq. (1) enter in two ways. First, the parton distribution functions in ultra-relativistic colliding nuclei differ characteristically from those in free protons, and hence, the PDFs are replaced by nuclear PDFs (nPDFs) [64][65][66][67][68]. Second, the partons leaving the high-momentum transfer vertex of a nucleus-nucleus collision enter a dense QCD medium that affects their par-ton shower. In the description of single inclusive hadron spectra, this is typically modeled by replacing the hard partonic vacuum cross section by a medium-modified differential parton cross section Here, P g/q ( ) denotes the probability for a gluon (quark) with momentum p ⊥ + to lose of its transverse momentum prior to being convoluted with the fragmentation function.
The nuclear modification of centrality averaged hadron spectra is expressed as the ratio of charged hadron cross sections in nucleus-nucleus (AA) collisions and pp collisions scaled by A 2 , where A is the total number of neutrons and protons in the nucleus The hadron nuclear modification factor is the main deliverable of our simple energy loss model. We work at mid-rapidity |y| < 1 and drop the explicit y-dependence in the following.
In the subsequent sections we describe in detail different model assumptions entering P g/q ( ) and how Eq. (3) is computed in the presence of medium modifications.

A. Medium induced gluon radiation
Inelastic processes provide the most efficient mechanism for degrading the energy of high momentum partons. In models of radiative parton energy loss, these are described by calculating the medium induced gluon emission rate dI g/q med /dω [21][22][23][24][25][26][27][28][29]. Following Ref. [44], the probability P g/q ( ) is given as a sum over the probability to emit n medium-induced bremsstrahlung gluons = n i=1 ω i , The factorial accounts for an arbitrary ordering of the emissions and the exponential normalizes the distribution to ∞ 0 d P ( ) = 1. Here, we use for the evaluation of the medium induced gluon emission rate a particularly clean and transparent reformulation of the BDMPS-Z formalism due to Arnold [43]. For a high-energy parton of species s with energy E moving through a medium, we write [43] ω where x is the momentum fraction carried by the emitted gluon, and s = g/q denotes the species of the emitting parton. In the vacuum, this gluon emission is dictated by the DGLAP vacuum splitting function P s→g . The factor ln |c(0)| determines to what extent the gluon emission rate dI s in the medium differs from that in the vacuum. The entire BDMPS-Z formalism can be reduced to the problem of determining |c(0)| from the function c(t), which satisfies the differential equation [43] with the boundary condition that c(t) → 1 and c (t) → 0 for t → ∞. Here, the complex frequency ω 0 (t) is given in the small x 1 limit by where ω = xE is the energy of the radiated gluon. For small x we have xP s→g (x) ≈ C s . All information about the interaction with the QCD medium enters the formalism via the quenching parameterq in Eq. (7). This parameter, multiplied by the Casimir C s of the corresponding color representation of the energetic parton, characterizes the average transverse momentum squaredq = C sq that is transferred due to soft interactions from the QCD medium to the energetic parton per unit path length. To leading order in the weak coupling expansion,q is independent of the particle species. It depends in general on the local density that the medium has at time t at position x(t), where x(t) is the trajectory of the hard parton through the medium. In this way, information about the density of the soft QCD medium and its time evolution enters the calculation of modified high-p ⊥ hadron spectra.

B. Background temperature parametrization
Many sophisticated hydrodynamic models exist for the evolution of the bulk QCD medium that have been validated phenomenologically against soft physics data in central and semi-peripheral collisions. In principle, any of these models could be interfaced with the present formalism via a simple prescription that determinesq(t, x(t)) from the soft bulk quantities evolved. However, in very peripheral collisions of 90% centrality and light-ion collisions (with number of participant nucleons N part ∼ 10) the assumptions about the fluid dynamic evolution of QCD matter may become more questionable.
Without entering a detailed discussion about the system size dependence of the soft physics modeling [69], we employ a particularly simple setup of the QCD medium evolution in which the system size dependence is given in terms of a few parameters. We will subsequently vary the background evolution to gain insight into the robustness of the parton energy loss signal. For background temperature evolution T (τ, x ⊥ ) we use a one-parameter (opacityγ) solution of a conformal kinetic theory in relaxation time approximation that interpolates between free-streamingγ = 0 and perfect fluidityγ = ∞ [70]. The spatio-temporial temperature profile is given by whereT is a scale invariant solution of the kinetic theory and dimensionful constants T * and R define the temperature normalization and radial size of the system. For different centrality classes and collision systems the radius R is calculated from the entropy density profile s(x ⊥ ), which we obtain from the TrENTo initial state model [71] Furthermore, we fixed the temperature normalization T * to reproduce the centrality dependence of the total entropy dS/dy = d 2 x ⊥ s(x ⊥ ), i.e., As a reference value, we choose to set the temperature at the origin in 0-10% PbPb collisions at time We note that none of the predictions of our models depend on the specific choice of T (τ ref , 0) as it can be reabsorbed in the quenching parameterq. The θ-function in Eq. (8) implements the model assumption that the medium modifications of hard partons cease at freeze-out at T F = 175 MeV. We include interactions between hard partons and the medium for τ > τ 0 = 0.05 fm/c. Kinetic solutionT is given for times τ 0.06R, so if needed the temperature is back-extrapolated to τ 0 = 0.05 fm/c using τ −1/3 scaling. The centrality dependencies of dS/dy and R are tabulated in the Appendix A. We choose the kinetic theory solution with an opacityγ = 16 which corresponds to an almost perfect (η/s ≈ 1/4π) fluid in central √ s N N = 5.02 TeV PbPb collisions [70]. We compare this fluid limit to the case of free-streaming (opacitŷ γ = 0). In addition to the azimuthally symmetric profile Eq. (8), we model the elliptical deformation of the background profile in off-central nucleus-nucleus collisions. This is achieved by adding a linearized kinetic theory solution of an elliptic background perturbation [70]. The magnitude of such deformation is fixed by the eccentricity in the initial conditions (see Appendix A).
The above formulation of background evolution clearly aims at simplicity rather than completeness. However, we checked by drastically changing the temperature evolution in Eq.  course, this does not mean that other observables are not sensitive to these details (see Sec. III C), but we leave a more refined description of the background evolution to future works.

C. Embedding hard partons in medium
The quenching parameterq is determined by the temperature profile along the trajectory x(t) of a particular particleq Here, the proportionality factor d is a model parameter that will be adjusted to reproduce the medium induced suppression of single inclusive hadron spectra in central PbPb collisions at p ⊥ ∼ 50 GeV (we keep α s = 0.3 constant in Eq. (5)). It is Eq. (11) that relates the modeling of the QCD evolution and the geometrical embedding of parton trajectories in that medium to the actual dynamics of parton energy loss. Hard partons are assumed to be produced in binary scatterings and to follow eikonal trajectories in the plane transverse to the beam For boost invariant medium evolution we can always find such a frame. The distribution of production vertices x 0 is set to reproduce the (hard) RMS radius R h of binary nucleus-nucleus collisions obtained from the product of the nuclear thickness functions of the two nuclei in the TrENTo model (see Appendix A). We discretize the velocity angle and initial radial location of the hard particles as shown in Fig. 1. A linear grid in radial coordinate ρ with leads to a Gaussian distribution of hard particles in the physical r coordinate. The values of R h and N coll are documented in the Appendix A. For each collision system and centrality, the nuclear modification factor Eq. (3) is obtained by averaging the energy loss of hard partons over the ensemble of starting locations and velocities shown in Fig. 1. We obtain minimum bias results by taking the N coll -weighted average over 10 centrality classes.

D. Vacuum parton and hadron spectra
In the absence of parton energy loss, the single inclusive hadron (parton) spectra can be calculated in collinearly factorized perturbative QCD according to Eq. (1). For the proton reference spectrum, we take PDFs provided by CT14 [72] and for oxygen and lead nuclei we use nPDFs derived from EPPS16 global fit [64]. We convolute the PDFs with LO QCD scattering matrix elements to produce the vacuum spectra dσ vac g/q of quarks and gluons (for the nuclear modification factor, the difference between LO and NLO results is negligible [63]). The charged hadron cross section is obtained from the partonic one by the convolution with the quark and gluon fragmentation functions D g/q h provided by NNFF1.1h [73] dσ h,vac g/q dp 2 where z is the momentum fraction of the parton that is carried by the leading hadron. We use LHAPDF6 interpolator for evaluating PDFs and FFs [74]. Details of the computation are summarized in the Appendix B.
In Fig. 2 we show the ratio of quark and gluon fragmentation contributions to the inclusive charged hadron (parton) cross section at √ s N N = 5.02 TeV for different collision systems, i.e.
Although gluons dominate the partonic spectra at momenta up to p ⊥ ≈ 300 GeV, they fragment to softer hadrons than quarks and therefore the hadron spectrum is dominated by quark fragmentation already at p ⊥ > 40 GeV. r(p ⊥ ) does not change significantly between pp and AA collisions (in the absence of energy loss), although the nPDF modifies the absolute yields. We computed such "vacuum" nuclear modification factor for hadrons and partons in OO and PbPb collisions, see Fig. 3. We emphasize that here we take the central values of nPDFs [64]. Within current nPDFs uncertainties, the modifications shown in Fig. 3 are consistent with zero for most of the kinematic range. Taking into account such uncertainties (and constraining them with further data) is crucial for disentangling the different sources of nuclear modification in comparison to experimental data. We address this question in detail in our companion paper [63], so we will not discuss nPDF uncertainties further here.
We see that nPDF effects become smaller with decreasing A. We find empirically that the nPDF contribution to the nuclear modification scales well with ( N part − 2) 1/4 , where N part is the average number of participant nucleons. As nPDF effects are expected to be smaller in peripheral collisions [75], we use our empirical scaling to estimate the nPDF effects in centrality selected events. For each centrality class we take this factor to be where k = 0.25 is a normalization such that for PbPb the N coll -weighted centrality average reproduces the minimum bias nuclear modification factor.

E. System size dependence of parton energy loss
For any generic quenching parameter Eq. (11) associated to a particular parton trajectory Eq. (12) through a QCD medium of given temperature profile Eq. (8), we can solve numerically the differential equation Eq. (6) and we can thus determine the medium-modified gluon energy distribution ω where The resulting emission rate is 2. For large (opaque) slowly varying systems with |ω 0 (t)| |ω 2 0 (t)|, Eq. (6) can be solved using adiabatic approximation The solution is for which Comparing the parametric estimates Eq. (21) and Eq. (22), one finds that the crossover between these two limiting cases occurs at a frequency ω kink In the upper panel of Fig. 4 we illustrate the characteristic interpolation between the non-abelian Landau-Pomeranchuk-Migdal (LPM) ω −1/2 -powerlaw of Eq. (22) in the limit of soft gluon energies, ω ω kink and the ω −2 powerlaw Eq. (21) of the opacity expansion for ω ω kink . As the integrals in Eq. (24) depend on the in-medium path length and the density of the system, ω kink depends on the QCD medium produced in the collision and is larger for systems of larger geometrical extent and/or larger density (see caption of Fig. 4 for numerical details).

F. Quenching of hadron spectrum
Having calculated for each trajectory Eq. (12) the medium-induced gluon rate ω dI med dω as illustrated in the upper panel of Fig. 4, we determine the corresponding probability P ( ) of parton energy loss in Eq. (4). To characterize the impact of parton energy loss, we consider the ratio of partonic medium modified and vacuum cross sections, i.e., the quenching factor [44] For 1 We employ an alternative rewrite of the Taylor series  (23) respectively. The frequency where the two asymptotic rates are equal defines ω kink that is for these systems approximately at ω PbPb kink ≈ 6.6 GeV, ω OO kink ≈ 1.0 GeV and ω pPb kink ≈ 0.7 GeV. (Here we have chosenq/T 3 = 2.46.). (botom panel) Integrand of the shift function Eq. (33). Area under the curves represents contributions to Sg by gluon emission at different energies scales for final hadron with p ⊥ = 100 GeV (thin line-p ⊥ = 50 GeV) and zn ≈ 3. The vertical lines correspond to ω kink , the shaded region to ω < 500 MeV.
where n g/q (p ⊥ ) is the spectral index Note that partonic spectra are falling steeply with n 5 in the kinetic regime 20 GeV < p ⊥ < 1000 GeV relevant for our study, see Fig. 5.
In close analogy to Eq. (25), we define also the suppression of charged hadrons due to parton energy loss by the ratio where σ vac h is the single inclusive charged hadron cross section in vacuum, and σ med h is the corresponding quan- tity with medium included modifications. Fragmented hadrons are produced at softer momenta, which results in the partonic cross section effectively changing momentum by a factor of z ≈ 0.5, i.e., σ vac g/q (p ⊥ )/dp 2 ⊥ ≈ σ vac h,g/q (0.5p ⊥ )/dp 2 ⊥ . Analogously we can write the result in exponential form Eq. (27) with reduced exponent 2 (30) where typically zn ≈ 3. In Fig. 5 we display the momentum dependence of zn g/q for hadrons produced by quark and gluon fragmentation.
The exponential form of Eq. (27) allows for a particularly simple evaluation of the integral over the probability distribution Eq. (4). For large hadron momentum the medium modification of the hadron spectra is pro- 2 The approximation amounts to assuming z n = z n . By doing the fragmentation of the quenched partonic spectra in Eq. (29) directly, we verified that this does not qualitatively alter the nuclear modification factor.
portional to the mean energy loss For generic p ⊥ , the result can be expressed with a shift function S g/q (u) as where S g/q (u) denotes the energy loss due to multiple medium-induced gluon emissions [44] S g/q (u) = 1 u log As discussed in Sec. II E, the characteristic emission energy ω dI g/q med dω has a UV cutoff at ω kink , Eq. (24), therefore if uω kink 1 (which is usually the case), the energy loss Eq. (33) becomes proportional to the integral over the gluon emission rate ω dI med dω . From Eq. (23) one finds for the quenching weight the parametric form In the following, the quenching factor will be calculated using the full integral Eq. (33). Finally, the hadron nuclear modification factor can be computed by multiplying the nPDF modification, Eq. (16), with appropriately weighted quark and gluon quenching factors for the hadron spectra where r(p ⊥ ) ratio is given by Eq. (15).

G. Model applicability in small collision systems
Parton energy loss models have been applied so far to relatively large collision systems. Here we ask whether the parametric range of applicability of the parton energy loss model extends to smaller systems like inclusive OO collisions or even pPb collisions.
The BDMPS-Z formalism was developed for the emission of sufficiently energetic gluons to which a perturbative reasoning applies. To establish to what extent this condition is met in our model calculations, we show in the lower panel of Fig. 4  depends only weakly on p ⊥ in the kinematical range of phenomenological interest, and the scale ω kink is seen to characterize the peak of the integrand for all collision systems. In calculations we consistently assumed ω p ⊥ , which is approximately fulfilled for p ⊥ > 50 GeV in the largest collision systems and holds for much lower momentum in smaller systems. The characteristic energy of medium induced gluon radiation ω kink decreases with decreasing density and geometric extent of the system, and the integral Eq. (33) receives an increasing contribution from very soft gluon emission for which the validity of our model becomes questionable. We note however that the extrapolation to small systems shown in Fig. 4 is smooth and roughly half of the computed energy loss can be attributed to radiation with ω 1 GeV for OO collisions. With these considerations we take a pragmatic approach of basing a first exploratory study of the systems size dependence of parton energy loss on a BDMPS-Z formalism that is not modified with additional assumptions for small systems.
We mention as an aside that we have performed other consistency checks of our model setup. In particular, the discussion above assumed x 1. We checked that relaxing this approximation has only mild effects on the results in Fig. 4 (data not shown). Within the model uncertainties quoted in the present work, these are negligible, and we do not discuss them further. We also checked that the phenomenological practice of mapping parton energy loss of a dynamically evolving QCD medium onto a parton energy loss calculation for a static brick of suitably chosen quenching parameter describes, over the entire ωrange, the energy loss curve in Fig. 4 within 5% accuracy. We do not employ this observation to simplify our calculation, but we note it here since it implies that our results could be reproduced in other existing approaches.

III. RESULTS
We now compare results of the parton energy loss model described above to the measured centrality and momentum dependence of the charged hadron nuclear modification factor R h AA (p ⊥ ) in PbPb and XeXe collisions at the LHC. We then predict the nuclear modification factors in minimum bias pPb, OO and ArAr collisions, and centrality selected OO collisions. We test the robustness of these results by varying model assumptions. Finally, we discuss to what extent parton energy loss can account for the observed azimuthal momentum anisotropy v 2 (p ⊥ ) at sufficiently high transverse momentum within our setup.

A. System size and momentum dependence
If the temperature profile of the QCD medium is fixed, the only remaining unconstrained parameter of the parton energy loss model of Section II is the proportionality factor d that sets the value of the quenching parameterq in units of T 3 in Eq. (11). We adjust d such that the model reproduces the measured centrality averaged hadron nuclear modification factor Once the overall normalization ofq is thus fixed, the p ⊥ -dependence of R h AA (p ⊥ ), its dependence on centrality, and its dependence on the nucleon number A in centrality averaged collisions are model predictions. Fig. 6 shows that the model describes well the observed p ⊥dependence in centrality averaged PbPb and XeXe collisions. Here the error bands account only for the above In Fig. 6 we also compare the same model to measurements of the nuclear modification factor in pPb collisions. At p ⊥ ∼ O(100) GeV the model predicts a slight enhancement of R h pPb indicating that the nuclear modification of the PDFs in the anti-shadowing region is numerically more important than the small parton energy loss [64]. We note that within current theoretical and experimental uncertainties no firm statement about the discrepancy between data and model predictions for pPb shown in Fig. 6 can be made.
Up to now we followed the standard assumption that parton energy loss is negligible in pp collisions. To check the internal consistency of our model we estimated the expected energy loss in pp collisions. The yellow band in Fig. 6 shows the ratio of hadron spectra with and without the medium effects. In light of other model uncertainties, this assumption seems justified.
In Fig. 7 we show how the nuclear modification factor in centrality averaged OO collisions evolves from √ s N N = 5.02 TeV to √ s NN = 7 TeV-the projected center-of-mass energy of the upcoming OO run at the LHC [53]. The effect of changing collision energy is twofold. First, an increase in √ s N N shifts the nPDF effects to higher transverse momentum. Second, the soft medium produced in the collision also depends on the collision energy. Here, we model this by assuming T * ∝ s 0.05 N N in Eq. (8), which is motivated by the charged particle multiplicity dependence on center-of-mass energy [78]. In Fig. 8 (Fig. 9) we compare the p ⊥ -and centrality dependence of the charged hadron nuclear modification fac-tor in our model and measured data at √ s N N = 5.02 TeV PbPb ( √ s N N = 5.44 TeV XeXe) collisions. The p ⊥ dependence of R h AA (p ⊥ ) mainly stems from the steeply falling particle spectra, while the centrality dependence is driven by the in-medium path-length, see Eq. (34). As seen in Fig. 8 and Fig. 9, the model reproduces without any parameter adjustment both the p ⊥ -and centrality dependence of R h AA between 0% and 70%. At very high p ⊥ the fractional energy lost by the parton is small and R h AA is dominated by nPDF effects. We note that systematic normalization uncertainties in the experimental data are shown by blue (green) boxes, which increase to ∼ 15% (∼ 30%) in the most peripheral bin. If these are taken into account, the tension between data and model results visible in the 70-90% (70-80%) centrality bin lies within the 2-σ uncertainty band. We note however that no parton energy loss model of BDMPS-Z type contains physics that could account for a stagnation or an increase of the suppression as the system size and the energy density reduces from the 50-70% to the 70-90% (70-80%) centrality bin.
We note that our model predictions of minimum bias inclusive nuclear modification factors in OO collisions addresses the same N part ∼ 10 range as 70-90% (70-80%) peripheral PbPb (XeXe) collisions. Measuring R h AA in OO collisions is a much wanted independent test of the expected system size dependence of parton energy loss, that is free of assumptions about the modeling of the soft physics that enter the baseline of peripheral R h AA measurements. We scrutinize the potential of discovering energy loss in small systems in our companion paper [63].

B. Robustness of model predictions
In the previous section we showed the results of a simple energy loss model based on the BDMPS-Z energy loss formula of hard partons in a kinetically evolved background. The system size dependence was modeled by TrENTo initial state model and we included nPDF and fragmentation effects. Although this setup is well motivated, many of the model assumptions have not been independently constrained. Therefore we now stress-test the robustness of model predictions by varying different model assumptions in the same framework.

Summary of models considered
First, to understand the relative importance of nPDF, system size modeling and fragmentation effects on our results, we consider four unphysical setups with some of these model components switched off.
Minimal: In this minimal implementation, the isotropic background geometry is scaled according to optical Glauber R and N part . Energy loss is modeled only for gluons and no nPDF or fragmentation effects are included. In essence, the gluon quenching factor Eq. (27) with n g = 6 is used as a proxy for R h AA .
Anisotropic: The same as Minimal, but the system size dependence of R and N part is now modeled using the TrENTo initial state model and we include the average elliptic deformation of the background.
nPDF: The same as Minimal, but nPDF effects are included. That is the (partonic) gluon quenching factor is multiplied by "vacuum" (partonic) R vac AA shown in Fig. 3. The centrality dependence of nPDF effects is scaled with ∼ N part − 2 1/4 , see Eq. (17). No fragmentation is included.
Fragmentation: The same as Minimal, but gluons are fragmented into hadrons, i.e., the hadronic quenching factor Eq. (32) is compared to R h AA . No nPDF effects are included.
Next we study how our results depend on the assumed background medium evolution. As explained in Sec. II B, by default we use a particular simple parametrized temperature profile. Here, we test to what extent our predictions dependent on this evolution. In all cases we include both nPDF and fragmentation effects.
Simple: This is our default model described in Sec. II and with the results shown in Sec. III A. It includes geometry scaling based on TrENTo, nPDF effects, fragmentation and energy loss for both quarks and gluons.
Simple τ 0 = 0.5 fm/c: The same as Simple, but the energy loss is calculated from the later starting time of τ 0 = 0.5 fm/c instead of τ 0 = 0.05 fm/c.
Simple T F = 120 MeV: The same as Simple, but energy loss is computed until a later time, namely until the temperature falls below T F = 120 MeV instead of T F = 175 MeV.
Lattice EOS: The same as Simple, but the temperature profile is determined using lattice equation of state T lat (e) [79], where e ≈ 15T 4 is the energy density in our (conformal) kinetic simulation. The freeze-out temperature is again set to T F = 120 MeV.
Bjorken: The same as Simple, but the kinetic temperature evolution Eq. (8) is replaced by Bjorken scaling Free streaming: The same as Simple, but with the free streaming (γ = 0) solution of kinetic theory for an azimuthally symmetric initial profile.
All model variations above used the parton energy loss formula derived by Arnold [43] in BDMPS-Z formalism. Here we use our simple framework to compare three characteristically different parametrizations of parton energy loss inspired by recent phenomenological studies [48][49][50][51]. We calculate the shift function S g/q for these parametrizations with free normalization constant κ.
A: Energy loss with weak path length and temperature dependence dE/dL ∼ −L 0.4 T 1.2 , leading to B: Energy loss with linear path length dependence and strong temperature dependence, dE/dL ∼ −LT 3 , leading to C: Energy loss implementing stopping with where τ stop = 1 2(κ/5) p

Discussion
In Table I we summarize the model variations introduced above. For each model we adjust the free parameter d =q/T 3 or κ to reproduce the centrality averaged R h AA at p ⊥ = 54.4 GeV in √ s N N = 5.02 TeV PbPb collisions as it was done in Sec. III A. We then compare these models in Fig. 11 to the centrality dependence of charged hadron nuclear modification factors measured in PbPb and XeXe collisions, and we extrapolate to OO collisions at √ s N N = 5.02 TeV.
Before entering a more detailed discussion, let us note that despite the dramatic approximations implemented in the different models in Table I, most of the models reproduce the p ⊥ dependence of R h AA in central and semicentral PbPb and XeXe collisions. They do so with values ofq/T 3 or κ that vary significantly with model assumptions. However the aim of the present paper is only to estimate the expected signal of parton energy loss in light-ion collisions. We can do this extrapolation without judging the completeness of the different model scenarios or the numerical value of the extracted medium parameterq/T 3 or κ.
Now we discuss individual model variations listed in Table I. In the first four, Minimal, Anisotropic, nPDF and Fragmentation, some of the model components were switched off. The spread of model predictions in Fig. 11 (dotted lines) informs us to what extent the detailed modeling of these components is important for system size extrapolations. Moreover, as fragmentation converts partons to much softer hadrons, the same R h AA is achieved with three times larger value ofq/T 3 . We note in addition that doing fragmentation directly of the quenched parton spectra in Eq. (29) instead of using Eq. (30) increasesq/T 3 by ∼ 20%.
Next we considered the parton energy loss dependence on the variations of the background temperature evolution (dashed lines). Starting energy loss at 0.5 fm/c requires a twice larger value ofq/T 3 than any other model variations in Table I. This model scenario shows also a more pronounced tension with experimental data in the mid-central PbPb and XeXe data. This suggests that data favors early onset of energy loss. Other variations of the temperature evolution-such as varying from Bjorken to free-streaming, extending the interaction down to T F = 120 MeV, or switching to lattice EOS, see Fig. 10-seem to have only a mild effect on R h AA . We finally consider parton energy loss formulas that differ significantly from BDMPS-Z (solid lines). Here, the formula assuming full stopping (C) is arguably the  Fig. 6).  Table I  most extreme choice, and it is the one that shows the most significant tension with the observed centrality dependence in PbPb and XeXe collisions. We therefore do not include it in our extrapolation to OO. The other two parametrizations (A and B) are comparable to our Simple model.
Given the range of model assumptions explored, we regard the envelope of the different predictions in Fig. 11 as a realistic theory uncertainty for R h AA in OO collisions at √ s NN = 5.02 TeV. In our companion paper [63], we  Table I for minimum bias pPb collisions. Data points are the same as in Fig. 6.
use the same range of model scenarios to compute the expected parton energy loss signal and its uncertainty for the proposed √ s NN = 7 TeV OO collisions.
For completeness we show in Fig. 12 results for the same set of model variations applied to minimum bias pPb. As there are no mechanisms in the considered models (other than nPDF effects) to produce larger than unity nuclear modification, none of the models go through the experimental data points.

C. High momentum hadron anisotropy
A more differential probe of parton energy loss is the high momentum anisotropy of the final particles. In a peripheral collision with elliptical shape parton energy loss is expected to depend on the orientation of the hard parton trajectory. This dependence can be parametrized as cos(2φ) modulation of the nuclear modification factor where φ is the azimuthal momentum angle and φ 2 characterizes the event-plane. Experimentally, v 2 (p ⊥ ) is obtained from the correlation between high-p ⊥ hadrons and soft particles.
It has been long a challenge to simultaneously describe the nuclear modification factor and the sizable high momentum anisotropy within the same model. Models that do not include early time parton energy loss typically fare better [47], because they concentrate the energy loss at later times where the background anisotropy is more relevant. Moreover, it has been shown that including eventby-event fluctuations of the underlying medium can increase the high-p T elliptic flow [48].
Our simple framework does not model event-by-event fluctuations of soft particle production and therefore we do not expect it to accurately reproduce the experimentally measured v 2 (p ⊥ ). Nevertheless, it is interesting to check how different model assumptions in Sec. III B affect the elliptic flow of high-p ⊥ . We determine v 2 (p ⊥ ) from the energy loss modulation of Eq. (39) with respect to the inputted background deformation. In Fig. 13 we compare our model predictions of v 2 (p ⊥ ) to data in different centrality bins of √ s N N = 5.02 TeV PbPb collisions.
We see that for the most of model scenarios v 2 (p ⊥ ) is underpredicted by a factor of ∼ 2. Possible exceptions are the scenario with Bjorken temperature profile and energy loss model with stopping C. The slow temperature evolution in Bjorken and the concentration of energy loss towards the end of the evolution in model C presumably allow for stronger correlation between initial state geometry and high-p ⊥ energy loss. Finally in Fig. 14 we show the Simple model predictions for v 2 (p ⊥ ) in small collisions systems, i.e., centrality selected OO collisions and minimum bias pp and pPb. The tendency of our model to underpredict the experimental data prevents us from making quantitative conclusions about high-p ⊥ elliptic flow in small systems. However we can make the following qualitative observations. First the large initial eccentricity in central OO collisions [54,81] results in monotonically decreasing v 2 (p ⊥ ) with centrality. Secondly, we find a small elliptic flow in minimum bias pPb collisions and even smaller in pp. Making more quantitative statements about the elliptic flow magnitude in small systems is outside the scope of the current manuscript.

IV. CONCLUSIONS
In the present paper we document a model for calculating the high momentum charged hadron spectra modifications due to the medium induced parton energy loss in small collision systems. Our baseline calculation of hadron spectra consists of the leading order QCD partonic cross sections convoluted with (nuclear modified) parton distribution functions and fragmentation functions. The parton energy loss is modeled by small-x gluon emission and the dynamical temperature profile is scaled to match the expected system size and entropy.
After tuning a single model parameter to a single data point of the charged hadron nuclear modification factor at p ⊥ ∼ 50 GeV in minimum bias PbPb collisions, we demonstrated that our model is consistent-up to a 2-σ tension in the most peripheral bin-with the p ⊥ and centrality dependence of R h AA in √ s N N = 5.02 TeV PbPb and √ s N N = 5.44 TeV XeXe collisions. Validated against these data, the model provides well motivated predictions for the charged hadron nuclear modification factors in the minimum bias pPb, OO, ArAr and in centrality selected OO collisions.
To ascertain the systematic uncertainties we varied the different model components, medium evolution and the energy loss formula. All modeling scenarios provide rather comparable momentum and system size dependencies of R h AA once fitted to the same point in the minimum bias PbPb collisions. These model variations predict up to ∼ 15% modification of hadron spectra in minimum bias OO collisions at p ⊥ ∼ 50 GeV. Such small nuclear modification could not be resolved within the systematic experimental uncertainties present in the comparable size peripheral PbPb or XeXe collisions. However, a measurement of R h AA in an inclusive OO collisions is free of model dependent uncertainties entering the centrality selected nuclear modification factor.
The ability to identify parton energy loss also depends on the accuracy with which the baseline without the medium effects can be calculated. In the companion pa-per [63], we show that the accuracy of the baseline R h,vac AA in an inclusive OO collisions is known with sufficient precision that the discovery of the medium induced parton energy loss in small systems with N part ∼ 10 is possible.
Parton energy loss is sensitive to the spatio-temporal extension of the QCD medium and its density profile. In the main text, we have described the physical picture underlying our modelling of the collision geometry. For completeness, we provide in this appendix quantitative information.
The simplest way to determine the initial geometry for PbPb, XeXe and OO collisions is to use standard optical Glauber model [82]. In Table II we present, for each collision system, the computed number of participants N part , the number of binary collisions N coll , the radius of the profile R, and the RMS radius R h as a function of centrality defined by the impact parameter b. The more sophisticated way to determine the initial geometry for each collision system is to use the TrENTo initial condition framework [71]. In the TrENTo model, the initial transverse entropy density profile is computed from where the parameter p controls the mixing of fluctuating thickness functions T A and T B . In this work, we use the following parameter values [83] to obtain the entropy density for each collision system: • reduced thickness parameter p = 0.013 • fluctuation parameter k = 0.93 • nucleon width σ = 0.6 • inelastic nucleon-nucleon cross section σ N N = 64 mb For all elements we used the standard settings in TrENTo, except for oxygen where for the nucleon positions we used the tables from [84] (see also [81]), as provided in [85]. We take an ensemble of 20 000 events and for each centrality (defined as a class of events ±5% from the midpoint value) we obtain an average of all values used in the main text, which is the radius of the entropy density (Table III), the average entropy density (Table IV, used in Eq. (10)), the radius of the hard parton scattering centers R h (Table V), the number of participating nucleons N part (Table VI, used for nPDF corrections), the eccentricity 2 (Table VII) and finally N coll (Table VIII, used for weighting centrality classes).

Appendix B: Parton and hadron production
In this section we summarize the leading order (LO) computations of inclusive parton and hadron cross sections including the discussion of parton distribution and fragmentation functions.

Single inclusive parton cross section
At LO in the strong coupling α s , the production of jets in the collisions of hadrons A and B with momentum P A and P B is given by the partonic 2 → 2 QCD scattering process where two incoming partons a and b are sampled from parton distribution functions and the scattered partons c and d can be identified as final state jets. According to the factorization theorem, the total twojet cross section may be written as where we factored out the coupling constant |M| 2 = (4πα s (µ 2 R )) 2 |M| 2 . Here, the coupling constant is evaluated at the renormalization scale µ R (for partonic cross section we take µ R = µ F = p ⊥ , where p ⊥ is the transverse parton momentum).
The single inclusive jet cross section at LO is then given as dσ AB j dy inc p ⊥ dp ⊥ = 1 8πs 2 where cosh y inc < √ s 2p ⊥ and the integration limits for y min < y < y max are given by, − log √ s p ⊥ − e −yinc < y < log √ s p ⊥ − e yinc . (B10)

Single inclusive hadron spectra
The single inclusive hadronic cross section at LO in the absence of medium modifications are given by the convolution of the jet spectrum Eq. (B9) with the fragmentation function D k h : where the cross section for producing q ⊥ momentum parton k is convoluted with the probability to fragment to momentum p ⊥ = zq ⊥ charged hadron.
Performing the integration over q ⊥ and inserting the partonic cross section formula, the invariant hadron spectra may be rewritten as In the expression above, the momentum fractions x A and x B appearing in Eq. (B9) are evaluated at the rescaled momentum p ⊥ → p ⊥ /z and z min = 2p ⊥ √ s cosh y inc . The gluon and the (averaged) quark fragmentation functions are taken from NNFF1.1h [73]. The renormalization and factorization scales are taken to be µ R = µ F = p ⊥ , where p ⊥ is the transverse hadron momentum.