Inclusive jet and hadron suppression in a multistage approach

We present a new study of jet interactions in the quark-gluon plasma created in high-energy heavy-ion collisions, using a multistage event generator within the JETSCAPE framework. We focus on medium-induced modifications in the rate of inclusive jets and high transverse momentum (high-$p_{\mathrm{T}}$) hadrons. Scattering-induced jet energy loss is calculated in two stages: A high virtuality stage based on the MATTER model, in which scattering of highly virtual partons modifies the vacuum radiation pattern, and a second stage at lower jet virtuality based on the LBT model, in which leading partons gain and lose virtuality by scattering and radiation. Coherence effects that reduce the medium-induced emission rate in the MATTER phase are also included. The TRENTo model is used for initial conditions, and the (2+1)dimensional VISHNU model is used for viscous hydrodynamic evolution. Jet interactions with the medium are modeled via 2-to-2 scattering with Debye screened potentials, in which the recoiling partons are tracked, hadronized, and included in the jet clustering. Holes left in the medium are also tracked and subtracted to conserve transverse momentum. Calculations of the nuclear modification factor ($R_{\mathrm{AA}}$) for inclusive jets and high-$p_{\mathrm{T}}$ hadrons are compared to experimental measurements at the BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider (LHC). Within this framework, we find that with one extra parameter which codifies the transition between stages of jet modification -- along with the typical parameters such as the coupling in the medium, the start and stop criteria etc. -- we can describe these data at all energies for central and semicentral collisions without a rescaling of the jet transport coefficient $\hat{q}$.

We present a new study of jet interactions in the quark-gluon plasma created in high-energy heavy-ion collisions, using a multistage event generator within the jetscape framework. We focus on medium-induced modifications in the rate of inclusive jets and high transverse momentum (high-pT) hadrons. Scattering-induced jet energy loss is calculated in two stages: A high virtuality stage based on the matter model, in which scattering of highly virtual partons modifies the vacuum radiation pattern, and a second stage at lower jet virtuality based on the lbt model, in which leading partons gain and lose virtuality by scattering and radiation. Coherence effects that reduce the mediuminduced emission rate in the matter phase are also included. The trento model is used for initial conditions, and the (2+1)dimensional vishnu model is used for viscous hydrodynamic evolution. Jet interactions with the medium are modeled via 2-to-2 scattering with Debye screened potentials, in which the recoiling partons are tracked, hadronized, and included in the jet clustering. Holes left in the medium are also tracked and subtracted to conserve transverse momentum. Calculations of the nuclear modification factor (RAA) for inclusive jets and high-pT hadrons are compared to experimental measurements at the BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider (LHC). Within this framework, we find that with one extra parameter which codifies the transition between stages of jet modification-along with the typical parameters such as the coupling in the medium, the start and stop criteria etc.-we can describe these data at all energies for central and semicentral collisions without a rescaling of the jet transport coefficientq.

I. INTRODUCTION
Jet modification in high-energy heavy-ion collisions [1,2] is currently one of the leading mechanisms to study the properties of the quark-gluon plasma (QGP) [3,4] created at the Relativistic Heavy-Ion Collider (RHIC) in Brookhaven National Laboratory (BNL) and the Large Hadron Collider (LHC) in CERN. Due to the much larger momentum scales associated with partons in a jet, these partons typically exchange momenta with the medium that are much larger than the thermal momentum scale. As a result, they probe the medium at much shorter distance scales than the thermal scale [5][6][7][8][9]. The primary observable for studying jet energy loss is the nuclear modification factor, R AA , defined as the ratio of the yield in heavy-ion collisions (typically in bins of transverse momentum p T ) to the corresponding yield in protonproton collisions, scaled by the number of binary nucleonnucleon collisions for a specified class of heavy-ion events.
Early experimental jet-modification results at RHIC were restricted to single-hadron spectra [10][11][12][13][14], dihadron correlations [15,16], and γ-hadron correlation [17,18]. Theoretical approaches at the time were likewise restricted to the calculation of energy loss of the leading parton in a jet [19]. Already at that time there existed several different approaches that described the nuclear modification of the single-hadron spectrum within error bars [20][21][22][23]. The differences in formalism between the approaches applied at RHIC manifested in the widely different values of the jet transport coefficient q, that was extracted by these different approaches when compared to the same data [24,25]. The jet transport coefficientq is the mean squared momentum exchanged between a jet parton and the medium, per unit length traversed by the jet parton, in a direction transverse to the momentum of the jet parton: The equation above sums over heavy-ion events where jet partons encounter varying momentum exchanges with the medium. The meaning of the expression above is that in event i we consider the propagation of a jet parton a distance L i , shorter than its lifetime τ i (L i < τ i ). It may encounter several scatterings in this length yielding a net transverse momentum k i ⊥ (different in each event). Under the assumption of short Debye screening lengths, multiple scatterings can be factorized into uncorrelated single scatterings. One can then reduce the length L i , such that the hard parton will at most engender one single scattering. In this limit,q becomes a local property of the QGP.
In the equation above, Γ el is the scattering rate between a jet parton and constituents from the QGP. In the limit of single scattering, it will include 2-to-2 matrix elements that describe the scattering off a single constituent in the medium. In principle, this rate, which encompasses the energy-momentum exchange between the jet and the medium, is not known a priori. It may be nonperturbative [26][27][28] or perturbative [29] in nature, or a combination of both [9]: nonperturbative for softer exchanges and perturbative for harder exchanges. One of the central goals of the study of jet quenching is the determination of this rate or distribution, and by extension, to determine the dynamical behavior of the medium constituents, off which the hard partons scatter. Observables that strongly depend on the soft component of the jet are sensitive to physics beyond the scattering rate, e.g., the energy-momentum deposition and thermalization in the medium [30][31][32]. However, these effects should be separable by comparing a sufficiently comprehensive simulator with a wide range of data.
Subsequent measurements of hadron production at the LHC, extending the transverse momentum (p T ) reach by an order of magnitude [33,34], revealed a reduction in the nuclear modification at the LHC, even accounting for the change in the shape of the hard spectrum, suggesting a weakening of the interaction strength between the medium and the hard parton, typically indicated by the ratioq/T 3 (where T is the ambient temperature). This effective reduction in the suppression at LHC, compared to RHIC, was established by the JET Collaboration via a comparison with the nuclear modification factor for high-p T single-hadron spectra for the most-central (0-5%) events at RHIC and LHC [35]. A wide range of approaches [36][37][38][39][40] to jet modification were constrained to calculate energy loss while propagating through an identical fluid medium, constructed using a realistic equation of state and by comparison with bulk observables. All these approaches were compared to the nuclear modification data, and all comparisons required a reduction in the overall normalization ofq at LHC compared to RHIC.
Hard sector measurements at the LHC and RHIC have since been extended to cover a variety of observables related to jets over a range of collision energies and centralities, presenting an opportunity to further constrain and refine the theoretical approach to modeling jet quenching in heavy-ion collisions. To systematically compare theoretical models with this growing assortment of observables requires a comprehensive and extendable simulation framework. The JETSCAPE Collaboration has developed such a framework for a multistage event generator to study and interpret bulk medium, jet-quenching, and heavy-flavor measurements in heavyion collisions [41]. The jetscape framework has been benchmarked against p + p collisions [42] and used for Bayesian parameter estimation of the bulk properties of the QGP [43][44][45]. An earlier, simplified multistage generator was used to carry out a Bayesian evaluation of the jet transport coefficientq [46] [comparing to central and semicentral R AA at RHIC (π 0 at √ s NN = 200 GeV) and the LHC (h ± at √ s NN = 2.76 TeV and 5.02 TeV)].
In this paper, we present results from a new calculation for nuclear modification factors for inclusive jets and charged particles, calculated for a range of centralities and collision energies using the publicly available jetscape 3 series [47]. This version incorporates modifications of a hard thermal-loop (HTL)q for fixed coupling, running coupling, and with a virtuality dependent factor that effectively modulatesq to account for a reduced medium-induced emission in the high virtuality phase, due to coherence effects. For the description of the medium response in this study, the energy and momentum exchanged between jet partons and the medium are tracked using a recoil and hole scheme in both the high-virtuality and the transport stages.
Results are presented in the form of a sensitivity study in which we vary the parameters governing these new features, along with parameters for formation time, hadronization temperatures, and switching virtuality. These results will be compared to nuclear modification factors for inclusive jets and hadrons over a range of collision centralities for √ s NN = 2.76 and 5.02 TeV, measured by the ALICE, ATLAS, and CMS experiments, and for √ s NN = 200 GeV measured by the PHENIX and STAR experiments. This effort will demonstrate that our multistage framework, with an in-medium coupling strength and a transition scale between the stages, set by comparison to data at one energy and centrality along with parameters typical of energy loss in a fluid medium, such as the energy-loss start and end time, is sufficient to describe R AA data for inclusive jets and hadrons, simultaneously at all centralities and energies. This work will set the stage for a future Bayesian parameter estimation over a wider range of parameters. The remainder of the paper is organized as follows: In Sec. II, we will describe the various components of the simulation framework: the evolution of the bulk medium (Sec. II A) that provides a substrate for the propagation of jets, carried out by a combination of the highvirtuality stage (matter, described in Sec. II B), and a lower-virtuality stage (lbt, described in Sec. II C). Jet partons and partons scattered by jets are then hadronized using a variant of the Lund pythia model described in Sec. II D. The parametrized interaction of the jet partons with the medium is described in Sec. III. This involves both a description of the theory behind the jet transport coefficientq in Sec. III A, and the parametric dependencies of coherence effects, as well as a description of the recoil formalism used in Sec. III B. A multistage simulator will engender multiple parameters, these are reca-pitulated and discussed in Sec. IV along with various technical details of the simulation. Results of the simulation are presented in Sec. V. A summary of our findings is presented in Sec. VI, followed by a discussion of alternate parameter choices, different background subtraction mechanisms, and other considerations in the appendices.

