Searching for QGP droplets with high-$p_T$ hadrons and heavy flavor

The search for the smallest quark-gluon plasma (QGP) droplets in nature has motivated recent small collisions system programs at RHIC and LHC. Unambiguous identification of jet quenching due to final-state interactions is key to confirming QGP formation in these reactions. We compute the nuclear modification factors $R_{AA}$ and $R_{p(d)A}$ of charged hadrons and heavy flavor mesons in large (Au-Au, Xe-Xe, Pb-Pb) and small ($d$-Au, $p$-Pb, O-O) colliding systems, respectively. Our results include the Cronin effect and initial-state parton energy loss in cold nuclear matter. In the final state, hard partons undergo collisional energy loss and branching that was recently derived using Soft-Collinear-Effective-Theory with Glauber Gluon (SCET$_{\rm G}$). In large colliding systems, medium-modified QCD evolution of the fragmentation functions dominates the nuclear correction. As the system size decreases, we find that cold nuclear matter effects, collisional energy loss, and QGP-induced radiations can become equally important. A systematic scan over the medium size and mass/flavor dependence of $R_{AA}$ provides the opportunity to separate these individual contributions and identify QGP signatures in small systems. Predictions for $R_{AA}^{h}$, $R_{AA}^{D}$, $R_{AA}^{B}$ in O-O collisions at $\sqrt{s}=7$ TeV are presented with and without the formation of a QGP and contrasted with the corresponding $R_{p(d)A}$ calculations. Upcoming single-hadron measurements at the LHC will not only test the O-O predictions for both light and heavy flavor production, but will shed light on the possibly very different dynamics of $p$-A and A-A reactions at similar soft particle production multiplicities.

Recently, puzzles have emerged in small colliding systems such as d-Au and p-Pb collisions. Similar "collective" behavior in the pattern of soft particle production that is attributed to QGP evolution in large systems has been observed [44][45][46][47][48][49], suggesting the possibility of final-state effects. However, clear evidence of jet quenching has not been observed [50][51][52][53][54][55][56]. This puzzle has led to intense discussion of the origin of the apparent collectivity [57][58][59][60] and the nature of the medium produced in small systems. Sensitive experiments have been designed to provide further insight into this problem, such as the geometry scan at RHIC using p-Au, d-Au, and 3 He-Au [61][62][63], and the upcoming O-O and p-O collisions program at the LHC [64,65]. Theoretical predictions of jet quenching in O-O collisions already exist [66][67][68][69] for light and heavy-flavor quenching, using various frameworks, including transport equations and energy loss calculations. The QGP effects are modeled by medium-induced gluon radiations using either higher-twist or the BDMPS-Z formula, some studies also included collisional processes. Nevertheless, what is still missing from the theory side is an analysis that combines QGP and full cold nuclear matter (CNM) effects and a comparison and contrast of results in p(d)-A and A-A reactions. This will be important to better understand the baseline without QGP formation and identify possibly different dynamics in symmetric vs asymmetric small systems. Finally, we use the in-medium QCD evolution formalism to consistently treat final-state parton shower effects for light and heavy flavor and understand the interplay with collisional energy loss 1 .
This paper aims to analyze the nuclear modification factors in various large and small systems. Calculations will be performed for both light hadrons and heavy-flavor mesons effects including both cold nuclear matter effects and hot QGP final-state effects, collisional energy loss and medium-induced radiation corrections to the baseline QCD factorization formalism. By fixing the jet-medium interaction parameter g s in large colliding systems, we make predictions for light and heavy-flavor productions in d-Au, p-Pb, O-O, with and without the assumption of the existence of a QGP.
The motivation for this comprehensive analysis is threefold. First, although initialstate cold nuclear matter effects are overwhelmed by the QGP-induced modifications in large systems, they can be important in small systems and provide an alternative explanation of the structures of modifications without always resorting to the existence of a QGP. Second, the radiative correction can be strongly reduced relative to the collisional effect in a QGP of decreasing size, as we will demonstrate in this paper. Therefore, a complete treatment of hot medium effects has to include elastic collisions. Using the distinct mass dependence of radiative and collisional processes, the difference between light and heavy flavor modifications in small systems can provide a handle on the relative contribution of the two. Third, a notable source of uncertainty that has hindered theoretical analyses in Hard Hard Figure 1. Left: illustration of cold nuclear matter effects. Vertical gluon lines represent multiple collisions between initial-state partons with the other nucleus. Horizontal lines represent the modified gluon radiation in the initial state. Right: illustration of final-state effects of jets in the quark-gluon plasma.
small systems is the decorrelation of charged particle production, which defines centrality classes, and the nuclear overlap function T pA . The large model-dependence in T pA makes it hard to interpret the present jet p-Pb modification data. By looking into various collision geometries, especially the O-O program at the LHC, the normalization uncertainty is expected to be significantly reduced.
The paper is organized as follows. In section 2 we review the factorization calculation in nuclear collisions and illustrate final-state and initial-state effects that will be considered. In section 3 we describe the dynamical approach for the cold nuclear matter effect. In section 4, we review in-medium QCD splitting functions obtained in SCET G , the modified QCD evolution for in-medium fragmentation function, and collisional energy loss in a thermal QGP. The model of QGP medium simulation is introduced in section 5. Results for R AA in large and small systems are presented in section 6. Finally, we present our conclusions in section 7. How to constrain the parameters in the hydrodynamic simulation is discussed in appendix A. We comment on the differences between dynamically calculated CNM effects and nPDF parametrization in appendix B.

