QGP modification to single inclusive jets in a calibrated transport model

We study inclusive jet suppression and modifications in the quark-gluon plasma (QGP) with a transport-based model. The model includes vacuum-like parton shower evolution at high-virtuality, a linearized transport for jet-medium interactions, and a simple ansatz for the jet-induced hydrodynamic response of the medium. Model parameters are calibrated to nuclear modification factors for inclusive hadron $R_{AA}^{h}$ and single inclusive jets $R_{AA}^{j}$ with cone size $R=0.4$ in 0-10% central Au-Au and Pb-Pb collisions measured at the RHIC and LHC. The calibrated model consistently describes the cone-size dependent $R_{AA}^{j}(R)$, modifications to inclusive jet fragmentation functions and jet shape. We discuss the origin of these modifications by analyzing the medium-induced jet energy flow in this model and elucidate the interplay of hard parton evolution and jet-induced medium response. In particular, we demonstrate that the excess of soft hadrons at $p_T\sim 2$ GeV/$c$ in jet fragmentation function and jet shape at large $r=\sqrt{\Delta \eta^2+\Delta \phi^2}$ are consequences of both soft medium-induced gluon radiation and jet-induced medium excitation.


Introduction
Jets produced in nuclear collisions probe the properties of the color-deconfined quark-gluon plasma (QGP). Strong final-state jet-medium interactions suppress the yield of jets as well as high-p T hadrons. This phenomenon, known as the "jet quenching", is among the key signatures of the formation of the QGP [1][2][3][4] and is also used to determine the stopping power of QGP to colored particles [5][6][7][8][9]. Besides yield suppression, measurements of modifications to jet constituents and internal structures become available in recent years at both the RHIC and the LHC energies [10]. It is therefore important to develop a consistent theoretical picture to understand the suppression patterns of jets and hadrons and the modified jet internal structure [11][12][13][14][15][16][17][18][19], which requires the modeling of microscopic parton-medium interactions.
Conceptually, interactions between energetic jet partons and medium constituents fall into two categories: elastic collisions between jet and medium partons and the mediuminduced parton splittings. The formula for a single medium-induced parton splitting [20][21][22][23][24][25][26][27] have been known for decades, and energy loss due to elastic collisions has also been studied in both weakly-coupled [28][29][30][31][32][33][34][35] and strongly-coupled field theory [36]. Applying these theoretical inputs to jet phenomenology in the dense QGP medium, one needs to include multiple interactions, which has led to the use of different evolution equations. For example, many studies use DGLAP-like QCD evolution equations with medium-modified evolution kernels or equations derived from high energy effective theory including medium effects [37][38][39][40][41]. In this work, we take another widely used approach based on transport equations from the kinetic theory [42][43][44][45][46][47][48][49]. In this approach, multiple radiations and collisions are included as the time-evolution of a rate equation as governed by the Boltzmann-type transport equation. The first goal of this work is to develop a simulation procedure around a previously developed linearized partonic transport model LIDO [50]. This simulation procedure will be referred to as "a transport-based model" in this paper.
Given the premise of using a transport equation for the jet-medium interaction, we are confronted with two problems in the modeling.
The first problem concerns with the scale evolution of the parton showers from the hard scale Q ∼ p hard T to Q 0 Λ QCD . Because jets have scale-dependent partonic contents, one has to choose a matching scale where the transport equation in medium starts to apply. Many studies, such as those used the Linearized Boltzmann Transport (LBT) model [47][48][49]51], included the full vacuum shower before applying the transport method, i.e., the transition was made at the same scale as the minimum evolution scale in p-p collisions. The study using the MARTINI model explored two scenarios [45]. One includes the full evolution of the vacuum parton shower before the onset of medium effects and the other uses a formation-time argument to motivate the choice of matching scale Q 0 ∼ p T /τ 0 , where τ 0 is the initial proper time of hydrodynamic evolution of the medium. It was found that different matching procedures affect the phenomenological choice of the effective jet-medium coupling constant. In a series of recent works [18,52], the authors carefully explored the phase-space boundaries between vacuum-like and medium-induced emissions in a finite and static medium. The matching scale between high-virtually vacuum-like emissions and the medium-induced emissions are motivated by the momentum broadening Q 0 ∼ (qE) 1/4 obtained using the multiple-soft collision approximation. The aforementioned studies assume the scale evolution of the parton is either unmodified by the medium or only restricted in emission phase-space. Other models, such as the MATTER event generator [53], include medium-modified splitting functions into the DGLAP-type evolution equation according to a higher twist approach [54][55][56]. The MATTER model has been combined with either LBT or MARTINI transport models in the JETSCAPE framework [57][58][59]. The dependence of model predictions on choices of Q 0 was studied [57,60,61].
In this work, we use the simpler method that the parton evolution above a certain scale Q 0 is unmodified by the medium and the evolution below Q 0 is modeled a time evolution governed by the transport equation. Despite Q 0 is treated as a free parameter in our model, we expect the physical value of Q 0 should be comparable to the momentum broadening k 2 ⊥ acquired during the formation time of an emission in the QGP medium 1 . We will argue that in a fast expanding medium, typical values of k 2 ⊥ are relatively insensitive to the parton energy and are mainly determined by medium properties. Therefore, it is reasonable to use a single matching scale Q 0 for all partons.
For partons going out of the QGP fireball, QGP modifications vanishes and we apply the vacuum-like evolution again, neglecting jet interactions with the hadronic matter. However, one should note that this procedure is contradictory to the current cross-over picture of the phase transition at small baryon chemical potential, where medium properties are continuous across T c . Moreover, a recent study indicates substantial interactions between high-p T pions and the hadronic matter near T c [62]. We will leave the study of jet modifications in the hadronic phase to future works.
The second problem is that a transport equation is only defined for well-defined quasiparticles. In the QGP medium, only hard particles with energy much greater than the temperature of the medium (E T ) can be identified as good quasiparticles, including both hard partons from partons showers as well as recoiled medium partons that acquires a large momentum in the collision with a jet parton. However, the change in the dynamics of the soft medium in response to the passing of a jet is beyond the scope of the linearized partonic transport equation with weakly-coupled interactions. The perturbed medium consists of particles with thermal momentum p ∼ O(T ) and its dynamics are better modeled in a hydrodynamic approach. Models that couple a transport model to hydrodynamic [51,63], models based on instant thermalization [64] and linearized hydrodynamic response [65], and methods of including medium recoil particles [47,66] have been used to address such a problem. In this work, we will develop a simple ansatz to model the hydrodynamic-like response of medium, including the freeze-out of particles. We use a single term "jet-induced medium response" to refer to both the hydrodynamic-like response of soft medium dynamics and the hard recoiled medium partons.
We will introduce this transport-based model in detail in section 2, including a review of the LIDO transport model and our approach for addressing the above-mentioned problems.
For reliable predictions with quantified uncertainty, we will perform a Bayesian model calibration of the model parameters to the inclusive jet and hadron suppressions in central nuclear collisions at both the RHIC and the LHC. This also tests the model's ability to describe both jet and hadron suppression simultaneously.
In section 4, using the calibrated model, we analyze the energy flow induced by jet-medium interactions in the two-dimensional space (z jet , r), where z jet is the longitudinal momentum fraction of particles relative to the momentum of the jet, and r is the radial distances between particles and jets Mapping out this energy flow and differentiating effects from different transport mechanisms help to interpret our prediction on jet modifications. Qualitatively, medium-induced radiations predominately transfer energy from high-z jet to low-p T particles, while elastic collisions and jet-induced medium response carry energy to large angles.
In section 5, we predict modifications to jet fragmentation functions and jet shape as well as the cone-size dependence of jet suppression using the calibrated model. We try to identify different energy transport mechanisms in these modifications, which provides a deeper understanding of jet-medium interaction.
Finally, we summarize in section 6 our results and comment on the jet transport parameter determined in this work.

A transport-based model for jet evolution in nuclear collisions
In this section, we introduce the jet evolution model that has been outlined in the introduction. It involves three stages of evolution: 1. The initial hard processes of jet production and vacuum-like evolution of parton shower from Q ∼ p hard T to Q = Q 0 is modeled by Pythia8.
2. The initial jet shower parton production is followed by partons-medium interactions, which are treated in a previously developed linearized partonic transport model LIDO.
Meanwhile, jet partons excite the medium in two ways: medium partons that acquire a hard recoiled momentum in collisions with jet partons; the perturbed soft medium partons seed a linearized hydrodynamic-like response. The subsequent evolution of the hard recoil partons is also treated by the transport equation. Together, hard recoils and the hydrodynamic response are termed as the jet-induced medium excitation.
3. A stage of vacuum evolution of parton shower from Q 0 to Λ QCD occurs outside the QGP medium. It is followed by hadronization and particle decays. This stage is again modeled by Pythia8.
We will discuss models for each stage in detail and the matching procedures between them. In addition, the soft QGP medium evolution is described by a 2+1D viscous hydrodynamics with event-averaged initial conditions. Due to the lack of event-by-event fluctuations and longitudinal evolution, we will only focus on midrapidity observables in most central nuclear collisions at the RHIC and the LHC in this paper. Details of the medium evolution model can be found in appendix A.