II. OVERVIEW OF JETSCAPE FRAMEWORK
To explore the multiscale dynamics of jets within the framework of jetscape, we embed the space-time information of the bulk medium in the parton shower and set up an effective parton evolution within this background. Thus, the fluid dynamical simulation is run first, the space-time profile of the energy-momentum tensor, along with local fluid velocities and temperature, are stored and then recalled in a second run of the framework, which simulates hard parton formation, energy loss, and fragmentation.
The high-virtuality parton evolution is handled by the matter event generator (Sec. II B) and the lowvirtuality parton evolution is handled by lbt event generator (Sec. II C). The transfer of a parton from one energy loss model to another is performed on a partonby-parton basis. In this first attempt to simultaneously describe the nuclear modification factor for inclusive jets and leading hadrons we use a single switching virtuality Q sw . Partons with a virtuality Q > Q sw are handled by the matter generator, while those below Q sw are handled by the lbt generator. Partons that escape the medium are transferred back to the matter generator to continue showering in the vacuum. For p + p simulations, the entire parton vacuum evolution is carried out in the matter generator by switching off the medium effect.
A very straightforward means to understand the cause of the transition from high virtuality to a lower virtuality gen-erator is given in Ref. [48], where one can see the medium modified contribution to the radiation spectrum steadily grow and surpass the vacuum contribution as the virtuality is reduced. As a parton propagates in a medium it undergoes multiple scattering, which iteratively adds virtuality up to the medium generated scale where τ is the lifetime of the parton andq(ξ) is the local value of the transport coefficient. Besides the energy E of the parton, emissions (or splits) depend on the virtuality of the parton Q 2 : they have a transverse momentum l 2 ⊥ y(1 − y)Q 2 (where y is the forward momentum fraction of one of the daughter partons in a split), and occur in a time, If Q is large, τ is small, and the accumulated virtuality from scattering is small compared to the actual virtuality of the parton Q Med Q. In this regime, scattering from the medium is at most a perturbative correction to the process of vacuum emission. This portion is simulated with the matter generator.
As the virtuality Q continues to drop with each successive emission, the formation time of splits increases and the medium generated virtuality becomes comparable to the virtuality of the parton. The time of onset of this stage can be estimated by setting Q = Q Med : By this point, the medium modified kernel far exceeds the vacuum emission kernel [48]. Thus, generators such as martini and lbt which simulate this phase typically ignore the vacuum contribution. Since the T and E dependent transition scale (Q Med ), between the high and low virtuality regimes, is only known approximately, we replace it with a parameter Q sw , which is then tuned in comparison with data.
As the shower proceeds through the dynamical medium, the ambient temperature will eventually begin to drop and this causesq and the medium generated scale to also drop. If this drop is sudden, e.g., in the case of the jet passing the parton-hadron transition surface, the jet parton may once again be in a regime where its virtuality is much larger than that generated by multiple scatterings in the medium. To mimic this effect, partons that cross the phase transition surface with Q 2 > 1 GeV 2 are fed back to the matter generator.
In the following subsections, we describe the basic features of the simulation of the bulk medium. This is followed by a discussion of the salient features of the matter and lbt event generators. Both of these generators have been extensively discussed in the literature. Hence, the descriptions that follow will be terse. The last subsection (II D) concludes with a discussion of the hard hadronization module, which fragments the jet partons, recoil partons, and hole partons into regular and hole hadrons respectively.
The simulations described in this paper have been carried out using the publicly available jetscape 3.0 [47] version of the event generator framework. Unlike prior versions, the current framework can simultaneously account for light and heavy flavor energy loss. It also contains modules that can nonperturbatively deal with the energy and momentum deposited from a jet to the medium. There are three separate modules for the hadronization of the hard sector and one for the soft sector. In the main body of this paper, we use the simple recoil description of the medium response (described in the next section), with other approaches discussed in the appendices.

A. Dynamics of the soft sector in A+A collisions
The QGP medium evolution is modeled by relativistic viscous hydrodynamics. We assume longitudinal boost invariance for heavy-ion collisions at the top RHIC and LHC energies. Event-by-event simulations are set up using the trento initial conditions [49] followed by a (2+1)-dimensional [(2+1)D] free-streaming preequilibrium evolution up to a proper time of τ hydro (=1.2 fm/c at LHC, and 0.5 fm/c at top RHIC energy [50]). After matching the system's energy-momentum tensor between the preequilibrium and fluid phase, the QGP medium evolution is described by the vishnu (2+1)D hydrodynamics [51,52]. As the system evolves to dilute densities, individual fluid cells are converted to hadrons via the Cooper-Frye prescription [53,54]. This conversion is performed on an isothermal hypersurface at T sw = 151 MeV [55]. The produced hadrons are transferred to a hadronic transport model urqmd for scatterings and decay [56,57]. We point out that while the default jetscape settings use the music generator [58] for the fluid dynamical simulation and the smash generator [59] for the hadronic cascade, the framework is devised to work with a variety of interchangeable generators (detailed comparisons between vishnu and music and between urqmd and smash are presented in Appendix H of Ref. [43]).
The causal relativistic hydrodynamical equation of motion is given by the second-order Israel-Stewart theory [60,61]. In addition to conservation of energy and momentum, second-order hydrodynamical equations also include relaxation-type equations for six independent viscous degrees of freedom, namely five in the shear stress tensor π µν with the remaining being the bulk viscous pressure Π. Energy-momentum conservation is expressed as: with the energy-momentum tensor where is the energy density, u µ is the flow four velocity, P is the thermodynamic pressure related to by the lattice QCD equation of state P ( ) at vanishing net baryon density [54,62]. We define the spatial projection tensor is the metric tensor. The dissipative degrees of freedom are evolved according to τ ππ µν + π µν = 2ησ µν − δ ππ π µν θ + λ πΠ Πσ µν − τ ππ π µ α σ ν α + φ 7 π µ α π ν α ,  (7) and (8).
The second-order transport coefficients present in Eqs. (7) and (8) were computed in Refs. [63,64] assuming a single-component gas of constituent particles in the limit m/T 1, where m is their mass and T is the temperature. Table I summarizes these transport coefficients, where c 2 s = ∂P/∂ is the speed of sound squared. The temperature dependent specific shear viscosity η/s(T ) -where s is the entropy density -and the specific bulk viscosity ζ/s(T ) are taken from a recent Bayesian model-to-data comparison [65].
The time evolution of hydrodynamic fields, such as local energy density, temperature, and flow velocity, are stored on disk event-by-event for the second stage for jet showering. During the preequilibrium stage τ < τ hydro , the Landau matching procedure [50] is performed at every time step to compute the local energy density and flow velocity from the system's energy-momentum tensor. Then the local temperature is estimated by the ideal massless quark-gluon gas equation of state with N f = 3 and N c = 3.
Individual jet partons start to interact with the QCD medium at a longitudinal proper time τ 0 = 0.6 fm/c. When τ 0 is smaller than the proper time of the fluid dynamical simulation τ hydro , which is the case for the LHC simulation in this work, the energy loss calculation is performed using the local temperature and flow velocity obtained via the Landau matching in the preequilibrium phase. We will vary this jet energy loss starting time τ 0 to quantify its effects on the R AA observables in Fig. 9 below.
As the jet develops its shower inside a dynamically evolving QCD medium,q is sampled according to the local temperature information for the jet partons boosted to the local rest frame of the fluid cell. We stop the jetmedium interactions at the energy loss termination temperature T c = 160 MeV, below which the partons propagate only with vacuum emissions in matter followed by fragmentation. In the jetscape two-stage approach, wwe note that neither the jet partons nor the fragmented hadrons interact with hadrons from the soft sector in the hadronic phase. To quantify the uncertainty from this approximation on R AA observables, we will vary the value of the energy loss termination temperature between 150 MeV and 170 MeV (see Fig. 10 below). We remind the reader that the transition temperature at which the fluid simulation undergoes Cooper-Frye particlization is T = 151 MeV. We do not take into account the parton energy loss in the phase space between T = 150 − 151 MeV after the entire fireball is frozen out.