Factorization approach with initial and final-state effects
Calculations of hadron and jet production in reactions with nuclei are based on incorporating medium corrections into the QCD factorization approach. The baseline p-p transverse momentum (p T ) and rapidity (y) differential hadron production cross-section is is the fragmentation function of parton k into hadron h carrying momentum fraction z. Final-state effects in the QGP modify D h/k (z, µ F ), and the fragmentation will generally depend on the parton energy E in the rest frame of the medium in addition to the medium transport properties. dσ k dq 2 T dy is the production cross-section of the hard parton and accounts for many-body scattering effects in large nuclei. It can be expressed as Equation 2.2 accounts for the fact that the initial parton can acquire a finite transverse momentum k i , k j from multiple collisions with the other nucleus, which we treat in the Gaussian approximation [20,71] that often used to describe the Cronin effect [72]. x i and x j are the longitudinal momentum that we find from Another impact of cold nuclear matter is that initial-state partons can lose fractions (∆x i /x i and ∆x j /x j ) of their energy due to CNM-induced gluon emissions [73]. The physical picture of equation 2.2 and 2.1 are illustrated in the left and right panel of figure 1. Finally, the coherent scatterings also lead to dynamical shadowing effect that shifts the x i by an amount proportional to the nuclear thickness function [74]. It only contributes at small Bjorken-x and low transverse momentum We will elaborate upon these CNM effects in section 3. Finally, cos θ c.m. = tanh η c.m. is the polar angle in the center-of-mass frame. The integration has been transformed to the rapidity of the produced particle in the partonic center-of-mass frame η c.m. . We will set the factorization and renormalization scale to the transverse momentum of the hard parton µ R = µ F = |q| = |p|/z.

The dynamical approach for cold nuclear matter effects
Despite the complicated nature of parton interaction in the nuclear environment, it is possible to model the initial-state effect from QCD interactions. In reference [73], one considers coherent multiple scatterings (vertical gluons lines in the left panel of figure 1) and induced soft gluon radiation (horizontal gluons lines in the left panel of figure 1) from initial-state partons as they traverse in the nuclear matter before the hard interaction. The collisions push particle production to slightly larger transverse momentum, resulting in an enhancement at p T around a few GeV, known as the Cronin effect [75]. The interaction potential between partons and cold nuclear matter (CNM) is chosen to be a screened Coulomb potential with typical transverse momentum transfer squared µ 2 = 0.12 GeV 2 and the mean free paths λ g = (C F /C A )λ q = 1.5 fm. As mentioned earlier, we employ the Gaussian parametrization of the Cronin effect and account for the power-law tails of the Moliere multiple scattering with a numerical factor ξ ∼ few. The Cronin effect is sensitive to the shape of the particle spectra and decreases at large colliding energies. For phenomenological applications we also consider 50% shorter mean free paths to estimate its uncertainty. At small transverse momenta and small values of Bjorken-x coherent multiple scattering also leads to dynamical shadowing [74,76], which we include in the calculation. Coherent power corrections scale as ∆x i /x i ∼ µ 2 A 1/3 /(−u) and ∆x j /x j ∼ µ 2 B 1/3 /(−t) in A-B reactions and t and u are the partonic Mandelstam variables for the hard scattering process. The Cronin and dynamical shadowing effects disappear at high p T .
On the contrary, the induced gluon emissions that causes initial-state parton energy loss in CNM continues to be important at high energy [73] x dN IS as can be seen from its weak dependence on the lightcone momentum p +2 . x is the momentum fraction carried away by the radiated gluon. L is the path length that parton propagates in the cold nuclear matter before the hard collision. The same transversemomentum transfer squared µ 2 and mean free path λ g are used in equation 3.1 as those for the Cronin effect. The CNM energy loss also causes a shift in the momentum fraction x of the initial parton in equation 2.2 At high energy, the CNM energy loss contribution dominates ∆x and is proportional to L, introducing the centrality dependence. Fluctuations due to multiple gluon emissions reduce the effect of the mean fractional energy loss, which we account for with fl < 1. For steeply falling final-state spectra fl can be as small as 0.4 [77]. The energy dependence of the initial parton flux is much more moderate and we use fl = 0.7. For phenomenological applications, we also consider the scenario without initial-state energy loss. The calculation that includes the Cronin effect, dynamical shadowing, and the CNM energy loss, hereafter referred to as the "dynamical CNM calculation" or "Cronin+eloss", will be used as the primary model for how the nuclear environment affects the initial parton density f (x, k). The advantage of the dynamical approach is that one can use only two parameters that control the magnitude of broadening and initial-state energy loss to systematically study the energy and A dependence predicted by QCD. In the large x region, the current calculation only implements the isospin effect without parametrizing the antishadowing, EMC, and Fermi motion regions. We expect these effects to be small for the moderate p T regions of hadron production considered in this paper.
In figure 2 we compute the spectra of partons produced in the hard interaction in Au-Au collisions relative scaled p-p baseline at √ s = 200 GeV. The blue bands are dynamical model calculations with Cronin and dynamical shadowing effects only, and shaded bands are results that further include the CNM energy loss. The spread of the bands represents 50% variation in the magnitude of transverse momentum broadening. Scattering in nuclear matter introduces a nuclear size-dependent enhancement at moderate p T , while depleting particle production below 2 GeV. Dynamical shadowing also contributes to this low-p T suppression [74]. The CNM energy loss suppresses the spectra at large p T . The dynamical CNM calculations are compared to results using collinear nuclear parton distribution functions (nPDF) from the (n)NNPDF Collaboration [78] (black dash-dotted lines). There is no momentum broadening in the collinear nPDF. The small low-p T suppression comes from the parametrized shadowing effect. At large p T , modifications results from the anti-shadowing, EMC, and the Fermi motion effects included in the nPDF. Unlike the dynamical CNM model, nNNPDF does not depend on the impact parameter, which is essential to include to study the centrality dependence of R AA and R pA in small colliding systems. Therefore, we primarily use the dynamical CNM model in this paper. We will compare the impact of using nPDF to the final results in appendix B.