Jet evolution in vacuum and initialization of transport equation
We use the Pythia8 event generator [67] with ATLAS A14 central tune [68] to sample initial hard scatterings and generate the initial parton showers. The proton partron distribution function (PDF) in Pythia8 takes the leading order CTEQL1 PDF [69], and the nuclear PDF takes the leading-order EPS09 parametrization [70].
The initial parton radiations are ordered in decreasing virtuality, which is just the transverse momentum of the radiated particle in this case. In proton-proton collisions, the scale Q starts from that of the initial hard collision p hard T down to a cut-off scale, which is close to Λ QCD . In nuclear collisions, we have argued that the high-virtuality vacuum-like evolution should be stopped scale Q 0 , comparable to the virtuality of partons acquired in the medium. The latter can be estimated in a parton splitting process including medium effects of momentum broadening: a parton with energy E radiates an almost collinear parton with energy xE and transverse momentum k ⊥ ∼ Q relative to the original parton. Within the radiation formation time collisions with medium broaden the transverse momentum to, whereq = d ∆k 2 ⊥ /dt is the jet transport parameter. Medium impacts are small when Q 2 ∆k 2 ⊥ , where we apply the vacuum-like evolution. For smaller Q 2 , the typical transverse momentum saturates at ∆k 2 ⊥ and the vacuum-like parton shower gives way to the inmedium evolution governed by the transport equation.
In the introduction, we have summarized different choices of Q 0 in the past studies, including the use of fixed Q 0 and parton-energy-dependent choices, such as Q 0 ∝ E 1/2 obtained using formation time argument and Q 0 ∝ E 1/4 motivated by the momentum broadening in a static medium due to multiple-soft collisions. Here, we argue that because of the fast dropping of medium temperature due to expansion, Q 0 ∼ ∆k 2 ⊥ is relatively independent of the parton energy. In section 3.3, We will verify this argument by running the model in a medium simulated by the hydrodynamics and show dependence of the distribution of ∆k 2 ⊥ on the parton energy. Right now, we first demonstrate this argument using a medium under Bjorken expansion [71]. In a medium with Bjorken flow solution q ∝ T 3 ∝ τ −1 orq(τ )τ ≈q(τ 0 )τ 0 . The momentum broadening is estimated from the self-consistent relation Approximately, k 2 ⊥ ≈q(τ 0 )τ 0 ln(E/q(τ 0 )τ 2 0 ) and it only increases logarithmically with the parton energy. Therefore, at given beam energy, collision system, and centrality, it is reasonable and also practical to use a single Q 0 scale to transit the entire parton shower from vacuum-like scale evolution to parton transport in the medium. Finally, to initialize the transport equation, which requires both the momentum and spatial information, we assign space-time coordinates for partons created in the initial parton shower. First, a transverse position x ⊥ is sampled according to the density of binary collisions in the transverse overlap region. The coordinates x µ = (0, x ⊥ , 0) are assigned to the hardest partons directly created from the hard collision. Consider parton a with energy E and coordinate x a splits into partons b and c, and parton b has a larger energy xE than parton c, x > 1/2. We assume that the radiation of the softer parton is delayed by its formation time τ f and is unmodified by the medium. The production space-time coordinate of parton c is

