Inconsistencies in, and short pathlength correction to, R AA ( p T ) in A + A and p + A collisions

. We present the first leading hadron suppression predictions in Pb + Pb and p + Pb collisions from a convolved radiative and collisional energy loss model in which partons propagate through a realistic background and in which the radiative energy loss receives a short pathlength correction. We find that the short pathlength correction is small for D and B meson R AA ( p T ) in both Pb + Pb and p + Pb collisions. However the short pathlength correction leads to a surprisingly large reduction in suppression for π mesons in p + Pb and even Pb + Pb collisions. We systematically check the consistency of the assumptions used in the radiative energy loss derivation—such as collinearity, softness, and large formation time—with the final numerical model. While collinearity and softness are self-consistently satisfied in the final numerics, we find that the large formation time approximation breaks down at modest to high momenta p T ≳ 30 GeV. We find that both the size of the small pathlength correction to R AA ( p T ) and the p T at which the large formation time assumption breaks down are acutely sensitive to the chosen distribution of scattering centers in the plasma.


Introduction
The modification of the spectrum of high transverse momentum (high-p T ) particles is one of the key observables used to understand the non-trivial, emergent, many-body dynamics of quantum chromodynamics (QCD) in highenergy collisions [1][2][3][4].One of the most important findings of the Relativistic Heavy Ion Collider (RHIC) was a roughly factor of five suppression of leading light hadrons with p T ≳ 5 GeV/c in central Au + Au collisions [5,6].This suppression, equal for pions and eta mesons [7], along with null controls of qualitatively no suppression of the weakly-coupled photons in Au + Au collisions [8] as well as of leading hadrons in d + Au collisions [6,9], clearly demonstrated that the suppression of leading hadrons in central collisions is due to final state energy loss of the high-p T partons interacting with the quark-gluon plasma (QGP) generated in the heavy ion collisions (HIC).Models of leading hadron suppression based on final state energy loss derivations using perturbative QCD (pQCD) methods qualitatively describe a wealth of these high-p T RHIC data [10][11][12].
One of the other major findings of RHIC was the near perfect fluidity of the strongly-coupled low momentum modes of the QGP formed in semi-central nucleusnucleus collisions as inferred by the remarkable agreement of sophisticated, relativistic, viscous hydrodynamics models with the spectra of measured hadrons with p T ≲ 2 GeV [13][14][15].
The data from the Large Hadron Collider (LHC) has been no less impressive.Of extraordinary importance have been the signs that the non-trivial, emergent, many-body QCD behavior associated with QGP formation in central and semi-central heavy ion collisions at RHIC and LHC are also observed in small collision systems such as p + p and p + A for large final measured multiplicity.For example, strangeness enhancement [16,17] and quarkonium suppression [18] appear to depend only on multiplicity but not collision system.And the same sophisticated, relativistic, viscous hydrodynamics models [19] also describe the spectra of measured hadrons in high-multiplicity p+p and p + A collisions [20,21].
One may thus conclude that small droplets of QGP form even in these smallest of collision systems at the LHC.If QGP is formed in high-multiplicity collisions of small systems, then high-p T partons should suffer some final state energy loss, as has been observed in large collision systems at RHIC and LHC.(Models already demonstrate the importance of final state energy loss in forward hadron production in cold nuclear matter [22,23].)Experimentally, there are tantalizing signs of the non-trivial modification of high-p T particles in small collision systems [24][25][26].However, there are likely non-trivial correlations between the multiplicity in small collision systems and the presence of high-p T , high-multiplicity jets.For example, these correlations likely impact the initial spectrum of high-p T partons [27] that enters the numerator of the nuclear modification factor, while the minimum bias spectrum in the denominator is unchanged; thus one should be cautious in interpreting a standard R AB measurement.
From the experimental side, it is likely very interesting to consider the ratio of the spectrum of a stronglyinteracting particle with a weakly-interacting particle, each from the same multiplicity class.
From the theoretical side, we may potentially make progress by considering the small collision system predictions of the energy loss models used to describe so well qualitatively the large collision system high-p T particle suppression.One obvious challenge for directly comparing these energy loss models to small collision system data is the assumption made in the energy loss derivations that the high-p T particles travel a large pathlength in the QGP medium.For example, energy loss models built on BDMPS-Z-type energy loss [28][29][30][31][32][33] utilize the central limit theorem [34], and so assume a very large number of collisions occur between the high-p T probe and the QGP medium.Even in large collision systems of characteristic size ∼ 5 fm, since the mean free path for these high-p T particles given by weakly-coupled pQCD is ∼ 1-2 fm [34], the application of the central limit theorem is dubious.
Even for the thin plasma approach of DGLV that naively seems the best suited for modelling the radiative energy loss processes in systems of phenomenologically relevant size, there is an explicit assumption that the partonic pathlength in the QGP medium, L, is large compared to the natural scale set by the Debye screening mass µ, L ≫ 1/µ.In the original derivation of the induced gluon radiation spectrum, contributions from the Gyulassy-Wang potential [35] were dropped as they are O(e −Lµ ).For µ ∼ gT ∼ 0.5 GeV [36], the characteristic size L ∼ 1 fm in high-multiplicity p + p and p + A collisions is not particularly large compared to 1/µ ∼ 0.4 fm.Thus to create an energy loss model to compare to data in these small systems, one needed the small pathlength corrections to the DGLV opacity expansion.
This small pathlength correction to the first order in opacity DGLV radiative energy loss were derived for the first time in [37].For later noting, the derivation in [37] benefited significantly from a simplification due to the large formation time assumption, τ form ≫ 1/µ, an assumption made also in the original DGLV derivations.
The small pathlength correction to the usual DGLV energy loss contained four surprises: due to the LPM effect (interference between the induced radiation and the usual vacuum emissions due to the hard initial scattering), the small pathlength correction reduces the energy lost; the reduction in energy loss is seen in all pathlengths (although the relative importance of the correction decreases with pathlength, as expected); the correction grows linearly with partonic energy (as opposed to the logarithmic growth of the usual DGLV energy loss); and the correction breaks color triviality, with the correction for gluons ∼ 10 times the size of the correction for quarks.
Having derived the correction to the radiative energy loss due to small pathlengths, it is of considerable interest to determine quantitatively the importance of the correction in phenomenological models of high-p T particle suppression.It is the goal of this manuscript to provide just such predictions.We are particularly interested in seeing the quantitative importance of the reduction in energy loss as a function of energy and of collision system: the reduction in energy loss from the short pathlength correction might provide a natural explanation for the surprisingly fast rise in the nuclear modification factor with p T for leading hadrons at LHC [38] and for the enhancement of the nuclear modification factor above unity seen in p + A collisions [39,40].
What we will see is that for light flavor final states at larger energies, p T ≳ 30 GeV/c, the small system "correction" becomes of order of the energy loss itself, which leads us to consider systematically the extent to which energy loss model energy losses are consistent with the various approximations used in the opacity expansion derivation.We will consider in detail the approximations used in the opacity expansion derivation of radiative energy loss [37,41,42].We will find that for the radiated gluons that dominate the energy weighted single inclusive emission distribution, the large formation time assumption is violated for E ≳ 30 GeV/c in large collision systems and for E ≳ 10 GeV/c in small systems, where E is the energy of the radiating parton; which implies the need for yet another derivation of radiative energy loss in the thin plasma limit but with the large formation time assumption relaxed.We also see that the usual WHDG treatment of the average elastic energy loss, with fluctuations given by the fluctuation-dissipation theorem, is not appropriate for small collision systems where the number of scatterings is not large; thus in future work one must implement an elastic energy loss appropriate for small and large collision systems.

Energy loss framework
We wish to make contact with known energy loss model results.We will in particular attempt to, in as reasonable way as possible, make an apples-to-apples comparison of the Wicks-Horowitz-Djordjevic-Gyulassy (WHDG) convolved radiative and collisional energy loss model [43], which has seen such success in describing a breadth of leading hadron suppression data, with an energy loss model with the same elastic energy loss but with a radiative energy loss that includes the short pathlength correction as derived in [37].Let us briefly review the radiative and collisional energy loss setups.

Radiative energy loss
The Djordjevic-Gyulassy-Levai-Vitev (DGLV) opacity expansion [42,44] gives the inclusive differential distribution of radiated gluons from a high-p T parent parton moving through a smooth brick of QGP.The expansion is in the expected number of scatterings or the opacity L/λ g , where L is the length of the QGP brick and λ g is the mean free path of a gluon in the QGP.
The 4-momenta of the radiated gluon, the final hard parton, and the exchanged Debye medium quasiparticle are given respectively in lightfront coordinates (using the same conventions as in [37]) by where M is the mass of the hard parton, m g is the gluon mass, P + is the initial hard parton momentum in the + direction, and x is the radiated momentum fraction.
The DGLV approach makes a number of assumptions related to the physical setup of the problem: -The large pathlength assumption, that L ≫ µ −1 .
-The well separated scattering centers assumption, that λ g ≫ µ −1 .-The eikonal assumption, that P + = E + ≃ 2E is the largest scale in the problem.-The soft radiation assumption, that x ≪ 1.
-The collinear radiation assumption, that k + ≫ k − .-The large formation time assumption, that k 2 /xE + ≪ µ and (k − q 1 ) 2 /xE + ≪ µ 2 + q 2 1 .The DGLV single inclusive gluon radiation spectrum is then [42,45] In Eq. 2 we have made use of the shorthand ω i , and ωm ≡ (m 2 g + M 2 x 2 )/2ω following [37,42].Additionally q i is the transverse momentum of the i th gluon exchanged with the medium; k is the transverse momentum of the radiated gluon; ∆z is the distance between production of the hard parton, and scattering; C R (C A ) is the quadratic Casimir of the hard parton (adjoint) representation (C F = 4/3 [quarks], and C A = 3 [gluons]); and α s is the strong coupling.
The quantity ρ(∆z) is the distribution of scattering centers in ∆z and is defined in terms of the density of scattering centers ρ(∆z) in a static brick, where ∆z is in the direction of propagation, N is the number of scattering centers, A ⊥ is the perpendicular area of the brick, and dz ρ(∆z) = 1.The analysis of realistic collision geometries adds complexity to the scenario, as detailed in Sec.2.7.