Final-state QGP effects
4.1 Medium-modified splitting functions from SCET G An effective theory of QCD ideally suited to studying jet physics is Soft-Collinear-Effective-Theory (SCET) [79,80]. The power counting parameter λ = |k|/p + can be thought of as the typical transverse momentum in the jet divided by its large lightcone component. In reactions with nuclei, the SCET G theory [26,27] was developed to couple the collinear fields to the background nuclear medium via Glauber gluon exchanges. Thus, hadron and jet production, and jet substructure can be described in different strongly-interacting environments without loss of generality. This framework has been applied to study both the jet broadening [81] and derive the medium-modified QCD splitting functions [30,82]. For phenomenological application, the Glauber gluon field is often approximated by a sum of screened color potential over the scattering centers in the medium Here, R denotes the color representation of the collinear parton while i is the representation of the color charge of the medium quasi-particle. g s is the jet-medium coupling constant.
is the Deybe screening mass of the plasma. There is no momentum exchange in the p + direction that scales as λ or stronger. The phase factor contains the position (z, z + )information of the medium color charge. The differential collision rate between a collinear quark (or a gluon) and the QGP is and we choose an effective N f = 2 for the QGP. The modified QCD splitting functions have been obtained to first order in opacity explicitly [27,30] and an iterative approach has been developed to generalize them to higher opacity orders [83]. For light partons, we take the full splitting function to first order, which we quote for completeness from Ref. [30]. The double differential spectrum of the quark to quark+gluon branching is where the gluon carries momentum fraction x. Note that this differs from the standard high energy notation and is done to make contact with the much studies energy loss soft gluon emission limit when x → 0. The quark to gluon+quark splitting function is obtained by . For gluon to quark+quark and gluon to gluon+gluon splittings, the distributions are In the above expressions, we have defined vectors and frequencies (which are inverse formation times) P qq , P gq , P gg , P qg are the standard splitting functions in the vacuum at leading order, The coupling constants associated with the vacuum splitting functions are evaluated at µ 2 = k 2 in this study, while the jet-medium coupling parameter g s that goes into the collisions rates dR/dq 2 is taken as a free parameter. The ∆z integration in equations 4.
The x-dependence can be significantly modified compared to light flavors when k 2 M 2 . The formulas for the medium corrections to equations 4.8-4.10 have been derived in reference [28]. Though we will not write down the full in-medium expressions, we emphasize that the heavy quark mass not only introduces corrections to the propagators and the interference phases but also generates many new terms proportional to the mass. We demonstrate the impact of mass corrections in the medium with numerical results, as shown in figure 3. In figure 3, we present the medium-correction to q → q + g splitting functions in 0-5% central Pb-Pb collisions at 5.02 TeV and in 0-5% O-O at 7 TeV. We have factored out the the vacuum-like P qq and P bb kernels given by equation 4.7 and equation 4.8. Medium corrections to light parton branching are suppressed by the Landau-Pomeranchuk-Migdal interference effect at large x. For the bottom quark, the vacuum splitting functions P bb further shows the so-called dead-cone effect that suppress radiations in the phase-space region x > |k|/M relative to P qq . However, after factoring out the vacuum factor, the ratio still displays strong mass modifications [43]. Only ratios in the energy loss region x → 0 are comparable to that of the light quark. Looking at the |k|-dependence in figure 3 or simply at equations 4.3, medium-induced branchings are suppressed at large k by at least another power of 1/k 2 as compared to the vacuum radiation. So, in principle, they do not contribute additional ln Q 2 enhancement upon integration for asymptotic energies. Its contribution peaks when |k| are comparable to the typical size of |q|. In central Pb-Pb collisions, the correction at its peak can be much larger than the vacuum contribution. In central O-O collisions at 7 TeV, it is estimated to be about a factor of three smaller than that in central Pb-Pb at 5 TeV, assuming QGP effects do exist in O-O.