B. MATTER event generator
The Modular All Twist Transverse-scattering Elasticdrag and Radiation (matter) is a higher-twist formalism-based event generator that simulates the parton modification both in a vacuum and within a medium. In this instance, a parton in matter will engender an arbitrary number of emissions, where stimulated emissions are calculated in the one-rescattering or twist-4 approximation. It is primarily applicable to the high-virtuality, the high-energy epoch of the parton shower, where the virtuality of the parton Q 2 √q E. In this phase, the medium-modified radiative processes are not dominant, and the successive emissions from the parton are ordered in virtuality.
In matter, the distribution of the medium-modified radiated gluon from a single scattering with the medium is given as where α s is the strong coupling constant and P a (y, Q 2 ) = P vac a (y) 1 + Here a = (g, q,q) is the parent parton species, P vac a (y) is the standard Altarelli-Parisi vacuum splitting function, y is the momentum fraction carried away by the emitted daughter parton, p + = (p 0 + p 3 )/ √ 2 is the lightcone momentum for the parent parton traveling along z-direction, and τ + = 2p + /Q 2 is the formation time of the radiated gluon. The parent parton started at ξ + o and completes the split at ξ + which lies between ξ + o and ξ + o + τ + . The quantity K a (ξ + , ξ + o , y, p + , Q 2 ) is singleemission-single-scattering kernel given as [9,[66][67][68] where, In above equation, δ a,q and δ a,q are Kronecker delta functions, while χ a = (δ a,q + δ a,q )y 2 m 2 a /(y(1 − y)Q 2 − y 2 m 2 a ) with m a being the mass of the parent parton a. The jet transport coefficientq a measures the average squared transverse momentum broadening per unit length of the medium. If the value ofq is zero, the distribution of the emitted gluon in Eq. (9) reduces to a vacuumlike distribution. The factor f in Eq. (11) accounts for the transverse size of the jet parton and is discussed in detail in the next section.
The virtuality-ordered shower is generated based on the Sudakov formalism where we solve the in-medium Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation using Monte Carlo (MC) method. The shower is initiated by a single hard parton produced at a space-time location r carrying the forward light-cone momentum p + = (p 0 +n · p)/ √ 2, where p ≡ (p 0 , p) is 4-vector momentum of the parton andn = p/| p| denotes the direction of the jet. Then, given a maximum allowed virtuality t max and minimum virtuality t min , one determines the virtuality of the parent parton a by sampling the Sudakov form factor given as where the Sudakov form factor represents the probability for a parton to transition from virtuality t max to t via "unresolvable" emissions. Here the path length integration in the in-medium splitting function of Eq. (10) is performed along the direction of the jetn. The virtuality of the parent parton is determined by generating a random number R from a uniform distribution between 0 and 1. If S a (t max , t min ) > R, then the parton is assigned t = t min and propagates freely to next energy loss routine. However, if S a (t max , t min ) < R, then the virtuality t is obtained by solving S a (t max , t) = R, and the splitting may happen. With the determined virtuality t, the splitting functionP a (y, t) is sampled to determine the momentum fraction y shared by the two daughter. This sets the daughters' momenta to be (1 − y)p + (daughter 1) and yp + (daughter 2). Now, the daughter parton's virtuality is determined using the Sudakov factor with (1 − y) 2 t as maximum virtuality for first daughter and y 2 t for the second daughter. Once the actual virtuality t 1 and t 2 of the daughters are known, their transverse momentum with respect to the parent is calculated from Next, the l − component is determined using energymomentum relation: l 2 1 = t 1 and l 2 2 = t 2 . Finally, the location ( r +nξ) of medium-induced splitting is determined by sampling a Gaussian distribution given as where τ + f is the mean life time given as τ + f = 2p + /t. The above procedure is repeated iteratively for each shower initiating parton until the virtuality reaches or goes below a switching scale Q sw . At this scale, the parton is transferred to the lbt event generator, discussed in the next subsection. The minimum virtuality t min in the Sudakov sampling is always fixed to 1 GeV 2 . If the parton exits the medium and the lower virtuality generator (in this case, lbt), it will return to matter and continue vacuumlike showering until the virtuality reaches or below 1 GeV 2 .

C. LBT event generator
The linear Boltzmann transport (lbt) model is a parton transport generator that is used to simulate the propagation and interaction of both the jet shower and recoil partons with elastic and inelastic collisions in the QGP medium [69,70]. It is primarily applicable to the lowvirtuality, high-energy epoch of the parton shower. In this phase, multiple scattering-induced emission is the dominant mechanism of parton energy loss. Vacuumlike emission is ignored in this stage.
The phase space distribution f a (x a , p a ) of the parton with momentum p µ a (E a , p a ) is determined by solving the linear Boltzmann equation: where C el a and C inel a are the collision integrals for elastic and inelastic scatterings. In this formalism, the total scattering probability is expressed as P tot a = P el a +P inel a − P el a · P inel a , where P el a and P inel a are elastic and inelastic scattering probability, respectively. These probabilities are sampled using Monte Carlo techniques to determine the type of scattering channel. The probability for a parton to undergo elastic scattering (2 → 2) in the given time step ∆t is given by P el a = Γ el a ∆t, where the elastic rate is given as where g b represents spin-color degeneracy, f b is the ther- 3 ], and |M ab→cd | 2 is the amplitude square of the process a + b → c + d. In the interaction kernel, is imposed to regularize the divergence in the matrix element |M ab→cd | 2 arising from small angle, u, t → 0. The Debye screening mass is given as where N f is the active quark flavors in the QGP.
Currently, the lbt model is set up to simulate inelastic channels via single scattering (2 → 2 + n) causing a multiple gluon emission (n) processes, where the mediuminduced gluons are independent. At each time step ∆t, the number of gluon emissions is sampled using the Poisson distribution, where the mean number of gluon N a g = Γ inel a ∆t. Thus, the probability for total inelastic scattering process to occur is given as P inel a = 1 − P (0) = 1 − e −<N a g > . The inelastic rate for medium-induced gluon radiation is given by where δ a g is the correction term for double counting of the process g → gg.
The medium-induced gluon spectrum in Eq. (20) is derived using the higher-twist energy loss formalism and given as where y is the momentum fraction, l ⊥ is the transverse momentum of the radiated gluon, t i is the initial time of the parent parton, and τ f is the formation time of the radiated gluon. Based on the probabilities P tot a , P el a , and P inel a , we first determine whether scattering occurs and whether the scattering is elastic or inelastic. Once these are known, the differential spectra given in Eq. (17) and Eq. (21) are sampled to determine the energies and momenta of the outgoing partons. The lbt model has one free parameter, the jet-medium coupling constant α s that controls both elastic and inelastic parton energy loss.
In the low virtuality transport stage of a heavy-ion collision, one expects multiple scattering per emission. The lbt generator, however, only includes one re-scattering (or two scatterings) per emission. In contrast to this, the martini generator [39] includes an arbitrary number of scatterings per emission. In Ref. [71] we studied multistage energy loss in a static medium using a combination of matter and lbt, as well as matter and martini. In these studies it became clear that for static media with lengths that lie in the range 2 L 8 fm, there is a negligible difference in the final results if martini is replaced by lbt. Further studies in fluid dynamical media were reported in Ref. [72], where a combination of mat-ter+lbt was compared with matter+martini for jets and leading hadrons at 2.76 TeV for two different centralities. The differences between these two implementations were less than 5%. Due to the order of magnitude longer compute time required by the martini generator, we carry out this first study of the nuclear modification of jet and leading hadrons in jetscape using the mat-ter+lbt generator. D. String hadronization jetscape 3.0 has three different hadronization modules: colored hadronization [41,42], colorless hadronization [41,42], and hybrid hadronization [73]. Both colored and colorless hadronization use the default Lund string fragmentation from pythia 8. The hybrid hadronization model is a combination of Lund string fragmentation and recombination. Since colorless hadronization is the only hadronization used in this study, we provide a brief explanation of colorless hadronization in this subsection.
Even though colorless hadronization uses string hadronization through pythia, it removes all color information prior to the hadronization process. All the partons generated from the collection of shower-initiating partons-the radiated partons, the recoils, etc.-are collected in one list. The module then reconstructs one to several strings depending on the number of quarks and antiquarks in the combined set of showers in that event.
Depending on whether the total number of quarks or antiquarks are even or odd, extra quarks or antiquarks are added at beam rapidities to make the net quark number of all the showers be zero. Gluons are assigned to a string with a quark and an antiquark at either end. Once all partons have been assigned to strings, color tags are generated for all partons, such that each string remains a color singlet. These are then hadronized.
The collection of all hole partons, which are particles introduced for the description of the medium response explained in the next section, is then combined and treated similarly to generate hole hadrons, which can then be subtracted from jet cones within which they appear. As will be discussed further in Appendix C, forming strings out of a large number of partons, especially where there are a lot of soft partons, may run into issues with pythia string-breaking algorithms. In cases where two partons have a |δ p| 4Λ QCD , the p z component of the parton at larger rapidity is increased until the bound is reached.