Short pathlength correction to DGLV radiative energy loss
The derivation of the modification to the radiative energy loss in the DGLV [41,42] opacity expansion approach with the relaxation of the large pathlength assumption L ≫ µ −1 was considered in [37,46].In the derivation of the short pathlength correction, all assumptions and approximations made in the original GLV and DGLV derivations were kept, except that the short pathlength approximation L ≫ µ −1 was relaxed.The single inclusive radiative gluon distribution, including both the original DGLV contribution as well as the short pathlength correction, is where the first two lines of the above equation are the original DGLV result [42,45], Eq. 2, while the last two lines are the short pathlength correction.We emphasize that contributions from all diagrams which are not suppressed under the relevant assumptions are included.Of particular importance is the large formation time assumption, which allows one to systematically neglect a significant number of diagrams in both the original DGLV derivation [41,42] and in the short pathlength correction [37,46].
Since dN/dx includes an integration over all ∆z, the correction is present for the energy loss of a parton going through plasma of any length; however, the relative contribution of the correction term does go to zero as the pathlength goes to infinity.
The finite pathlength correction originates from not neglecting the q z = iµ 1 pole in the Gyulassy-Wang potential, as was originally done [42,45], which leads to the overall exp(−µ 1 ∆z) scaling of the correction term in Eq. 4 [37,46].
There is a significant literature of energy loss derivations and corrections to earlier energy loss derivations.Even though the focus of this work is the numerical implementation of Eq. 4 and the examination of its underlying assumptions, it is worth taking some time to contextualize the short pathlength correction in Eq. 4 within the literature.In particular, there is currently no other derivation of short pathlength corrections to any energy loss formalism in the literature.
The original BDMPS [28][29][30][31] energy loss derivation explicitly neglects the q z = iµ pole in the Gyulassy-Wang potential.In principle, then, one could derive a short pathlength correction in the original BDMPS-Z formalism analogous to the one derived in [37,46].Subsequent work within the BDMPS-Z formalism [28][29][30][31][32][33] considered the saddle point approximation of the path integral, that in the limit of a large number of scatterings one could make a simple harmonic oscillator approximation (via the Central Limit Theorem).This SHO approximation explicitly requires a large opacity L/λ ≫ 1.For a perturbative calculation, one requires the scattering centers are wellseparated, λ ≫ 1/µ, and so a large opacity implies a large pathlength; one therefore cannot determine a short pathlength correction to the SHO approximated BDMPS-Z approach.If one assumes that the system is strongly coupled (see, e.g., [47,48]) and λ ≪ µ −1 , then all paths are long and there is no short pathlength correction.
In the Improved Opacity Expansion (IOE) [49,50], the starting point is already the z-integrated path integral; i.e. the IOE starts with the equation of motion for the propagator in transverse position space, and is completely insensitive to questions about the interplay between any longitudinal scales such as the pathlength, mean free path, and distance between scattering centers.Within the IOE formalism the "Gyulassy-Wang" potential is taken to be V q 2 ∼ 1/ q 2 + µ 2 2 where, importantly, q is the transverse momentum exchanged with the scattering center; i.e. any potential pole from the z component of the fully three dimensional Gyulassy-Wang potential is already ne-glected.Thus the IOE formalism is unable to compute any short pathlength correction to the energy loss.
Similar to the IOE approach, the finite sizeimprovement [51] to the AMY formalism [52][53][54] begins with Zakharov's path integral formalism and considers only the transverse momentum transfer q from the inmedium scattering cross section.Thus, like in the IOE approach, unless the scattering center cross section is considered in full three dimensions, the finite size-improvement to AMY cannot capture the short pathlength corrections to energy loss.
In [67] the authors couple the opacity expansion to the collective flow of the quark gluon plasma.This work explicitly makes the assumption µ∆z ≫ 1; thus short pathlength corrections to this derivation are possible, but not included in the original derivation.
In [68][69][70] the Gyulassy-Wang static scattering center potential in canonical opacity expansion calculations was replaced by HTL propagators communicating between the high-energy parent parton and the in-medium thermal parton.In the first derivation [68,69], the authors work purely in momentum space to compute the interaction rate; since the authors do not Fourier transform into position space, their result is completely insensitive to the particulars of the path.This work was improved in [70] where a Fourier transform into position space was done and, as is done in the opacity expansion approach, the phases were kept.However the limit L → ∞ is explicitly taken.In principle one could derive a short pathlength correction to this derivation by relaxing the assumption of L → ∞.
The Higher Twist (HT) approach [71,72] in general only keeps the most length-enhanced contributions [73].In principle, one may include less enhanced contributions from the assumed factorized nuclear expectation values of the various four point functions.
Gradient jet tomography [74,75] couples a high momentum parton to an asymmetric medium.The practical implementation of this procedure utilizes the dipole approximation in the path integral approach where the entire nuclear medium is treated as a sheet [76], and so any sensitivity to longitudinal physics is lost.

Numerical implementation of the radiative energy loss
For all numerical calculations we neglect the running of the strong coupling constant, and use α s = 0.3, consistent with [37,42,43].Additionally we use charm and bottom quark masses of m c = 1.2 GeV/c 2 , and m b = 4.75 GeV/c 2 , respectively.The effective light quark and gluon masses are set to the asymptotic one-loop medium induced thermal masses, of m light = µ/2 and m g = µ/ √ 2 respectively [77].The upper bounds on the |k| and |q| integrals are given by k max = 2x(1 − x)E and q max = √ 3Eµ, following [43].This choice of k max guarantees that the momentum of the radiated gluon; and the initial and final momenta of the parent parton are all collinear.
In order to maintain consistency with the WHDG model, we assume the distribution of scattering centers is exponential: The exponential distribution serves to make the integral in ∆z analytically simple.The physical motivation for this is that an exponentially decaying distribution of scattering centers captures the decreasing density of the QGP due to Bjorken expansion [78].It is likely that using an exponential rather than, say, a unit step or power law decay distribution is an overestimate of the effect of the expansion, since Bjorken expansion obeys a power law decay along the incident partons path, not exponential.It turns out that, for the uncorrected DGLV result, the distribution of scattering centers ρ(∆z) affects the characteristic shape of dN g /dx; however we will see that the radiative energy loss is largely insensitive to ρ(∆z) [34].Once the short pathlength correction is included, however, the energy loss becomes far more sensitive to the distribution of scattering centers; particularly at small ∆z [37].
The integral in Eq. 4 can be dramatically simplified if one assumes that the kinematic bound on q, q max = √ 3Eµ, can be taken to infinity.Assuming that q max → ∞ allows one to perform the angular and k integrals analytically, after which one may return the kinematic limit q max = √ 3Eµ.This procedure was done for all previous WHDG predictions [38,43,79] in order to make the numerics simpler; however we performed the full numerical calculation without using this approximation.A full description of this kinematic approximation is provided in [42].
To gain some familiarity with the effect of the short pathlength correction, we calculate the fractional radiative energy loss numerically according to Eq. 4. The mean radiative energy loss is calculated from dN g /dx using The results of evaluating Eq. 5 as a function of the medium length L are shown in Fig. 1.The four surprises mentioned in the Introduction are in evidence here: the short pathlength correction reduces the energy loss; the short pathlength correction does not disappear for L ≫ µ −1 , since all possible distances of scattering ∆z are integrated over in Eq. 5; the short pathlength correction grows linearly with energy; and the effect of the reduction in vacuum radiation is particularly strong for gluons because of the breaking of color triviality.One may also consider the relative size of the short pathlength correction to the radiative energy loss in comparison to the uncorrected radiative energy loss, by calculating the ratio ∆E corrected /∆E DGLV , shown in Fig. 2.
To understand the energy E, length L, and incident Casimir C R dependence of the relative size of the short pathlength correction to the radiative energy loss, we examine the asymptotic dependence of both the DGLV and short pathlength corrected DGLV fractional radiative energy loss.For asymptotic energies the short pathlength corrected fractional energy loss is given by [37] ∆E which was calculated in this asymptotic analysis for simplicity with k max = 2xE and, also for simplicity, an exponential distribution of scattering centers.The equivalent asymptotic result without the short pathlength correction is given by [45] ∆E The relative size of the correction is then given by keeping only leading terms in E. From Eq. 8, the relative size of the short pathlength correction: increases linearly in energy, is C A /C F = 9/4 times larger for gluons in comparison to quarks, and is ∼ 15 times larger for a system with µ = 0.5 GeV and L = 1 fm as opposed to a system with L = 5 fm.One may see in Fig. 2 that the detailed numerics display these three behaviors.

Multi-gluon emission
The DGLV energy loss kernel Eq. 4 gives the inclusive spectrum of emitted gluons.Thus the expected number of gluons can be greater than 1.In fact, one sees that for hard partons emerging from the center of a central heavy ion collision, the expected number of emitted gluons is ∼ 3 [80].To take into account multi-gluon emission we assume that the multiple gluon emissions are independent, following [80].This assumption of independent emissions allows us to convolve the single inclusive gluon emission kernel given by dN g /dx into a Poisson distribution.Explicitly we can write where the P n are found via the convolution and we have P 0 (ϵ) ≡ e −⟨N g ⟩ δ(ϵ).Here, and for the rest of the paper, we define 1 − ϵ as the fraction of initial momentum kept by the parton, that the final energy of the parton in terms of the initial energy of the parton is The Poissonian form of Eq. 9 guarantees the distribution is normalized to one, and the expected number of emitted gluons is n dϵ n P n (ϵ, E) = ⟨N g ⟩.
The bounds on the x n integral are max(0, ϵ − 1) ≤ x n ≤ min(ϵ, 1), which are determined by ensuring that no functions are evaluated outside of their domains.The support of P rad (ϵ) ends when the bounds of the x n integral are equal, i.e. ϵ ∈ (0, 2) is the region of support.The support of P rad (ϵ) past ϵ = 1 is unphysical, and we interpret this as the probability for the parton to lose all of its energy before exiting the plasma.Under this interpretation we put the excess weight We note that the random variable N g should rigorously be thought of as where N g med is the number of radiated gluons occurring due to medium interactions, and N g vac is the number of DGLAP vacuum radiation gluons.With this understanding, the independent gluon emission assumption means that (N g vac + N g med ) and N g vac should each be modeled by a Poisson distribution.Then N g is actually given by a Skellam distribution [81], the difference between two Poisson distributions.In the current model, energy gain relative to the vacuum at some x * corresponds to dN g (x * )/dx < 0; whereas it should rigorously correspond to dN g (−x * )/dx > 0. Following previous work [38,43,80], we will simply model P (ϵ) as a Poisson distribution; the effect of this simplification is not obvious and requires future work.