Collisional energy loss
In obtaining equations 4.3 and 4.4, one assumes that p + is conserved. This is a good approximation for the computation of radiative correction because a mismatch in p + due to collisional energy losses is expected to be subleading in powers of λ. However, collisional energy loss should be taken into account to compute R AA at intermediate and small p T [10,15,84]. It was found that it can be comparable to the radiative energy loss for partons up to p T = 10 to 20 GeV/c. For heavy flavor particles, due to the reduced phase space for radiation when p T is only a few times the heavy quark mass M , the dead cone effect renders collisional energy loss even more important to describe heavy meson suppression at the intermediate p T . Finally, as we will see immediately, the different path length dependence of induced radiation and collisions energy dissipation makes the latter an indispensable component in the analysis of small collision systems.
We take the collisional energy loss obtained in hard-thermal loop calculations [7]. The energy loss of a quark per unit length in a weakly-coupled thermal plasma is given by, with v = p/E being the velocity of the parton in the rest frame of the QGP, applied to both heavy and light quarks. Gluon energy loss is related to that of the quark via the C A /C F quadratic Casimir ratio. The running coupling value at one loop is used in the above expression [84], where the maximum α s is cut-off at g 2 s /(4π). We remark that the running coupling effect cancels the ln(ET /m 2 D ) enhancement from phase-space integration and results in an approximately energy-independent collisional energy loss.
In this study, we use averaged collisional energy loss ∆E = x ⊥,0 +∆zmaxφ x ⊥,0 dE coll d∆z d∆z obtained by averaging over the production location x ⊥,0 and orientationφ of a parton with initial energy E. We will apply the averaged energy loss to the partons created in the hard process before invoking the in-medium splitting function modification of the fragmentation function. The justifications of this approximation are that 1) collisional energy loss is almost E-independent, and 2) a ∆E el mismatch will only cause a small difference in the QCD evolution so long as ∆E el Q = |q|.

Medium size dependence of radiative and collisional energy loss
While we employs the full expressions (equations 4.3 and 4.4) in the calculation, it is instructive to look at the radiative parton energy loss obtained in the x → 0 limit. The energy loss fraction ∆E/E = x dN med dx dx in this approximation is with u t = µ(∆z) 2 ∆z 2xE and the second line for the asymptotic behavior at high energy. Compared to the scaling of the elastic energy loss fraction one sees that radiative energy loss only dominates over the elastic one by α s ln E. Furthermore, they scale differently with medium size. For example, in a QGP that undergoes Bjorken expansion such that T 3 τ = T 3 0 τ 0 and ∆z = τ − τ 0 , the typical momentum transfer and the mean free path evolve with proper time as µ 2 = µ 2 0 (τ 0 /τ ) 2/3 and λ g = λ g,0 (τ /τ 0 ) 1/3 . Therefore, the radiative energy loss fraction only scales linearly with size; while the elastic energy loss fraction changes much slower with L. Consequently, one expects that collisional processes becomes increasingly important in a small-sized QGP.