The LIDO transport model
The jet shower partons at scale Q 0 provide the initial condition for the transport equation.
We use a previously developed linearized partonic transport model LIDO [50] for this stage. This model relies on several assumptions: • The first assumption is the weakly-coupled interactions in which the coupling constant α s or g s = √ 4πα s are parametrically small. In this limit, jet partons interact with medium quasiparticles via screened gluon exchanges. We use the three-flavor leadingorder running coupling constant with Λ QCD = 0.2 GeV. At finite temperature, we assume that the running scale has to be greater than a medium scale µ min ∝ T . This results in the following ansatz of running coupling at finite temperature, (2.5) • The linearization assumption: the equation only treats hard particles with p T , where the population is small compared to thermal partons in the medium. We also neglect collisions between hard partons and their interactions with the jet-induced medium excitation. Therefore, the transport equation is linearized with respect to hard partons.
• The final assumption requires the background medium is composed of massless partons in local thermal equilibrium. Thermal partons are then sampled from Boltzmanntype distribution function f 0 (p) = e −p·u/T , where u, T are the flow four-velocity and temperature of the local medium obtained, respectively, from the hydrodynamics simulation.
Given these premises, a Boltzmann-type transport equation evolves the hard-parton distribution function f (t, x; p) The left-hand side of the above equation is the ballistic transport term. We have inserted step functions to emphasize that the transport equation only handles hard partons with energy greater than E min in the co-moving frame. Throughout this work, E min = 4T , which is slightly larger than the average kinetic energy of a thermal parton.
Elastic jet-medium interactions Elastic interactions refer to the processes that conserve the number of hard partons. In the LIDO model, elastic collisions that only involve small-momentum transfer q ⊥ < Q cut transverse to the hard parton are included in drag and diffusion operator D, Here, theq S denote the transport parameter due to soft interactions [46] q a,s (p, where a represents quarks or gluons. The respective Casimir factors are C A = 3 for gluons (a = "g") and C F = 4/3 for quarks (a = "q"). m D = 3/2g(µ min )T is the Debye screening mass of the QGP with three flavors of light quarks and η D is the momentum drag coefficient determined by the fluctuation-dissipation relation η D =q S /(4ET ). We set Q cut = min{2m D , √ 6ET } throughout this work. Here √ 6ET estimates the average center-of-mass energy between jet partons and medium partons 2 .
Elastic collisions that involve large momenta transfer q ⊥ > Q cut are included explicitly in a two-body collision term C 2↔2 for particle specie a, where p b , p c and p d , p e are four momenta of the initial-state partons and final-state partons of the collision, respectively. q = p b − p d is the four-momentum change between the hardest initial-state and final-state particle in the co-moving frame. |M (q)| 2 bc,de is the squared amplitude of the two-body collision in the vacuum 3 , averaging over initial-state and summing over final-state spins and degeneracies. The running coupling constant is evaluated 2 Separating small-q collisions into a diffusion equation seem to be an unnecessary step if we only want to include leading order weakly-coupled physics. But a diffusion formulation has a different range of applicability than matrix-element based collisions. For example, one can easily implement diffusion coefficients from strongly-coupled theories for the soft-momentum exchange process. 3 The difference between screened propagators and propagators in the vacuum is assumed to be small when Qcut mD [46], rationalizing the use of vacuum matrix-elements.
is the distribution function of the initial-state hard parton, and f 0 (p c ) is the equilibrium distribution function of the initial-state medium parton. We first integrate over the entire phase space of the initial-state {p i } and final-state {p f } particles (2.10) Then, we insert in the second line of equation 2.9 a set of δ-functions in momentum and particle specie projectors to project out contributions to the distribution function of hard parton species a. The first negative term in the projector corresponds to the loss term, and the second summation that goes over the final-state particles corresponds to the gain term. Finally, the expression sums all reaction channels, and ν c is the degeneracy factor of the medium parton c, including spins, colors, and the three light flavors: up, down, and strange.
Because of the separation in q ⊥ , the jet transport parameter for parton specie a is also composed of a soft part and a hard part, whereq a,s has been defined in Equation 2.8, and we getq a,h by computing the average momentum broadening due to large-q scatterings per unit timê Medium-induced parton radiation Medium-induced parton radiations change the number of hard partons in collisions. Similar to the separation for elastic collisions, inducedradiations are separated into a diffusion-induced radiation term C 1↔2 and a radiation term C 2↔3 induced by large-q collisions with medium partons. The formulae that follow immediately only applies to the Bethe-Heitler limit of incoherent radiation, where the formation time of the radiation τ f is much shorter than the collisional mean-free-path λ el . We will then introduce corrections to these formulae in the deep-Landau-Pomeranchuk-Migdal regime of τ f λ el . The large-q collision-induced radiation term C 2↔3 is defined similarly as in equation 2.9, except that the radiation introduces another hard parton in the final state, labeled by four-momentum k The first term on the right includes two-to-three-body collisions for parton radiation. Its squared amplitudes are listed in reference [50] in the appendix. The coupling at the radiation vertices is evaluated at scale k ⊥ -the radiated parton's transverse momentum relative to the original parton. Rigorously speaking, there are also detailed balancing processes corresponding to medium parton absorption. However, the balancing term only becomes important when the distribution f (p) is close to kinetic equilibrium. Because we only consider particles with p · u > 4T in the transport equation, the absorption term is neglected. The diffusion-induced radiation collision integral is obtained by taking the |q| k ⊥ limit of two-to-three-body squared amplitudes and integrating the phase-space over |q| < Q cut . In such an approximation, the collision integral reduces to an effective one-to-twobody integral, which is proportional to the soft transport parameterq S in equation 2.8, where P b cd (x) is the vacuum splitting function of parton b to parton c and d and we does not write down the explicit temperature dependence of the jet transport parameters.

Radiations in the deep-LPM region
The diffusion-and collision-induced radiations discussed above only apply to those parton splittings whose formation time (τ f ) is short compared to the collisional mean-free-path (λ el ). However, collinear splittings in a dense medium often have τ f λ el in the so-called deep-LPM region. We have to correct the above formulae to apply the transport equation to these splittings in the deep-LPM region. In reference [50], we developed a simple correction method inspired by the next-to-leading-log induced radiation rate in [72]. First, one samples the incoherent radiation collision integral in equations 2.13 and 2.14 at time t 1 . Momentum broadening to the radiated partons and the formation time are determined self-consistently in the simulation, where k ⊥,t 2 is the broadened transverse momentum at time t 2 = t 1 +τ f . Then, the radiation probability is suppressed by the factor .
is the effective mean-free-path. Q max is an estimation of the maximum momentum transfer in an elastic collision, which is taken to be the average center-of-mass energy between jet partons and medium partons Q max = √ 6ET . This method has been tested in [50] and we refer to the detailed discussion there. It was demonstrated to agree well with theoretical calculations in the deep-LPM region for induced radiations in an infinite static medium [72]. It also qualitatively captures the path-length [73] and expansion rate dependence [74] of single radiation probability in a finite and expanding medium. To better demons rate that this model achieves quantitative agreement with our theoretical knowledge, we have included a detailed comparison between the theoretical single gluon emission rate and the LIDO model simulations in a finite and expanding medium in appendix B.
To demonstrate the importance of the LPM effect and provide a benchmark result for this transport model, we plot the simulated single-parton energy loss in 0-10% Pb-Pb collisions in figure 1. Initial partons with energy E at midrapidity with randomized azimuthal orientation. The transverse positions are sampled according to the transverse density of binary nucleon collisions. The energy loss is defined by the difference between the initial-state and final-state energy of the parton. We show the elastic (lines) and radiative (bands) energy loss of gluons (red solid), light quarks (green dashed), and heavy quarks including charm (blue dash-dotted) and bottom (black dotted). The parton energy loss has been rescaled by the color factor C R so the results are comparable for quarks and gluons. In a finite and expanding medium, the elastic energy loss is a slowly varying function of parton energy and the radiative energy loss grows as ln E at high energy. Elastic parton energy loss dominates at low energy E < 10-20 GeV, while radiative energy loss dominates at high energy. Moreover, the radiative parton energy loss of quark shows a clear mass dependence ∆E q > ∆E c > ∆E b . Finally, we also include a calculation for light quark energy loss without the LPM effect (green dashed line) for comparison. Without the LPM effect, radiative energy loss would increase faster than ln(E) at high energy.

Medium excitations induced by hard partons: an approximate treatment
The linearized partonic transport model only handles hard particles with p · u > 4T . The energy-momentum transfer to the soft sector (p · u < 4T ) is considered in this section, and its impact on the final observables-transverse-momentum density dp T /dη/dφ-is modeled in a hydrodynamic-motivated ansatz. This simple model consists of three parts: instant thermalization of energy-momentum transfer, angular redistribution of energy-momentum, and freeze-out.
Instant thermalization We assume that, at each instant, the medium only responses to the four-momentum exchanges with jet partons 4 . Within each time step ∆τ and for each collision labeled by "i" at coordinate (τ i , x ⊥,i , η s,i ), we compute the differences of the total four-momentum of hard particles in the initial state (IS) and in the final state (FS), This gives the four-momentum deposited to the soft medium.
Angular redistribution of the energy-momentum perturbation If ∆p µ i only causes a small perturbation to the hydrodynamic evolution of the bulk medium, we can compute the linearized response of this energy-momentum deposition in energy-momentum density. However, even the linearized response on top of an arbitrary hydrodynamic background is highly non-trivial. We reduce the problem to a manageable level using the following simplifications: 1. Assume the typical temporal and spatial length scales of the perturbation are much shorter than those of space-time evolution of the background, including that of the transverse flow velocity, are neglected 2. The perturbations are localized in the space-time rapidity.

Assume a constant speed-of-sound c s .
Under such assumptions, energy-momentum perturbations e and g satisfy simple linearized ideal hydrodynamics in a static medium in the Bjorken frame. Take one of the sources in equation 2.17 and make the space-time Fourier transformation from τ, x to ω, k. The resultant equations are 5 , where k andk are the magnitude and the unit direction vector of k and ∆p 0 and ∆p i are the energy-momentum deposition. The solutions are e = 1 ∆t (2.21) Now, we can avoid the transformation back to space-time coordinates by focusing on the angular distribution of energy (dG 0 /dΩ) and momentum (dG/dΩ)that are observed at a point far from the source. Such distributions share the same angular structure as the numerators of equations 2.20 and 2.21; therefore, where G µ = drg µ . One can directly check the conservation of energy and momentum by integrating 2.22 overk and verify that dG 0 = ∆p 0 and dG = ∆p.
The angular structure in equation 2.22 is comprised of an isotropic term induced by energy deposition and a cos φ p,k modulated term induced by momentum deposition. Taking a light-like energy-momentum deposition ∆p 0 = |∆p| and choosing typical c 2 s ≈ 0.25, this ansatz propagates an excess of energy-momentum to the direction parallel to p, but it propagates a small energy-momentum deficit to the direction anti-parallel to p. This qualitative feature is similar to those observed in the coupled transport-hydrodynamic simulation [51].
Naïve freeze-out The formula in equation 2.22 distributes energy-momentum as a function of the spatial solid anglek. To calculate the contribution of perturbed transverse momentum density to jet energy as a function of azimuthal angle φ and pseudorapidity η, we convolve dG µ /dk with momentum-space information using a simple freeze-out procedure.
Though transverse flow is neglected in deriving spatial equation 2.22, flow effects are essential for particle distributions in the momentum space. We assume, in central nuclear collisions, an averaged fluid velocity of v r = 0.6c at freeze-out and approximate the direction of the flow byk, so where δT = c 2 s δe/w and δu µ = δg ν · (g µν − u µ u ν )/w, where w is the specific enthalpy of the background medium. We approximate the system at freeze-out by a non-interacting gas of massless particles so that w = 4T 4 /π 2 per degree-of-freedom. Putting these pieces together and expanding the result to the first order in δT and δu, we arrive at the perturbation to the transverse momentum density . We remark that details of the medium, such as the number of species and degeneracy factors, cancel in the final expression. This is because we are computing the transverse-momentum density, which is less sensitive to the composition of the gas. On the contrary, corrections to multiplicity distribution are sensitive to the chemical composition of the medium. Nevertheless, we will use this model to estimate corrections to charged-particle multiplicity in section 5.1 when computing jet fragmentation function in nuclear collisions.

Evolution outside of the medium, hadronization, and jet defintion
We assume the jet-medium interaction should become negligible when temperature drops to T f that is around or below the pseudo-critical temperature T c = 154 MeV [76]. We call T f the termination temperature of jet-medium interaction. For partons entering regions with T < T f , they stop interaction with the medium according to the transport equation and we again use Pyhtia8 to evolve the parton shower in the vacuum. The parton shower hadronizes through the Lund string fragmentation mechanism and is followed by hadronic decays. The Lund string model requires the system consists of only a color singlet string of partons. Colors of hard partons in the transport equation are initialized by the color generated by Pythia8 and the system, including beam remnants, is color neutral. However, after the interactions with the QGP medium, the hard parton shower is not color-neutral in general. Currently, it is intractable to trace the color exchange and conservation between hard partons and the entire medium. This is partly because the derivation of the in-medium scattering cross-sections already takes the average over colors of both hard partons and the medium. More importantly, even if we can keep the color information in the matrix-element of hard-soft collisions, it is impossible to model the transport of the soft carries of color charges in a strongly-coupled QGP. In the present study, we use the following prescription to model color exchange between hard parton and the medium and impose color neutrality before applying the Lund string fragmentation model. During the evolution, we assume the color of a hard parton is not altered in the small-angle diffusion process but is changed instantaneously in large-angle elastic collisions with a randomly sampled color exchange with the medium. At the end of the transport evolution in QGP, we connect all possible final-state hard partons into strings. The endpoints of some strings may carry net color charges due to interactions with the medium. Then, thermal quarks or anti-quarks with the desired color or anti-color are sampled and attached to these endpoints so that the entire string is color neutral. This procedure is similar to the one used in a study with the MARTINI model [45]. To guarantee energy-momentum conservation, the four-momentum of the sampled thermal partons will be treated as a negative source term to the medium using the ansatz for the hydrodynamic-like medium response that is introduced in section 2.3. Finally, each hard parton is assigned a virtually using the momentum broadening accumulated during the formation time of the latest medium-induced radiation k ⊥,final .
Pythia8 generates vacuum-like emissions out of the hot QGP medium from k ⊥,final to the non-perturbative scale and performs hadronization.
We define jet at the hadronic level. In each event, a grid of energy towers E T,ij is defined using the hadronic final state from hard partons and the contribution from linearized hydrodynamic response to transverse-momentum distribution, Jets are defined from these "energy towers" using the anti-k T jet finding algorithm implemented in the FastJet package [77,78]. Because contributions to E T,ij from the linearized hydrodynamic response can be either positive or negative, in the first round of jet finding, we use only positive towers (E T,ij > 0) to define jets. Then, the jet momentum is recalculated using all E T,ij contributions within the jet cone. This method implicitly subtracts the unperturbed medium as the background from the jet energy. Finally, Grid sizes ∆φ and ∆η should be small compared to jet radii, but not too fine that increases jet clustering complexity significantly. In practice, we used a 300 by 300 grid within the acceptance −3 < y < 3, −π < φ < π with resolution ∆y = 0.020 and ∆φ = 0.021, respectively.

Model calibration
The model introduced in the last section has multiple parameters and its predictive power relies on the ability to constrain these parameters from limited data. In this section, we calibrate a subset of model parameters on the suppression of inclusive jet and hadron spectra data measured at the RHIC and the LHC using Bayesian inference techniques [79]. We will show that with reasonable choices of parameters, the model consistently describes the inclusive jet and high-p T hadron suppression within ±10%-20% uncertainty. More importantly, we report sets of representative parameters for predictions with uncertainty quantification in the next section.

Summary of parameters
As the first work applying this model to jet quenching phenomenology, we only vary a subset of all possible parameters in the model. In particular, the jet-medium interaction is completely determined by the perturbative jet-medium coupling at scale µ min . Next, we list the parameters that are tuned in the calibration with ranges of variation and explain why we fixed other parameters in this analysis.
1. Parameter µ min in the jet-medium coupling, chosen to be proportional to the local QGP temperature. It serves as a cut-off for the running scale at finite temperature α s (Q) = α s (max{Q, µ min }). The actual parameter varied is ln(µ min /πT ), ranging from 0.7 to 4.0. The appearance of π is not necessary and one can simply translate it into 2.2T < µ min < 12.6T , which is around the typical energy of a thermal parton.
2. The scale parameter Q 0 at which one transits jet shower evolution from high-virtuality vacuum parton shower evolution to transport evolution in the QGP. Q 0 is assumed to only depend on the collision system, beam energy, and centrality. Practically, this is equivalent to two independent Q 0 parameters for central Au-Au collisions at the RHIC and central Pb-Pb collisions at the LHC respectively: . A current limitation is that the Pythia8 parameter that controls Q 0 is only allowed to be varied from 0.4 GeV to 2.0 GeV, which precludes the exploration of higher Q 0 values. We will consider other event generators in the future that allows a large range of variation.
3. The termination temperature T f of QGP effects. Near the pseudo-critical temperature, the QCD equation of state deviates significantly from the Stefan-Boltzmann limit. To account for this gradual onset of confinement, we stop the transport evolution in hot QGP at T f that can be varied from 0.15 GeV to 0.17 GeV. Because the temperature drops slowly at late times of the hydrodynamic expansion, a 0.02 GeV change in T f causes a large variation in the time duration that the system stays in the transport stage. This is demonstrated in appendix A.
It is necessary to explain why we do not vary the following parameters in this work.
1. We do not vary the Q cut parameter introduced in section 2.2. Because given the perturbative input to the diffusion coefficients and large-q collision matrix-elements, the Q cut dependence of elastic energy loss due to the diffusion process is canceled by that of the large-q collisions. This parameter is of interest if one intends to compare detailed effects due to diffusion dynamics versus large-q collisions dynamics. The default choice in this study is Q cut = 2m D .
2. We do not vary "hard" particle cut-off whose default value is set at 4T . This parameter does not affect the energy loss of hard partons and therefore, high-p T inclusive hadron R AA spectra. Also, because both elastic recoil and hydrodynamic response transport energy to large angles, the jet suppression with R = 0.4, for which we calibrated the model, is also insensitive to its choices.

Experimental dataset
We calibrate parameters to inclusive jet and hadron suppression in central nuclear collisions. Jet data include R = 0.4 jet R AA in 0-10% central Pb-Pb collisions at √ s N N = 5.02 TeV measured by the ATLAS Collaboration [80] and by the ALICE Collaboration [81], R = 0.4 charged jet R AA measured by the STAR Collaboration [82] in 0-10% Au-Au collisions at Figure 2. Emulated calculations of inclusive jet and single hadron suppression after model parameters are systematically calibrated to these observables. In each set of plots, the R AA is shown in the upper panel, and the ratio to experiment is shown in the lower panel. The blue bands correspond to 95% credible limit given by the model emulator. The gray bands in the ratio plot denote ±10% and ±20% level of discrepancy respectively.
√ s N N = 200 GeV. There is also a recent CMS measurement at this beam energy [83], but we do not include it in the calibration because its p T coverage is similar to those highest p jet T ATLAS data points. We will compare the prediction to the CMS measurement after the calibration in section 5.2, especially focusing on the jet cone-size dependence. The p jet T coverage at LHC ranges from 60 GeV to 1 TeV, while STAR data at RHIC extend the reach of p chg jet T down to 15 GeV/c 6 . Data on single inclusive hadron spectra include CMS measurement of charged particle R AA [84] and D 0 meson R AA [85] in 0-10% Pb-Pb 6 Note that the STAR Collaboration measurements charged jet, and uses a leading p ch T > 5 GeV trigger in nuclear collisions. The triggering effect can bias the low-pT jet suppression. The STAR Collaboration has estimated the triggering bias to be small for p chg jet earlier study of open heavy-flavor using the LIDO transport model Only data on hadron R AA at p T > 10 GeV/c are used in the calibration for two reasons. First, the soft-diffusion coefficient and the approximate scattering cross-sections for single-parton evolution are under better theoretical control in the high energy limit. Therefore, we limit our comparison to the single-hadron production at high p T . Second, the model only implements the Lund string hadronization but neglects the processes of parton recombination in the medium. The latter is important for hadron production at low and intermediate p T (p T 8 GeV/c ) [87]. The choice of p T > 10 GeV/c does have a certain degree of arbitrariness. In appendix C.2 we checked the effect using different p T cut of 5 GeVc and 15 GeV/c to see its impact on the parameter extraction.

Parameter calibration
In this section, we focus on the results of the model calibration and provide detailed procedures of the Bayes parameter calibration in appendix C. In figure 2, we demonstrate the level of agreement of the model (or model emulator) with data, after systematic parameter calibration. Within each of the four plots in figure 2, the upper panel shows the direct comparison to the data used for calibration and the lower panel shows the model-to-data ratios. At the LHC energy, the calibrated model describes jet R AA within ±10% and charged particle suppression at high p T within ±20%. The shape of the calibrated D-meson R AA is slightly different from the CMS measurement and is closer to the ALICE measurement [88], though the latter is not included in the calibration. At RHIC energy, the jet and pion R AA are consistent with data within experimental uncertainty. The jet R AA at RHIC in the considered p T range decreases with p T while the STAR measurement increases.
The extracted probability distribution of model parameters (posterior distribution) is shown in figure 3.3. The logarithmic of µ min /πT has a well-constrained distribution around 0.3 with 95% credible limits from -0.07 to 0.62. This range translates into a 95% credible limits of 2.9T < µ min < 5.8T . This is a reasonable range for the in-medium scale in the coupling constant at finite temperature, considering typical thermal momentum is about 3T .
The medium-scale parameter for parton shower at LHC Q LHC 0 is found to be large, with a 95% credible limits from 1.39 GeV to 1.99 GeV. At the moment, we are not sure whether a larger Q 0 will be preferred if it is allowed to vary above 2 GeV. This scale parameter at RHIC energy Q RHIC First, we find that values of ∆k T increase rather slowly with p T and occupy a relatively compact region as compared to the range of virtuality from p hard T to Q min = 0.5 GeV during the evolution. Both features are in favor of approximating the separation scale between vacuum and medium-induced shower by a single Q 0 value. At the LHC energy, typical values of ∆k ⊥ are slightly lower than Q LHC 0 , but still around the same order of magnitude for partons with p T > 10 GeV/c . At RHIC energy, the large uncertainty of Q 0 covers the same range of ∆k ⊥ distribution.
Finally, the termination temperature of jet-medium interaction T f prefers large values towards 0.17 GeV, though it is in general not well constrained. Note that T f strongly correlates with the in-medium scale parameter in the coupling constant as shown in the bottom-left corner plot: a largerµ min -smaller coupling strength-is correlated with a smaller termination temperature, which means the jet-medium interaction lasts longer to compensate for the reduced coupling in the calibration.

Parameter sets for central prediction and error estimation
In principle, to make predictions with quantified uncertainty from the calibrated model, we should pull random parameter samples from the posterior distribution and make an ensemble of predictions using the sampled parameters; then, central values and uncertainties can be defined accordingly from the ensemble prediction. However, it is unnecessary and impractical to propagate arbitrary many parameter samples to predictions. Instead, we define a small number of representative parameter sets to characterize parameter variations in the high-likelihood region of the posterior, including a central parameter set and a few error estimation parameters sets for the prediction. There is not a unique way to define the representative parameter sets. The method adopted here relies on the posterior showing a small degree of non-linear correlations among different parameters. Neglecting the non-linear correlations, we define linear combinations of parameters q i = c ij p j so that new parameters q i are linearly independent from each other. Then in the q-space, the central parameter set is defined to be the collection of medians of parameters, q central = {median(q 0 ), · · · , median(q 3 )}, (3.1) and error sets are obtained by changing one of the four values in the central set to the 2.5% and 97.5% quantile numbers of that parameter, q error,k± = {median(q 0 ), · · · , quantile 97.5% 2.5% (q k ), · · · median(q 3 )}, k = 0, 1, 2, 3. (3.2) This introduces eight error parameter sets in total. Finally, q central and q error,k± are transformed into the original model parameter space p with values listed in table 1.
As a validation, we compute the single inclusive jet and hadron suppression using the representative inputs in table 1 and show the results in figure 5. Uncertainty bands are drawn as the envelope (blue bands) of these nine different calculations. The uncertainty is at about ±15% level, compatible with the emulator predictions that use random draws of samples from the full posterior distribution. These representative parameter sets will be used for predictions and analysis for jet suppression and modification.

Energy flow induced by jet-medium interaction
Before using the calibrated model to predict jet modifications in section 5, we would like to understand how energy-momentum is transported within the parton shower. We studied a simple case-back-to-back quark or gluon pairs produced with fixed p hard T , using the central parameter set in table 1. We focus on the energy distribution E(z H , r) and E(p T , r) of partons in the two-dimensional spaces (z H , r) and (p T , r), where z H is the longitudinal momentum fraction relative scale of the hard process p hard T and r is the radial distance of a particle relative to initial hard parton, or in transverse momentum p T and r, The subscript "H" reminds us that the reference scale of z H is the momentum of the initial hard process. The definition of z H if different from that of z jet used for inclusive jet fragmentation function measurement, whose reference scale is the transverse momentum of the reconstructed jet. We will comment on the difference between z H and z jet at the end of this section and emphasize again in section 5.1. This simple example elucidates the respective roles of elastic collisions, induced radiation, and medium excitations in modifying E(z H , r) and E(p T , r), which helps to interpret the prediction in section 5. We initialize pairs of back-to-back quarks or gluons with initial hard momentum p hard T at mid-rapidity. Following the steps described in section 2, we use Pythia8 to generate vacuum parton showers and initialize the transport equation in the medium at scale Q 0 . Toward the end of transport in the medium at T f , the jet shower continues to evolve in vacuum down to Q 0 = 0.5 GeV. We omit hadronization in this simple example.

Energy transport via vacuum and induced radiations
We first focus on energy transport induced by radiative process and momentum broadening due to elastic collisions, while elastic energy loss and medium excitation are turned off 7 . Energy cascades due to radiative processes in the vacuum and in a static medium were studied in [89][90][91][92][93] under certain approximations. Here is one of its key observations. Compared to the vacuum-like jet shower governed by DGLAP-type evolution with vacuum splitting kernels P (0) (x), parton splittings in the medium has an enhanced soft splitting rate because its typical inverse formation time increases for softer parton energy. For example, in a static medium, dR/dx ∼ α s P (0) (x)τ −1 f ∼ α s P (0) (x) q/xE, where x is the momentum fraction of the softer daughter parton relative to the momentum of the mother parton. Such behavior builds up a power-law like energy spectrum dE/dz H ∼ z −1/2 H for z H 1, which causes a finite amount of energy to flow into z H = 0 within any arbitrary time duration [93]. In practice, the "point-like" energy sink at z H = 0 is replaced by a finite region z z BH . z BH is the energy fraction corresponding to the Bethe-Heitler energy scale z BH = ω BH /E with ω BH ≡ 1 2q λ 2 el . The flow of energy ceases below z BH , because the Bethe-Heitler radiation rate dR BH /dx ∼ α s P (0) (x)λ −1 el lacks additional enhancement at small-x. While splittings from higher z H regions continue to pump energy into this region, energy starts to accumulate among particles with typical energy ω BH .
Though some detailed assumptions employed in [93] become invalid in a finite and expanding plasma, one still preserves the most important qualitative feature that the inverse formation time of induced radiations increases for soft splittings. One also notices that as a consequence of the formation of the thermal medium such thatq ∼ α 2 s T 3 and λ el ∼ 1/(α s T ), the Bethe-Heitler energy ω BH = O(T ) coincides with the thermal scale 8 . This suggests that induced-radiation processes alone can build up energy storage near the thermal scale p ∼ T , no matter how energetic the initial hard jets are. Such behavior has also been confirmed in recent numerical simulations including the effective kinetic theory with resummed radiation 7 This is achieved by rescaling the magnitude of a particle's momentum to its value before the elastic collision or a diffusion step. 8 Neglecting the logarithmic factor, the jet transport parameter for a gluon isq ≈ αsCAT m 2 D . Effective mean-free-path has been defined as m 2 D /q. These lead to ωBH ≈ πT .  vertex [94] and the LIDO transport model in the large-medium L med τ f limit, which is discussed in section 3.5.3 of [95]. We will see immediately that this feature indeed emerges from our simulation.
Because ω BH coincides with the thermal scale, the distribution of the accumulated energy-momentum around ω BH will be further modified by other transport mechanisms including elastic collisions and jet-induced medium responses that are important for particles with p ∼ T . This will be discussed in section 4.2.
In figure 6, we plot the energy distribution E(z H , r) from the vacuum-like parton shower evolution of a jet produced at p hard  Vacuum radiations transfer the energy of the initial hard quark/gluon into a broad range of z H and r, while most of the energy remains in the small-r and high-z H region. The distribution from an initial quark is harder and narrower than that from an initial gluon, as gluons radiate more frequently.
Next, we evolve the jet shower in the QGP medium with induced radiations and the elastic broadening, following it by vacuum radiations outside the QGP. The modified E(z H , r) for a initial gluon is shown on the left panel of figure 7. Compared to figure 6 (right), the energy distribution is softened and widened. More importantly, an "energy bump" build up in the large-r and low-z H (and small p T ) region. To confirm that this bump is indeed the consequence of radiative energy transfer to the thermal scale, we show the p T -dependent energy distributions E(p T , r) with different p hard T : 100, 400, and 1600 GeV/c on the right panel of figure 7. Evidently, changing p hard T by an order of magnitude does not affect the typical momentum scale p T ∼ 2 GeV/c at which a bump develops. Therefore the bump indeed locates at a soft medium scale that is independent of the initial hard scale.

Energy transport with full effects: radiation, elastic energy loss, and medium excitation
We have seen in figure 1 that the elastic energy loss is almost independent of the momentum of a parton. In the meantime, the parton multiplicity at low-z H is much higher than that at high-z H during jet evolution in the QGP. Combining these two arguments, elastic energy loss is most important for low-z H region while becomes insignificant at high-z H . The lost energy-momentum due to elastic scatterings are then redistributed by either hydrodynamic responses or hard medium recoiled partons. The former induces near-equilibrium perturbation that modifies particle distribution around p T ∼ T at large angles. Momentum of the recoil particles in a t-channel process range from gT to the average center-of-mass energy √ 6ET . Therefore, the inclusion of elastic energy loss plus medium excitation mainly affects particles within O(T ) < p T < O( √ ET ). On the left of figure 8, we show the resulting E(z H , r) with full medium effects: inducedradiations, elastic broadening, and energy loss, and the medium excitation. Comparing to the left of figure 7, elastic energy loss plus medium excitation deform the low-z H bump into a sharper peak around p T ∼ 1-2 GeV/c (z H = 0.005-0.01).
On the right of figure 8, we show the ratio of the modified energy distribution to the theoretical baseline E vac (z H , r) due to parton shower in vacuum down to Q 0 = 0.5 GeV. There is an energy deficiency at intermediate z H compared to the theoretical baseline in the vacuum, and energy excess develops at the small z H with soft momentum. This qualitatively agrees with the finding of authors of [18,52], where both vacuum and BDMPS-Z type medium-induced emissions are included. However, we emphasize that in our simulation the energy excess is a combined result of both medium-induced radiation and jet-induced medium excitations: induced radiation builds up the energy excess at the soft-scale and jet-induced medium excitation reshapes the bump into a sharper peak closer to the thermal momentum.
When z H approaches one, the ratio can be larger than unity. This is a consequence of the way we interface the vacuum-like evolution equation and the transport equation and we explain as follows. In our calculation for jets in medium, the scale evolution from Q 0 to k ⊥,final , where k ⊥,final is the momentum broadening accumulated during the formation of the last medium-induced radiation, is replaced by the evolution in QGP using the transport equation. In figure 1, we have already seen that the high-energy parton energy loss in the medium grows logarithmically with parton energy: ∆E med ∼ ln E. In simulations for jets the vacuum, soft vacuum radiation between Q 0 and k ⊥,final causes energy shifts of the leading parton ∆E vac ( Q 0 k ⊥,final ). Because ln Q 0 k ⊥,final ln Q 0 Q min ≈ 1.2 is not a large number, we can estimate ∆E vac as αs . As a result, at large energy, the ∆E vac that we removed from scale evolution is larger than the energy loss in the medium caused by jet-medium interactions, which explains why the ratio is larger than unity for z H ∼ 1 for this simple test where z H = p parton T /p hard T . This phenomenon is also discussed in [18] However, such an excess is only visible when we focus on very energetic particles so that ∆E vac > ∆E med and analyze the result using the variable z H = p parton T /p hard T . In fact, the peak near z H = 1 can be smaller than unity for a lower-p hard T or be smeared out by hadronization effects. Furthermore, in inclusive measurement, the momentum fraction is defined using z jet = p T /p jet T . Unlike using p hard T , jets in a sample with a fixed p jet T in nuclear collisions starts are produced with a higher scale and favors jets with harder constituents due to jet energy loss. Therefore, the jet fragmentation function ratio between AA and pp collisions being large than one at large z jet is mainly the result of the trigger bias [18,96]. To remove the trigger bias and test this effect due to a switching scale between vacuum-like and medium-induced radiations, one would like to use γ-tagged and Z-tagged jets samples [19]. Using z γ/Z = p T /p γ T or p T /p Z T , one should focus on the ratio of fragmentation function large-z γ/Z and its dependence on the transverse momentum of the trigger-boson. It would be interesting to see how this effect manifests in the fragmentation functions of γ-jet and Z-jet. It can also help to test the current modeling of the switching between vacuum-like evolution and in-medium evolution. This will be pursued in future studies.

Remarks on trigger bias
The energy flow picture discussed in this section helps understand the theoretical origin of certain features that we are going to see in data. But experimentally, it is impossible to determine the exact initial hard scale p hard T for jet production in both nuclear collisions and proton-proton collisions. Experiments usually compare single inclusive jets in nuclear and p-p collisions at the same final p jet T . One has to be careful when interpreting the difference between jet properties in nuclear and the p-p baseline. For example, jets in central Pb-Pb collisions are produced at a higher scale on average than the reference jets in p-p collisions due to energy loss, which means ratios deviating from unity are not always direct consequent of energy flow.

Prediction of jet modifications
In this section, we predict modifications to jet fragmentation function and jet shape, as well as the cone-size dependence of jet R AA . These predictions use full-model calculation with representative parameter sets from the calibration in section 3.3. Jets are defined at the hadron level with the same kinematics cuts in the measurements.

Modified jet fragmentation function in central Pb-Pb collisions
For the ATLAS measurement [97], the jet fragmentation function is either defined as charged particle distributions in longitudinal momentum fraction z jet or in transverse momentum p T , normalized by the number of jets, Note that z jet uses the transverse momentum of the reconstructed jet as the reference scale, which is different from that for z H in section 4.
In figure 9, we compare the Pythia8 simulation as a baseline to the ATLAS measurement of D(z) in p-p collisions [97] for 126 < p jet T < 158 GeV/c , 200 < p jet T < 251 GeV/c , and 316 < p jet T < 398 GeV/c . There is a reasonable agreement for a wide range of z jet except in the very large-z jet bins, where the simulation drops faster than data toward z jet = 1.
Our predictions for ratios between jet fragmentation function D(z) in Pb-Pb and p-p collisions are shown in figure 10. The predictions include bands of uncertainty using the representative parameter sets in table 1. The model constrained by inclusive jet and hadron suppression consistently describes the modification of jet fragmentation function observed Figure 9. Comparison of charged particle fragmentation function of inclusive jets obtained by Pythia8 simulation to data from the ATLAS experiment [97]. Ratio of experiment to simulation is shown in the lower panel. Figure 10. Prediction of modification to the jet fragmentation function (D(z)) and comparison to ATLAS measurements for three different p jet T bins [97]. Bands corresponds to the envelop of calculations using the parameter sets listed in table 1. Dashed lines are calculations using the central parameter set but excluding contributions from medium excitations when calculating fragmentation functions. in experimental data. To illustrate the respective roles of radiative processes and medium excitations, we compare the full model calculation to the calculation that excludes contributions from medium excitations (dashed lines). We emphasize that jets are still defined with contributions from medium excitations so that the jet samples are unbiased. It is clear that even without contributions from medium excitations, low-z enhancements of the ratio D AA (p)/D pp (z) are already sizeable. In the previous section, we have discussed that the medium-induced radiative process effectively transfers energy towards small-z and causes the accumulation of energy at thermal scales. Because the overall impact of jet-medium interactions is to excite medium partons to higher energy, contributions from medium excitations are negative and become positive at high transverse momenta. As a result, medium excitations develop a notable peak in the modification of fragmentation functions around p T ∼ 0.7 GeV for all three p jet T bins. Modifications of D(z jet ) are known to demonstrate different scaling behaviors in the hard and the soft regions [97]: ratios for D(z jet ) at high-z jet for different p jet T collapse on the same curve; while ratios for D(p T ) at low-p T are approximately independent of the p jet T . In figure 11, the central predictions for D(z) (left) and D(p T ) (right) ratios for three different jet p T bins are compared with experimental data. The high-z jet scaling of D(z) is well reproduced. In the model, the energies of high-z jet particles are sufficiently higher than other scales such as temperature and elastic recoil energy from O(T ) to O( √ ET ); therefore, it is natural that its modification, dominated by induced radiation, is solely a function of z jet . For the modifications of D(p T ), the model captures the low-p T enhancement of D(p T ) at p T 4 GeV/c . From discussions in section 4, we see that the origin of low-p T scaling is Figure 12. Predictions of cone-size dependence of inclusive jet nuclear modification factor in Pb-Pb collisions at √ s N N = 5.02 TeV, as compared to data from ALICE [81] and ATLAS [80] experiments, and preliminary data from the CMS experiment [83].
a combination of both induced-radiation and medium excitation: energy first builds up at soft-momentum scale due to induced radiation, whose distribution is further modified by the jet-induced medium response. Comparing the magnitude of the enhancement for different p jet T , the p jet T dependence in the simulation is stronger than that observed in experimental data. This could be caused by several reasons. For example, it is possible that the transition scale parameter Q 0 should be slightly higher for higher p jet T -triggered jets. This would decrease the multiplicity of the parton shower at the initialization of the transport equation and reduce the amount of energy transfer to lower energy partons.

Cone-size dependence of jet suppression and jet-shape modifications
We also study the cone-size dependence of jet R AA and jet shape modifications that are closely related to the radial transport of jet energy.
Cone-size dependence of jet suppression Jet reconstruction with a larger cone-size should recover more energy in the hard process. The cone-size dependence of jet R AA comes from the competition of two effects: the increase of jet yield with R in nuclear collisions and that in p-p collisions. Recent measurements demonstrate that jet R AA has a weak R dependence. At the RHIC energy, the STAR experiment observes that the charged jet R AA Figure 13. Predictions of cone-size dependence of inclusive jet nuclear modification factor in Au-Au collisions at √ s N N = 200 GeV, compared to charged-jet measurement from the STAR experiment [82]. The experimental data colored in red (the three highest p chg jet T data points) are considered to be "unbiased", while the rest data points shown in back are subject to strong leading particle triggering bias.
for R = 0.2, 0.3, and 0.4 are the same within the experimental uncertainty [82]. At the LHC energy, the ALICE measurement [81] shows similar R AA at R = 0.2 and R = 0.4. the CMS measurement [83] has extended jet measurement to R = 1.0 for p jet T > 400 GeV/c , and the R = 1.0 R AA shows less than 10% increases compared to R = 0.2 results. This suggests energy lost at the core of the jet is transported mostly to very large angles.
In figures 12 and 13, the model calculations are compared to both the RHIC and the LHC data on jet R AA with different jet cone-sizes. We notice that the CMS jet R AA at R = 0.4 is about 0.15 higher than the ATLAS measurement in the similar p T and rapidity window, though the discrepancy is only at 1-2 σ level. Because we calibrate the model to ATLAS data at R = 0.4, the predictions are systematically lower than the CMS measurements. The model, however, explains the weak cone-size dependence at both the LHC and RHIC. A weak cone-size dependence is also found in the strongly-coupled hybrid approach to jet modification [17]. Another weakly-coupled transport equation calculation predicts a stronger cone-size dependence [49].
The STAR experiment uses a leading particle trigger of p trigger T > 5 GeV to define jets in Au-Au collisions to suppress "fake jets" from the background. This trigger inevitably biases the reconstructed inclusive jet spectra. Especially, the jet yield vanishes below the momentum of the trigger particle p chg jet T < p trigger T . Jet with higher transverse momentum is less affected by the trigger and the STAR experiment has estimated that the measured R AA with p chg jet T > 16 GeV, denoted by red symbols, are effectively free from the trigger bias. For this reason, we only include the "unbiased" data points in the calibration in section C. Nevertheless, we can estimate the effect of trigger bias by comparing the inclusive jet R AA calculated without triggering (open solid box) to that with leading-particle triggering (gray bands). With the triggering, the charged-jet R AA is almost flat as a function of p chg jet T . Figure 14. Pythia8 simulations of inclusive jet shape in proton-proton collision (red-solid line), compared to data from the CMS experiment [98]. The red dashed line corresponds to a "corrected" simulation to better agree with the CMS data at large r. The "corrected" results has been subtracted by a 2.2 GeV per radian per rapidity background.
Our calculation with and without triggering agrees well above p chg jet T ∼ 15-20 GeV and diverges for lower momentum jet, which is in agreement with the estimation from the STAR experiment that R AA values above 16 GeV are effectively unbiased. The calculation with trigger explains the dropping of the charged-jet R AA below 20 GeV. The calculation also shows that jets with smaller cone-size are less sensitive to the trigger bias because hard particles satisfying the momentum cut are more likely to be contained in jets with smaller cone-size at fixed p chg jet T . Jet shape Jet shape directly measures the radial distribution of energy around the jet. Following the CMS measurement [98], we define the radial distribution of the transverse momentum of charged particles associated with R = 0.4 jets P (r) = 1 N jets jets |r −r|<∆r/2 p T ∆r .
The nested summation goes over charged particles for each jet. Then, the jet shape is defined as a self-normalized transverse momentum distribution within r < 1.
Note that there is a small but finite kinematic cut p track T > 0.7 GeV/c on charged particles in the CMS measurement.
In figure 14, the p T profile P (r) in p-p collisions simulated by Pythia8 is compared to CMS data for p jet T > 120 GeV/c . Inside the jet cone of r < 0.4, the simulation (red solid line) describes the data well. For r > 0.4, Pythia8 simulations overshoot the data points, with a discrepancy about 50% at r = 1.0. Such a deviation in the baseline can affect the interpretation of the ratio between jet-shape in Pb-Pb and p-p. For a fairer comparison in Pb-Pb, we try to use an ad hoc method to correct this discrepancy at large r. Because the measured transverse momentum density is already small at large r: P (r)/(2πr) ≈ 0.55 GeV per squared radian, we attribute the discrepancy at large r to a small mismatch in background estimation between experiment and simulation 10 , P calc. (r) − P exp. (r) = 2πr 0 . The fitted 0 is around 0.37 GeV per squared radian. We then correct the simulated jet shape by subtracting this constant background. The corrected Pythia8 result is plotted as the red dashed line in figure 14, which agrees better with the measurement at a larger r. Note that 0 ≈ 0.37 GeV hardly affects the jet energy determination as 0 πR 2 ≈ 0.2 GeV is much less than the jet energy and energy loss.
In figure 15, the ratio between jet shape ρ(r) in Pb-Pb and that in p-p collisions is compared to CMS data. The uncorrected prediction (blue band) describes well the dip of the ratio around r = 0.1 and the fast-rising with increasing r. At large radius, the prediction starts to deviate from the measurement where ρ AA /ρ pp ≈ 1.7. As we have assumed, this stems from a mismatch in the background definition. We apply the ad hoc correction to the prediction where the estimated mismatch in background 0 is subtracted from both the p-p and Pb-Pb calculations, and N AA and N pp are the self-normalization factor of the corrected momentum profile P (r) − 2πr 0 for Pb-Pb and p-p collisions respectively. The correction increases the calculated extended jet shape for r > 0.4 to the measured level of "energy excess", shown as the open blue bands. Finally, the dashed line shows the jet shape calculation where we exclude the contribution of medium excitations. Comparing to the full model calculation, we conclude that medium excitation is negligible at small r. At r ∼ 1, its contribution is essential to explain the large-r excess of jet shape modification.
Modifications of jet fragmentation function and jet shape only contain information of the one-dimensional projection of the two-dimensional energy distributions dE/(dp T dr) and dE(dz jet dr), and a direct comparison to the latter impose unprecedentedly detailed constraints to the modeling of jet-medium interaction. For this purpose, the CMS experiment has measured the jet shape with different cuts on the transverse momentum of charged particles.
In figure 16, we show the ratios of the transverse momentum profile P (r) for charged particles with p T > 0.7, 2.0, and 4.0 GeV/c between Pb-Pb and p-p collisions (filled bands) and compare to data. The open bands are the "corrected" ratio as calculated in equation 5.5 where a constant transverse momentum background is subtracted so that the radial profile of transverse momentum density in proton-proton collisions agree with the data outside of the jet cone. The residual background levels that we subtracted from the transversemomentum distribution are 0 = 0.37, 0.29, and 0.12 GeV for p T cuts equals 0.7, 2.0, and  4.0 GeV/c , respectively. After the correction, we see a better description of the extended jet shape at large r. Outside the R = 0.4, increasing the p T cut from 0.7 GeV/c to 2 GeV/c reduces the transverse-momentum outside the jet cone by about 50%. This can be understood from the simple test shown on the left of figure 8 that the large-r enhancement is dominated by soft p T ∼ 2 GeV/c particles and an increased p T cut effectively removes such contributions. This effect is more prominent with p T cut increased to 4 GeV/c .

Summary and outlook
In this work, we have performed a transport-model based study of jet production and modifications in relativistic heavy-ion collisions. Initial hard processes and vacuum-like parton showers are generated by Pythia8, jet-medium interactions are treated in a previously developed linearized partonic transport model LIDO, and the jet-induced medium excitation contains both hard parton recoils and a linearized hydrodynamic response.
Using data in single inclusive jet and single hadron suppression in central nuclear collisions at the LHC and RHIC, we systematically calibrated four model parameters: scale parameter µ min for the in-medium coupling, separation scale parameter Q 0 between vacuumlike and in-medium parton shower (separately tuned at the RHIC and the LHC), and a termination temperature T f for jet-medium interactions. The calibrated model consistently describes inclusive jet and hadron suppression at high-p T . Parameter sets for central prediction and error estimation are reported.
The calibration also extracts of the effective medium-coupling strength g med = g(µ min ) (left) and the jet transport parameterq(p, T ) (right) from both jet and hadron suppression data as shown in figure 6. The large calibrated values of the effective coupling at the medium scale g med ∼ 2-3 stands out as a potential problem because it challenges the use of weak-coupling assumption for jet-medium interaction. We will revisit this problem in future studies with non-perturbative modeling of the medium constituents. The extracted 95% credible region of the scaled jet transport parameterq(p, T )/T 3 of a light quark is consistent to earlier estimation by the JET Collaboration [6], but is larger than a recent extraction (preliminary) by the JETSCAPE Collaboration [60]. A possible explanation is that the JETSCAPE multi-stage approach includes medium effects in both the high-virtuality evolution and the transport equations of parton showers, but the present study accounts for all the jet-medium interactions in the transport equation. As a result, the present study may have over-estimated the magnitude of the jet-medium coupling to compensate for the neglected medium effects in the evolution of high-virtuality partons. In the future, we seek to include the LIDO model as an alternative transport model in the JETSCAPE framework to benefit from the more advanced implementation of high-virtuality parton evolution in the medium using the MATTER event generator.
In the second half of this paper, we first analyzed the energy flow induced by jet-medium interactions to help understand the observed jet modifications. We found that a mediuminduced radiative process transfers energy from the large-z jet partons down to partons with momentum p ∼ ω BH ≈ πT . Elastic energy loss and jet-induced medium excitation further redistribute the energy that has been accumulated around ω BH .
Using the calibrated model, we made predictions for modifications of jet fragmentation function, cone-size dependence of jet R AA , and modification of jet shape. The calibrated model produces a good agreement with the measured fragmentation function of inclusive Figure 17. The resultant 95% credible region of medium coupling strength g(µ min ) as a function of temperature T (left) and jet transport coefficient for light quarkq q scaled by T 3 (right). Red and green bands correspond to the value ofq/T 3 for a light quark with 10 GeV/c and 100 GeV/c momentum respectively. Black symbols with error bars are values estimated by the JET Collaboration [6] and the hatched band is one of the preliminary results extracted by the JETSCAPE Collaboration [60] . jets at the LHC in central Pb-Pb collisions. The excess of low-p T fragmentation function and jet shape at large-r are found to be consequences of both medium-induced radiation and jet-induced medium excitation.
In conclusion, the transport-based model can simultaneously explain inclusive jet and hadron suppression at high-p T . The same model with calibrated parameters yields comparable jet-transport coefficient to earlier extractions and provides consistent descriptions of inclusive jet modifications including jet fragmentation function and jet shape. These observations support the perturbative-based transport modeling of jet-medium interactions.

A Medium evolution
We simulate the medium evolution using a hydrodynamic-based hic-eventgen package that is described in detail in [99]. The package assumes boost-invariance and includes a parametric initial condition [100], a free-streaming model for pre-equilibrium dynamics [101], the 2+1D relativistic viscous hydrodynamics [102,103] 11 . For simplicity, we use event-averaged initial conditions to generate events for 0-10% central nuclear collisions. This is an acceptable approximation considering the observables that we consider, inclusive hadron and jet suppression, are most sensitive to long-wavelength structures of the medium.
We only extract the equilibrium information from the medium for the simulations of jet evolution. This means that during the hydrodynamic stage, jet partons interact with a medium in local equilibrium, which is comprised of massless quarks and gluons, whose energy density and flow velocity u µ match the viscous hydrodynamic simulation. The temperature is defined using the lattice equation of state T = T ( ) [76].
In figure 18, we plot the temperature profiles (left) and the evolution of the maximum temperature with proper time (right). From the right plot, we see that a 0.02 GeV decrease in T f leads a longer duration of jet-medium interaction: ∆τ ≈ 2.5 fm/c. Figure 19. Comparisons of the full theory calculations (dashed lines and symbols) of the single splitting rate of the process g → g + g to the LIDO simulations (bands). Left: differential rate dR/dω in a static medium as a function of patron evolution time t + , plotted at different energy of the daughter gluon: ω/T = 3, 6, 11, 20, 37, 68. Right: differential rate dR/dω in ax expanding medium as a function of patron evolution time t + , plotted at different energy of the daughter gluon: ω/T = 6, 20, 60. The medium expands according to Bjorken flow where T 3 ∝ τ −1 .

B Comparisons of single-gluon emission rates
We shall demonstrate in this appendix that the LIDO model is able to quantitatively reproduce the single parton emission rate as calculated in the full theory [22,23,73], where multiple collisions to the single parton emissions are resummed to all orders with full leadingorder weakly-coupled jet-medium interaction. In [50], we had compared our simulation to the published numerical results in reference [73] for a quark at relatively low energy E = 16 GeV. In this work, we use a new implementation [95] of the approach introduced in [73] to solve the single parton radiation rate from the full theory numerically. It allows us to compare our simulation to the theory calculation at higher parton energy that is relevant for jets at the LHC energies and also extend to the case of an expanding medium. In this appendix, we choose to demonstrate this comparison using the g → g + g splitting channel with the initial gluon energy E = 100 GeV. Both the simulation and the numerical solution of the theory use a fixed coupling constant α s = 0.3.
On the left panel of figure 19, we compare the differential rate (bands) dR/dω of gluon splitting g → g + g as a function of time at which the radiation takes place in a finite and static medium. omega is the energy of one of the radiated gluons. Bands and dashed lines correspond to LIDO simulations and theoretical calculations, respectively. Different colors label the differential rates at different gluon energies. We can see that simulations quantitatively agree with the theoretical calculations at late times. The model also captures the radiation rate in the transition region at early times, which is critical for the model to 11 It also includes a particlization model that samples hadrons from the fluid and a hadronic afterburner.
We only include jet-medium interaction in the QGP phase and neglects jet interaction with hadronic matter.
reproduce the non-linear path-length dependent radiative energy loss of leading partons. This agreement is better for energetic daughter parton with ω/T 1. At lower energy ω 3T , one observes notable deviations at early times. On the right panel of figure 19, we make the comparison between theory (symbols) and simulations (bands) in an expanding medium. The temperature of the medium drops according to that of the Bjorken flow [104], where T 3 ∝ τ −1 with an initial temperature T 0 = 0.5 GeV at the initial time of τ 0 = 0.5 fm/c. We observe that simulations capture the time dependence of the differential radiation rate very well for high energy daughter gluons. However, the current model produces more low-energy daughter gluons at late times than what is expected from the theory.
One concludes that in both cases, the agreement between theory and simulation is satisfactory for those emitting gluons with thermal energies ω ∼ 3T . We do not worry too much about this discrepancy at this moment for two reasons. First, for a moderate coupling α s = 0.3, the screening mass is of the same order as the kinetic energy of a gluon with thermal energy. The collinearity between the final-state and initial-state partons as the theoretical formula assumed breaks down. Therefore, the theory calculation already approaches its limits of applicability. Second, emissions of gluons with thermal momenta do not strongly affect the radiative energy loss of the leading parton, and their impact on soft partons are overwhelmed by those caused by elastic collisions. Therefore, despite the current discrepancy between our model and the full theory around ω ∼ 3T , we do not expect it will strongly affect the prediction on physical observables.

C.1 Procedures of Bayesian parameter calibration
We refer the readers to the reference [99] for a detailed review of the application of the Bayesian method to model calibration in the context of relativistic heavy-ion physics. Specific to our application, we perform the Bayesian model calibration in the following steps.
Model computation at the design parameter points Fifty parameter sets are sampled using the Latin-hypercube sampling procedure from the four-dimensional parameter space with the range introduced in section 3.1. The Latin-hypercube sampling algorithm is implemented in the R package [105,106]. Each of these fifty parameter points, p i , predicts a set of inclusive jet and hadron R AA through the full model simulation. The list of R AA for jet and hadron at different colliding energies is concatenated into a single training observable vector y i for the parameter point labeled by i.
Dimensional reduction of the training data We use principal component analysis, as implemented in the scikit-learn [107] package, to compress the information of the observable vectors y i into vectors of lower dimensional z i where dimz dimy. First, the reduction of the dimension of the observable is possible because many outputs vary in a highly correlated manner when parameters are varied. For example, in our study, an increase in coupling strength decreases all high-p T R AA . Second, the outputs are dominated by linear correlations. For example, the amount of change in jet R AA approximately scales linearly with the change in inclusive hadron R AA , when we vary the model parameters. In this case, the change in the full observable vectors can be approximately decomposed into linear combinations of just a few bases vectors b j labeled by j y i = Npc j=1 z ij b j + (Truncation errors). (C.2) Now, each row of z ij becomes the transformed observable vector for parameter set "i", and different columns of z are guaranteed to be linearly independent. When the dimension of z-N pc -is chosen to be smaller than dim y, the linear combination does not fully agree with the original vector and introduces finite truncation errors. For a well-behaved model with a few parameters, a small number of N pc is enough to reproduce the original vector with small truncation errors. We choose N pc = 5 in this work, which accounts for 97.5% of the variation of the original data. Variances in the truncated data are propagated to the interpolation uncertainty as white noises.
Emulate model prediction by interpolating training data. Next, we define an interpolation z = f (p) to obtain model prediction z at a general point p in the parameter space. The interpolation is realized by training the Gaussian processes on the fifty training (design) parameter sets and the corresponding predicted observable vectors. A Gaussianprocess emulator is a procedure to sample a random n-dimensional scalar function whose outputs at different inputs follow a multivariate Gaussian distribution. The training procedure adapts the hyper-parameters in the Gaussian-process emulator so that the distribution of the generated random functions is conditioned to the training output z i at training inputs p i . At a novel input point x i , the ensemble of random function values defines the central prediction and uncertainty of the interpolation. Finally, we transform the predicted z back to predictions of physical observables y. The prediction uncertainty that is represented by a diagonal covariance matrix in the principal component space is also transformed back to the physical space Σ GP . All these functionalities of building and training Gaussian process emulators are provided by the scikit-learn package [107].
Model discrepancy It is hardly true that a model, even with optimal parameters, can accurately describe every detail of an observable. This is the so-called model discrepancy.
Without accounting for the model discrepancy in the uncertainty estimation, a calibration potentially leads to over-confident constraints on parameters. In reality, for a complex model as in this work, it is impossible to quantify the discrepancy between the model and the true underlying theory. We adopt a strategy of including model discrepancy used in [99]. In this method, a relative uncertainty σ model is assigned to the prediction uncertainty of each principal component in addition to the uncertainty of the Gaussian-process emulators. Because we do not know the exact value of σ model , we treat it as another parameter to be marginalized in the Bayesian parameter calibration. Moreover, we give σ model an informative prior, This results in a 15% average model discrepancy prior to the model-to-data comparison.
Model-to-data comparison from Bayes' theorem Assisted by a fast model emulator, we can infer model calculations with uncertainty at any point in the considered range of parameters. Next, we apply Bayes' theorem to combine the model and experimental data to define the posterior probability distribution of parameters P (p).
where δy p = y(p) − y exp . Σ is the uncertainty matrix that combines interpolation uncertainty and experimental uncertainty. where Σ GP is the covariance matrix of the interpolation procedure. δy stat exp + δy sys exp is the experimental statistical uncertainty plus systematic uncertainty, which only contributes to the diagonal term of Σ. In reality, systematic uncertainties usually have a certain degree of correlation. For simplicity, they are treated as uncorrelated in this analysis except for the normalization uncertainty. The normalization uncertainty δy norm exp comes from binary collision T AA uncertainty and luminosity uncertainty in R AA measurement; therefore, they are fully-correlated uncertainty over different transverse momentum bins in one measurement. The ⊗ denotes the outer product of the uncertainty array for one of the experiments, and the ⊕ denotes the direct sum of normalization uncertainty covariance matrices among different experiments. The normalization uncertainty are included as percentage of the central reported experimental values δy norm exp = norm×δy exp . The normalization factor for different experiments in 0-10% Pb-Pb or 0-10% Au-Au collisions are listed in table 2. To obtain the maginalized posterior distribution that has been discussed in section 3.3, we use Markov chain Monte Carlo (MCMC), as implemented in the emcee package [108], to pull random samples from the posterior distribution defined in equations C.4 and C.5 and project onto the one-or two-parameter subspace.
ATLAS jet ALICE jet CMS chg. part. CMS D 0 STAR chg. jet PHENIX π 0 0.0085 0.028 0.025 0.025 0.07 0.12 Table 2. Normalization uncertainty level of different measurements in 0-10% Pb-Pb collision or 0-10% Au-Au collisions. GeV/c (blue dash-dotted lines). Including more low-p hadron T data slightly decreases the preferred in-medium coupling scale, which increases the coupling strength. It also drives the preferred termination temperature for jetmedium interaction toward lower values, increasing the effect of hot-medium on jet partons. Both of these changes try to enhance the medium modifications on a single particle. The matching scale parameter Q 0 at LHC energy favors higher values when we include more low-p hadron T data. A higher Q 0 effectively reduces vacuum showers, which compensates for the increased suppression on high-p T particles. The posterior of Q 0 parameters at RHIC is consistent when changing p hadron T due to large calibration uncertainty.
The posterior's dependence on p hadron T cut-off on the calibration data set suggests that our model is not entirely consistent with the p T dependence of inclusive hadron data. As a result, lowering p hadron T cut-off, the calibration has to balance the losses and gains with the low-p hadron T data. This is particular prominent for the T f parameter when we switch p hadron T cut-off from 10 GeV/c to 5 GeV/c . Despite this dependence on p hadron T , the three posterior distributions still largely overlap with each other. In particular, the 95% credible interval are consistent as listed in table 3.