Elastic energy loss
Elastic energy loss is taken into account using the result derived by Braaten and Thoma (BT) [82], wherein analytic calculations are done in two asymptotic energy regimes.The elastic energy loss of a quark is calculated in the regions E ≪ M 2 /T and E ≫ M 2 /T , where M is the mass of the incident quark, and T is the temperature of the medium.For E ≪ M 2 /T the differential energy loss per unit distance is where B(v) is a smooth function satisfying constraints listed in [82], v is the velocity of the hard parton, and n f is the number of active quark flavors in the plasma (taken to be n f = 2 throughout).For E ≫ M 2 /T the differential energy loss per unit distance is The energy loss at arbitrary incident energy is taken to be the connection of these two asymptotic results such that dE/dz is continuous (determined numerically).It is assumed that: for incident hard gluons the energy loss scales simply by a factor of C A /C F = 9/4; and that there are enough elastic collisions such that the central limit theorem is applicable, following the WHDG model [43].The latter assumption implies that the distribution of elastic energy loss is Gaussian with mean provided by the BT energy loss formula, and width by the fluctuation dissipation theorem [83] where z integrates along the path of the parton, and T (z) is the temperature along the path.Thus where ∆p is found by integrating Eqs.11 and 12 over z, and thus where the additional p i is the Jacobian resulting from changing variables from p f to ϵ.

Total Energy Loss
As done in [43], we convolve the radiative and elastic energy loss probabilities to yield a total probability of energy loss, Note that P rad (ϵ) contains Dirac delta functions at both ϵ = 0 and ϵ = 1 while P el (ϵ) is only a Gaussian.We will see below that the lack of a probability of nothing happening in the elastic energy loss probability is a major shortcoming in modelling the suppression in small collision systems.

Geometry
For large systems it is standard to use the Glauber model for the collision geometry, with Wood-Saxon distributions for the nucleon density inside the heavy ions [84].For p+A collisions the Glauber model cannot be applied in its most simple form, since subnucleonic features of the proton are expected to be important [85].Additionally, the subsequent evolution of the medium can be treated in a more sophisticated way than simply assuming a Bjorken expansion of the initial Glauber model geometry as was done, e.g., in [86].In this work we will use collision profiles generated with [85], and sourced from [87].In these calculations, initial conditions are given by the IP-Glasma model [88,89], which are then evolved with the MUSIC [15,90,91] viscous relativistic (2+1)D hydrodynamics code, followed by UrQMD microscopic hadronic transport [92,93].In this first comparison, we wish to make as few changes as possible from the original WHDG model and so we will use the initial temperature profile T (τ = τ 0 = 0.4 fm) for our collision geometry, where τ 0 is the turn-on time for hydrodynamics.This means that we are effectively using the IP-Glasma model [88,89] as the initial condition coupled with Bjorken expansion time dependence for all presented phenomenological results (unless otherwise stated.)This method for generating the initial conditions also allows for fluctuating initial conditions, which will be useful for future calculations of azimuthal momentum anisotropy of hard partons v n .
In addition, future work investigating the use of the complete and realistic time-dependent collision profiles from [85,87] will be of interest.
The QGP is treated as an ultrarelativistic mixture of a Fermi and Bose gas, following [43,94].This results in standard expressions, for various thermodynamics quantities [94] where ζ is the Riemann zeta function, T is the temperature, ρ q (ρ g ) is the density of quarks (gluons), n f is the number of active quark flavors (taken to be n f = 2 throughout), σ qg (σ gg ) is the gluon-gluon (quark-gluon) cross section, and we have used N c = 3.The cross section weighted density is denoted as ρ, which we will subsequently refer to as the density for simplicity.
The radiative (Eq.4) and elastic (Eqs.11 and 12) energy loss results were derived using a "brick" model, which represents a medium with a fixed length L and constant temperature T .In order to capture fluctuations in temperature and density, we need a mapping from the path that a parton takes through the plasma, to a brick with an effective length L eff and effective temperature T eff .
We follow WHDG [43] and define the effective pathlength as and the effective density as ) Here, the effective pathlength L eff includes all (x i , ϕ) dependence, and ρ eff is a constant for all paths that a parton takes through the plasma for a fixed centrality class.In principle, one can allow both the effective density and effective pathlength to depend on the specific path taken through the plasma.However, such a numerically intensive model is beyond the scope and objective of this work.Nonetheless, we will discuss some implications of the details of the geometry modelling in Sec. 6.
In WHDG [43] the prescription ρ ≡ ρ part was made 1 , where ρ part is the participant density-the density of nucleons which participate in at least one binary collision.This prescription is not necessary in our case, since we have access to the temperature profile [85,87].The temperature is extracted from the hydrodynamics output [85,87], and for L eff one evaluates the temperature at the initial time set by the hydrodynamics simulation, τ 0 = 0.4 fm.There is no unique mapping from realistic collision geometries to simple brick geometries, and more options are explored in [86].
Bjorken expansion [78] is then taken into account by approximating where in the last step we have evaluated T (τ, x) at the average time τ = L/2, following what was done in [43,95,96].In [43] this was found to be a good approximation to the full integration through the Bjorken expanding medium.
For a given collision system we can then calculate the distribution of effective pathlengths that a hard parton will travel in the plasma.We assume, as is standard and consistent with WHDG [43], that the hard partons have starting positions weighted by the density of binary nucleon-nucleon collisions, provided by IP-Glasma [87]. Figure 3 shows the distribution of effective pathlengths P L for central p + Pb (red), semi-central Pb + Pb (blue), and central Pb + Pb collisions (black) for √ s = 5.02 TeV.We also indicate the average effective pathlength in these three systems by single vertical lines.One can see that all three distributions are broad, the p + Pb system perhaps surprisingly so.One may further see that the average effective pathlength for central p + Pb collisions is not particularly small, ∼ 1 fm, which is not too different from the average effective pathlength for semi-central Pb + Pb collisions, ∼ 2 fm. Figure 4 compares the temperature of the plasma as a function of proper time in the rest frame of the plasma calculated via hydrodynamics (solid lines) versus the temperature from the Bjorken expansion formula (dashed lines).The effective temperature using hydrodynamics is calculated using Eq.19 and the Bjorken expansion approximation to the time dependence of the effective temperature is given by Eq. 20.
Calculations are performed for the same three collision systems as in Fig. 3. Due to the fluctuations of the initial conditions of the plasma in our model-because of both nucleonic and subnucleonic fluctuations [85]-we obtain a distribution of effective temperatures at each point in proper time.We show the mean and the 2σ width of the Bjorken temperature estimates in Fig. 4 .The width in the Bjorken result arises solely from the variation of the initial hydrodynamics temperature profile at τ = τ 0 .We computed the mean and standard deviation of the temperature distribution as a function of τ for the full hydrodynamics simulation; however the widths of these distributions are not plotted in Fig. 4 as for τ > τ 0 they are negligible compared to the width in the Bjorken expansion results.(Bjorken) [78] to the τ dependence of temperature, T (τ ) ≈ T eff (τ0)[τ0/τ ] 1/3 , is plotted along with the effective temperature T eff (τ ) calculated as a function of time via hydrodynamics (Hydrodynamics) according to Eq. 19.Uncertainty bands represent a 2σ (95% CI), and are shown only for the Bjorken estimate only as the uncertainty on the Hydrodynamics result is negligible.Both curves are plotted for the collision systems Pb + Pb at 0-5%, and 50-60% centrality; as well as p + Pb at 0-5% centrality.The freeze-out temperature T f.o.= 0.155 GeV is shown as a horizontal gray line, which is the temperature at which it is assumed that hadronic degrees of freedom take over in the plasma.Hydrodynamic temperature profiles are taken from [85].
One can see that for the A + A collisions, the Bjorken formula does a very good job of approximating the temporal evolution of the temperature from the full hydrodynamics simulations.For the central p+A system, however, the Bjorken formula significantly overestimates the temperature as a function of time.We comment below on the effect of this overestimation on the energy loss model.
The overestimation of the temperature in p + A collisions from the Bjorken expansion formula implies that the effective temperature used in the energy loss model is overestimated; as a result, the effects of energy loss are all overestimated in the small collision system.However, we expect this overestimation to be a small effect; using ⟨L eff ⟩ ∼ 1 fm, one can see in Fig. 4 that the difference between the full hydrodynamics temperature and the Bjorken estimate at ⟨τ ⟩ ∼ 0.5 fm is very small.Exploring the effect of dynamical effective pathlengths and time-dependent temperatures, which take medium expansion into account using the full hydrodynamic temperature profiles as a function of time, is left for future work.
The average pathlength in 0-5% most central p + Pb collisions is L ∼ 1 fm, which has an average temperature of T ≈ 0.35 GeV, and correspondingly a mean free path of λ g ≈ 0.75 fm.As one can see from Fig. 3, the large pathlength assumption thus breaks down for the vast majority of hard scatterings in central p + A collisions.Further, in these small collision systems L/λ g ∼ 1 for most of the distribution of effective pathlengths, which implies that approaches that assume many soft scatterings are likely inapplicable.
We note that although the hydrodynamics simulation turns off the QGP phase at the freeze-out temperature T f.o.= 0.155 GeV, the Bjorken expansion formula Eq. 20 has no such turn off.High-p T particles can of course interact with matter in the hadronic phase, which is, in part, captured by using the Bjorken expansion formula for determining the effective temperature as an input into our energy loss model.It is worth emphasizing that the effective pathlengths are determined at the initial time τ 0 = 0.4 fm, well before much of the plasma has had a chance to cool down; thus we do not lose any contribution to the effective pathlengths from hadronization.