Modified QCD evolution equations for in-medium fragmentation
To compute hadron production in a nuclear environment, we take the modified DGLAP approach [85] to evolve the vacuum fragmentation function from an initial non-perturbative scale Q 0 to Q = p T + ∆E el with medium modified QCD splitting functions described in section 4.1.
Evolution in the vacuum. The evolution equation for the fragmentation function D 0 h/i of hadron specie h from parton i produced in the vacuum is The above equations solve to c HH (r) = 1 1 + r 2 + 2r 2 + 1 2(1 + r 2 ) 2 + 2r 2 1 + r 2 − 2 ln 1 1 + r 2 . (4.25) Finally, the relation between Q 2 and x, k and the allowed kinematic ranges are summarized in table 1 for each channel.
Channel Table 1. Q 2 in terms of splitting kinematics and the kinematic ranges of Q 2 (or of x).
Evolution in the medium. For the QCD evolution in the medium, the splitting functions and virtual corrections in equation 4.16 are replaced by the medium-modified ones, and The medium-induced correction to the splitting function does not introduce a ln Q 2 dependence as large as the vacuum one because it decays as 1/k 4 or faster at large k 2 and is screened by the Debye mass m D at small k 2 . In the low-Q 2 region, the medium modifications become comparable to or can even dominate over the vacuum contribution. When this happens, one has to consider how to implement such corrections. For example, in [86], the authors only applies the modified DGLAP equation to the region Q > 1 GeV, and use the medium-modified QCD splitting function below 1 GeV to build an in-medium  initial condition of the evolution. Other methods treat the low-Q 2 region in a transport approach [87][88][89][90] where multiple medium-induced emissions are generated sequentially in a time-ordered fashion because the number of soft emissions is enhanced by the medium size. In our calculation, we notice that with the choice Q 2 = k 2 /(x(1 − x)), Q 2 is inversely proportional to the formation time of the branching τ f = p + /Q 2 . For soft emissions that do not significantly change p + , the Q 2 -ordered evolution is the same as a formation-time ordered approach to compute radiative parton energy loss. However, they are not the same for energetic splittings that take a large fraction of the parton energy. We further remark that the QCD evolution approach can be applied to regions where the branching fraction x is large and regions where the medium-correction is negative due to interference, which are beyond the scope of the transport equation. This may be important for small collisions systems, as the interference effects are very sensitive to a small path length. Therefore, we consider the QCD evolution approach to be a better choice for the system-size scan down to small collisions systems such as O-O and p-Pb.
Both the vacuum and in-medium evolution use vacuum fragmentation functions at Q 0 = 0.4 GeV as the initial condition. The QCD evolution approach, similar to the traditional energy loss approach, assumes that hadronization happens outside of the medium. However, for heavy mesons, the formation time can be significantly shortened by the large mass [94]. This may have additional phenomenally consequences and we will come back to  this point in section 6.2. We take the charged pion fragmentation functions as parametrized in reference [95] at Q = 1 GeV. They are evolved backward from Q = 1 GeV to Q 0 = 0.4 GeV to provide the common initial condition for evolution in the vacuum and medium. For heavy mesons, we use the Lund-Bowler function [96] as the initial condition, with parameters a = 0.68, b = 0.98 taken from Pythia8 [97] and M 2 ⊥ ≈ M 2 h + (0.7 GeV) 2 . We have verified that the evolved heavy meson fragmentation functions with Lund-Bowler type initial conditions provide a reasonable description to heavy-meson fragmentation measurements by CLEO [91] Collaboration at Q = M Υ(4S) /2 and ALEPH [92,93] Collaboration at Q = M Z /2. In figure 4, the red lines are the initial conditions and the blue lines are results evolved to Q = p T = M z /2 for ALEPH experiment and Q = p T = M Υ(4S) /2 for CLEO experiment. The blue bands denote variation p T /2 < Q < 2p T .
In figure 5, we compare the fragmentation functions in p-p (blue solid lines) and in 0-5% central Pb-Pb collisions (red dashed lines and bands) for four channels evolved from Q = Q 0 to Q = p T = 50 GeV. The blue solid lines are the evolved results in the vacuum, and the red dashed lines with bands are results evolved in the medium with g s = 1.6, 1.8, 2.0. Compared to D(z) in the vacuum, medium-modified DGLAP evolution "red-shifts" D(z). In the calculation of inclusive hadron spectra, D(z) is always folded with a steep falling partonic cross-section dσ/dq ∼ 1/q n , n 1. The resulting spectra depend on the integral 1 0 z n−1 D(z)dz. Because n 1, R AA of both light and heavy mesons are mostly sensitive to the modification in the region z 1, even though the pion fragmentation functions (top row) are much softer than those of heavy mesons (bottom row).