III. JET TRANSPORT COEFFICIENT AND MEDIUM RESPONSE
In the previous section, we divided the history of a jet shower into the production, the initial propagation of high virtuality partons (in matter), the evolution of lower virtuality partons (in lbt), and their fragmentation into hadrons. In both cases of matter and lbt, the scattering in the medium identifies and correlates medium partons which were nudged forward in the direction of the jet. Along with the partons generated by showering, the entire collection of jet-correlated partons now includes the incoming partons from the medium, referred to as holes and their scattered versions called recoils.
Given the differences in virtuality, there is some difference between stimulated emission in the matter phase versus the lbt phase. Due to the small transverse size of the emission antenna in matter, the effect of scattering in the medium is diminished. This leads to a reduction of the effective value of the transport coefficientq and the distribution of recoils and holes. These details, along with the method of subtraction of the holes, are discussed in the following subsections.
A. Jet transport coefficient and coherence effects at high virtuality The dominant mechanism of jet energy loss in the plasma is due to bremsstrahlung radiation, triggered by soft gluon exchanges with the medium. The jet transport coefficientq defined in Eq. (1) effectively characterizes the momentum broadening of a single parton due to the in-medium scattering, leading to induced emission.
In the last decade, several attempts have been made to computeq using first principles and model-to-data comparisons. In the limit of high temperature and weakcoupling approximation, the hard-thermal-loop (HTL) based calculation yields aq given by [29]: In the equation above, C a is the representation specific Casimir, the constant C arises from different choices of the upper limit of the k ⊥ integral in Eq. (1), which leads to the logarithmic term. The weak coupling calculation has been extended to next-to-leading order (NLO) by the author of Ref. [74], where a large additive contribution to the above equation was found. Quantitative determinations ofq based on lattice gauge theory have also been formulated [75][76][77][78]. Recently,q has been computed for (2+1)-flavor QCD plasma on 4D lattices [28]. These results are similar to and somewhat lower than those extracted from phenomenology. The calculations of transport coefficient q based on N =4 Super-Yang-Mills theory have also been carried out, however, these yield an order of magnitude higher results compared to phenomenology-based extractions [79][80][81][82].
A first systematic extraction ofq based on phenomenology was carried out by the JET Collaboration [35]. Extractions were based on a comparison of jet quenching model calculations to the experimental measurement of only the hadron R AA , in only the mostcentral collisions at RHIC and LHC energies. These were performed independently, for five different parton energy loss approaches: GLV-CUJET [20,83], HT-M [36], HT-BW [38], martini [39], and McGill-AMY [84]. These jet energy loss calculations were run on identical (2+1)D viscous hydrodynamical media [51,[85][86][87]. The striking result of this work was that the interaction strengthq/T 3 for the QGP at RHIC energy appeared to be higher compared to that at LHC energy. In other words, one has to re-adjust the fit parameter inq when comparing with data from LHC collision energies versus RHIC energies.
Even within the work of the JET Collaboration, it was clear that different energy loss models had made slightly different assumptions regarding the temperature and parton energy dependence ofq. While some used a variant of Eq. (22), other models assumed a scaling with some density profiles of the medium, e.g., entropy density. As a result of the JET effort, it becomes necessary to generalize the functional dependence ofq on T and E and use a more data-driven approach to isolate this dependence. This was carried out recently by our collaboration in Ref. [46]. In that effort, the formula forq was generalized to allow for an additive dependence on logarithms of energy along with the generic HTL form [Eq. (22)]. Comparisons were carried out with the hadron R AA at two centralities over three different collision energies ( √ s NN = 0.2 TeV, 2.76 TeV and 5 TeV). The outcome of the effort in Ref. [46] was that an additive dependence ofq on logarithms of energy and temperature did allow for a simultaneous description of the hadron R AA at RHIC and LHC, without the need for refitting. However, there was no noticeable improvement in the fit from a multistage versus a single-stage model. Up to this point, almost all attempts to extract the transport coefficientq have at most assumed dependence on E and T , which are the only possibilities for an onshell hard parton propagating through the plasma. This may not be the case for a highly virtual hard parton. Recently, several authors have argued that medium-induced radiation should depend on the resolution scale of the medium [6,8,88,89]. Early in the history of a jet propagating in a medium, partons are very virtual, and splits engender large transverse-momentum scales. As a result, the transverse sizes of the QCD antennae are very small, compared to the resolution scale in the medium, Q 2 med ≈qτ (τ is the formation length of the parton). The inability of partons emanating from the medium to resolve such small antennae are often referred to as "coherence effects" in jet propagation.
While the authors of Refs. [6,8,88,89] advocated the use of a sudden approximation-neglecting any mediuminduced emission in the high-virtuality phase-the authors of Ref. [9] derived a more gradual reduction of medium-induced emission as a function of the virtuality of the parton. In Ref. [9], the reduction in mediuminduced emission is cast as a reduction in the effective value ofq as a function of the parton virtuality Q 2 . The latter formalism will be used in the current effort as it is a better approximation to the reduced energy loss at high virtuality. Also, with minor modifications, as outlined below, we will be able to study the onset of coherence effects at high virtuality. The reduction in the effectiveq will only take place in the matter phase of the simulation. The lbt phase will include the fullq with running coupling, described below. As part of our analysis, we will vary the transition scale between matter and lbt as well.
For a high energy on-shell parton (E m D ), the hardthermal-loop (HTL) calculation yields the following form of transport coefficientq, given as [69], In the above calculation, Bose and Fermi distributions for plasma constituents are invoked. Comparing with Eq. (22), we have set C = ln (2) and N c = N f = 3 in Eq. (18) for the Debye mass, and ζ(3) = 1.202 is Apéry's constant. In the above formula, only the α s within the logarithm is designated as α fix s , while such a qualification is suppressed for the factors of α s outside the logarithm. We will present results with these factors of α s either fixed at the Debye scale or running with the scale of the exchanged k ⊥ in an improved perturbation theory calculation of Eq. (1).
Incorporating the reduction in the medium induced emission described in Ref. [9], we propose a virtuality (Q 2 ) dependent modulation factor f (Q 2 ) as the f in Eq. (11). One could incorporate the factor f withinq, yielding a virtuality dependentq(T, E, This f factor effectively decreasesq as virtuality increases, in the high virtuality stage, and reduces to the HTL result (f = 1) in the low virtuality stage of the parton shower. In this effort, we shall explore the following three formulations and carry out a sensitivity study of the in-medium coupling constant (α fix s ) and switching scale (Q sw ) parameters in the simultaneous description of inclusive jet and charged-particle R AA : (1) Type 1: HTLq with fixed coupling (f = 1 applies for any Q 2 ), where m 2 D = 6πT 2 α fix s is the Debye mass for N f = 3 flavors.
(2) Type 2: HTLq with running coupling (f = 1 applies for any Q 2 ), where Q 2 max = 2ET (3) Type 3: HTLq with a virtuality (Q 2 ) dependence factor where E is the energy of the hard parton, T is the local temperature of the medium, and Q 2 is the running virtuality of the hard parton. The form for f (Q 2 ) is a simplified form of the formula derived in Ref. [9]. We point out here that theq is the same between type-2 and type-3. The extra factor of f (Q 2 ) in type-3 diminishes the interaction between the jet and the medium due to coherence effects [6,7,9,88]. Theq in this case is the same as in type-2.  To visualize the new functional dependence (Q 2 ) proposed above, we plotq/T 3 for Type 2 andqf (Q 2 )/T 3 for Type 3 in Fig. 1. In the case of Type 3 (solid line), the virtuality dependence in the matter phase is shown, exhibiting a reduction as Q 2 increases. In addition to this, qf (Q 2 ) of Type 3 reduces to traditional HTLq in the lbt phase where the parton's virtuality becomes Q 2 ≤ Q 2 sw . In Sec. V, we shall show that the experimental data for the charged-hadron R AA strongly favors the one with Q 2 dependence.

B. Medium response in a weakly-coupled approach by recoils
Jets exchange energy and momentum with the soft medium, during which they excite medium constituents. Some of these excited partons are clustered with the jet and modify the structure of reconstructed jets. In the current effort, the medium responses are incorporated using perturbation theory (nonperturbative methods are further described in Appendix C 3), with nonperturbative effects incorporated solely during hadronization. In this medium response process, some portion of jet energy and momentum are transported through further successive interactions among the medium constituents far be-yond the jet cone. Their contribution appears as jetcorrelated particles in the final state and affects the jet medium correlations (these will be described in detail in a future effort).
In this study, the medium response is described as the propagation of recoil partons and their successive interactions in the jetscape framework [90][91][92][93][94][95][96]. In both matter and lbt, the energy-momentum transfer between jets and the medium is executed via 2-to-2 partonic scatterings. For each scattering, a medium parton is sampled from a thermal bath of three-flavor ideal QGP gas. The sampled thermal parton after being scattered by the jet parton, referred to as a recoil parton, is assumed to be on shell. Since it is not virtual, its subsequent in-medium evolution will be carried out by lbt, assuming its weak coupling with the QGP. Note, in matter energy loss phase, the elastic scattering probability is also reduced by f (Q 2 ) to be consistent with the modified kernel for transverse momentum broadening.
The jet shower partons, including these recoil partons, are hadronized together via the Colorless string hadronization routine. On the other hand, the recoil parton leaves an energy-momentum deficit (hole) in the medium. We also keep track of the hole partons and subtract their contribution to ensure energy-momentum conservation. The hole partons are assumed to freestream in the medium and are hadronized separately from other regular jet-shower partons. The subtraction of the hole contribution for the final reconstructed jet momentum is performed as Here p µ shower is the four-momentum of the jet reconstructed from all particles from the hadronization of jet showering partons including the recoil contribution by the anti-k T algorithm [97] implemented in the fastjet package [98,99] with a jet cone size R. In the second term on the right hand side, the sum of four momenta p µ i is taken over particles hadronized from hole partons inside the jet cone This recoil prescription gives a reasonable description of the medium response as long as jet shower partons have sufficiently large energy and are far from the thermalization, where their mean free paths are long enough to apply the kinetic theory. This recoil approximation breaks down when the showering partons' energy approaches the typical scale for the thermalized medium constituents [2,100,101]. To extend this region of applicability, one needs to incorporate the hydrodynamic description for the soft modes of jets as presented in Refs. [96,[102][103][104][105][106][107][108][109][110][111][112][113]. In this study, we do not include this hydrodynamic description for the medium response to jets, which exacts a huge computational cost for the systematic studies of jets presented here. Although they are essential for a more precise description of jet-correlated particle distribution, the recoil prescription is still a good approximation for the estimation of jet transverse momentum with typical jet cone sizes R ≈ 0.4 [96]. We leave systematic jet-quenching studies including more detailed modeling of medium response for future work (however, see Appendix C 3 for a discussion of alternative approaches).

IV. SIMULATION WITH MULTISTAGE ENERGY LOSS APPROACH IN JETSCAPE
The jetscape framework is a general-purpose numerical framework to simulate the complete evolution of heavy-ion collisions. Currently, it provides several alternative implementations of physics models to simulate both the QGP evolution and the jet-medium interactions at different epochs of the parton shower. In this paper, we carry out the calculation for both p + p and A+A collisions using the publicly available jetscape 3.0 [47]. In this section, we discuss the choice of modules and parameters that will be explored in the current effort.
The p + p baseline is simulated using the jetscape PP19 tune [42] 1 . This tune generates the hard scattering in a p + p collision using the pythia module where initial state radiation (ISR) and multiparton interaction (MPI) switches are enabled, but the final state radiation (FSR) is turned off. Then, the produced partons from pythia hard scattering are transferred to the matter energy loss module for final-state radiation. As matter embeds the space-time structure of the bulk medium in the parton shower and is capable of performing both the vacuum (q = 0) and in-medium (q = 0) energy loss, it is the desired choice for final state radiation module to ensure the consistency between p + p and A+A collisions.
The soft products in A+A collisions are generated using fluctuating trento initial conditions [49], evolved hydrodynamically using the (2+1)D vishnu code package [52] (the underlying physics setup is discussed in Sec. II A). Events at LHC were simulated using the maximum a posteriori (MAP) parameters obtained in Ref. [65] using Bayesian calibration, while hand-tuned parameters were used for top RHIC energy. To obtain events within different centrality classes, we made those centrality tables for each beam energy first and then simulated these events by inputting entropy ranges corresponding to these centrality bins for trento. To make the centrality tables, we first ran 10 5 minimum-bias trento, events for each energy using the MAP parameters and then sorted them by entropy for binning. Virtuality to switch parton energy loss from matter phase to lbt phase 2 GeV τ0 Starting longitudinal proper time for in-medium jet energy loss 0.6 fm/c Tc Temperature below which the jet-medium interaction is turned off 160 MeV Once the bulk simulations are complete, the space-time distribution of the energy-momentum tensor is saved. This allows simulations of in-medium jet evolution for several hard scattering events on a single bulk event. 2 In the next step, the initial hard scatterings are generated using pythia. The trento initial-state module generates the initial energy density profile for the following free-streaming evolution. It also provides the binary collision profile to sample the transverse positions for the hard scatterings. To incorporate multiscale dynamics of the parton energy loss, the jet evolution is carried out as follows. First, all the partons from pythia hard scatterings are transferred to the matter module, which assigns an initial virtuality limited by the maximum Q 2 init ≡ p 2 T /4, based on the parton's initial transverse momentum and t max in Eq. (13) for the first Sudakov sampling is set to Q 2 init [42]. In both p + p and A+A collisions, the partons generated during initial state pythia hard scattering with p T < 2 GeV/c are discarded [42]. This treatment is because the matter module can not create a DGLAP based parton shower if the maximum initial virtuality of the parton Q 2 init < 1 GeV 2 . Second, the matter module generates virtuality-ordered showers for individual partons from pythia. It includes both vacuumlike and medium-induced radiation. As parton's virtuality drops below a switching scale Q sw , it is switched to propagate with the lbt module. Once the partons fly outside the medium where the local temperature is below the energy loss termination temperature T c , they are transferred back to matter for vacuumlike radiation if the partons have virtuality larger than the minimum virtuality Q 2 0 = 1 GeV 2 . After all the partons are outside the QGP medium and have virtuality Q 2 ≤ Q 2 0 , they are hadronized by the Colorless string fragmentation (see Sec. II D for details).
The parameter set for this multistage jet evolution model in A+A collisions is summarized in Tab. II. The exact functional form of the jet transport coefficientq of the first parameter is unknown from theory. Given prior efforts in Refs. [35,46], it is clear thatq depends on more than the ambient temperature T . The choice of an additive dependence on T and parton energy E in Ref. [46] showed only modest improvement on describing the charged hadron R AA measurements. Based on the work of Ref. [9], we have invoked a multiplicative dependence on the virtuality of a parton Q via the modulation factor f (Q 2 ). In Appendix A, we will show conclusively that the multiplicative dependence on the parton virtuality Q is essential for a simultaneous description of the modification of the inclusive jet and charged particle spectra.
The next parameter is the in-medium coupling α fix s at the Debye scale. It also controls the overall normalization ofq and the strength of recoil scattering of jet partons in the medium. The value of this nonperturbative parameter is unknown and needs to be calibrated by the R AA data. The switching scale Q sw controls the relative phase spaces between the matter and lbt parton shower. The value of Q sw should be close to the medium scale Q sw ≈qτ with τ being the parton formation time.
In this work, we choose a constant switching scale Q sw between the matter and the lbt modules. In principle, Q sw is a dynamical quantity depending on localq, energy and virtuality of the parton. Thus, the constant value of Q sw introduced in our current work should be interpreted as the typical transition scale averaged over those quantities of all partons. We leave studying the effects of a dynamical switching scale Q sw = Cqτ [71] on observables for future work. The last two parameters τ 0 and T c specify the start time and stop temperature conditions for jet-medium interactions. We choose τ 0 to be smaller than the starting time of hydrodynamics τ hydro to take into account parton energy loss in the preequilibrium stage. The R AA 's sensitivity to the parameter τ 0 will be investigate in Fig. 9. We stop the jet-medium interactions below the energy loss termination temperature T c . This approximation neglects energy loss in the dilute hadronic phase. The uncertainty from this approximation will be quantified by varying the termination temperature T c in Fig. 10.

V. RESULTS
In this work, we cover three collision energies: √ s NN = 5.02 TeV, 2.76 TeV, and 200 GeV, and show comparisons with selected experimental data available from ALICE, ATLAS, CMS, PHENIX, and STAR Collaborations. Table III shows the list of references for the data used in the comparisons. We will first tune the model parameters in Table II to achieve a simultaneous description of the inclusive jet and charged-particle R AA in 0-10% Pb+Pb collisions at √ s NN = 5.02 TeV. With the optimized model parameters, we will present jetscape postdictions for these observables in semi-peripheral centralities of Pb+Pb collisions at √ s NN = 5.02 TeV and in central Pb+Pb and Au+Au collisions at 2.76 TeV and 200 GeV, respectively. We will also present predictions for inclusive jets at 200 GeV.
Once the R AA for the most central collisions at 2.76 TeV has been obtained (for inclusive jets or leading hadrons) and compares favorably with experimental data, the variation with centrality is easily obtained. This is discussed in Refs. [46,72], and will not be repeated here. Also semicentral R AA for jets at √ s NN = 200 GeV are not currently available and as a result we do not present centrality dependence of the R AA here. Nuclear modification data for neutral pions in semicentral events is available. However, once the R AA for most central collisions compares favorably with the experimental data, the centrality dependence is easily obtained, almost independent of the choice of multistage model used. Such comparisons were already presented in Ref. [46].

A. The p + p baseline
To quantify the medium effects in inclusive jet and charged-particle spectra as the nuclear modification factor R AA , calculations for p + p collisions are necessary to obtain the baseline for A+A collisions. A systematic study of inclusive jet, intra jet, and charged-particle observables has been extensively carried out using the jetscape PP19 tune and presented in Ref. [42]. Thus, we shall only present a selection of plots for inclusive jets and charged-particle yield for p + p collision energies √ s =5.02 TeV, 2.76 TeV, and 200 GeV. The presented results are further restricted to those based on the jetscape Colorless hadronization, as it is the hadronization module employed in the A+A sector. The uncertainty in the final observable spectra from the two different prescriptions, colored and colorless hadronization, are roughly of the order of 10%, and we refer readers to Ref. [42] for more details. Figure 2 shows our p + p simulation results for inclusive jet spectra around midrapidity at √ s = 5.02 TeV, compared to experimental data from ATLAS [114] and ALICE [115]. The jets are reconstructed with jet cone size R = 0.4 using the anti-k T algorithm in fastjet. The comparison with the results from pythia 8 with default parameters are shown in Fig. 3.
Results from jetscape PP19 describe the experimental data very well for jets with large transverse momentum or small rapidity. Our results for |y jet | < 0.3 are compatible within uncertainties with data from ATLAS throughout the entire p jet T range (up to 1 TeV). The results for the wider rapidity range |y jet | < 2.8 agree with ATLAS data in the region with p jet T 300 GeV but deviate at low-p jet T . In comparison to the ALICE measurements for |η jet | < 0.3, jetscape PP19 overestimates the data and tends to be similar to pythia 8. We also show the ratios of differential cross section for inclusive jet at midrapidity in p + p collisions at √ s = 2.76 TeV and 200 GeV in Fig. 4. Both jetscape and pythia 8 tend to overestimate the measured jet cross-section at those collision energies.
In Fig. 5, we show the ratio of charged particle crosssections at mid-rapidity in p + p collisions at various collision energies. Throughout all the collision energies and p T range, jetscape results match with pythia 8 results within the statistical errors and describe the experimental data.

B. Comparison between differentq · f formulations
In this section, we perform a multistage based jet energy loss calculation and explore three different forms ofq · f presented in Eqs. (24)- (26). As described in Sec. III A, in the Type 3 formulation, the reduction in the medium induced emission due to coherence effects [7,9] can be incorporated with a scale-dependent reduction  [122] factor f (Q 2 ) included with the HTL expression forq. We will demonstrate that this form is essential for the simultaneous description of inclusive jet and charged-particle R AA . In Fig. 6, we present the nuclear modification factor for inclusive jets and charged-particles at most-central (0-10%) Pb+Pb collisions at √ s NN = 5.02 TeV. The calculations are carried out using the multistage jet quenching model (matter+lbt) for three differentq · f formulations and compared with the experimental data from ATLAS [114], CMS [120], and ALICE [115] for inclusive jets R AA and from CMS [118] for charged-particles R AA . In the above calculation, the switching virtuality parameter is set to Q sw = 2 GeV. The best fit to inclusive jet R AA and charged-particle R AA yields α fix s = 0.3 for Type 3 and α fix s = 0.25 for Type 1 and Type 2. All three formulations produce similar results for inclusive jet suppression given one needs to readjust the fixed coupling constant α fix s when going from Type 1/Type 2 to Type 3. Results for Type 1/Type 2 for other values of α fix s = 0.3 with Q sw = 2 GeV are shown in the Appendices. (see Fig. 16 in Appendix A for Type 1 and Fig. 18 in Appendix B for Type 2). We note the inclusive jet R AA presented in the top panel of Fig. 6 shows good agreement with the ATLAS data for jet p T < 250 GeV. However, for jet p T > 300 GeV, inclusive jet R AA deviates ( 10%) from the ATLAS data and strongly favors the CMS inclusive jet data. In addition to this, the inclusive jet R AA for eachq formulation are in agreement with the ALICE experimental data for all measured jet p T .
In contrast to the reconstructed jet, the suppression of high-p T charged-particle yields exhibits strong sensitivity to the virtuality dependence of coherence effects. The results from matter+lbt with a virtuality dependent modulation factor f (Q 2 ) (Type 3) reproduces the experimental data from the CMS Collaboration quite well throughout the entire p T range. The other two formulations without the Q 2 dependence (Type 1 and 2) agree with the data and Type 3 at low p T , but show significant over suppression at high p T .
It should be noted that this behavior of Type 1 and 2 formulations cannot be improved by changing the parameters (see Appendices A and B). This strongly indicates that the virtuality dependence in energy loss (the gradual onset of coherence effects) is essential to describe the p T dependence of the R AA of charged particles.
The insensitivity in jet R AA and sensitivity in charged particle R AA to the virtuality dependence in parton energy loss can be interpreted as follows. The charged particle distribution is dominated by the contribution of hadrons from leading partons in jet shower. Leading partons with large p T are more likely to have large virtuality during the early stage of the energy loss. When partons undergo energy loss based on the virtuality dependent formulation in Eq. (26), the strength of interaction with the medium is diminished for the high-p T leading partons due to their large virtuality. This results in the weaker suppression of larger-p T charged particles.
On the other hand, the jet energy loss is mainly brought by the medium effects on partons at larger angles, which causes energy outflow from the jet cone. Results for yjet < 0.3, compared to ATLAS data [114]. Middle panel: Results for yjet < 2.8, compared to ATLAS data [114]. Bottom panel: Results for ηjet < 0.3 with p lead, ch T > 7 GeV, compared to ALICE data [115]. TeV, compared to CMS data [116]. Bottom panel: Jets with R = 0.6 and |ηjet| < 1.0 at √ s = 200 GeV, compared with the STAR data [117]. Additionally, we show STAR measurements for inclusive jets based on Midpoint-cone algorithm with R = 0.4 [123].
Through successive interactions with the medium, the large-angle region of a jet is dominated by soft daughter partons in the low-virtuality phase and becomes less relevant to the inner core structure directly radiated from the leading parton. Thus, jet suppression has small sensitivity to the details of the leading-parton energy loss, in particular for large jet cone sizes.
Although the difference between the formulations is not visible in the p jet T dependence of jet R AA , it can be seen in the inner structures of jets, particularly in the core region dominated by the leading parton; Jets are more likely to have particles with a larger-p T fraction in their core for the case with the virtuality-dependent formulation. Thus, one may see more details in the energy loss of high-virtuality partons by studying the medium modification of jet substructure observables, e.g. jet fragmentation function, which will be discussed in an upcoming effort. FIG. 5. Ratio of differential cross-section for inclusive charged-particle at mid-rapidity in p + p collisions. The ratio is taken w.r.t. the default pythia 8. The solid red lines and dashed blue lines show the results from jetscape and pythia 8, respectively. Statistical errors (black error bars) and systematic uncertainties (grey bands) are plotted with the experimental data. Top panel: Results for inclusive charged particle with η < 1.0 at √ s = 5.02 TeV, compared to CMS data [118]. Middle panel: Results for inclusive charged particle with η < 1.0 at √ s = 2.76 TeV, compared to CMS data [34]. Bottom panel: Results for charged pion with y < 0.35 at √ s = 200 GeV, compared to PHENIX data [119]. FIG. 6. Nuclear modification factor for inclusive jets and charged-particle in most central (0-10%) Pb+Pb collisions at √ sNN = 5.02 TeV from matter+lbt simulations for three differentq · f formulations within the jetscape framework. Top panel: Results for inclusive jet RAA with R = 0.4 and yjet < 2.8, compared to ATLAS data [114] (black circles) and CMS data for ηjet < 2.0 [120] (magenta squares). Middle panel: Results for inclusive jet RAA with R = 0.4, yjet < 0.3 and p lead, ch T > 7 GeV, compared to ALICE data [115]. Bottom panel: Results for charged-particle RAA with η < 1.0, compared to CMS data [118].

C. Exploring A+A model parameters
In this section, we explore the sensitivity of the free parameters in the multistage jet quenching model (matter+lbt). We shall employ the virtuality dependent formulation (Type 3) as it gives the best simultaneous description of the data. The free parameters in the jet quenching model are summarized in Tab. II. By changing the parameters from the default values, we present the inclusive jet R AA and charged-particle R AA at most central (0-10%) collision at √ s NN = 5.02 TeV and show the parameter dependence in our model.
1. Sensitivity to jet-medium coupling α fix s In Fig. 7, we present the nuclear modification factor for inclusive jets and charged-particles for a jet-medium coupling parameter α fix s = 0.25, 0.3 and 0.35. The jetscape results are obtained using the multistage jet energy loss approach (matter+lbt) where a virtuality dependent factor f (Q 2 ) given in Eq. (27) is employed. Increasing α fix s from 0.25 to 0.35 leads to suppression (≈ 10%) in the inclusive jet R AA for all jet p T . A similar trend is observed in the charged-particle R AA as well. Since the variation of Q sw from 1 to 3 GeV increases the effective length of the LBT energy loss stage, the parton at lowp T undergoes significant energy loss. This leads to an overall suppression of R AA for inclusive jets and charged particles.
Comparison of both inclusive jet R AA with AT-LAS/CMS data and charged-particle R AA with CMS data favors α fix s to be between 0.

Sensitivity to switching virtuality parameter Qsw
In Fig. 8, we present inclusive jet R AA and chargedparticle R AA for the switching virtuality parameter Q sw = 1 GeV, 2 GeV and 3 GeV. The jetscape results are obtained using the multistage jet energy loss approach (matter+lbt) where a virtuality dependent factor f (Q 2 ) given in Eq. (27) is employed. Increasing Q sw from 1 GeV to 3 GeV leads to suppression (≈ 10%) in the inclusive jet R AA for all jet p T . Similar trend is observed in the charged-particle R AA as well. Since the variation of Q sw from 1 GeV to 3 GeV increases the effective length of lbt energy loss stage, where there is no virtuality-driven suppression of the medium effect. This leads to an overall suppression of R AA for inclusive jets and charged particles.
Comparison of inclusive jet R AA with ATLAS and CMS data indicate Q sw to be between 2 GeV and 3 GeV, whereas ALICE data favors Q sw 2 GeV. Moreover, CMS data for the charged-particle R AA disfavors Q sw = 1 GeV in the low-p T region and Q sw = 3 GeV in the highp T region. Overall, the optimized value of Q sw from the above comparison comes out to be about 2 GeV.

Sensitivity to start time of jet energy loss parameter τ0
In Fig. 9, we present inclusive jet R AA and chargedparticle R AA for the jet quenching start time parameter τ 0 = 0.3 fm/c, 0.6 fm/c and 0.9 fm/c. The jetscape results are obtained using the multistage jet energy loss Increasing τ 0 from 0.3 fm/c to 0.9 fm/c does not seem to affect the inclusive jet R AA or the charged-particle R AA . The effect of the starting time of jet-medium interaction as studied in Refs. [124,125] is not visible in the modifications of inclusive jet and charged-particle spectra within our model. This is primarily due to the fact that the jet initiating parton is highly virtual, and the virtuality-dependence in Eq. (26) highly suppresses the medium effect in the early stage. It is still possible, though unlikely that the choice of τ 0 will affect the azimuthal anisotropy. This will be explored in a future effort.

Sensitivity to temperature parameter Tc in jet energy loss
We vary the temperature cut-off parameter for jet energy loss T c = 150 MeV, 160 MeV and 170 MeV, and present results for inclusive jet R AA and charged-particle R AA in Fig 10. The jetscape results are obtained using the multistage jet energy loss approach (matter+lbt) where a virtuality dependent factor f (Q 2 ) given in Eq. (27) is employed. All other parameters are set to their default value. Increasing T c from 150 MeV to 170 MeV decreases the in-medium portion of the jet energy loss. Since, T ∈ [150, 170] MeV corresponds to the late time dynamics of jet energy loss, the relevant energy loss stage is lbt. The decrease in the effective length of the lbt energy loss stage enhances the low p T region of the charged-particle R AA significantly, compared to highp T region. This leads to an overall enhancement (≈ 10%) in the inclusive jet R AA at all jet p T .
Comparison of inclusive jet R AA with ATLAS and CMS data indicate T c to be between 150 MeV and 160 MeV, whereas ALICE data favors T c 160 MeV. Moreover, CMS data for the charged-particle R AA favors T c ∈ [160, 170] MeV. Overall, the optimized value of T c from the above comparison comes out to be 160 MeV.

D. Inclusive jet and hadron suppression at semi-peripheral collisions
In this section, we use the optimized value of the free parameters listed in Table II and present comparisons for centrality dependence of inclusive-jet R AA and charged-particle R AA at √ s NN =5.02 TeV. The jetscape results are obtained using the multistage jet energy loss (matter+lbt), where we employ a virtuality dependent factor f (Q 2 ) along with theq in Eq. (26), to account for a reduction in the rate of medium induced emission in the high virtuality phase. Figure 11 shows our inclusive-jet R AA results for different centrality classes at √ s NN =5.02 TeV, compared with the experimental data from ATLAS [114]. Our full results (solid red lines) are typically consistent within the uncertainties of the ATLAS data for centralities 20-30%, 40-50%, and 50-60% for all jet p T , while we observe a deviation of about 5% in highest jet p T bin at 400-600 GeV for centrality bins 10-20% and 30-40%. For the case of charged-particle R AA shown in Fig. 12, our multistage jet quenching model can describe CMS data [118] in the high-p T region very well with a deviation 10% in the low-p T region.
The deviation from the data in charged-particle R AA at p T < 30 GeV, for more peripheral events, is mainly due to the absence of the jet energy loss in the hadronic phase. As one moves away from central collisions, the expectation is that the hadronic phase will be a larger fraction of the entire system and hence play a nonnegligible role in the quenching of jets in semi-peripheral and peripheral collisions. We emphasize that the calculation presented here employed the same free parameters as the calculation for most-central (0-10%) collisions at √ s NN =5.02 TeV, and no further re-tuning of the free parameters has been performed.

E. Effects of medium response and hole subtraction on inclusive jets
In this section, we highlight the importance of the recoil-hole formalism and demonstrate their effect on inclusive jet R AA at √ s NN =5.02 TeV. We use the optimized value of the free parameters listed in Tab. II and present the effect of hole subtraction in inclusive jets at different centralities (dashed blue lines) in Fig. 11. The jetscape results are obtained using the multistage jet energy loss (matter+lbt), where we employ a virtuality dependent factor f (Q 2 ) that modulates the effective value ofq in the high virtuality (matter) stage to account for a reduction in the medium induced gluon radiation rate due to coherence effects [Eq. (26)].
We remind the reader that the hole hadrons are the product of hadronization of the thermal partons sampled from the medium during the jet-medium interaction. The jetscape framework keeps track of such partons, which are then used to determine the background correlated to the jet. For inclusive jets, we subtract the hole hadrons according to the criteria discussed in Eq. (28).
The effects of holes subtraction are roughly 5% at jet p T > 400 GeV and becomes gradually larger as one goes to lower value of jet p T . The calculation demonstrates that, for inclusive jets, hole subtraction correctly reproduces the jet p T dependence for all centralities, except the deviation of roughly 5% at the highest jet p T bins for centrality 10-20% and 30-40%. The deviation from the data at low jet p T can be attributed to the fact that we do not have jet energy loss in the hadronic phase. In this section, we present the comparisons for nuclear modification factor for inclusive jets and charged-particle R AA at lower collision energies and demonstrate that the multistage jet quenching model, with a recoil-hole formalism and a virtuality dependent factor f (Q 2 ) that modulates the effective value ofq in the high virtuality phase, captures the essential aspects of the parton energy loss in the QGP.
First, we present the inclusive jet R AA for most central (0-10%) PbPb collisions at √ s NN =2.76 TeV in the top panel of Fig. 13. The inclusive jets are constructed using anti-k T algorithm with cone size R = 0.4 and |η jet | < 2, and compared with CMS data [116]. The theory calculation shows a very good agreement with the experimental data for all jet p T . In the bottom panel of Fig. 13, we present the charged-particle R AA for most central (0-5%) collision at √ s NN =2.76 TeV. The comparison of the theory calculation with the CMS data [34] is quite remarkable. Second, we present in Fig. 14 the comparisons for the charged-particle jet R AA and charged-pion R AA at RHIC collision energy √ s NN =200 GeV. The charged-particle jets are constructed using anti-k T algorithm for cone sizes R =0.2, 0.3 and 0.4 with kinematic cut |η jet | < (1.0−R). We impose the leading charged-particle trigger bias p lead,ch T > 5 GeV to select true jets in the same way as done in the experiment and compare the results with STAR data [121] at most central (0-10%) Au+Au collisions. The charged-jet R AA for all three jet cone sizes show a good agreement with the STAR measurements. In the bottom-right panel (Fig. 14), we show comparisons of the charged-pion R AA at most central (0-10%) RHIC collision energy √ s NN =200 GeV. The theoretical calculations are compared to PHENIX data [122], where we only compare with data available at top RHIC energy that covers hadron p T up to 20 GeV. It should be pointed out that such comparisons are meaningful due to the iso-spin symmetry between charged pions and neutral pions. Here again, the theoretical calculations are in good agreement with the experimental data at all hadron p T . Next, we present the prediction of inclusive jet R AA in Fig. 15 for most-central (0-10%) Au+Au collisions at √ s NN = 200 GeV. The calculation is performed using the multistage jet quenching model (matter+lbt) with virtuality dependence (Type 3) for the default values of free parameters presented in Tab. II. We have shown inclusive jet R AA with the jet cone size R = 0.4 for kinematic cuts |η jet | < 0.5 (solid red line) and |η jet | < 1.0 (dashed blue line). The jets show significant suppression, but a weak jet p T dependence.

VI. SUMMARY AND DISCUSSION
In this manuscript, we have presented a multistage [matter+lbt (recoils on) + colorless hadronization] jet quenching model within the jetscape framework and demonstrated, for the first time, a simultaneous description of the nuclear modification factor for inclusive jets and single hadrons from the top RHIC to the top LHC collision energies. We covered three collision energies √ s NN =5.02 TeV, 2.76 TeV, and 200 GeV and performed model-to-data comparison for selected data sets from ALICE, ATLAS, CMS, PHENIX, and STAR experiments. Event-by-event bulk medium simulations, without jets, were carried out first and calibrated to data [65]. Binary collision profiles extracted from these individual simulations are sampled to yield locations of hard scatter- performed using a multistage jet quenching model (matter+lbt). Jet transport coefficient multiplied by virtuality dependent factor [q run HTL f (Q 2 )] is used. The free parameters employed in the jet quenching model are extracted from simultaneous fit to inclusive jet RAA and charged-particle RAA at most central (0-10%, √ sNN=5.02 TeV) Pb+Pb collisions (top left plot for jets) and no further re-tuning has been performed. Also shown in the dashed blue lines, is the effect of not subtracting the holes.
Results are compared to ATLAS data [114] (black circles) in all centrality cases and CMS data for ηjet < 2.0 [120] (magenta square) in only the 0-10% case. [q run HTL f (Q 2 )] is used. The free parameters employed in the jet quenching model are extracted from simultaneous fit to inclusive jet RAA and charged-particle RAA at most central (0-10%, √ sNN=5.02 TeV) Pb+Pb collisions and no further re-tuning has been performed. Results are compared to CMS data [118] for both centralities.
ing. The pythia generator with ISR and MPI turned on is used to simulate hard scatterings that produce final state partons, without any final state shower. These are transferred to the matter generator, where they are imbued with a timelike virtuality Q, which depends on their transverse momentum p T .
To incorporate the multiscale dynamics of jet energy loss within the jetscape framework an effective parton evolution was set up, in which we encoded the space-time profile of the QGP, obtained from the bulk simulations, within the parton energy loss process. The initial high virtuality stage is modeled by the matter event generator, followed by the low virtuality stage, modeled by the lbt event generator. The switching of jet energy loss stage from high virtuality to low virtuality is carried out on a parton-by-parton level, depending on the virtuality of the parton: Those with Q > Q sw remain in matter, while those partons, whose virtuality drops below Q sw , in the process of multiple emission, are transferred to lbt. Further medium-induced emission within lbt is assumed to maintain the virtuality at or below Q sw . Partons that escape the medium and still have virtuality larger than Q 0 = 1 GeV undergo vacuum evolution using the matter module. A weakly-coupled description of the medium in terms of thermal partons (recoils/holes) is used to include the medium response to the jet.
We systematically explored the three functional forms of the parameter dominating the jet-medium interaction strength. We demonstrated that the inclusion of a virtuality dependent factor [f (Q 2 )] which modulates the effective value ofq to account for a reduced medium-induced emission in the high virtuality phase, due to coherence effects [7,9], is essential for a simultaneous description of inclusive jet R AA and charged-particle R AA at LHC and RHIC collision energies.
The success of this approach, in comparison with inclusive jet and hadron data, at all energies and centralities (at both RHIC and LHC), may imply thatq/T 3 does not have a cusplike behavior at 300 MeV< T < 400 MeV. Moreover, the effective reduction inq/T 3 , at LHC collision energies compared to RHIC, is mainly due to the fact that the energy of the jets produced at the LHC is an order of magnitude higher, compared to those produced at RHIC energies. The jets emanating from the hard parton, at LHC collision energies, start out with significantly higher virtuality, and hence, experience a significantly smaller stimulated emission rate, compared to RHIC energies. Alternatively, as the virtuality of the hard parton increases, the transverse size of the dipole formed by the hard parton and the emitted gluon decreases, due to which the plasma starts to appear more dilute (one often states that the partons in the medium cannot resolve the different parts of the dipole for stimulated emission). While coherence effects have been included in the high virtuality phase, effects such as the rise in the effectiveq in radiative processes compared to purely scattering processes, at lower virtualities, have not [126,127]. However, these effects can be approximately applied after extraction ofq has been made. We explored A+A model parameters: α fix s , Q sw , τ 0 , and T c , and studied their effects on inclusive jets and charged-particle R AA . These explorations were carried out using a virtuality dependent factor f (Q 2 ) effectively modulatingq in the matter phase, as this gives the best simultaneous description of the inclusive jet and charged-particle R AA data. The study suggests that the τ 0 parameter does not seem to have much effect on inclusive jets and charged-particle spectra, indicating that FIG. 13. The inclusive jet RAA and charged-particle RAA at most central (0-5%) Pb+Pb collisions at √ sNN=2.76 TeV. The calculation is performed using the multistage jet quenching model (matter+lbt) with virtuality dependence (Type 3). The free parameters employed in the jet quenching model are extracted from simultaneous fit to inclusive jet RAA and charged-particle RAA at most central (0-10%, √ sNN=5.02 TeV) Pb+Pb collisions and no further re-tuning has been performed. Top panel: Results for inclusive jet RAA with R = 0.4 and ηjet < 2, compared to CMS data [116]. Bottom panel: Results for inclusive charged-particle RAA with η < 1.0, compared to CMS data [34].
the medium effects on the jet energy loss may be highly suppressed during the early stage. The jet interaction termination temperature T c seems to have a noticeable effect on both jet and hadron spectra, particularly at low p T . The neglect of jet energy loss in the hadronic phase is the most probable reason for the reduced suppression for p T < 30 GeV, particularly at more peripheral collisions, which have a proportionately larger hadronic phase. In spite of these minor offsets, the model proposed in this paper presents a good comparison with a wide swath of experimental data.
Most of the results presented in Sec. V included simulations with a modulating factor f (Q 2 ) that reduces the effectiveq in the matter phase to account for the re-duction in the medium induced radiation. We have also carried out similar simulations without this modulating factor. These results are presented in Appendix A and Appendix B. As the reader will note, neither of these formulations can describe the inclusive charged particle R AA and the inclusive jet R AA simultaneously. The emerging picture from these simulations is that the modification of jets is dominated by the migration of the softer widerangle components beyond the jet cone, with only minor modifications of the core region. Thus a considerable portion of jet energy loss is possible in the nonperturbative region. This is discussed in Appendix C, and will be further explored in the upcoming effort on jet medium correlations. Along with these studies, other efforts will focus on the heavy-quark sector, jet substructure, and azimuthal anisotropies.
Our simulations have increased the number of parameters typically invoked in jet simulations. Some of this increase is simply due to the increased sophistication of a multistage jet modification framework. There is currently no well-established theory for how to transition from the higher virtuality phase simulated by matter to a lower virtuality phase simulated by lbt. While we have used an energy scale Q sw , in an effort to get a direct handle on the average scale of the transition, a more apt method may use Q sw = C √ 2qE. The case for C = 1 was already explored in Ref. [71], however, no comparison to data was carried out. The parameters in f (Q 2 ) are a fit to the form derived in Ref. [9], however, this could also be further parametrized as more observables are included in a more global fit. We point out that another user of the jetscape framework may easily choose to replace either one or both of the energy loss modules that we used with any number of their own modules, with different criteria for transition. This may lead to a different collection of parameters.
More parameters are expected to arise in the introduction and simulation of the nonperturbative wake of the jet [96], which is required for the simulation of the jet medium interactions and jets at large angles [120]. The physics of soft energy-momentum spreading away from a jet is inherently nonperturbative, and thus, the inclusion of new parameters is unavoidable. The current effort evades the need for many of these parameters by limiting its scope to only jets and single hadrons. As each new tranche of parameters is introduced, a new Bayesian analysis should be carried out to constrain and correlate these against the existing parameters. FIG. 14. The charged-particle jet RAA and charged-pion RAA at most central (0-10%) Au+Au collisions at √ sNN=200 GeV. The top two-panels and the bottom left panel are charged-jet RAA for cone sizes R=0.2 (ηjet < 0.8), 0.3 (ηjet < 0.7), and 0.4 (ηjet < 0.6), compared to STAR data [121]. The bottom right panel is charged-pion RAA compared to PHENIX data of neutral pion RAA [122]. The calculation is performed using the multistage jet quenching model (matter+lbt) with virtuality dependence (Type 3). The free parameters employed in the jet quenching model are extracted from simultaneous fit to inclusive jet RAA and charged-particle RAA at most central (0-10%, √ sNN=5.02 TeV) Pb+Pb collisions and no further re-tuning has been performed.  is visible for both jets and charged particles throughout almost the entire p T range. Most notably, the α fix s dependence in jet suppression for the case with the fixedcouplingq is much larger than that for the case with virtuality-dependent formulation (top and middle panels in Fig. 7). Since a larger number of daughter partons, whose interactions with the medium give the main contribution to the jet energy loss, are branched off the leading partons, the stronger sensitivity is seen in the fixedcouplingq case. This strongly implies that the modification pattern of the inner jet structure can be very different between the formulations, Type 1 and 3, even when we tune the coupling strength to make their jet energy loss the same. Figure 17 shows results of the nuclear modification factor of inclusive jet with R = 0.4 and charged particle in Pb+Pb collisions at √ s NN = 5.02 TeV from jetscape (matter+lbt) calculations with the fixedcoupling HTLq (Type 1) for different switching virtuality Q sw = 1, 2, and 3 GeV. The other free parameters are set to their default values shown in Tab. II. The same trend as that from the virtuality-dependent formulation (Fig. 8) is shown here: stronger suppression with increasing Q sw in both inclusive jet and charged particle R AA .
Here we emphasize that the rising behavior of the data in the charged particle R AA as a function of p T cannot be reproduced with any values of the main free parameters α s and Q sw when we employ the formulation of the fixedcoupling HTLq (Type 1).
Appendix B: Parameter dependence for Type 2: HTLq with running coupling Figure 18 presents the nuclear modification factor of inclusive jet with R = 0.4 and charged particle in Pb+Pb collisions at √ s NN = 5.02 TeV calculated using jetscape (matter+lbt) with the running-coupling HTLq (Type 2) for different coupling strengths α fix s = 0.2, 0.25, and 0.3.
The same trend as that from the other formulations can be seen also here: stronger suppression with increasing α fix s in both inclusive jet and charged particle R AA . Here the strength of the dependence on α fix s in jet suppression is closer to that for the fixed coupling case (Type 1) than for the virtuality-dependent case (Type 3). This means that the effect of the running coupling does not affect the jet inner structure drastically.
In Fig. 19, we show the nuclear modification factor of inclusive jet with R = 0.4 and charged particle in Pb+Pb collisions at √ s NN = 5.02 TeV obtained from jetscape (matter+lbt) calculations with the runningcoupling HTLq (Type 2) for different switching virtuality Q sw = 1, 2, and 3 GeV. The other free parameters are set to their default values shown in Tab. II. The trend of the stronger suppression with increasing Q sw is also shown here in both inclusive jet and charged particle R AA . Again we would like to note that the data for the p T dependence of the charged particle R AA cannot be de- scribed with any values of the main free parameters α s and Q sw unless the virtuality-dependence via f (Q 2 ) is introduced. In this appendix, we discuss how nonperturbative effects can modify the reconstructed jet spectra. An analysis of the p T distribution of partons in jets, that have propagated through a dense medium, shows a large population of soft partons within quenched jets. This leads to noticeable shifts in the jet spectrum during hadronization carried out by string-based models. As outlined below, one can highlight these shifts by comparing distributions of jets constructed using final partons versus those constructed using final hadrons. In this section, we study these differences and discuss alternative approaches.

Comparison with Partonic Jet
The contribution of nonperturbative effects in our model is limited mostly in the hadronization process. To quantify the effect, we conduct an analysis also for the jets reconstructed from the final state partons, just before being passed to the hadronization module. In Fig. 20, the R AA of the inclusive partonic jet for most-central (0-10%) Pb+Pb collisions at √ s NN = 5.02 TeV is compared with the full calculation results with hadronization, for the case of running coupling, with (Type 3) and without (Type 2) the virtuality dependent modulation factor f (Q 2 ). Note that the denominator of R AA , for partonic jets, is the jet spectrum in p + p collisions calculated by turning off the hadronization module. The ratio between the spectra for parton jets and hadron jets in p + p collisions is shown in Fig. 21. Thus, we find that hadronization effects are not significant in p + p collisions and gives almost no modification for jets with p jet T 200 GeV. In both cases, with and without the virtuality dependence, additional suppression from the colorless string hadronization can be seen. This additional suppression comes from the modification of the soft parton spectrum prior to hadronization. As jets propagate through a dense medium, a considerable number of low-energy partons are branched off collinearly in the jet shower. In string hadronization, strings connecting such collinear soft partons do not have the minimum mass necessary to produce hadrons, and cannot be included in the hadronization process as is.
In this work, we follow the methodology devised in our prior work on jets in p + p collisions [42]: Consider a pair of partons that possesses a |δ p| < 4Λ QCD . The parton with larger |p z | in this pair has its p z shifted to p z ± 4Λ QCD . The sign of the added momenta is set in a way to increase the relative momentum |δ p|. As a result, no net z-component of momentum is added to jets on average. The energy of the modified parton is changed accordingly to ensure that it remains on-shell. One continuously iterates through the entire parton list until no pairs meet the condition |δ p| < 4Λ QCD . While cast differently, in terms of the three-momentum, this condition is consistent with the minimum mass condition in Ref. [132]. The reader will note that the p T of the jet is not affected by this procedure, though the energy and mass of the jet may be affected. The predominant effect of this procedure is that a large fraction of the soft partons are pushed to larger rapidities and outside the jet cone, leading to extra jet suppression (energy loss). This procedure does not affect the hard partons in a jet and thus has no effect on the single hadron spectrum or its nuclear modification. In Fig. 20, the effect of this prescription shows up as a noticeable difference between the spectra of jets clustered with partons and those clustered with hadrons. Whereas this prescription shows almost no difference between partonic jets and hadronic jets in p + p collisions, as shown in Fig. 21. These shifts for jets modified in heavy-ion environments are indicative of the limitations of applying Lund string hadronization to systems with several soft partons.
A large fraction of these soft partons are included with the hard jet in the process of scattering and recoil. In this sense, they represent extra energy and momentum from the medium which has become correlated with the jet shower. As a result, these partons are subject to background subtraction. The exact method in which this background subtraction is carried out affects the difference between the partonic and hadronic spectra but has minimal effect on the final hadronic jet spectrum. While the conventional method of background subtraction was described in Sec. III B, we describe alternate mechanisms and their effect on the parton-hadron offset in jet spectra in the subsequent subsection.

Background Subtraction at Source with Ecut
In this sub-section, we study the effects of background subtraction of jet partons that are at the thermal scale. Throughout this paper, the method of background subtraction has been that of Eq. (28). All partons generated by the jet shower, and those from the medium which scatter with a jet parton, and are included in the modified shower, are retained all the way down to p T → 0. All these partons are assumed to be weakly coupled with the medium and can escape the medium. The four momentum of the incoming partons from the medium, which scatter with the jet partons, are also retained (referred to as holes). These holes are assumed to be free from medium interactions. Once jets have been clustered, the four momenta of these incoming partons (holes) which fall within a jet's area are subtracted from the full four momentum of the jet.
A more approximate algorithm, which allows for faster simulations is to remove all holes from the event record, as they appear in the simulation, along with the softest partons (with energy E in the rest frame of the fluid cell) which range from 0 < E < E cut . The remaining partons with E > E cut are clustered to form jets. The upper limit of E cut is varied until the jet spectrum at the parton level is almost identical to that obtained by the subtraction of holes within clustered jets, as outlined in the paragraph above [i.e., using Eq. (28)].
The physical picture underlying this methodology is that holes, partons that arise from the medium and are scattered by the jet to form recoil partons, are removed from the medium, thus constituting "holes". Also, one would expect the softest partons in the jet to thermalize within the strongly interacting medium. Thus, we are assuming that as the energy lost from the medium to the jet, thermalizes, it balances this negative contribution with the positive contribution of the softest partons correlated with the jet. The upper limit E cut , of the soft parton spectrum which is balanced by the hole contribution is varied to ensure that both methods for background subtraction yield identical results. Our simulations indicate that an E cut = 3.2T , where T is the ambient temperature in the rest frame of the unit cell where the scattering takes place, yields the same jet spectrum as that obtained from Eq. (28), across all energies and centralities.
While the presence of a single such value may be surprising, it should be pointed out that the mean momentum of the hole parsons emanating from the medium, for a given thermal distribution, is of the order of ≈ 2.5T -3.5T . We also point out that since E cut is determined by the comparison between two methods of background subtraction, it does not constitute a new parameter in jet modification.
Along with the increased speed in the simulation, brought on by neglect of holes and soft partons with E < E cut , this second method of background subtraction has one further advantage. Since a large number of the soft partons have been removed from the parton showers, the offset between the partonic jet spectrum and the hadronic jet spectrum is greatly reduced. Note that as discussed in the preceding subsection, the major contribution to the offset is the presence of a large number of soft and collinear partons within the jet shower. In Fig. 22, we present inclusive jet R AA for two different choices of parton energy cuts: E cut = 3.2T and E cut = 2 GeV. The partons with energy E ≤ E cut (in the rest frame of fluid cell) have been removed from the parton shower. The results in Fig. 22 shows that the offset between the parton-level and hadron-level jets goes away if one removes the thermal partons from the parton shower. Moreover, the results also demonstrate that the contribution of such nonperturbative effects is essential to describe the experimental data of R AA for reconstructed jet.

Energy Loss due to Medium Response
As shown in Fig. 22, the removal of holes and partons with E < E cut = 3.2T , yields coinciding hadronic and partonic jet spectra. Also with this value of E cut the two methods of background subtraction yield identical partonic jet spectra. However, with only an E cut = 3.2T the suppression observed in the jet spectrum is not consistent with the data on jet R AA . Comparison with the hadronic jet spectrum in Fig. 20, indicates the presence of further jet modification by medium response.
In a complete simulation of jets in a heavy-ion collision (via the jetscape package), partons with E ≤ E med ≈ 10T would be considered soft enough to be thermalized within the medium [96,109,110,112,133], the fourmomentum of these soft partons would then become part of an energy-momentum source term for a subsequent bulk medium simulation, which would start with the exact initial state that generated both the primary bulk simulation and the distribution of hard scattering, that led to jet production. The space-time-dependent energymomentum source term enacts a boundary between the portion of jet modification that can be carried out using perturbation theory and the part that should be carried out nonperturbatively. This source term takes the four momentum of the partons with E < E med and diffuses these out in space-time using a causal diffusion equation [96]. The second hydro simulation, while starting out with the same initial state as the prior simulation, is modified by the presence of the source term. The hadrons produced in the freeze-out of this new bulk simulation are then combined with the hadrons from the fragmentation of the jet. Jet reconstruction algorithms will have to be run on the "full" event. Background subtraction and unfolding can be carried out by statistical subtraction of the jet distribution clustered from the hadrons produced using the primary bulk simulation (without jets). What we outline above is a computationally challenging problem. A secondary bulk simulation would have to be run for every hard scattering event. One would not be able to avail of the standard methodology of running up to 1000 hard scattering events per a single bulk event. However, due to the diffusion of the source term and the hydrodynamical response of the bulk medium to supersonic energy deposition, one would expect the formation of a Mach cone at an R 1 from the jet axis [134]. The exact Mach angle depends on the equation of state and thus varies as the jet passes through the plasma. However, energy and momentum from partons with E < E med would, over time, propagate away from the jet, at the Mach angle. The loss of these partons from the clustered jet would constitute an additional source of jet energy loss.
In the last part of this appendix we model this process of energy loss via medium response by instituting an additional correlated broadening on partons with E cut < E < E med . Partons with energy in the local fluid rest frame that lie within this range receive a transverse momentum kick such that over a small length δl, k 2 ⊥ =qδl, whereq is the local jet transport coefficient. The difference from the regular process of transverse broadening is that subsequent kicks are always radially away from the jet cone and always in the same direction. These correlated kicks continue until the parton is at or beyond the Mach angle from the jet axis, or it exits the medium, whichever occurs first. After crossing the Mach angle it continues to undergo transverse kicks from the medium, but with subsequent kick directions randomized (regular 2-dimensional diffusion).
This approximate method yields both a loss of partons from within the jet cone and also their reappearance at the Mach angle. The result of this shift on the partonic and eventual hadronic R AA is shown in Fig. 23. As one would notice, the effect of additional medium response with changing α fix s = 0.35 allows us to obtain a simultaneous fit to both the inclusive jet and hadron suppression, and leads to consistent partonic and hadronic R AA .
As the reader will note, the effect of energy loss via medium response seems to have little effect on the hadronic jet R AA ; comparing to the simulations in the earlier sections, its primary effect is to make the partonic R AA consistent with the hadronic one which has been presented in the earlier sections. As a result, for the current work, which only focuses on the inclusive jet and single hadron R AA , the new at-source method of background subtraction, followed by the jet modification due to medium response, has little import. Hence we discuss these in the appendix. The exact nature of how energy is transferred out of the jet cone is relevant to substructure observables such as the jet shape, which will be discussed in our future efforts.