Nuclear modification factor
The observable which we will be computing is the nuclear modification factor R AB (p T ) for a collision system A + B, defined experimentally by where dN AB/pp→h /dp T is the differential number of measured h hadrons in A + B/p + p collisions, and ⟨N coll ⟩ is the expected number of binary collisions (usually calculated according to the Glauber model [84].)To access this observable theoretically we make several assumptions about the underlying quark and gluon partons.In the following we will only refer to quarks, but all assumptions and formulae apply equally well for gluons.We first assume, following [43,94], that the spectrum of produced quarks q in the initial state of the plasma (before energy loss) is dN q prod /dp i = N coll × dN q pp /dp i where dN q pp /dp i is the quark production spectrum in p + p collisions.We further assume that the parton production spectra can be approximated by a power law, where n(p i ) is slowly varying, and A is a proportionality constant.The n(p i ) function is fitted using the initial parton spectra according to Eq. 22.For charm and bottom quarks, the initial parton spectra are computed using FONLL 2 at next-to-leading order [97]; and for gluons and light quarks, production spectra are computed 3 to leading order [98] as in [38,41].From Sec. 2.6, we have the total probability density function P tot (ϵ|p i ), which is the probability for a parton with initial transverse momentum p i to lose a fraction ϵ of it's energy such that the final momentum is p T = (1−ϵ)p i .Assuming that particle spectra are modified primarily due to energy loss implies, where dN q (p i ) (dN q (p T )) is the differential number of quarks in the initial (final) state with momentum p i (p T ).
Finally we assume that P tot (ϵ|p i ) varies slowly with p i , leading to the following expression for the partonic R q AB (neglecting hadronization for now) [94] In the above we have used P (ϵ|p T /[1 − ϵ]) ≈ P (ϵ|p T ) and n[p T /(1 − ϵ)] ≈ n(p T ) which follows from the slowly varying assumptions about n(p T ) and P (ϵ|p T ) and the soft assumption ϵ ≪ 1.The assumption that n(p T ) varies slowly, can easily be verified when fitting Eq. 22 to spectra.Since radiative energy loss is dominant, P (ϵ|p T ) will be peaked at ∆E/E determined with Eq. 5. Asymptotically this grows as log(E)/E for the (D)GLV result [45], and as log(E) for the correction [37].It is safe to assume this logarithmic growth is slow, however at high energies where the short pathlength correction dominates this assumption may become worse.
To incorporate our model for the collision geometry, we expand Eq. 26 to average over the effective pathlengths according to the length distribution described in Sec.2.7, where P L (L eff ) is the normalized distribution of effective pathlengths weighted weighted by the binary collision density (Fig. 3) and T (L eff /2) is the effective temperature determined according to Eq. 20.
The spectrum dN h /dp h for a hadron h is related to the spectrum dN q /dp q for a parton q via [94] where z = p q /p h ∈ (0, 1], p h is the observed hadron momentum, D h q (z, Q) is the fragmentation function for the process q → h, and Q is the hard scale of the problem taken to be Q = p h /z.The hadronic R h AB is then found in terms of the partonic R q AB (Eq.27) as [94] For details of the derivation needed for both Eqs.26 and 29, refer to Appendix B of [94].The fragmentation functions for D and B mesons are taken from [99], and the fragmentation functions for π mesons are taken from [100].Note that the fragmentation functions for π mesons were extrapolated outside their domain in Q 2 , as we found that the extrapolation was smooth.There are uncertainties in both the fragmentation function fits and the generation of initial parton spectra's; however these uncertainties are very small compared to others in this problem (running coupling, first order in opacity, etc.), and so we will not take the fragmentation function uncertainties into account.

Initial results
As a first exploration of the effect of the short pathlength correction to the DGLV radiative energy loss model, we consider Pb+Pb and p+Pb collision systems.We emphasize that in this work we are not trying to create the best possible fit to data (which could be done by tuning various parameters), but are instead focused on the impact of the short pathlength correction.For this reason all numerical values are used in consistency with the original WHDG predictions [38,43,79], and agreement with data is not the primary focus of this work.