Dynamical simulations of QGP and its existence in small systems
Finally, we discuss the model for the medium evolution. The dynamical simulation of the quark-gluon plasma produced in nuclear collisions is performed using the Duke hic-eventgen code package [98]. In this calculation, the TRENTo initial condition model of the collision geometry provides the energy deposition profiles at the proper time τ = 0 + . The model of quark-gluon plasma dynamics consists of a pre-equilibrium stage modeled by free-streaming [99], followed by the 2+1D boost-invariant relativistic viscous hydrodynamics [100,101]. Finally, the hydrodynamic fields are particlized into hadrons at transition temperate T sw slightly below the pseudo-critical temperature of the QGP equation of state [102], and the hadronic interactions are handled by the Ultra-relativistic-Quantum-Molecular-Dynamics (UrQMD) [103,104]. The model parameters have been tuned to the experimental measurement of particle production, flows, and correlations in previous studies [98].
In the current study, it will be very computationally intensive to obtain the full splitting functions and perform DGLAP evolution on an event-by-event basis. Therefore, we simulate events with centrality-averaged initial conditions. For this reason, the scale parameter (normalization) of the TRENTo energy deposition model is re-tuned for each system to reproduce the centrality-dependent charged particle yield and transverse energy. Because no data is available for O-O collisions at 7 TeV, we interpolate the normalization tuned at RHIC and LHC energies using a third-degree polynomial in ln √ s to predict the normalization at 7 TeV. The details can be found in appendix A. Now, we make an important remark on the QGP effects in small systems. This may naively seem to be an unnecessary discussion because the hydrodynamic simulations already provide the time evolution of the temperature profile. One can, in principle, use this to compare the medium temperature and the QGP pseudo-critical temperature T c to determine if the jet propagates in the QGP phase or in a hadronic phase. In fact, for those high-multiplicity events in small system collisions, the simulated medium temperature starts from a point well above T c . However, the definition of temperature in the hydrodynamics-based simulation bears certain ambiguity when the system is far from equilibrium, for example, in small colliding systems. In the hydrodynamic picture, energy density is converted into temperature using the lattice QCD EoS, which means that the number of scattering centers defined in this manner would approach the thermal limit. If the system is far from equilibrium, the density of scattering centers can significantly deviate from this expectation. In this study, we will therefore investigate two extreme limits of small colliding systems.
• Calculations with cold nuclear matter effect only.
• Calculations with cold and hot medium effects that assume the QGP is described by the hydrodynamic-based model. We will let future experiments falsify either scenario.

Results and discussion
In the result section, we first illustrate the type of modification cold nuclear matter effect and QGP effects may induce in the cross-section ratio R AA . We then fix a range for the jet-medium coupling parameter g s in large colliding systems Au-Au, Pb-Pb, and Xe-Xe.
With the same set of parameters, we then present predictions for R h AA , R D AA , and R B AA in d-Au, p-Pb, and O-O collisions. These calculations are performed with the dynamical cold nuclear matter effects. We discuss the impact of using the dynamical approach and the nuclear PDF in appendix B.

Interplay of cold (initial-state) and hot (final-state) medium effects
In figure 6, we sequentially include the contribution from Cronin effect and coherent power corrections, CNM energy loss, elastic and radiative effect in the QGP in 0-10% central Pb-Pb and O-O collisions (rows). The parameter used in this demonstration is g s = 1.8.
Columns from left to right show the modifications for charged hadrons, D mesons, and B mesons. The dashed lines include only Cronin momentum broadening and the peak for light hardons is around 3 GeV. For heavier mesons, it moves to slightly higher p T . The inclusion of CNM energy loss (black dash-dotted lines) results in an overall suppression as the CNM energy loss fraction is almost independent of energy (see equation 3.1). CNM effects are much smaller in O-O collisions than those in Pb-Pb collisions, as expected from the A 1/3 scaling.
In Pb-Pb collisions, the modified QCD splitting functions (calculations shown as blue dashed lines) lead to significant suppression of light hadron production. The further inclusion of collisional energy loss (solid blue lines) is a sub-leading effect in ln(E). For heavy mesons, the radiative correction in the region p T < 5M is strongly suppressed, and the modifications are largely attributed to collisions energy loss. One should be careful, however, when interpreting the heavy-flavor results at p T M . In this region, the "jet approximation" E M completely breaks down, and the heavy quark's orientation can change randomly as it "diffuses" in the QGP. The low p T regime is better modeled by Langevin or other transport approaches that fully evolve the phase-space density of the heavy quark [105][106][107]. In O-O collisions, there is a notable change in the relative importance of radiative correction and collisional energy loss. This can be understood from the different medium-size scaling of the two processes as discussed in section 4.3. Especially for heavy meson R AA , collisional processes are responsible for at least 50% of the QGP modifications in 0-10% O-O events.

Nuclear modifications in Pb-Pb, Xe-Xe, and Au-Au
For simplicity of the uncertainty estimation we vary the jet-medium coupling g s = 1.8 ± 0.2. In principle, g s may run with the medium scale and this has been investigated in other studies [90]. In figure 7, calculations with bands that show the sensitivity to the interaction strength are given for R h AA (top row), R D AA (middle row), and R B AA (bottom row) in Pb-Pb collisions at 5.02 TeV (left column), Xe-Xe collisions at 5.44 TeV (middle column), and Au-Au collisions at 200 GeV (right column). We include the Cronin effect, cold nuclear matter energy loss, and coherent power corrections. Within each panel, the nuclear modification in 0-10% and 30-50% central collisions are shown and compared to available data from the ALICE [108], ATLAS [109], CMS [54,110], PHENIX [111,112], and STAR Collaboration [113]. For comparisons to b-decay electrons and non-prompt J/Ψ, smearing functions extracted from Pythia8 [97] simulations are applied to the B meson spectra.
For light hadron suppression, the range 1.6 < g s < 1.8 provides a good description of the LHC data in Pb-Pb and Xe-Xe collisions. At RHIC energy the data suggests a larger coupling 1.8 < g s < 2.0. This trend is consistent with many other findings that the effect of GeV for centrality classes 0-10% and 40-50%. The calculations are compared to measurements by the ALICE [108], ATLAS [109], CMS [54,110], PHENIX [111,112], and STAR [113] experiments. For bottom flavor, the B-decayed J/Ψ is computed for the LHC energies and B-decayed electron is presented at RHIC energy. Note that the centrality classes for the STAR measurements of D meson is 0-10% and 40-80%.
jet-medium interactions is larger at lower temperatures relevant for collisions at the RHIC beam energy [36,114].
Switching to the flavor/mass dependence of the suppression, the calculation agrees well with D-meson suppression at high transverse momentum but slightly overestimates R D AA at low p T . Data tend to lie on the lower edge of the band. The overestimation (not enough suppression) is evident for the bottom quarks. However, we are not going to tune a separate set of parameters for the heavy sector in this paper. Instead, tension with data is a useful indicator of physics that might be missing in the calculation. One possible explanation for the systematic deviation at low p T with increasing quark mass is the collisional dissociation of heavy mesons in the QGP [94,115]. In section 4.4, we have argued that after the evolution the light parton fragmentation should take place outside of the QGP medium. However, the formation time of low-p T heavy mesons is considerably shorter so that they can be produced inside the nuclear medium. The calculation in reference [94] considers the collisional broadening and break-up of the D and B in the QGP that further suppresses R D AA below 10 GeV and R B AA below 30 GeV. This effect is not included in the present study. However, it should be much less important in small colliding systems to be discussed in subsection 6.3. Other works consider collisions between D and π, ρ meson in the hadronic phase [116,117], which is important in the low-p T region. In addition, reference [118] studied M/E-type drag-induced radiations of heavy quarks. These two additional effects qualitatively push the calculation in the right direction, but their overall magnitudes are too small to explain the large difference that we saw in bottom-flavor R AA .

Small systems
With the same range of parameters, we now turn to systematic predictions for small collision systems. At LHC energies there have been extensive measurements of jet production in p-A collisions. However, the interpretation of the results suffers from the ambiguity of the geometric model of nuclear collisions in the presence of large fluctuation. This situation is illustrated in figure 8 obtained using the TRENTo initial condition model used in this study. From the left to the right panel, we plot the histograms of the self-normalized "multiplicity" at the initial condition level versus the number of binary collisions (N coll ≡ T AB /σ inel pp ) for p-Pb, O-O, and Pb-Pb collisions. In large colliding systems, there is a strong correlation Figure 9. Nuclear modification factor R p-Pb compared to ATLAS data [53] scaled by the overlap functions from three different calculations of nuclear collision geometry. The results for 0-1% and 60-90% centralities are colored in blue and red, respectively. The calculations only include cold nuclear matter effects. The shaded bands include dynamical shadowing and Cronin effect, while the filled bands further reflect consideration of CNM energy loss. between the nuclear geometry and the final-state multiplicity, and the determination of N coll or T AB that normalizes R AA is less sensitive to subnucleonic modeling and fluctuations. This relation strongly decorrelates in p-Pb collisions, making the determination of T AB extremely sensitive to proton shape, fluctuations, and particle production mechanisms. One of the motivations of O-O program is to partly recover the correlation between collision geometry and the multiplicity to provide unambiguous signatures of nuclear modification in small systems.
First, we study R pA in p-Pb collisions. In figure 9 we compare theoretical predictions with only CNM effects to the ATLAS data. In theoretical calculations, we always know the correct normalization for R pA such that R pA = 1 in the absence of nuclear effects. On the experimental side, the published R pA data are strongly model-dependent: the ATLAS Collaboration obtains the normalization T pA in the conventional Glauber model (left), and the improved Glauber-Gribov model with two choices of a parameter that controls the proton fluctuation (middle and right panels) [53]. Here we label them as ATLAS T pA #1, #2, and, #3. The resulting Q pP b = dN pA→h / T pA /dσ pp→h is shown to be extremely sensitivity to the experimental choice of the nuclear geometry models.
This model dependence clearly cannot be controlled within this study. Nevertheless, we argue that it is unlikely that medium corrections can truly cause 50% enhancement of hadron production at p T = 20 GeV as suggested by geometric model # 1. We consider the geometric models # 2 and # 3 to be much more realistic from the point of view that R AA ≈ 1 at large p T . Focusing on scenarios #2 and 3, the cold nuclear matter calculation nicely explains the peak at low p T in the top 1% high-multiplicity events and its disappearance in 60-90% p-Pb collisions, though small residual enhancement remains in peripheral collisions. Also, the peak in the data is at a slightly higher p T . This "centrality" dependence comes from the nuclear-thickness (impact-parameter) dependence of the Cronin effect. Either calculations with or without CNM energy loss could be consistent with the current data in scenarios #2 and 3, though geometry #2 favors no CNM energy loss. Clearly, a better understanding of nuclear geometry in p-A is needed to further constrain cold nuclear matter effects.
Despite the large model-dependent uncertainty, the current measurements in p-Pb leave little room for the hot QGP effects discussed in section 5. In figure 10, the calculations include both cold and hot nuclear effects. QGP effects introduce a strong centrality dependent suppression of R pP b at intermediate and large p T . This is not consistent with data in either scenario. Again, we emphasize that in our calculation, the temperature of the hot QGP is determined by matching the initial condition to hydrodynamic equations using lattice EoS of state. This cannot exclude models of large non-equilibrium corrections to the density of collision centers in small systems.