Suppression of heavy flavor mesons
In Fig. 5, we show R AA (p T ) for D mesons in √ s = 5.02 TeV Pb + Pb collisions at 0-10% and 30-50% centrality from our convolved radiative and collisional energy loss model, with and without the short pathlength correction to DGLV radiative energy loss, compared to data from CMS [101], and ALICE [40].In Fig. 6, we show R AA (p T ) for B mesons in √ s = 5.02 TeV Pb + Pb collisions at 0-10% centrality from our convolved radiative and collisional energy loss model, with and without the short pathlength correction to DGLV radiative energy loss, compared to minimum bias data from CMS [102].From the figures we conclude that for heavy flavor final states, the short pathlength correction to the R AA is small up to p T ∼ 100 GeV, and that the difference grows with p T .This small difference was expected from the numerical energy loss calculations performed in [37] and reproduced in Fig. 1 that showed the small pathlength corrections are a relatively small correction for quarks.Agreement with data for D 0 mesons (Fig. 5) is especially good for all p T ≳ 5 GeV, given that the calculation is leading order.The D 0 meson suppression prediction is underestimated for moderate p T ≲ 20 GeV in comparison to a previous prediction with the WHDG model in [79].We believe the main cause of this difference with past results is from not using the approximation q max → ∞ (see Sec.   [103].In the same figure, we also show the prediction of R pA of D mesons from our energy loss model with the collisional engery loss turned off.One can see that the energy loss model that includes both collisional and radiative energy loss, both the small pathlength corrected and the uncorrected version, dramatically overpredicts the suppression in this small system.At the same time, the predictions from the model with the collisional energy loss turned off are significantly less oversuppressed compared to the data.The surprising sensitivity to the presence of the collisional energy loss process is due to our naive implementation of the collisional energy loss as discussed in Sec.2.6: the WHDG [43] collisional energy loss model assumes an average collisional energy loss with Gaussian fluctuations, which is inappropriate for a small system with very few elastic scatterings.

Suppression of light flavor mesons
Figure 8 shows the R AA (p T ) for π mesons as a function of p T in 0-5% most central Pb + Pb collisions at √ s = 5.02 TeV from both our convolved radiative and collisional energy loss model compared to the R AA (p T ) for charged hadrons measured by ATLAS [104], CMS [105], and AL-ICE [106].In our energy loss model, we used the fraction of π mesons originating from light quarks versus gluons as in [38].Figure 8 shows that the short pathlength correction does, in fact, lead to a stronger p T dependence in R AA (p T ) than the prediction without the correction.In fact, the correction leads to a large change in predicted R AA (p T ), even at moderate momenta 20 GeV ≲ p T ≲ 100 GeV, when compared to the size of the correction for heavy flavor mesons.This large change in R AA (p T ) is consistent with the large change in the average energy loss as calculated numerically in [37] and reproduced in Figs. 1 and 2; the correction is almost a factor of 10 times larger for gluons compared to quarks due to the specific way in which color triviality is broken by the short pathlength correction to the energy loss.For p T ≲ 200 GeV the corrected result is tantalizingly consistent with data, however for p T ≳ 200 GeV the corrected result predicts anomalously large enhancement up to R AA = 1.5 at p T ≃ 450 GeV, a shocking ∼200% increase over the uncorrected result.Figure 9 shows predictions for the R pA (p T ) of π mesons in 0-10% most central p+Pb collisions from our convolved radiative and collisional energy loss model, with and without the short pathlength correction to the radiative energy loss, as well as predictions from our model with the collisional energy loss turned off.Fig. 9 also shows AT-LAS charged hadron R pA (p T ) data [39].In this figure the breakdown of the elastic energy loss is even more obvious than for the heavy flavor mesons, since the elastic energy loss for gluons is enhanced by a factor C A /C F ≈ 2. The uncorrected radiative-only result (no collisional energy loss) predicts mild suppression of R AA ≈ 0.9 for all p T ; and the corrected radiative-only result predicts enhancement that grows monotonically in p T , reaching R AA ≈ 1.5 at p T = 150 GeV.The predicted enhancement for p T ≳ 100 GeV is in excess of the measured enhancement [39], however for p T ≲ 100 GeV the presence of enhancement is qualitatively consistent with data.
The short pathlength "correction" leads to changes in R AA (p T ) of 100% or greater in the light hadron sector, with predictions of R AA (p T ) well in excess of 1 for extremely high momenta, in both small and large collision systems leads us to question whether or not the energy loss calculation is breaking down in some fundamental way.In particular, are we applying the energy loss formulae in our energy loss model in some regimes where the assumptions made in the derivation of those energy loss formulae no longer apply?

Consistency of assumptions in DGLV
The prediction of significant enhancement of high-p T light flavor mesons shown in Figs. 8 and 9 stems from the asymptotic dependence of the short pathlength correction on energy [37].We see from Eqs. 6b and 7 that for asymptotically large values of energy, ∆E corr./E ∼ E 0 while ∆E DGLV /E ∼ ln E/E.Thus, inevitably, the correction becomes larger than the uncorrected result in the large E limit.Presumably, then, there's some intermediate value of the energy at which the assumptions that went into either the DGLV derivation, the derivation of the correction, or both are violated.As noted in the Introduction, the derivations of both DGLV energy loss and its small pathlength correction assumed: (1) the Eikonal approximation which assumes that the energy of the hard parton is the largest scale in the problem; (2) the soft radiation approximation which assumes x ≪ 1; (3) collinearity which assumes k + ≫ k − ; (4) the impact parameter varies over a large transverse area; and (5) the large formation time assumption which assumes For the large formation time assumption we found that in the original calculation (details in [46]), ω 0 ≪ µ was only used in the weaker form ω 0 ≪ µ 1 ; and so we will be considering this weaker assumption instead.
In this section we numerically check the consistency of these assumptions with the final radiative energy loss result.In particular, the analytic properties of the matrix element mean that it may have non-zero support for momenta that are unphysical (even complex).The relevant question for us is: does the matrix element (modulus squared) give a significant contribution to the energy loss in kinematic regions that are integrated over but for which the derivation of the matrix element is not under control?
In an attempt to (partially) answer this question, we are motivated to calculate expectation values of ratios assumed small under the various assumptions, weighted by the absolute value of the mean radiative energy loss distribution, Eq. 5 4 .Explicitly the procedure to calculate the expectation value of a function R({X i }), depending on the set of variables {X i }, is ⟨R⟩ ≡ where {X i } can be any of {k, q, x, ∆z} and d{X i } ≡ i dX i .Also note that R can depend on quantities that are not integrated over, such as {L, E, µ}.It is important to note that this expectation value is not an expectation value in the usual sense, where the distribution is the distribution of radiated gluons, because we are weighting by the radiative energy loss and not radiated gluon number.It is also important to note that even if a particular assumption is violated in the sense of this weighted average, that violation does not necessarily mean that the correction computed by relaxing the assumption is large; rather, we only know that the correction is not necessarily small.
1 , where the DGLV and correction derivations assume that ω 1 /µ 1 ≪ 1, using L = 5 fm, λ g = 1 fm, and µ = 0.5 GeV.The large formation time assumption is explicitly violated for the energy loss with and without the short pathlength correction.For the energy loss without the correction, the large formation time assumption is violated for E ≳ 100 GeV for both charm quarks and gluons, while for the energy loss with the correction the large formation time assumption is violated for E ≳ 35 GeV for gluons and for E ≳ 50 GeV for charm quarks.This breakdown calls into question the validity of the large formation time assumption in DGLV radiative energy loss for p T ≳ 100 GeV, regardless of whether the energy loss receives a short pathlength correction.The increased rate and magnitude of large formation time assumption violation once the short pathlength correction is included in the energy loss indicates that the short pathlength corrected energy loss model predictions may be breaking down at moderate to high p T in central A + A collisions.One sees that the large formation time assumption, always breaks down before enhancement, ∆E/E < 0, is predicted.
One possibility for the increased rate and degree of large formation time breakdown once the short pathlength is included -and the erroneously large correction at high p T for π mesons-is the emphasis placed on short pathlengths by the exponential distribution of scattering centers.The exponential distribution was originally chosen to make the analytics simple, with the physical motivation that it captures the Bjorken expansion in the medium [42].
Bjorken expansion leads to a power law decay of the plasma density in time [78], and so an exponential distribution likely overestimates the amount of expansion.This biases scatterings to occur at smaller ∆z than is physical, and likely overestimates the contribution from the short pathlength correction.Additionally it is not obvious how to model the time dependence of the collision geometry before thermalization τ ≲ τ 0 , as in principle the medium should be thermalizing during this time.Furthermore the treatment of scatters that occur for times τ ≲ τ 0 is not obvious as it is possible that the well-separated scattering centers assumption λ ≫ µ −1 breaks down in this phase of the plasma.It was found numerically that DGLV energy loss results are insensitive to the exact distribution of scattering centers used [34].Perhaps not surprisingly, the small pathlength correction has a large sensitivity to the exact distribution of scattering centers used [37].
We are thus motivated to consider an alternative distribution of scattering centers as we consider whether or not the various assumptions made in the energy loss derivations are consistent with our final energy loss numerics.In this paper we will consider, in addition to the usual exponential distribution, the truncated step distribution from [37].The truncated step distribution is given by ρstep (∆z) ≡ (L − a) −1 Θ(∆z − a)Θ(L − ∆z), where a is a small distance cut off.The truncated step function attempts to capture the effect of a "turn on" of the QGP, before which no energy loss takes place, with subsequent equal probability for a scattering to occur until the end of the pathlength.We think of the exponential distribution and truncated step distributions as limiting cases for what the real distribution of scattering centers may be.The exponential distribution maximally emphasizes the effect of early-time physics, while the truncated step distribution completely neglects early-time physics.A more realistic distribution is likely somewhere in between these two extremes.
One choice for a is a = µ −1 , since for ∆z ≲ µ −1 , the production point and scattering center are too close to be individually resolved.Another choice is a = τ 0 , the hydrodynamics turn-on time; τ 0 is a particularly reasonable choice since we already only consider the medium density evolution for times τ > τ 0 in computing our pathlengths, Eq. 18.For a typical value of µ = 0.5 GeV in central Pb + Pb collisions, τ 0 ≃ µ −1 ≃ 0.4 fm; however in central p + Pb collisions µ −1 ≃ 0.25 fm < τ 0 = 0.4 fm, and so the distinction between the two options for a might be important.In this report we have chosen to use a = τ 0 throughout for simplicity.
We now check all of the assumptions made in the computation of the DGLV radiative energy loss for: the corrected and uncorrected results; the exponential and truncated step scattering center distributions; charm quarks and gluons; and for large (L = 5 fm) and small (L = 1 fm) systems (Figs.11-16).All calculations use constant µ = 0.5 GeV, and λ g = 1 fm which are approximate averages in A + A collisions, and standard benchmark choices [34].Figs.11 and 12 show the expectation values of ω 1 /µ 1 = (k − q 1 ) 2 /2xE µ 2 + q 2  1 and ω 0 /µ 1 = k 2 /2xE µ 2 + q 2 1 , respectively.The large formation time assumption is equivalent to both ratios much less than one.For an ex-ponential distribution of scattering centers with L = 5 fm (top panes of Figs.11 and 12), the large formation time assumption: breaks down for E ≳ 40 GeV for the corrected result, and breaks down for E ≳ 100 GeV for the uncorrected result.For the truncated step distribution with L = 5 fm, the large formation time assumption breaks down for both the corrected and uncorrected results for E ≳ 100 GeV.We see in the plots the known numerical insensitivity to the distribution of scattering centers in the uncorrected DGLV radiative energy loss result.For L = 1 fm (bottom panes of Figs.11 and 12), the shape of all curves are approximately the same, but scaled by a factor ∼ 1 fm/5 fm in E. Thus the breakdown in the large formation time assumptions occurs roughly five times earlier in E for the L = 1 fm pathlengths compared to the L = 5 fm pathlengths.The reason for this simple approximate scaling is that all of the nontrivial dependence of ∆E/E on E and L in the distribution of emitted gluons Eq. 4 comes from terms ∼ ω α ∆z where α ∈ {0, 1, m}.Once integrated these terms become ω α L ∼ L/E; assuming k and q have negligible dependence on E. If the kinematic cutoffs on the k, and q integrals are important then this scaling breaks down; any deviation from this simple scaling must be due to the effects of the cutoff.
We found that finite kinematic bounds do not significantly affect the consistency of the assumptions.The above scaling argument holds for all expectation values ⟨R⟩ so long as R does not depend on ∆z.Thus, in order to keep the number of plots shown to a manageable number, most assumption consistency plots from now on will be shown for only L = 5 fm.
The collinear and soft assumptions are tested for consistency in Figs. 13 and 14, respectively.We find that both of these assumptions are consistently satisfied for both the short pathlength corrected and uncorrected DGLV results, for both the exponential and the truncated step distributions of scattering centers.
Figure 14 shows ⟨x⟩ as a function of parent parton energy E, where ⟨x⟩ ≪ 1 is assumed under the soft approximation.The expectation value ⟨x⟩ decreases monotonically in energy for the uncorrected result with an exponential distribution of scattering centers; all other expectation values appear to converge numerically to a constant nonzero value.For the DGLV result with an exponential distribution, one can calculate ⟨x⟩ analytically for asymptotically high energies.For asymptotically high energies we take m g → 0, M → 0, k max → ∞, and q max → ∞.These simplifications allow the angular, k, and q integrals to be done analytically, as described in [42,45].Proceeding in this way we obtain the following asymptotic expression for ⟨x⟩: =⇒ ⟨x⟩ → 0 as E → ∞.
In a similar way we can derive the same result for the short pathlength correction with an exponential distribution, using the asymptotic result from [37] (see Eq. 6a) Note that numerical investigation of ⟨x⟩ indicates that the convergence to the asymptotic values is slow for both the corrected and uncorrected energy loss.For the truncated step distribution it is more difficult to perform asymptotic calculations, but numerical investigation shows that the uncorrected result with truncated step distribution converges to ⟨x⟩ ≈ 1 3 , and the corrected result with truncated step distribution converges to ⟨x⟩ ≈ 1 2 .Figure 15 tests the consistency of the large pathlength assumption with the DGLV result as a function of parent parton energy E for charm quarks.For L = 1 fm, we see that the large pathlength assumption is not a good approximation for either the uncorrected or small pathlength corrected DGLV calculation for the exponential scattering center distribution; with ⟨1/∆z µ⟩ ∼ 0.6, the large pathlength assumption is not a particularly good approximation for the truncated step distribution, either.Even for L = 5 fm, the large pathlength assumption breaks down for the short pathlength corrected energy loss when the exponential distribution of scattering centers is used.While not shown, results for gluons are essentially identical.We see that, as expected, one must quantitatively determine the importance of the short pathlength correction terms to the radiative energy loss derivation for short pathlengths, L ∼ 1/µ for all values of E and for all parton types.We note that there was one final assumption implicitly made in the derivation of the short pathlength correction to the DGLV radiated gluon distribution, the short formation time (with respect to scattering centers) assumption, ∆z ω 0 ≫ 1 [37,46].This in conjunction with the large formation time assumption furnishes a separation of scales ∆z −1 ≪ ω 0 ≪ µ 1 .Alternatively one can view this as a large pathlength assumption, wherein we have replaced ∆z µ ≫ 1 with ∆z ω 0 ≫ 1; which is guaranteed to be a weaker assumption, according to the large formation time assumption.Fig. 16 shows ⟨∆z ω 0 ⟩ as a function of parent parton energy E for charm quarks with the short pathlength corrected energy loss.We find that for both large and small systems, and exponential and truncated step distributions this assumption holds self-consistently.While not shown, results for gluons are essentially identical.The short formation time (with respect to scattering centers) assumption is intimately tied to the large formation time assumption, and so if future work aims to remove the large formation time assumption, then this assumption should also be removed.Note that we computed ⟨∆z ω 0 ⟩ (assumed ≫ 1) instead of ⟨(∆z ω 0 ) −1 ⟩ (assumed ≪ 1), due to numerical convergence issues with the latter.

Suppression with Exponential and Step Distributions
We saw in the previous section that for p T ≳ 40 GeV the large formation time approximation breaks down for energy loss calculations using the exponential distribution of scattering centers (which bias the scattering to shorter pathlengths) while the large formation time approximation breaks down only for p T ≳ 100 GeV for an energy loss calculation using the truncated step distribution of scattering centers.Additionally, we claimed that the exponential distribution of scattering centers and the truncated step distribution of scattering centers represent two extreme possibilities for what a more realistic distribution of scattering centers likely will be.We are thus motivated to explore the sensitivity of our suppression predictions to the choice of distribution of scattering centers.For pure DGLV energy loss without the short pathlength correction, one saw an insensitivity to this choice of scattering center distributions [34].We will see that when including the short pathlength correction to the radiative energy loss, the heavy flavor observables are still relatively insensitive to the distribution of scattering centers.However, we find that hadron observables that include a contribution from gluons are very sensitive to the choice of scattering center distribution, with the use of the truncated step function dramatically reducing the effect of the short pathlength correction to pion R AA (p T ).This dramatic reduction in the effect of the short pathlength correction to the phenomenologically accessible pion suppression observable is expected from the dramatic reduction in the effect of the short pathlength correction to the average radiative energy loss, as seen in [37].
Figure 17 shows R AA (p T ) for D mesons in central 0-10% and semi-central 30-50% Pb + Pb collisions at √ s = 5.02 TeV for both the exponential and truncated step distributions of scattering centers.Figure 18 shows R AA (p T ) for B mesons in central 0-10% Pb+Pb collisions at √ s = 5.02 TeV for both the exponential and truncated step distributions of scattering centers.For both heavy meson cases the truncated step distribution decreases the R AA by up to 10% (20%) for the original (corrected) WHDG results.For both D and B mesons with a truncated step distribution, the short pathlength correction is negligible for p T ≲ 100.As seen before, with an exponential distribution of scattering centers the effect of the short pathlength correction is ≲ 10%, which is small compared to other theoretical uncertainties in the model (e.g. higher orders in α s , treatment of the early times, etc.).Agreement with data is good for all predictions (corrected/uncorrected, and exp./trunc.step distributions), except for p T ≲ 10 where bulk effects may be important and the eikonal approximation is likely breaking down.
The radiative-only nuclear modification factor R pA for D 0 mesons in 0-10% most central p+A collisions is shown in Fig. 19, with and without the short pathlength correction, and with both an exponential and a truncated step distribution of scattering centers.We only show the radiative-only R pA since in Sec. 3 we determined that the elastic contribution was erroneously large in p + A collisions due to the inapplicability of the average elastic energy loss in previous WHDG calculations in small collision systems.The R pA for the exponential distribution without the correction and the truncated step distribution with and without the correction all agree with each other to within 5% and predict mild suppression of R pA ≈ 0.9 for all p T .The corrected R pA with an exponential distribution predicts mild suppression at low p T ≲ 20 GeV and consistency with unity at moderate p T ≳ 20 GeV.Measured R pA is shown for D 0 mesons produced in 0-5% most central p + Pb collisions from ALICE [103].Data predicts mild enhancement for all p T , not inconsistent with our results shown here.
Figure 20 shows R AA (p T ) for pions in 0-5% central Pb + Pb collisions at √ s = 5.02 TeV.The convolved collisional and radiative energy loss model predictions both include and exclude the short pathlength correction to radiative energy loss, and use either the exponential or the truncated step distribution for the scattering centers.The predictions are compared to data from ATLAS [104], CMS [105], and ALICE [106].While the predictions excluding the short pathlength correction are insensitive to the particular scattering center distribution chosen, one can see the tremendous sensitivity of the predictions to the scattering center distribution when the short pathlength correction to the radiative energy loss is included; when the truncated step distribution is used, the effect on R AA (p T ) of the short pathlength correction to the DGLV radiative energy loss is dramatically reduced.For our theoretical predictions of pion suppression we only include radiative energy loss, as the average elastic energy loss of the WHDG model is inappropriate to use here.We show predictions with (solid) and without (dashed) the short pathlength correction to the radiative energy loss, and we show predictions when using either the exponential (red) or truncated step (blue) distribution of scattering centers.Charged hadron suppression data is from ATLAS [39].
The difference between the R π pA (p T ) predictions from the two scattering center distributions is small when excluding the short pathlength correction.We see again that the effect on nuclear modification factor from the short pathlength correction to the radiative energy loss is large when using the exponential distribution of scattering centers and relatively small when using the truncated step distribution of scattering centers.The prediction of enhancement by the corrected R pA with an exponential distribution is qualitatively similar to the observed enhancement for moderate momenta p T ≲ 60 GeV.If we are provocative, we may thus suggest that the experimentally measured excess in R pA (p T ) is actually due to final state effects rather than initial state or normalization effects.

Discussion
The primary goal of this work, as outlined in Secs.2-5, is to implement the WHDG convolved energy loss model [38,43] with the novel inclusion of the short pathlength correction to the radiative energy loss [37].
Due to the complexity of this model, there are many points at which one must decide on the level of approximation to proceed with.In line with the motivation of this work, we have always chosen to maintain consistency with previous work such as WHDG and DGLV.This has occasionally led us to overlook a more physically reasonable prescription (in our view) for a component of the energy loss model.One such instance concerns the treatment of the realistic collision geometry (see Sec. 2.7) and how it is mapped to the brick geometry.We will now present a derivation of a more realistic effective pathlength and effective temperature, which could be implemented in future work.
An integral part of radiative energy loss is the distribution of scattering centers, ρ(∆z), normalized to a single hard scatter at first order in opacity (see Eq. 3).The quantity ρ(∆z) is a model for the shape of the plasma, as the parton propagates through it.In theory, it is possible to improve the realism of the model by integrating through the plasma, in some sense replacing ρ(∆z) with a realistic plasma density ρ(⃗ x, τ ) and correspondingly T with a realistic plasma temperature T (⃗ x, τ ).In this case the pathlength L no longer needs to be specified a priori.This approach would allow for a more accurate simulation of the hard partons propagation through the plasma, and is implemented similarly in the CUJET model [107,108].Unfortunately implementing this approach is exceedingly computationally expensive, even for the original DGLV result, and would present significant computational challenges for the short pathlength corrected results.
To account for a realistic collision geometry, we instead need effective temperatures, densities, and lengths as inputs to the simple brick models for elastic and radiative energy loss.Essentially, we establish a brick with characteristic {L eff , T eff , ρ} for each parton that propagates through the plasma.The relationship between the plasma density ρ(x, τ ) and the distribution of scattering centers ρ(∆z), is a separation of the shape of the plasma density from its magnitude.Schematically, we can write this separation as where d∆z ρ(∆z) = 1.In the above x i is the hard parton production point, φ is the direction of propagation, and ρ is the density of the QGP.Equation 33 is equivalent to what is done in CUJET [107,108] while Eq.34 is the approximation made in GLV [45], DGLV [42], WHDG [38,43] and the short pathlength correction to DGLV [37].Note that the step from Eq. 33 to Eq. 34 is exact for a brick of plasma.The power of the approximation made in Eq. 34 lies in separating the dependence of the path taken by the parton through the plasma from the rest of the energy loss calculation.This approach allows us to prescribe ρ(∆z), for instance exponential decay or truncated step, and perform the ∆z integral analytically.A more realistic approach to this, but less numerically intensive than integrating through the realistic plasma, would be to fit a trial ρ(∆z|x i , ϕ) to the realistic plasma density ρ(x i , ϕ) for each path taken by a parton through the plasma.
The magnitude of the density is definitionally related to the opacity n via [45] n where ρ is the density from Eq. 17d, and we have used Eq.36 to define the effective length L eff and effective gluon mean free path λ eff .Equations 17e and 36 yield where we have the freedom to prescribe ρ eff .Note that Eq. 38 differs from Eq. 18 in the fact that the density is evaluated at τ = z, which follows from the definition of the opacity in Eq. 35.
Breaking apart the opacity n as n = L eff /λ eff serves to: make contact with the effective pathlength prescription in WHDG [43], with a more rigorous derivation; -obtain a prescription for the effective density ρ eff which can then be used to calculate other thermodynamic quantities in Eq. 17, most importantly the Debye mass µ; -and obtain a length scale L eff in the problem, which is important for prescribing the distribution of scattering centers ρ(∆z).
In this manuscript, we have followed WHDG in using Eq.19 as the effective density, which prescribes a single density to the entirety of the plasma.A more natural approach can be motivated by considering the step form Eq. 33 to Eq. 34.In this step, we have approximated µ(z) ≈ µ eff , which leads to a natural definition for ρ eff as the average temperature along the path through the plasma, weighted by the plasma density.
This means both L eff and ρ eff depend on the specific path that the parton takes.This dramatically increases the numerical complexity as we must now evaluate the energy loss distribution for a distribution in (L eff , ρ eff ).Note that Bjorken expansion is naturally taken into account with this prescription for the effective density ρ eff and so we do not need to use the approximation in Eq. 20.This approach could be implemented in future work, offering the advantage of a more realistic collision geometry compared to the implementation in this work and WHDG [43]; while still being significantly less computationally expensive than integrating through the realistic plasma as in CUJET [108].

Conclusions
In this article we presented the first predictions for the suppression of leading high-p T hadrons from an energy loss model with explicit short pathlength corrections to the radiative energy loss.We included collisional energy loss in the model, as well as averages over realistic production spectra for light and heavy flavor partons that propagate through a realistic QGP medium geometry generated by second order viscous hydrodynamics.Thus our calculations here are, to a very good approximation, those of the WHDG energy loss model [86], but with short pathlength corrections [37] to the DGLV opacity expansion [42,45].Predictions were presented for central and semicentral Pb+Pb collisions and central p+Pb collisions and compared to data from the LHC.
We saw that the inclusion of the short pathlength correction to the radiative energy loss led to a reduction of the suppression of leading hadrons.This reduction is well understood as a result of the short pathlength correction enhancing the effect of the destructive LPM interference between the zeroth order in opacity DGLAP-like production radiation and the radiation induced by the subsequent collisions of the leading parton with the medium quanta.The reduction in the suppression also increases as a function of p T , which is a result of the different asymptotic energy scalings of DGLV energy loss (∆E ∼ log E) compared to the short pathlength correction (∆E ∼ E).
For heavy flavor observables, the inclusion of the short pathlength corrections leads to only a modest ∼ 10% enhancement of R AA (p T ) in Pb + Pb and p + Pb collisions.Even though the relative R pA enhancement of ∼ 10% is similar to that of R AA , as one can see in Fig. 2 the influence of the short pathlength correction on the energy loss is significantly larger for shorter pathlengths and therefore also in the smaller collision system.The reason that the R pA and R AA have similar relative enhancements is due to the scaling R AA ∼ (1 − ϵ) n−1 , where ϵ ≡ ∆E/E is the fractional energy lost by the leading parton.One can determine that the effective short pathlength correction, averaged over the Poisson convolution and geometry, is about 100% stronger in p + Pb compared to central Pb+Pb.(Note that the Poisson convolution, with its large probability of no interaction or energy loss for short pathlengths, is crucial for in fact reducing the enormous short pathlength correction influence seen in Fig. 2 in the energy loss model.)Thus, since R pA ∼ 1 and n ∼ 6, even though the short pathlength correction is about 100% stronger in p + Pb, the influence on R pA is similar to that in R AA .
When we assume that the distribution of scattering centers that stimulate the emission of gluon radiation from high-p T partons is given by an exponential, which biases the leading partons to scatter at shorter distances, the short pathlength correction to R π AA (p T ) grows dramatically with p T .This very fast growth in p T is due to the very large short pathlength correction to the gluonic radiative energy loss: the short pathlength correction to radiative energy loss breaks color triviality, and the correction to gluonic radiative energy loss is about ten times that of quark radiative energy loss (instead of the approximately factor of two that one would expect from color triviality) [37].When a truncated step function is used as the distribution of scattering centers that stimulate the emission of gluon radiation, the short pathlength correction to R π AA (p T ) becomes a much more modest ∼ 10%.It is interesting, but perhaps not totally surprising, that the short pathlength correction to the energy loss introduces an enhanced sensitivity to the precise distribution of scattering centers used in the energy loss model.In either case of distributions of scattering centers, the more rapid growth in the short pathlength correction as a function of p T than that of the uncorrected DGLV energy loss [38] suggests that at least part of the faster-than-expected growth of the measured R π AA (p T ) as a function of p T may be due to the influence of the short distance corrections to the energy loss of hard partons; cf. the reduction of suppression due to running coupling [109].
We found that the average collisional energy loss, with fluctuations of this energy loss given by a Gaussian whose width is dictated by the fluctuation-dissipation theorem, of the WHDG energy loss model [86] is inappropriate for small colliding systems.Considering radiative energy loss only, R pA (p T ) for heavy flavor hadrons was again only modestly, ≲ 10%, affected by the short pathlength correction to energy loss.R π pA (p T ) was significantly affected, ∼ 50%, by the short pathlength correction for an exponential distribution of scattering centers, but more modestly so, ∼ 10%, for a truncated step distribution of scattering centers.For both distributions of scattering centers, the short pathlength corrected pion nuclear modification factor sees a tantalizing enhancement above 1, similar to data [39,40].One may then provocatively suggest that the experimentally measured enhancement of R pA (p T ) > 1 may be due-at least in part-to final state effects.
We also investigated the self-consistency of the approximations used in the derivations of DGLV and short pathlength correction to DGLV single inclusive radiative gluon emission kernels for the phenomenologically relevant physical situations of RHIC and LHC.We constructed dimensionless quantities that represented the approximations and checked whether, when averaged with a weight given by the strength of the energy loss kernel determined by DGLV or the short pathlength corrected DGLV, those quantities were small (or large) as required by the approximations that went into deriving those same DGLV and short pathlength corrected DGLV single inclusive radiated gluon spectrum kernels.
We found that, when weighted by the energy loss kernel, the soft and collinear approximations were selfconsistently satisfied when computed with phenomenologically relevant parameters.We further found that both the original DGLV derivation and the DGLV derivation with the inclusion of the short pathlength correction are not consistent with the large formation time approximation for modest O(10-100 GeV) energies and pathlengths O(1-5 fm), independent of the choice of distribution of scattering centers.Finally, we see that the large pathlength approximation breaks down for small pathlengths ∼ 1 fm and even for large pathlengths ∼ 5 fm for large enough ≳ 100 GeV energies.We noted in Sec.2.1 one more assumption, that of large transverse area.This assumption is very difficult to assess using the methods of this article, especially as the utilization of the assumption occurs very early in the derivation.Qualitatively, even in p+A collisions, one has that the transverse size of the system will be ∼ (1 fm) 2 whereas the typical scale set by the scattering process itself is 1/µ 2 ∼ (0.5 fm) 2 .It thus seems likely that this large transverse size assumption holds even for small collision systems.
Instead of thinking of the self-consistency of the numerics with the assumptions that went into the derivation of the energy loss, one may rather formulate the issue as whether or not one is integrating the matrix element (modulus squared) beyond the region under which its derivation is under control.Thus one way of understanding that, e.g., ⟨ω 0 /µ 1 ⟩ > 1 as shown in Fig. 12 is that the matrix element is integrated over regions of kinematics under which the derivation is not under control.One may consider restricting the kinematics that are integrated over to those for which the derivation is under control.We find, for example, that the expectation of the collinear approximation ⟨k − /k + ⟩ is self-consistently less than 1, but note that the gluon kinematics are restricted in such a way as to enforce collinearity |k| max = 2x(1−x)E ⇔ k − < k + .As was shown in [34,110], the DGLV inclusive gluon emission kernel is not under good control near the kinematic bound k − ∼ k + .One may thus consider restricting the gluon kinematics such that the large formation time assumption is respected, for example, by taking |k| max = Min( √ 2xEµ 1 , 2x(1−x)E) ⇔ ω 0 < µ 1 and k − < k + .Then presumably one would find that ⟨ω 0 /µ 1 ⟩ ≲ 1.However, as was shown in [110], there would then be a significant sensitivity in the energy loss model predictions to the exact kinematic bound chosen.We thus conclude that in order to confidently compare leading hadron suppression predictions in A + A collisions at ≳ 100 GeV or in p + A collisions at ≳ 10 GeV from an energy loss model based on an opacity expansion of the single inclusive gluon emission kernel, future work is needed to re-derive the opacity expansion single inclusive gluon emission kernel with both the large pathlength and the large formation time approximations relaxed.Further work of numerically implementing finite pathlength effects in elastic energy loss [86,111] will also play an important role in any quantitative comparison of an energy loss model and leading hadron suppression in small systems such as p + A collisions.

Fig. 1 .
Fig. 1.The induced radiative fractional energy loss ∆E/E is plotted as a function of energy E for charm quarks (c), and gluons (g); and at E = 10 GeV (top panel), and E = 100 GeV (bottom panel).Note that ∆E/E < 0 is energy gain relative to the vacuum.Calculations were done with constant µ = 0.5 GeV and λg = 1 fm.

Fig. 2 .
Fig. 2. The ratio of the magnitude of the correction to the DGLV radiative energy loss −∆Ecorrected and the uncorrected DGLV radiative energy loss ∆E DGLV is plotted as a function of incident energy E. This ratio is plotted for charm quarks (c), and gluons (g); and at L = 1 fm (dashed) and L = 5 fm (solid).Calculations were done with constant µ = 0.5 GeV and λg = 1 fm.

Fig. 3 .
Fig. 3. Distribution of the effective pathlengths PL(L) in p+A and A + A collision systems, weighted by the binary collision density.The collision systems Pb + Pb at 0-5% and 50-60% centrality; as well as p+Pb at 0-5% centrality are shown, all at √ s = 5.02 TeV.Vertical lines indicates the average lengths for the respective collision systems, which numerically are 0.99 fm, 2.06 fm, and 4.68 fm left to right.

Fig. 4 .
Fig. 4. Plot of the temperature T as a function of the proper time τ (in the plasma rest frame).The Bjorken expansion approximation(Bjorken)  [78] to the τ dependence of temperature, T (τ ) ≈ T eff (τ0)[τ0/τ ] 1/3 , is plotted along with the effective temperature T eff (τ ) calculated as a function of time via hydrodynamics (Hydrodynamics) according to Eq. 19.Uncertainty bands represent a 2σ (95% CI), and are shown only for the Bjorken estimate only as the uncertainty on the Hydrodynamics result is negligible.Both curves are plotted for the collision systems Pb + Pb at 0-5%, and 50-60% centrality; as well as p + Pb at 0-5% centrality.The freeze-out temperature T f.o.= 0.155 GeV is shown as a horizontal gray line, which is the temperature at which it is assumed that hadronic degrees of freedom take over in the plasma.Hydrodynamic temperature profiles are taken from[85].

Fig. 5 .
Fig. 5.The nuclear modification factor RAA as a function of final transverse momentum pT is calculated for D 0 mesons, with and without the short pathlength correction to the radiative energy loss.Calculations were done for 0-10% centrality as well as 30-50% centrality.Data from CMS [101], and AL-ICE [40] are shown for comparison; where error bars (boxes) indicate statistical (systematic) uncertainties.The global normalisation uncertainty on the number of binary collisions is indicated by the solid boxes in the top left corner of the plot (left to right: 0-10% CMS, 0-10% ALICE, 30-50% ALICE).

Fig. 6 .
Fig. 6.The nuclear modification factor RAA as a function of final transverse momentum pT is calculated for B mesons with and without the short pathlength correction.Data for B ± mesons from CMS [102] is shown for comparison, where error bars (boxes) indicate statistical (systematic) uncertainties.The global normalisation uncertainty on the number of binary collisions is indicated by the solid boxes in the top left corner of the plot.

Fig. 7 .
Fig. 7.The nuclear modification factor RpA for D 0 mesons as a function of final transverse momentum pT is calculated with and without the short pathlength correction.The RpA is calculated both with collisional and radiative energy loss (el.+ rad.), and with radiative energy loss only (rad.only).Data from ALICE [103] for D 0 mesons is shown for comparison, where statistical (systematic) uncertainties are represented by error bars (boxes).The global normalisation uncertainty on the number of binary collisions is indicated by the solid box in the center left of the plot.

Figure 7
Figure7shows predictions for the R pA of D mesons as a function of the final transverse momentum p T for central p + Pb collisions at √ s = 5.02 TeV from our convolved radiative and collisional energy loss model compared to ALICE data[103].In the same figure, we also show the prediction of R pA of D mesons from our energy loss model with the collisional engery loss turned off.One can see that the energy loss model that includes both collisional and radiative energy loss, both the small pathlength corrected and the uncorrected version, dramatically overpredicts the suppression in this small system.At the same time, the predictions from the model with the collisional energy loss turned off are significantly less oversuppressed compared to the data.The surprising sensitivity to the presence of the collisional energy loss process is due to our naive implementation of the collisional energy loss as discussed in Sec.2.6: the WHDG[43] collisional energy

Fig. 8 .
Fig.8.The nuclear modification factor RAA for π mesons produced in Pb + Pb collisions at 0-5% centrality is calculated using WHDG, with and without the short pathlength correction to the radiative energy loss.Charged hadron suppression data from ATLAS[104], CMS[105], and ALICE[106] is plotted for comparison; with error bars (boxes) corresponding to statistical (systematic) uncertainties.The global normalisation uncertainty on the number of binary collisions is indicated by the solid boxes in the center left of the plot (left to right: CMS, ALICE).The normalisation uncertainty for the ATLAS data is included in the systematic uncertainties.

Fig. 9 .
Fig. 9.The nuclear modification factor RpA for π mesons as a function of final transverse momentum pT , is calculated with and without the short pathlength correction to the radiative energy loss.The RpA is calculated both with collisional and radiative energy loss (el.+ rad.), and with radiative energy loss only (rad.only).Data from ATLAS [39] for charged hadrons is shown for comparison, where statistical (systematic) uncertainties are represented by error bars (boxes).The global normalisation uncertainty on the number of binary collisions is indicated by the solid box in the top left corner of the plot.

Figure 10
Figure10investigates the large formation time assumption.The figure shows the expectation value of ω 1 /µ 1 = (k − q 1 ) 2 /2xE µ 2 + q2  1 , where the DGLV and correction derivations assume that ω 1 /µ 1 ≪ 1, using L = 5 fm, λ g = 1 fm, and µ = 0.5 GeV.The large formation time assumption is explicitly violated for the energy loss with and without the short pathlength correction.For the energy loss without the correction, the large formation time assumption is violated for E ≳ 100 GeV for both charm quarks and gluons, while for the energy loss with the correction the large formation time assumption is violated for E ≳ 35 GeV for gluons and for E ≳ 50 GeV for charm quarks.This breakdown calls into question the validity of the large formation time assumption in DGLV radiative energy loss for p T ≳ 100 GeV, regardless of whether the energy loss receives a short pathlength correction.The increased rate and magnitude of large formation time assumption violation once the short pathlength correction is included in the energy loss indicates that the short pathlength corrected energy loss model predictions may be breaking down at moderate to

Fig. 11 .
Fig. 11.Plot of ⟨R⟩ ≡ ⟨ω1/µ1⟩ as a function of parent parton energy E. ⟨R⟩ ≪ 1 implies consistency with the large formation time assumption.⟨R⟩ is computed without (dashed) and with (solid) the short pathlength correction for charm quarks [gluons] with scattering centers distributed according to the exponential distribution (dark red [orange]) and truncated step function (dark blue [light blue]).L = 5 fm in the top pane and L = 1 fm in the bottom pane.All curves use constant λg = 1 fm and µ = 0.5 GeV.

Fig. 12 .
Fig. 12. Plot of ⟨R⟩ ≡ ⟨ω0/µ1⟩ as a function of parent parton energy E. ⟨R⟩ ≪ 1 implies consistency with the large formation time assumption.⟨R⟩ is computed without (dashed) and with (solid) the short pathlength correction for charm quarks [gluons] with scattering centers distributed according to the exponential distribution (dark red [orange]) and truncated step function (dark blue [light blue]).L = 5 fm in the top pane and L = 1 fm in the bottom pane.All curves use constant λg = 1 fm and µ = 0.5 GeV.

Fig. 13 .
Fig. 13.Plot of ⟨R⟩ ≡ ⟨k − /k + ⟩ as a function of parent parton energy E. ⟨R⟩ ≪ 1 implies consistency with the collinear approximation.⟨R⟩ is computed without (dashed) and with (solid) the short pathlength correction for charm quarks [gluons] with scattering centers distributed according to the exponential distribution (dark red [orange]) and truncated step function (dark blue [light blue]).L = 5 fm in the top pane and L = 1 fm in the bottom pane.All curves use constant L = 5 fm, λg = 1 fm, and µ = 0.5 GeV.

Fig. 14 .
Fig. 14.Plot of ⟨R⟩ ≡ ⟨x⟩ as a function of parent parton energy E. ⟨R⟩ ≪ 1 implies consistency with the soft approximation.⟨R⟩ is computed without (dashed) and with (solid) the short pathlength correction for charm quarks [gluons] with scattering centers distributed according to the exponential distribution (dark red [orange]) and truncated step function (dark blue [light blue]).L = 5 fm in the top pane and L = 1 fm in the bottom pane.All curves use constant L = 5 fm, λg = 1 fm, and µ = 0.5 GeV.

Fig. 15 .
Fig. 15.Plot of ⟨R⟩ ≡ ⟨1/∆z µ⟩ as a function of parent par-ton energy E. ⟨R⟩ ≪ 1 implies consistency with the large pathlength assumption.⟨R⟩ is computed without (dashed) and with (solid) the short pathlength correction for charm quarks [gluons] with scattering centers distributed according to the exponential distribution (dark red [orange]) and truncated step function (dark blue [light blue]).L = 5 fm in the top pane and L = 1 fm in the bottom pane.All curves use constant λg = 1 fm and µ = 0.5 GeV.

Fig. 16 .
Fig. 16.Plot of ⟨R⟩ ≡ ⟨∆z ω0⟩ as a function of parent parton energy E. ⟨R⟩ ≫ 1 implies consistency with the large formation time (with respect to scattering centers) assumption used in the short pathlength correction derivation [37].⟨R⟩ is computed with the short pathlength correction for charm quarks for pathlengths of 5 fm [1 fm] with scattering centers distributed according to the exponential distribution (dark red [dark blue]) and truncated step function (orange [light blue]).All curves use constant λg = 1 fm and µ = 0.5 GeV.

Fig. 17 .
Fig. 17.Plot of the RAA for D mesons as a function of final transverse momentum pT in √ s = 5.02 TeV Pb+Pb central 0-10% (top) and semi-central 30-50% (bottom) collisions.Predictions with (solid) and without (dashed) the short pathlength correction to the radiative energy loss are shown using the exponential (red) and truncated step (blue) distributions for the scattering centers.Data are from ALICE [40] and CMS [102].The global normalization uncertainty on the number of binary collisions is indicated by solid boxes in the top left corner of the plot (left to right, top to bottom: CMS 0-10%, ALICE 0-10%, ALICE 30-50%).

Fig. 18 .
Fig. 18.The nuclear modification factor RAA as a function of final transverse momentum pT is calculated in Pb + Pb collisions at √ s = 5.02 TeV for B mesons.Predictions with (solid) and without (dashed) the short pathlength correction to the radiative energy loss are shown using the exponential (red) and truncated step (blue) distributions for the scattering centers.Data are from CMS [102].The experimental global normalization uncertainty on the number of binary collisions is indicated by the solid box in the top left corner of the plot.

Fig. 19 .
Fig. 19.The nuclear modification factor RpA for D mesons as a function of final transverse momentum pT in 0-10% central p + Pb collisions at √ s = 5.02 TeV.Only radiative energy loss is included; predictions with (solid) and without (dashed) the short pathlength correction are shown using the exponential (red) and truncated step (blue) distributions for the scattering centers.Data are from ALICE [103].The experimental global normalization uncertainty on the number of binary collisions is indicated by the solid box in the top left corner of the plot.

Fig. 20 .
Fig. 20.Plot of the RAA for π mesons produced in 0-5% most central Pb + Pb collisions at √ s = 5.02 TeV, as a function of the final transverse momentum pT .Predictions with (solid) and without (dashed) the short pathlength correction to the radiative energy loss are shown using the exponential (red) and truncated step (blue) distributions for the scattering centers.Data from ATLAS [104], CMS [105], and ALICE [106].The experimental global normalization uncertainty on the number of binary collisions is indicated by the solid boxes in the center left of the plot (left to right: CMS, ALICE).Note that in the ATLAS data the normalization uncertainty is included in the systematic uncertainty.

Fig. 21 .
Fig. 21.The nuclear modification factor RpA for π mesons as a function of final transverse momentum pT in 0-10% p + Pb collisions at √ s = 5.02 TeV.Only radiative energy loss is included; predictions with (solid) and without (dashed) the short pathlength correction are shown using the exponential (red) and truncated step (blue) distributions for the scattering centers.Data for charged hadrons are from ATLAS [39].The experimental global normalization uncertainty on the number of binary collisions is indicated by the solid box in the top left corner of the plot.

Figure 21 shows
Figure21shows R pA (p T ) for π mesons and charged hadrons in 0-10% central p + Pb collisions at √ s = 5.02 TeV.For our theoretical predictions of pion suppression we only include radiative energy loss, as the average elastic energy loss of the WHDG model is inappropriate to use here.We show predictions with (solid) and without (dashed) the short pathlength correction to the radiative energy loss, and we show predictions when using either the exponential (red) or truncated step (blue) distribution of scattering centers.Charged hadron suppression data is from ATLAS[39].The difference between the R π pA (p T ) predictions from the two scattering center distributions is small when excluding the short pathlength correction.We see again that the effect on nuclear modification factor from the short pathlength correction to the radiative energy loss is large when using the exponential distribution of scattering centers and relatively small when using the truncated step distribution of scattering centers.The prediction of enhancement by the corrected R pA with an exponential distribution is qualitatively similar to the observed enhancement for moderate momenta p T ≲ 60 GeV.If we are provocative, we may thus suggest that the experimentally measured excess in R pA (p T ) is actually due to final state effects rather than initial state or normalization effects.