Summary
In this paper, we investigated systematically the modification of light and heavy-flavor production in small and large colliding systems at moderate and high p T . Our goal was to differentiate the impact of cold nuclear matter and hot QGP effects. We performed the calculations by including in-medium corrections to the QCD factorization framework for the more elementary p-p collisions. Cronin effect, coherent power corrections and parton energy loss in the cold nuclear matter modify the initial-state parton densities, while HTL-type collisional energy loss and medium-induced radiative correction change the hadron fragmentation function in the medium. A modified DGLAP evolution approach handles the scale evolution of the fragmentation function in the medium. With jet-medium coupling g s = 1.8±0.2, the calculation results in a reasonable agreement with the light flavor suppression in A-A collisions at RHIC and LHC. This range can even accommodate the quenching of charm mesons, albeit with the largest of the studied couplings studied. However, the same range of g s underestimates the bottom meson suppression pointing to remaining tensions with the description of bottom quark dynamics even after the inclusion of collisional energy loss.
In small colliding systems, we found that the CNM effect only can already explain the basic patterns observed in p-Pb collisions scaled by the improved Glauber-Gribov model. Room for improvement in the description of such systems is available as with the CNM transport parameters used here the magnitude of the Cronin enhancement and/or cold nuclear matter energy loss can be overestimated. In order to place better constraints on parton transport in large nuclei, an improved understanding of centrality determination in p-A reactions will be greatly beneficial. In spite of the remaining uncertainties, we managed to establish that the current model of QGP formation in p-A as described by hydrodynamics, leads to quenching of hadron spectra that is inconsistent with the p-Pb data. The same cannot be said for d-Au data, but these two sets of measurements have very different and opposite high-p T behavior vs centrality. This once again points to the importance of understanding the centrality determination in highly asymmetric small-onlarge systems.
As for O-O collisions, we found that the CNM effects alone are only responsible for very small corrections, while the formation of a QGP can suppress charged particle spectra by more than a factor of two and bottom flavor up to 20%. Unlike the suppression in large systems that is dominated by induced radiation, collisional energy loss in O-O collisions leads to comparable modifications as the effect of medium-induced evolution. The predicted suppression in small systems at LHC energies with and without QGP formation is very distinct. We finally observed that if QGP quenching effects are identified in O-O, the enhanced contribution from collisional processes can be tested by simultaneously looking at the flavor dependence of R AA . with nuclear PDF. In figure 15, R AA using dynamical CNM approach (blue dotted bands) are compared to nPDF calculation (gray bands) for Pb-Pb and O-O collisions. For light hadrons and charm mesons, the major differences appear at high p T , since the dynamical approach does not include modifications in the valence region. B-meson displays a surprising sensitivity to the CNM models at low p T . Again, these differences show up in the region where M/p T = O(1), which is not a well-controlled region in the current framework. Nevertheless, this difference suggests that when hot medium effects are suppressed by the large mass of b quark, it is possible to be used to probe the details of the CNM calculation.