Azimuthal jet flavor tomography with CUJET2.0 of nuclear collisions at RHIC and LHC

A perturbative QCD based jet tomographic Monte Carlo model, CUJET2.0, is presented to predict jet quenching observables in relativistic heavy ion collisions at RHIC/BNL and LHC/CERN energies. This model generalizes the DGLV theory of flavor dependent radiative energy loss by including multi-scale running strong coupling effects. It generalizes CUJET1.0 by computing jet path integrations though more realistic 2+1D transverse and longitudinally expanding viscous hydrodynamical fields contrained by fits to low $p_T$ flow data. The CUJET2.0 output depends on three control parameters, $(\alpha_{max},f_E,f_M)$, corresponding to an assumed upper bound on the vacuum running coupling in the infrared and two chromo-electric and magnetic QGP screening mass scales $(f_E \mu(T), f_M \mu(T))$ where $\mu(T)$ is the 1-loop Debye mass. We compare numerical results as a function of $\alpha_{max}$ for pure and deformed HTL dynamically enhanced scattering cases corresponding to $(f_E=1,2, f_M=0)$ to data of the nuclear modification factor, $R^f_{AA}(p_T,\phi; \sqrt{s}, b)$ for jet fragment flavors $f=\pi,D, B, e$ at $\sqrt{s}=0.2-2.76$ ATeV c.m. energies per nucleon pair and with impact parameter $b=2.4, 7.5$ fm. A $\chi^2$ analysis is presented and shows that $R^\pi_{AA}$ data from RHIC and LHC are consistent with CUJET2.0 at the $\chi^2/d.o.f<2$ level for $\alpha_{max}=0.23-0.30$. The corresponding $\hat{q}(E_{jet}, T)/T^3$ effective jet transport coefficient field of this model is computed to facilitate comparison to other jet tomographic models in the literature. The predicted elliptic asymmetry, $v_2(p_T;\sqrt{s},b)$ is, however, found to significantly underestimated relative to RHIC and LHC data. We find the $\chi^2_{v_2}$ analysis shows that $v_2$ is very sensitive to allowing even as little as 10\% variations of the path averaged $\alpha_{max}$ along in and out of reaction plane paths.

In the pQCD framework, radiative energy loss is assumed to be the dominant dynamical mechanism along with elastic energy loss for jet-medium interactions. In the past two decades, a wide range of jet quenching models have been formulated and applied to explain or predict high transverse momentum p T measurements at RHIC and LHC. These models are based on different medium property assumptions. BDMPS-Z [15][16][17] and ASW [18][19][20][21][22] multiple soft scattering models are assumed to describe the medium as a series of static color scattering centers, the incoming parton is subject to Brownian motion due to multiple soft scatterings with the medium, and a constant jet transport coefficientq is presumed to characterize adequately the jet-medium interaction process; Higher twist (HT) [23][24][25] models formulate the medium in terms of matrix elements of gauge field operators, and the properties of the plasma are specified through the entropy density s; AMY [26][27][28] characterizes the medium as a thermally equilibrated plasma, and describe it in the context of finite temperature field theory using the Hard Thermal Loop (HTL) rate equation approximations with all properties of the plasma specified by its local temperature T and baryonchemical potential µ B fields.
CUJET2.0 is the most recent extension of the opacity series formalism of Gyulassy-Levai-Vitev (GLV) [29-31, 43, 44] for applications to jet tomography of the QGP produced at RHIC and LHC. Jet tomography assumes that initial production of hard jets occurs before QGP formation and is reliably predicted via collinear factorized pQCD. The depletion or quenching of the intial rates jet fragments as a function of (p T , y, A, b, √ s, M ) can be used to probe the dynamical properties of QGP at short wavelengths under the assumption that jet medium interactions can be calculated via perturbative QCD multiple collision theory. In GLV theory, high energy jet energy loss is formulated as an expansion in the number of parton-medium scatterings and is found to be dominated by the first hard contribution in a medium in kinematic regions involving coherence of the long formation time compared to the size of the medium. If the QGP medium is well described by a quasi-parton Hard Thermal Loop plasma then the density of scattering centers ρ and the Debye screening mass µ as well as the plasmon mass, m g , can all be computed as a function of only the temperature T and the effective thermal coupling α s (4T 2 ). By relaxing the assumptions of HTL approximation, different non-perturbative models of the QGP can be tested.
The GLV theory correctly predicted in 2002 the general form of the √ s evolution of the high p T pion nuclear modification factor R AA (p T , η = 0; √ s, b) = dN AA→π / (T AA (b)dN pp→π ) from SPS, RHIC to LHC energies. GLV was generalized to DGLV [32] to include the kinematic effects due to thermal masses and extend to charm and beauty quark flavors. However, it was found that DGLV radiative energy loss significantly underpredicted the quenching of non-photonic electrons from charm and bottom quark jets. This led to the WHDG [33] generalization of DGLV theory to include elastic scattering as well as more realistic jet path length fluctuations. Those effects were found to be insufficient to solve the "heavy quark puzzle". This led Djordjevic [45] to develop a dynamical generalization of DGLV, replacing the Gyulassy-Wang (GW) [29] static color electric scattering potential with dynamical color magnetic interaction included potential through the HTL weakly coupled QGP ansatz. Including the above generalizations of DGLV/WHDG approach, a powerful numerical code, CUJET1.0, was developed by Buzzatti [34,46], that finally solved the "heavy quark puzzle" and predicted a novel quark flavor inversion of the normal nuclear modification factor hierarchy as a unique signature of perturbative QCD based jet tomography models, which is not shared by the jet holography counterpart.
The CUJET1.0 [34,46] Monte Carlo code developed by Buzzatti implemented numerically the dynamical DGLV opacity series and featured (1) an interaction potential that could interpolate between the pure HTL dynamically screened color magnetic limit and static Debye color electric screening limit; (2) the ability to calculate high order opacity corrections to radiative energy loss up to 9th order; (3) explicit integration over jet path in diffuse nuclear geometries including Bjorken longitudinal expanding HTL QGP; (4) inclusion of fluctuating elastic energy loss; (5) the evaluation of the convolution over numerical tables of pQCD √ s dependent initial jet production spectra of all flavors; and (6) the final convolution over jet fragmentation functions and evaluation of semi-leptonic final decay into non-photonic electrons.
CUJET1.0 was found to explain for the first time [34] the anomalous high quenching of non-photonic electrons ("heavy quark puzzle") within a pure HTL QGP paradigm as due to the enhanced dynamical magnetic scattering effects proposed by Djordjevic. It also predicted a novel inversion of the π < D < e − < B flavor ordering of R f AA at high p T that has yet to be tested at RHIC and LHC. However, the "surprising transparency" [47] of QGP produced at LHC severely challenged the fixed coupling version of the CUJET1.0 inherited the failure of the fixed coupling approximation used in the WHDG extrapolation from RHIC to LHC.
This led to the present extension of the dynamic DGLV/WHDG theory called here CUJET2.0. This Monte-Carlo code generalizes CUJET1.0 to include multi-scale running coupling effects as well as full 2+1D transverse and longitudinal expanding medium which azimuthal tomography is very sensitive to. Both versions of CUJET were developed as part of the ongoing DOE Topical JET Collaboration Project [48] with the mission to develop more quantitative jet quenching codes coupled to state of the art bulk observable constrained viscous hydrodynamic fields. First results with this code for azimuthally averaged R π AA and the χ 2 determination the jet transport coefficientq/T 3 from fits to RHIC and LHC data, as well as comparison to other JET collaboration model approaches were recently reported in [49]. This long write up aims to document the physics and numerical details of the current CUJET2.0 model as well as present a systematic χ 2 comparison not only to azimuthal averaged R π AA but also to its elliptic azimuthal moment v 2 (p T ) that remains an open problem at this time. In addition, the quark flavor dependence of the above observables in this latest version is also presented.
Motivated by the high energy kinematic regions being probed presently at LHC and the complex reaction vertices involved in the jet-medium processes, in CUJET2.0 we introduce a physically motivated multi-scale running of strong coupling factors α s (Q 2 i (x, k, T )) in the DGLV opacity expansion, where α s is the vacuum QCD running coupling bounded by an upper limit α max coupling strength. In addition, we explore the sensitivity to nonperturbative color magnetic and color electric screening mass deformations by varying multiplicative screening parameters (f E , f M ). (The perturbative HTL QGP limit is (1,0).) The three control parameters (α max , f E , f M ) define the dynamical model in this running coupling extension of DGLV.
An especially important new feature of CUJET2.0 is its ability to adaptively read in a variety of 2+1D viscous hydrodynamic temperature T (x, y, τ ) evolution grids to be used to perform jet path opacity integrations. Thus far we have used only the VISH2+1 event averaged grids available from the JET Collaboration depository. Future applications to event by event fluctuating hydro grids are planned. This paper is organized as follows: first, we postpone many of the technical details and development of both CUJET 1.0 and 2.0 models to a series of appendices. See Appendix B for a review of the fundamental ingredients of CUJET. In this appendix, dynamical DGLV opacity expansion, elastic energy loss, path length and energy loss fluctuation, convolution of energy loss probability distribution over initial production spectra and fragmentation functions will be discussed.
In Section 2, we discuss the choice of scales in the multi-scale running strong coupling extension of DGLV. We quantify the CUJET running coupling effects on jet energy loss in terms of a phenomenological "abc" power law energy loss model [50][51][52][53][54][55]. We introduce an effective jet medium interaction potential which is able to interpolate color electric and magnetic screening effect, as well as extrapolate to non-HTL scenarios.
Section 3 presents the main numerical results obtained with CUJET2.0 including pion nuclear modification factor R π AA at RHIC and LHC central and semi-peripheral A+A collisions; χ 2 /d.o.f. analysis and discussions on the consistency of the model in various collision configurations; jet transport coefficientq/T 3 and its variation with temperature and jet energy; flavor dependent suppression pattern for pion, D meson, B meson and nonphotonic electron and the mass ordering in R AA ; jet quenching with respect to reaction planes and single particle azimuthal anisotropy v 2 .
As a further exploration with CUJET2.0, in Section 4, we study the thermalization time's effect on pion suppression factor, and non-HTL scenario's implication of nonperturbative near T c physics in the CUJET2.0 framework. Finally, we summarize our main results and conclusions, and discuss possible future works, improvements and tests on CU-JET2.0 in Section 5.

Running strong coupling effects
Recent data from the LHC showed a significantly steeper rise of R AA with p T in the range of 5 − 100 GeV/c [3,5] than predicted by fixed coupling extrapolation from RHIC via WHDG [47]. This indicated a "surprising transparency" of the QGP at LHC to high energy jets as compared to bulk multiplicity scaling (by a factor (dN LHC /dy)/(dN RHIC /dy) ≈ 2.2) of the QGP opacity assumed in WHDG. The fixed coupling version of CUJET1.0 also encountered the same difficulty [56]. An generic dE/dx model analysis [51,52] also found that effective jet medium κ ∝ α 3 s coupling required to fit the slope and magnitude of central R AA (p T ) at both RHIC and LHC is ∼ 30% less at LHC than at RHIC.
The above problem motivated us to study whether UV running coupling effect could account for the relatively greater transparency of the QGP at LHC with CUJET2.0.

Multi-scale running coupling for radiative and elastic energy loss
Earlier estimates of running coupling effects were made by Zakharov [57,58]. In CUJET2.0 model, we follow a similar scheme using the 1-loop pQCD running coupling that is cutoff in the infrared when the coupling reaches a certain maximum saturation value α max for Q ≤ Q min : where the saturation scale Q min is fixed by α max as Q min = Λ QCD exp 2π 9α max .
The choice of the α max parameter is not obvious and is regarded here as a key infrared control parameter of the CUJET2.0 model. In principle, it can also depend on the local temperature field. At T = 0, the α max ≈ 0.7 estimate arose historically from the analysis of the heavy quark production in the vacuum [59]. However, in a QGP, suppression of α s at scales ∼ T is also expected. Lattice QCD qq potential studies [60] found the effective thermal α s (T ) coupling decreases monotonically from 0.5 at T ∼ 175 MeV to 0.35 at T ∼ 400 MeV. More generally, the effective running coupling α s (Q, T ) depends on the relevant Q of an observable as well as on T .
In the dynamical case of CGC gluon production the generalization of fixed coupling rates in [61,62] led to intricate multi-scale running coupling modifications of the fixed coupling formulae. A corresponding generalization of DGLV to running coupling remains a challenging and open problem because virtualities in radiative amplitudes depend on multiple kinematic (x + , k, q) 1 as well as temperature dependent infrared screening scales and plasmon masses.
We explore in CUJET2.0 the running coupling effects using physics motivated ansatz for relevant virtuality scales. At leading opacity order, we can identify in the fixed coupling DGLV forumla (see Eq. B.17), three distinct scales Q i (i = 1, 2, 3) controlling the strength of different aspects of the physics: 1. Two powers α 2 s (Q 2 ) clearly originate from the jet-medium interaction vertices from the exchanged transverse momentum q, and so for these we simply take Q 2 1 = q 2 .
2. One power α s (Q 2 ) originates from the radiated gluon vertex. The off-shellness in the intermediate quark propagator for one of the three amplitudes where the gluon is emitted after the scattering is Where k is the transverse momentum of the radiated gluon, M is the mass of onshell quark and m g is the plasmon mass of gluon. An ambiguity arises from other amplitudes for example if the radiated glue scatters with q instead of the quark. In the limit when k q and mass effects are negligible , (2.3) as in DGLAP radiative splitting.
3. Running thermal couplings can arise from the Debye mass µ(α s (Q 2 ); T ) and plasmon mass. We allow these to run with scale Q 2 3 = (2T ) 2 . 2 Note in the above choices of running scales there is no explicit dependence on the jet energy, which comes instead from the kinematic limits of the q and k integrations.
4ET . With the above scheme, the running coupling DGLV inclusive radiative fractional energy loss distribution at first order in opacity is then given by (the notation is as defined in Appendix A and B): where C R is the quadratic Casimir of the jet (C F = 4/3 for quark jets, C A = 3 for gluon jets); z = (x 0 + τ cos φ, y 0 + τ sin φ; τ ) is the path of the jet created at (x 0 , y 0 ) in the production plane along azimuthal angle φ; ρ(z) and T (z) is the number density and temperature evolution profile of the medium; χ 2 (z) = M 2 x 2 + + m 2 g (z)(1 − x + ) controls the "dead cone" and Landau-Pomeranchuck-Migdal (LPM) destructive interference, squared gluon plasmon mass m 2 g (z) = f 2 E µ 2 (z)/2, HTL Debye mass µ(z) = g(z)T (z) 1 + n f /6, g(z) = 4πα (4T 2 (z)); integration limit 0 |q| min(|k|, 4ET(z)), 0 |k| We further include running coupling effects in the elastic portion of the energy loss following the work of Peigné and Peshier [65]: both powers of α s in Eq. (F.1) run witht, and when integrated over dt in Eq. (F.2), we obtain respectively. Note in the running coupling scenario, the calculation of average number of collisions involves recursively solving the E(z) integral equation.
The choice of running scales Q i is of course subject to significant uncertainties at present. To estimate the systematic uncertainties on the nuclear modification due to the variation of running scale we increase or decrease the running scales Q i by 25 or 50 percent respectively, and refit the fixed reference point at p T = 30 GeV by changing the free α max parameter. The results are shown in Fig. 1.
With a static Glauber transverse background in CUJET, to compensate 50% decreased running scales Q i , α max needs to be 25% lower; while to compensate 25% increased running scales Q i , α max needs to be 50% higher. Since the radiative energy loss depends on α max with approximately a cubic law according to Eq. (2.4), making the running scales Q i smaller by 50% or greater by 25% can result in approximately 5% decrease or 10% increase in R AA . The systematic error bar coming from varying running scales is therefore significantly smaller than the experimental errors at this time, indicating the relative insensitivity of CUJET to the precise choice of running scales.
where P (τ ) corresponds to the momentum of a massless jet passing though a plasma characterized by a local temperature T . The power of T is constrained by simple dimensional analysis, and the index a and b are set by the asymptotic LPM behavior of the GLV model. In the fixed coupling case, For the range of energies of interest, log (E/T )/E ∼ E a , with a ∼ 1/3 − 1/4. and LHC (right) in CUJET with Glauber static transverse plus Bjorken longitudinal expanding background. The gray opaque curves use a fixed coupling with α s = 0.3, while the black curves use a running coupling with α max = 0.4. The difference is notable, especially in the higher energy range available at the LHC, while RHIC results are left almost unchanged. The sensitivity to the variation of running scales Q i (cf. Eq. (2.1) and following) is measured by the red curves: on one side we decrease the value of all scales Q i by 50% and lower α max to 0.3 (red dashed), on the other we increase all scales Q i by 25% and increase at the same time α max to 0.6 (red dotted). α max is constrained to fit R π,LHC AA (p T ≈ 30 GeV) = 0.35.
In Fig. 2 we show the value of the index a as a function of the jet energy E, for five different cases: α s fixed, only α s (4T 2 ) running, only α 2 s (q 2 ) running, only α s (k 2 /(x(1 − x)) running and finally all couplings running. In this example, we use a non-uniform density profile generated from Glauber model with Bjorken expansion that approximates the thermal medium formed in a Pb+Pb central collision at LHC. The results are insightful: • As expected, the fixed case shows a ∼ 1/3 − 1/4.
• By introducing the thermal coupling, only the absolute value of the energy loss is affected and the energy dependence of the index remains unaltered. The scale at which the thermal α s is evaluated is in fact independent of E. Not noticeable in this plot, at very high temperatures the reduced thermal coupling causes a stronger quenching compared to the fixed coupling case, since the smaller Debye mass diminishes the screening in the plasma. This running effect is however small: for most of the temperature ranges, α s is in fact equal to the saturated value α max (with Q 0 ∼ 1 GeV, T needs to be greater than 0.6 GeV to start feeling the running effects).
• The couplings α 2 s (q 2 ) and α s (k 2 /(x(1 − x)) sensibly reduce the dependence of ∆E/E on E, and as a consequence the value of the index a gets smaller and closer to 0. The α 2 s (q 2 ) contribution is smaller since the q distribution is peaked at small values of q ⊥ , as opposed to the α s (k 2 /(x(1 − x)) contribution which is larger due to the high tails of the k distribution.

E GeV
Index a E Figure 2. Energy loss index a(E) (cf. Eq. (2.9)) for different assumptions of the running coupling in CUJET: fixed effective α s = 0.3 (black), only thermal coupling running (dashed red), only α 2 s (q 2 ) running (purple), only α 2 s (k 2 /(x(1 − x))) running (magenta), all couplings running (pink). The saturated α max value is chosen to be equal to 0.4, which corresponds to approximately Q 0 ∼ 1 GeV. The plot shows the energy loss of a light quark (M = 0.2 GeV) traveling from the origin of the transverse plane and through a gluonic plasma (n f = 0) of size L = 5 fm, whose density profile is generated from Glauber model and resembles the medium created in a Pb+Pb √ s N N = 2.76 TeV b = 0 fm collision.
• The all-running case shows almost no dependence of ∆E/E on E, and a(E) ≈ 0.
According to Fig. 2, the running coupling drastically alters the jet energy dependence of the energy loss, making ∆E/E approximately independent of E. This naively implies less quenching at high energies and an increase in the R AA slope. Fig. 1 in fact proves our assertions. What is remarkable is the fact that the change in the slope of R AA cannot be mimic by a rescaling of the fixed coupling α s : this measurement constitutes a potentially clear signature of running coupling effects.

Bulk evolution profile
The evolution profile of the bulk is encoded in CUJET as local density of scattering centers ρ(z) and local temperature T (z), in both the radiative energy loss Eq. (2.4) and elastic energy loss Eq. (B.27). To gain meaningful information about the parton medium interaction mechanism from comparing predicted jet quenching observables with experimental measurements, QGP evolution profile plays an essential role, and carefully constrain it with for example bulk low p T flow data is critical. In CUJET1.0, a static transverse geometry is generated from Glauber model, and longitudinally a Bjorken expansion is applied. In CUJET2.0, this picture is replaced with more realistic fluid fields generated from 2+1D viscous hydrodynamics, and the medium has dynamical expansion both transversely and longitudinally.

CUJET1.0: Glauber initial with 1+1D Bjorken expansion
In CUJET1.0, Glauber model [66] is used to generate the geometry of the collision, the plasma density profile and the jet production point distribution. Time evolution of the plasma density ρ QGP is given by the Bjorken picture. A simple analytical expression is obtained by making few assumptions: (1) the system expands only longitudinally, along the beam direction (1+1D expanding plasma); (2) the plasma is a perfect fluid; (3) the computation is carried out in a relativistic hydrodynamical framework.
In particular, the following steps are implemented in the model calculation: firstly, Density profile of nucleus A is generated with Woods-Saxon parametrization . (2.11) This density is normalized to mass number A. R is the nuclear radius and a represents the surface thickness. Subsequently, the thickness function of the nucleus A is defined as After acquired the corresponding thickness functions, the distribution of participants in a collision between two nuclei A and B that collide with impact parameter b is then given by where σ in is the inelastic nucleon-nucleon cross section. The A, R, a, σ in parameters chosen for Au+Au (RHIC) and Pb+Pb (LHC) collisions in Eq. (2.11)(2.13)(2.16) are listed in Table 1: Au 197 6.37 0.535 42 Pb 207 6.48 0.535 63 Table 1. Woods-Saxon parameters used in CUJET. A is the mass number of the nucleus, R is the nuclear radius, a is the surface thickness and σ in is the inelastic nucleon-nucleon cross section.
In determine the proper time τ (= √ t 2 − z 2 ) dependence of QGP density profile, which is characterized by longitudinal boost invariance, for practical applications, we concentrate on the mid-rapidity region of the collision (y = 1 2 log t+z t−z = 0). The QGP density field ρ QGP,0 (x, τ ; b) in the azimuthal direction is given by: where dN g /dy represents the gluon rapidity density, it is proportional to the measured charged hadrons rapidity distribution dN ch /dy, with dN g /dy = (3/2)dN ch /dy. Thermalization time τ 0 is chosen to be τ 0 = 1 fm/c in CUJET1.0. We introduce f (τ /τ 0 ) to characterize the pre-thermal stage and the evolution profile after the medium is fully thermalized. In absence of a clear theoretical answer to the way high energy jet couples to the medium before thermalization, we make such phenomenological assumptions: The density "seen" by the jet grows linearly until thermalization time τ 0 is reached, thereafter it decreases as 1/τ , converges to the Bjorken expansion picture. Note different parametrizations for the temporal evolution of the system exist, and by choosing the linear thermalization scheme systematic uncertainties are inevitably introduced. A detailed discussion about the choice of thermalizing schemes and associated errors can be found in Appendix I. The jet production points are distributed according to the binary collision distribution, which is given by (2.16) The ability of CUJET to perform a full jet path integration allows us to parametrize the evolution of the system in different ways, and we can perform comparisons with experimental measurements and draw insightful conclusions on the physics of the collision. See Appendix I for a comprehensive analysis of thermalization phase effects on jet quenching physics.

CUJET2.0: viscous hydrodynamics
In Appendix I's Fig. 28, we show a smooth temperature profile of a symmetric plasma. The shape of the region of interest takes the form of a perfect circle when the impact parameter is null, or an almond when b = 0.
The reality is however different. Nuclear matter is very granular at short distances, and the nucleons are not distributed in a perfectly symmetrical way. The naive picture of a circle or an almond is an idealization of the collision geometry in most situations, and the identification of a reaction plane determined by the orientation of the impact parameter is often a hard experimental task. The average over multiple collisions might lead to a smooth temperature profile, but this is not the case on an event-by-event basis, where fluctuations of initial conditions might lead to considerably different results.
Therefore, the full three-dimensional hydrodynamic expansion may differ substantially from the Bjorken ideal 1+1D hydrodynamic evolution in CUJET1.0. A complete description of the system must include not only transverse expansion, but also viscous corrections to the perfect fluid. And furthermore, effects such as initial condition fluctuation or jet energy deposition into the medium should be considered. We hence adapt the grid of QGP fluids in CUJET2.0 to encompass more realistic 2+1D viscous hydrodynamical fields, incorporating a dynamical medium expanding both transversely and longitudinally. This is for the first time the DGLV opacity series is fully coupled to viscous hydro fields, and the combination of strongly coupled bulk dynamics and weakly coupled pQCD energy loss theory generates rather indicative results of quenching observables computed from CUJET2.0 under RHIC and LHC conditions. The numerical analysis of CUJET2.0 output in RHIC and LHC nuclear collisions will be presented in Section 3.
In principle, the present CUJET2.0 framework can be coupled to virtually any complex geometries and plasma evolution profiles. A wide range of flow fields generated by external hydrodynamical code can be applied, e.g. transverse blast wave model [67,68], viscous RL hydro [69,70], VISH2+1 [71][72][73], etc. At present stage, we have used only the VISH2+1 event averaged evolution profile available from the JET Collaboration depository to study the azimuthal angle and transverse momentum dependence of high-p T light to heavy flavor quenching pattern, incorporating event by event fluctuating hydro fields is a work in progress.
VISH2+1 [71][72][73] utilizes viscous hydrodynamics to describe the fireball evolution. In the version adapted by CUJET2.0, MC-Glauber initial conditions are used to sidestep issues related with hypothetical early non-equilibrium evolution; using Cooper-Frye algorithm [74] along a hypersurface of constant temperature T f = 120MeV, a sharp transition from viscous fluid to free-streaming particles is generated to describe hadronic rescattering and kinetic freeze-out; the s95p-PCE (partial chemical equilibrium) equation of state (EOS) is constructed according to [75], which matches Lattice QCD data at high temperature and recovers the hadron resonance gas at low temperature. The various input parameters are adjusted to fit final hadron spectra and elliptic flow in low transverse momentum p T < 1.5(2.5) GeV/c region in [72]. In particular, experimental data of pion and proton spectra in 200AGeV Au+Au central collisions (0-5% centrality, b=2.33 fm), pion, proton and charged hadron elliptic flow v 2 (p T ) in semi-peripheral collisions (20-30% centrality, b=7.5 fm) are compared. With MC-Glauber initial conditions, s95p-PCE EOS and 120MeV freeze-out temperature, for a QGP with number of quarkonic flavor n f = 2.5, the starting time τ 0 at which the system is sufficiently close to local thermal equilibrium for viscous hydrodynamics to be applicable is calculated to be τ 0 = 0.6 fm/c 3 , and the key QGP transport parameter, shear viscosity η/s (the ratio between shear viscosity η and entropy density s), is phenomenologically extracted to be η/s = 0.08.
For initial jet production distributions, presumably binary distributions generated from corresponding viscous hydro should be used. However, in CUJET2.0 the binaries are given according to Eq. (2.16), because the spacing of the VISH2+1 x-y grid in CUJET2.0 is 0.5 fm while in computing final jet spectra the initial jet production points have a uniform 1.0 fm step in the radial direction. Systematic uncertainties resulting from binary distributions are therefore negligible. Nevertheless, this is a potentially interesting case if the quenching of jet occurs mostly at the edge of the medium, under which circumstance the fineness of the grid will play an important role. This topic will be explored in future studies with CUJET.

CUJET2.0 Results at RHIC and LHC
In this section, CUJET2.0 is applied to RHIC Au+Au 200AGeV and LHC Pb+Pb 2.76ATeV collisions. Impact parameter b=2.4fm and 7.5fm is used to simulate 0-10% and 10-30% centrality respectively. For radiative energy loss, the DGLV opacity series is calculated to first order because of the computational efficiency of the Monte Carlo algorithm. The convergence of the DGLV expansion is discussed in detail in Appendix C, from there we see, since CUJET2.0 concentrates on high energy jet suppression and averages over all possible path lengths in a realistic heavy ion collision, the first order in opacity can be regarded as a good approximation to the series. If to impose an artificial systematic uncertainty on higher order contributions, in terms of inclusive π, D, B, e − R AA , one can estimate the associated variation to be less than 15% (cf. Appendix C.2).
We implement CUJET2.0 to study jet quenching observables, in particular R AA and v 2 in the high transverse momentum p T > 5 GeV/c region where eikonal and soft approximation are applicable. RHIC inclusive neutral pion suppression factor R π AA (p T = 15GeV/c) = 0.3 at Au+Au 200AGeV 0-5% centrality is set as a reference point to fix the maximum coupling constant α max . We then extrapolate our calculation to p T = 5 ∼ 20 GeV/c region at RHIC, and p T = 5 ∼ 100 GeV/c region at LHC central to semi-peripheral collisions for pion, D meson, B meson and non-photonic electron R AA and v 2 . Rigorous χ 2 /d.o.f. study is conducted for inclusive pion spectra, and comprehensive azimuthal tomography is applied to pion single particle anisotropy v 2 in A+A collisions at both RHIC and LHC.
To elucidate our theoretical predication of quenching observables, we choose to plot R AA curves without adding error bands, and we will for completeness list the known systematic uncertainties which contribute to hadron spectra for CUJET at present stage in the next paragraph, besides the error from first order DGLV opacity calculation which has been discussed above. The uncertainties associated with elliptic flow v 2 can only be partially interpreted from errors in R AA because v 2 depends on more complex factors such as fluctuation, anisotropy, inhomogeneity, etc. Systematically induced variation of azimuthal anisotropy in our model is subject to future exploration.
A list of systematic uncertainties in CUJET 4 : • Running scale variation. Increase the running scale by 50% generates approximately 5% enhancement of R AA , while decrease the scale by 25% leads to about 10% lower R AA , cf. Section 2.1.
• Kinematic limits for k ⊥ integration in Eq. (B.17) and (2.4). Depending on the interpretation of energy fraction x E or light-cone x + in the DGLV formula, and the treatment of large angle emissions which break down the collinear approximation, fractional energy loss varies. The coupling constant needs to be tuned accordingly ±10% at most, cf. Appendix D.1.
• Number of quarkonic flavors n f for the medium. A rescaling of coupling constant by 6% perfectly reconciles the scenario of n f = 0 and n f = 2.5, cf Appendix E.
• Fragmentation temperature T f . In terms of R AA , lower T f to be below freeze-out temperature, typically around 100 MeV, by about 50% has no significant influence. But increase it by approximately 100% effectively generates a 20% weaker coupling, cf. Appendix E.
• Initial rapidity density dN/dy. In CUJET1.0, the variation of rapidity density significantly changes the magnitude of R AA , but the shape of suppression curves remains semi-stable, thus one can rescale the coupling constant to balance this effect, cf. Appendix E. However, the results presented in the following sections are calculated from CUJET2.0, within which framework the bulk is assumed to be properly modeled by the 2+1D viscous hydrodynamical fields, therefore variations coming from dN/dy are beyond the scope of CUJET2.0.
• Thermalization scheme. Fixing initial time τ 0 , divergent and free streaming prethermal scenario creates approximately an effective 10% larger and 7% smaller coupling constant comparing to linear scheme. For heavy flavor R AA , the scheme variation slightly affects the slope of the quenching pattern, however this change is minor when juxtaposed with experimental error bars, cf. Appendix I.
• Energy loss fluctuations. The assumed Poisson distribution for radiative energy loss has negligible effects on R AA , but it may significantly influence v 2 , cf. Appendix G.
• Partonic pp spectra variations. Error bands from NLO and FONLL initial cross sections span 5% at partonic R AA level, being insignificant to a certain extent. Absolute normalizations of the spectra drop out when calculating ratios such as R AA , but steepnesses matter. Relative steepness between production spectra has influence on the calculation of pion and non-photonic electron spectrum, which is fragmented from gluons and light quarks, and charms and bottoms receptively, cf. Appendix H.

Pion nuclear modification factor
Nuclear modification factor R AA is a key observable at RHIC and LHC relativistic heavyion collisions, it describes the relative magnitude of jet quenching. R AA is defined as the ratio of the quenched A+A spectrum to the unquenched p+p spectrum, scaled according to the number of binary collisions N binary : N binary dσ pp dp T (p T ) . (3.1) We have suppressed the explicit dependence on rapidity y and c.m. energy √ s in Eq. (3.1). Here the R AA is understood as azimuthally averaged, and azimuthal angle φ is integrated out.

Pion suppression
We calculate in CUJET2.0 the inclusive pion R AA at RHIC Au+Au 200AGeV and LHC Pb+Pb 2.76ATeV, central (b = 2.4 fm) and semi-peripheral (b = 7.5 fm) collisions, the results are illustrated in Fig. 3. Theoretical R AA curves are compared side by side with corresponding experimental measurements of charged hadron suppression factor at ALICE [3] and CMS [5], and neutral pion suppression factor at PHENIX [10]. In the radiative energy loss sector, DGLV opacity series is calculated to first order, with f E = 1, f M = 0 in Eq. (2.4) interpolating the pure dynamical scattering potential in the hard thermal loop scenario. Maximum coupling constant α max is adjusted to α max = 0.26 to fit the R AA reference point at RHIC central collision.
After adjusting α max = 0.26 to match the calculated pion nuclear modification factor from CUJET2.0 with the p T = 15GeV/c reference point of experimentally measured inclusive neutral pion suppression factor from PHENIX 2012 [10] Au+Au 200AGeV central collisions, the rest of the R AA curve in the range of p T = 5 ∼ 20GeV/c shows reasonable compatibility with both PHENIX 2008 and 2012 data. More importantly, when move on to simulate RHIC 20-30% centrality collisions by changing solely the impact parameter to b = 7.5 fm and fixing all other parameters in CUJET2.0, the theoretical R AA result demonstrates even better agreement with experimental data.
When switch to LHC after constrained all CUJET model parameters with RHIC data, fixed coupling CUJET1.0 used to encounter difficulties explaining the surprising transparency of the QGP at LHC high p T region [34], though this problem is eased by running coupling CUJET1.0 which has effectively reduced coupling strength at high energies, the pion R AA 's steep rising and successive flattening pattern at LHC remains only partially explained [56]. This issue is fully solved in CUJET2.0 which has a more realistic bulk evolution profile. As shown in Fig. 3(c)(d), at both ALICE and CMS, both central and semi-peripheral collisions, the CUJET inclusive pion R AA curves seamlessly simulate both the low p T steep rising and high p T saturating behavior of π 0 or h ± nuclear suppression factor.
One of the most appreciable signatures of CUJET2.0 which has 2+1D viscous hydrodynamic background is the drastic reduction of the strong coupling constant compared to CUJET1.0 where transversely a static Glauber medium is assumed. In running coupling CUJET2.0 calculations, the maximum coupling strength α max is adjusted to α max = 0.26 to fit inclusive pion R AA 's at both RHIC and LHC, central and mid-central A+A collisions, this value is distinguishably smaller than running coupling CUJET1.0's α max = 0.4. Compared to CUJET1.0's static Glauber bulk, CUJET2.0's average medium density is reduced because of the transversely expanding hydro fluids, and one would intuitively expect less quenching in such a medium. However, the fact that the effective strong coupling constant demands a tremendous reduction in itself to generate the same hadron suppression factor in CUJET2.0 as in CUJET1.0 indicates that there is another factor contributes more    [3] and CMS [5] h ± R AA for Pb+Pb collisions at 2.76ATeV and (c) 0-5%, (d) 20-30%(ALICE)/10-30%(CMS) centrality. The R AA is calculated at leading n=1 order in opacity, and the maximum coupling constant α max is constrained by fitting to central RHIC Au+Au data at a reference point R π AA (p T = 15GeV/c) = 0.3, setting α max = 0.26. For the effective potential in Eq. (2.4), the parameters are set to be f E = 1, f M = 0, i.e. dynamical QCD medium with HTL approximation. Impact parameter b=2.4 and 7.5 fm is used in CUJET2.0 to simulate the central and semi-peripheral collisions respectively. The 2+1D viscous hydro grid is generated from VISH2+1 [71][72][73], with τ 0 = 0.6 fm/c, η/s = 0.08, 120 MeV freeze-out temperature, MC-Glauber initial conditions and Lattice QCD s95p-PCE EOS at both RHIC and LHC. Compared to CUJET1.0 [56] which has a static transverse profile, with an transversely expanding medium in CUJET2.0, pion R AA flattens out more clearly at high p T region at LHC. The strong coupling however decreases from α max = 0.4 in running coupling CUJET1.0 to α max = 0.26 in CUJET2.0, indicating longer path length of jet in a transversely expanding medium overrides the reduction of density, and contribute to overall enhanced quenching.
remarkably to jet quenching and in fact dominates the modification of parton shower in the expanding medium. This key contributor can only be the path length of propagating jet in QGP. The longer jet path length generated from 2+1D viscous hydro fields in CU-JET2.0 plays a decisive role on parton energy loss, it overrides the effect of less quenching resulted from diminished medium density, and induces an overall more strongly quenched jet spectra in relativistic heavy ion collisions.

Chi squre per degree of freedom
Calculating relative variance per degree of freedom χ 2 /d.o.f. is one of the best quantitative methods to test to what extent the theoretical R AA results are in agreement with experimental measurements. This quantity is the average of relative variance, which is the ratio of squared difference between experimental data point and its theoretical counterpart to the quadratic sum of all theoretical and experimental statistical and systematic errors associated with that point, over all selected data. It is defined as follows: Where V th is the theoretical value, V exp is the experimental value, t σ 2 t stands for the quadratic sum over all types of errors that one chosen point has, and N is the number of data points selected.
We vary the maximum coupling constant α max in CUJET2.0 from 0.20 to 0.35 with 0.01 steps, and maintain the dynamical HTL scenario by fixing f E = 1, f M = 0, to study the most compatible CUJET2.0 one parameter (α max ) fit at RHIC and LHC, and test the consistency of the model at different √ s N N 's. Fig. 4 shows the pion R AA curves with those different α max values at RHIC Au+Au 200AGeV and LHC Pb+Pb 2.76ATeV central (b = 2.4 fm) and semi-peripheral (b = 7.5 fm) collisions. The experimental data being compared with are PHENIX 2008 [7], 2012 [10] and STAR [12] π 0 R AA at RHIC; and ALICE [3] and CMS [5] h ± R AA at LHC. In all four panels of Fig. 4, focusing on p T < 40 GeV region, the magnitude of inclusive hadron suppression has near uniform increment with an uniformly increasing maximum coupling α max , with exception at relatively large α max 's where the spacing between R AA curves becomes smaller, but the monotonic lowering of R AA proceeds. Since the saturation scale Q min for the running coupling depends solely on the maximum coupling constant α max , i.e., Q min = Λ QCD Exp{2π/9α max }. Take α max = 0.20, 0.25, 0.30 and 0.35 for example, the saturation scale Q min = 6.56, 3.26, 2.05 and 1.47GeV respectively. At relatively low α max , because of the large saturation scale, the strong coupling recovers asymptotically the fixed coupling scenario up to a relatively high energy, this explains the near uniform increment in the panels. The influence of running coupling is substantial at relative high α max where the minimum running scale is low, in that situation the logarithmic decay of coupling strength resulted from vacuum running shrinks the spacing of R AA 's more effectively.
A significant phenomenon shows up in the bottom panels of Fig. 4 -the flattening pattern (slope) of R AA in high p T (p T > 50 GeV) region at LHC is almost independent of   the choice of α max , this implies the relative insensitivity of R AA saturation to the running coupling effect, and we can therefore exclude to a certain extent the influence of running on the saturation of R AA for ultra-high energy jet. Note in [56], the previous calculation of multi-scale running coupling combined CUJET1.0, whose medium assumes static Glauber transverse profile plus 1+1D Bjorken longitudinal expansion, did not exhibit a clear signature of R AA flattening. Therefore, evident R AA saturation comes largely from the kinematics in a medium with both transverse and longitudinal expansion, which feature distinguishes CUJET2.0 from running coupling CUJET1.0. A dynamically transverse expanding medium, for instance a 2+1D viscous hydro fluid, plays a very important role in the R AA flattening of high p T jet in ultra-relativistic heavy ion collisions, and shall receive more attention in predicting jet quenching observables in A+A collisions from pQCD energy loss models. The essential eikonal and soft approximation in dynamical DGLV opacity expansion may break down at low p T region, hence for the purpose of is an indicative signature of model consistency, and the constrained α max range should be independent of whether p T > 5 GeV or p T > 8 GeV is chosen as long as the minimum p T is sufficient for preserving basic assumptions of the CUJET2.0 model and number of points being selected at high p T is large enough. Hence for safer comparison we choose p T > 8 GeV, and   [7,10,12], at LHC, ALICE and CMS data [3,5] are compared respectively. The α max ranges forχ 2 < 1 andχ 2 < 2 ( in Fig. 5, we list detailed α max ranges withχ 2 < 1 andχ 2 < 2 (χ 2 ≡ χ 2 /d.o.f.) at RHIC and LHC in Table 2.  Table 2. The ranges of α max withχ 2 < 1 andχ 2 < 2 ( The combination of Fig. 5 and Table 2 provides quantitative information about the consistency of CUJET2.0 HTL model in various A+A collision configurations. We see that if strictly constrainχ 2 to be less than 1, for b = 2.4 fm central collisions, CUJET2.0 results at RHIC and LHC have 0.01 offset in α max , and for b = 7.5 fm semi-peripheral collisions, the results are in perfect agreements with RHIC and LHC at α max = 0.23 − 0.25 range. We also notice the averaged best fit α max value at RHIC and LHC in semi-peripheral collisions is around 0.03 lower than central collisions, and at either centrality the best fit LHC α max is approximately 0.03 lower than RHIC. These observations will trigger useful analysis in Section 3.3.
At present stage, in the CUJET2.0 HTL scenario, if restrict separate maximumχ 2 to be 2, we find that the intersecting α max range for all four collisions at RHIC and LHC, i.e. Au+Au 200AGeV and Pb+Pb 2. Based on all the above discussions, we conclude that from testing CUJET2.0 HTL scenario's agreement with centrality dependent neutral pion and charged hadron suppression factors at RHIC and LHC in the mid-rapidity region, the maximum coupling constant in the model is constrained to α max = 0.25 − 0.27, in which range the averaged χ 2 /d.o.f. is strictly less than 1.5; if allow average χ 2 /d.o.f. < 2, then α max = 0.23 − 0.30. As discussed in Section 3.1.1, the small α max value itself can be attributed to the dominating longer jet path length feature for the parton shower modification in a transversely expanding medium. Notice in the calculation of χ 2 /d.o.f., we did not include in it any intrinsic CUJET2.0 systematic uncertainties which were discussed at the beginning of Section 3, this suggests the best fit α max region can in fact be broadened after taking those factors into account and by the mean time maintain the stringent χ 2 /d.o.f. limit. Given the complexity of estimating the complete systematic errors in the model, we will use α max = 0.25 − 0.27 to extract the effective jet transport coefficient in Section 3.1.3, and stick to α max = 0.26 to extrapolate the suppression pattern of open heavy flavors and heavy flavor leptons in Section 3.2.

Jet transport coefficient
The suppression of hadrons at large p T is understood to be caused by scatterings of the leading parton with color charges in the near thermal QGP. This process can be characterized by the jet transport coefficientq, which is defined as the squared average transverse momentum exchange per unit path length. CUJET2.0 treats thermal excitations in the assumed homogeneous QCD medium as partonic quasi-particles, and the transport parameterq in CUJET2.0 is related to the effective partonic differential cross section by the relation:q where the energy E and temperature T dependence comes in naturally from the partonic kinematics and plasma density. In CUJET2.0,q depends also on the maximum strong coupling constant α max , as well as electric and magnetic screening mass deformation parameters (f E , f M ), all of which originate from the effective cross section of the quark-gluon process: with the Debye mass µ(T ) = T 4πα s (4T 2 )(1 + n f /6). Note here the temperature is nonlocal. Eq. (3.4) differs from the effective scattering potential in Eq. (2.4) where the same cross section form appears but varies with local temperature T (z). We assume ρ ∼ 2T 3 for an idealized uniform thermal equilibrated medium, and calculate the CUJET2.0 jet transport coefficientq according to Eq.  [49]. At various initial jet energies, the reduction ofq/T 3 with rising T invariantly follows an approximate logarithmic law. Since ρ ∼ 2T 3 , this logarithm indicatesq grows slightly slower than medium density with increasing T, but there is still more elastic transverse momentum transfer per mean free path at higher temperature, as intuitively expected.  0, which set of parameters generates consistent fits to neutral pion and charged hadron suppression factor R AA at both RHIC and LHC both central and semi-peripheral A+A collisions.q/T 3 versus QGP temperature T at fixed incoming jet energy E is plotted on the left panel; the right panel showsq/T 3 versus E at fixed T. When E is fixed, the decrease ofq/T 3 with the increasing T follows approximately a logarithmic law, indicatingq itself gains slightly slower than medium density with rising T in CUJET2.0. This feature however still suggests there is more transverse momentum transfer per mean free path at higher temperature, as intuitively expected. When T is fixed, the logarithmic E dependence of q/T 3 at high energy region comes naturally from the kinematic limit of the exchanged transverse momentum.
For an idealized static equilibrium QGP with fixed temperature T, as illustrated in the right panel of Fig. 6, the transverse momentum transfer between the quark jet and dynamical scattering centers shows logarithmic dependence on initial jet energy if the T 3 contribution from the medium density is factored out. This is expected from the kinematic limit of transverse momentum exchange, i.e. (q 2 ) max = 4ET in Eq. (3.3). On the other hand, for quark jet with fixed initial energy E, the absolute jet transport coefficientq/T 3 drops at an diminished rate when temperature grows and reaches high T region. This is expected from the Debye mass µ's temperature dependence, i.e. µ(T ) = T 4πα s (4T 2 )(1 + n f /6). The thermal coupling has negligible contribution until when T is high, at that time the logarithmic decay of the coupling strength will weaken the linear increase of the Debye mass with rising temperature.
The QGP produced in heavy ion collisions is interpreted as strongly interacting medium [13], whose collective flow is known to be well described by relativistic hydrodynamics with a negligible shear viscosity, and effective perturbation theory may not be applicable to study interactions in a medium which is not dominated by quasi-particles. In [78], the authors derived a general expression relating the jet quenching parameterq with the shear viscosity η of a weakly coupled QGP, and the deviation from this relation is conjectured to be a more broadly valid measure of "strong coupling" of the medium than considering solely the shear viscosity divide by entropy density η/s. The relation is expressed as follows: η s ≈ 1.25T 3 /q for weak coupling , 1.25T 3 /q for strong coupling .

(3.5)
In high energy region where the QCD coupling is supposed to be weak, the simplified CUJET2.0q calculation from Eq. (3.3) and (3.3) showsq/T 3 ≈ 3.7 for typical temperatures reached by RHIC and LHC, and 1.25T 3 /q ≈ 4.2/4π. This value is larger than the quantum limit η/s = 1/4π, or the MC-Glauber VISH2+1's η/s = 0.08 which is extracted from fitting to hadron spectra and harmonics at low p T .
One way to reconcile the discrepancy is proposed in [49], where the authors suggest lattice calculation indicates that the non-perturbative soft modes in the collision kernel can double the value of the NLO pQCD result for theq [79,80]. And there are also recent formal pQCD calculations showing that NLO corrections can result in more than 50% increase inq [81][82][83][84]. Another possibility is that the presentq calculation in CUJET2.0 is over-idealized by disregarding anisotropy and heterogeneity/inhomogeneity, both factors can influence the jet-medium interaction significantly and require more careful considerations in both the bulk evolution sector and the energy loss sector.
We briefly summarize Section 3.1 here: CUJET2.0 inclusive pion R AA calculated with maximum coupling constant α max = 0.25 − 0.27 in a dynamical QCD medium with HTL scenario is strictly consistent with both RHIC Au+Au 200AGeV π 0 R AA and LHC Pb+Pb 2.76ATeV h ± R AA in both central and semi-peripheral collisions. Rigorous χ 2 /d.o.f. calculation indicates this parameter fit has averaged χ 2 /d.o.f. being stringently less than 1.5. And if allow average χ 2 /d.o.f. < 2, then α max = 0.23 − 0.30. The combined effect of multi-scale running coupling and transverse expanding medium give rise to the steep rising and subsequent flattening of inclusive hadron R AA at LHC. The small value of α max itself implies that longer jet path length in a transverse expanding medium overrides the reduction of density and contribute to enhanced overall quenching. Idealized CUJET2.0 effective jet transport coefficientq/T 3 is consistent with not only LO pQCD estimates, but also theq/T 3 's extracted from HT-BW, HT-M, MARTINI and McGill-AMY models [49] by fitting to the same set of experimental hadron suppression factors at RHIC and LHC A+A central collisions.

D meson, B meson and non-photonic electron
We apply CUJET2.0 to study the suppression pattern of not only neutral pions or charged hadrons, but also D mesons, B mesons and non-photonic electrons at RHIC and LHC. The open heavy flavor and heavy flavor lepton nuclear modification factors calculated from CUJET2.0 (α max , f E , f M ) = (0.26, 1, 0) HTL model are shown in Fig. 7. Inclusive nonphotonic electron R AA 's in central and semi-peripheral A+A collisions are compared with experimental data at PHENIX [8] and STAR [11], and D meson R AA 's are compared with measurements at ALICE [1].
Both running coupling CUJET2.0 and fixed coupling CUJET1.0 [34] predict a novel crossing pattern of p T dependent π, D, B, e − nuclear suppression factors. Compare the left panels of Fig. 7 to Fig. 1 in [34], one finds that regardless of the inclusion of multi-scale running coupling and dynamical viscous hydro fields in CUJET2.0, and the crossings of π, D, B, e − R AA 's in CUJET2.0 constantly occur at about same p T as in CUJET1.0. For example, at RHIC Au+Au 200AGeV central collisions, pion R AA intersects D meson, nonphotonic electron, B meson at p T ≈ 9, 19, 24 GeV, and at LHC Pb+Pb 2.76ATeV central collisions, pion R AA intersects D meson, B meson at p T ≈ 23, 33 GeV. Note all these p T 's are within the range of p T < 40 GeV where, as discussed in Section 3.1.1 and 3.1.2, running coupling effect contributes significantly to a steeper rising R AA . The running coupling induced change of slope emerges from pion, D meson and B meson R AA in a similar manner, this is as expected, because when considering the gluon radiation vertex for running coupling, the mass scale is small comparing to kinetic terms hence being dropped, cf. Eq. (2.2) and (2.3). The robust p T interval of R AA crossings in CUJET also suggests the mass ordering of π, D, B, e − suppression pattern comes intrinsically from the DGLV gluon radiation spectrum and TG elastic energy loss formula, and the bulk evolution profile has limited effect on this mass hierarchy.
Notice also the quenching pattern of non-photonic electron at Au+Au 200AGeV collisions calculated from both CUJET1.0 [34] and CUJET2.0 Fig. 7(a)(b) are in agreement with the RHIC data for central and semi-peripheral centralities. This clearly indicates the solution to "heavy quark puzzle" 6 is built intrinsically in the structure of CUJET energy loss framework. As discussed in Appendix F.2, the combination of dynamical medium effect and elastic energy loss first significantly brings down the light to heavy quark energy loss ratio. Secondly, appropriately weighing path length fluctuations of initial jet production coordinates plays a pivotal role, in CUJET this is realized by cutting off the DGLV integral at a dynamical T (z)| τmax = T f hypersurface which is parametrized by fragmentation temperature T f 7 . The third factor is the Poisson distribution assumption for inelastic energy loss, since one has noticed Appendix C.2 the variation of opacity can alter the gluon radiation spectrum for light and heavy quark jet differently.
Concentrating on the flavor dependent suppression pattern prediction from CUJET2.0, it is not only a revolution of [34] with multi-scale running coupling and transverse expanding medium effect, but also a critical supplementary with comprehensive semi-peripheral A+A predications and decisive non-photonic R AA at LHC. We can make several observations about the flavor dependent quenching scenarios in Fig. 7: firstly, R AA for inclusive D meson and pion tangles together at low p T at both RHIC and LHC. We notice in Appendix G the radiative energy loss probability distribution for charm and light quark almost overlap when jet has reasonably low initial energy, the similar suppression pattern of D meson and pion in this region would suggest comparable elastic energy loss probability distribution for them, and the D meson A+A production spectrum is expected to have a steeper slope than pion at low p T (cf. also Appendix H).
Secondly, B meson is much less quenched than both pion and D meson in the low p T region at both RHIC and LHC. Hint that the B measurement will be a clean refinement of pQCD energy loss models. Finally, at LHC, the R AA B > e > D mass ordering 8 at low p T evolves into e > B > D at p T ≈ 23 GeV, and RHIC seems to have the same inversion at slightly larger p T but less discernible than LHC. This mass ordering comes from a complex interplay between total energy loss probability distribution, and initial production spectra for charm and bottom jets, and fragmentation functions in the hadronization processes. Since non-photonic electron spectrum is the combination of B → e, D → e, and B → D → e channels, the change in mass hierarchy can partially be attributed to a significant enhancement in the B → D → e channel in certain p T range, this range occurs at lower p T at LHC which has larger multiplicity density and higher temperature than RHIC, and a semi-exclusive measurement of non-photonic electron production in AA can thus be a crucial benchmark.
One final comment about the flavor dependent suppression pattern from CUJET2.0 calculations is the theoretical agreement with ALICE [1] average prompt D R AA . Note impact parameter b = 2.4 fm in CUJET typically simulates 0-10% centrality, and b = 7.5 fm is applicable for 10-30% centrality, the integrated 0-20% centrality experimental D R AA should be within the b = 2.4 fm and b = 7.5 fm CUJET2.0 D meson R AA curves, and this feature shows up explicitly in Fig. 7 (c) and (d).
This brings us to a mini-summary of Section 3.2. We conclude that the robust crossing pattern of π, D, B, e − R AA 's is rigorously encoded in the flavor dependent energy loss structure of DGLV opacity expansion combined with TG elastic, and a transverse expanding medium has minor effect on the mass hierarchy; solutions to the "heavy quark puzzle" are intrinsically integrated in the framework of CUJET; and CUJET2.0 predicts a decisively less quenched B meson R AA which is well above D meson and pion at 5 GeV < p T < 15 GeV, as well as a critical alternation of low p T R AA 's mass ordering from B > e > D to e > B > D at p T ∼ 25 GeV.

Azimuthal flow
Anisotropic collective flow is a key observable in relativistic heavy ion collisions, it relates directly to the formation of QGP. With a large impact parameter for A+A event, the region of interest gains an increasingly asymmetric shape. And in a strongly coupled medium, the pressure gradients due to this initial azimuthal anisotropy are effectively transferred into the collective flow of its components. The different types of collective flows are quantified in terms of Fourier components of the azimuthal angle distribution [89]: Here dN h /dydp T dφ represents the number of hadrons of species h observed at rapidity y, with transverse momentum p T and azimuthal angle φ. Both dN h dydp T and the Fourier coefficients v i (p T , y, h) depend on the initial rapidity density dN i /dy. And dN i /dy is a function of the energy √ s and centrality b of the collision. Generally speaking, in the transverse plane with respect to the beam axis, the collective flow generated from non-central A+A collisions is centrosymmetric if there is no fluctuation, and odd number Fourier components drop out. Thus among all collective flow harmonics, elliptic flow v 2 is of the most significance, and by the mean time being least sensitive to fluctuations. Therefore, a comprehensive study of single particle v 2 can provide critical information about the azimuthal anisotropy, as well as useful information about the jet medium interaction mechanism.
3.3.1 Pion production with respect to reaction plane Simultaneously fit particle suppression pattern and azimuthal flow asymmetry is a decisive benchmark for all jet tomography models. 9 To visualize this simultaneity more clearly, one can calculate the R AA with respect to reaction plane. The typical choice of azimuthal angle set is the in plane φ = 0 and out of plane φ = π/2, and consequent nuclear suppression factors R in AA and R out AA are defined as: The AA, pp superscript and N binary has the same meaning as in Eq. (3.1). Since the p+p collision is presumably central, and generally no azimuthal anisotropy is expected, we have where φ 0 is an arbitrary azimuthal angle. In terms of rapidity y, the region we are interested in has y ≈ 0, hence we short-write dN pp h dydp T | y=0 as dN pp h dp T . Neglecting fluctuations by setting odd number harmonics to zero, in the mid-rapidity region, we get (3.8) We calculate pion's R in AA and R out AA in mid-rapidity region for p T up to 18 GeV/c in CU-JET2.0, for Au+Au 200AGeV central and semi-peripheral collisions, and compare with corresponding PHENIX [10] data. The results are shown in Fig. 8.
In Section 3.1.2, we have constrained the maximum coupling constant α max in the CUJET2.0 HTL scenario to be 0.25 − 0.27, this range of α max renders the most consistent pion R AA at RHIC Au+Au 200AGeV and LHC Pb+Pb 2.76ATeV central and semi-peripheral collisions through stringent χ 2 /d.o.f. calculations. However, the R in AA and R out AA for α max = 0.26 in both panels of Fig. 8 have smaller gaps than the PHENIX measurements, indicating over-isotropized high p T single inclusive pion spectra in the model, despite mean values of them are in agreement.   Table 2.
Nevertheless, we note in Section 3.1.2 that, Fig. 5 and Table 2 suggest even if strictly limit χ 2 /d.o.f. to be less than 1, α max maintains a non-negligible range which varies for different collisions. Using this flexibility, CUJET2.0 may create reaction plane dependent pion quenching patterns which are compatible with experiment measurements. To be rigorous, we first plot the χ 2 /d.o.f. vs α max for the R in AA and R out AA in Figure 8 using p T > 8 GeV, which p T range matches the choice in Section 3.1.2, the results are shown in Figure 9.
We find that in the left panel of Figure 9, for b = 2.4 fm central collisions at RHIC, in the CUJET 2.0 HTL scenario R in AA is best fitted by α max = 0.26 − 0.27, while R out AA is best fitted by α max = 0.28 − 0.30. In the right panel of Figure 9, for b = 7.5 fm semi-peripheral collisions at RHIC, R in AA is best fitted by α max = 0.23 − 0.25, while R out AA is best fitted by α max = 0.26 − 0.27.
If we choose α max = 0.26 R in AA and α max = 0.29 R out AA for b = 2.4 fm, α max = 0.23 R in AA and α max = 0.26 R out AA for b = 7.5 fm, the CUJET2.0 results are able to be perfectly consistent with RHIC R in AA and R out AA date. Since we have R AA = (R in AA + R out AA )/2, this set HENIX π 0 0-10% ∆φ=75-90° HENIX π 0 20-30% ∆φ=75-90°F of α max 's effectively generates R AA with α max = 0.275 at b = 2.4 fm and α max = 0.245 at b = 7.5 fm. Based on Figure 5 and Table 2, we see that the for average R AA with these α max 's in the CUJET2.0 HTL scenario is strictly less than 1.5 in all four collisions. It means this modest variation in α max is intrinsically allowed by our model without jeopardizing its consistency for averaged hadron R AA at a variety of collision configurations.
Notice among this set of α max parameters, for both b = 2.4 fm and b = 7.5 fm, the difference in α max for R in AA and R out AA is 0.03. And for either R in AA or R out AA , the variance in α max for b = 2.4 fm and b = 7.5 fm is 0.03 (which surprisingly coincides with the RHIC and LHC averaged α max gap discussed in Section 3.1.2). The maximumly 10% α max deviations in R in AA and R out AA , b = 2.4 fm and b = 7.5 fm imply that the anisotropic path averaged effective coupling strengths can be originated from local effects. For instance, in principle the strong coupling runs with not only energy but also local temperature, since in Eq. (2.1) only momentum appears, the local temperature dependence will effectively lead to azimuthally different α max . Further exploration of this topic is a work in progress. 10 .
After constrained the azimuthal anisotropy in CUJET2.0 at RHIC, our next step is to test the model consistency with R in AA and R out AA for Pb+Pb 2.76ATeV at LHC central and semi-peripheral collisions. Because of the absence of published reaction plane depen- 10 There is also a noteworthy ordering in αmax: at τ0, the length of the medium in b = 2.4 fm and b = 7.5 fm collisions along φ = 0 • and φ = 90 • direction can be approximately ordered as 7.5fm+0 • < 7.5fm+90 • ≈ 2.4fm + 0 • < 2.4fm + 90 • , and the best fit αmax in corresponding situations is 0.23 < 0.26 = 0.26 < 0.29. dent neutral pion or charged hadron suppression data at LHC, we turn to compare the CUJET2.0 results of v 2 with experimental elliptic flow measurements, which comparison is more sensitive to fluctuations from the theoretical point of view. And we will discuss this in detail in Section 3.3.2.

Elliptic flow
If presuming no fluctuations in the azimuthal plane, recall Eq. (3.8) suggests that R in AA and R out AA depend solely on even harmonics. If further drop higher order components assuming they have much smaller magnitude comparing to v 2 , we get 9) in this limit. And the elliptic flow v 2 follows from R in AA and R out AA via . (3.10) We compute in CUJET2.0 pion v 2 's using Eq.  Fig. 10, and corresponding ALICE [2], ATLAS [4], CMS [6], and PHENIX [9] data are compared. Fig. 10 shows that if varying solely the maximum coupling constant α max from 0.20 to 0.35 with 0.01 steps in the CUJET2.0 HTL scenario, none of the theoretical curves matches the single pion v 2 at both Au+Au 200AGeV and Pb+Pb 2.76ATeV, central and semi-peripheral collisions. Nevertheless, as already noted in Section 3.3.1, due to the non-negligible influence that anisotropy and heterogeneity/inhomogeneity have on the jetmedium interaction, local effects can alter the CUJET2.0 framework significantly. Under present circumstances that the strong coupling running has no simultaneous energy and local temperature dependence, these effects can effectively generate azimuthally anisotropic path averaged α max .
By choosing α max = 0.26 R in AA and α max = 0.29 R out AA for b = 2.4 fm, α max = 0.23 R in AA and α max = 0.26 R out AA for b = 7.5 fm in the CUJET2.0 HTL scenario, we effectively constrained our model at RHIC with assumed azimuthal α max anisotropy caused by possible local temperature field effects. The top panels of Fig. 10 show the consequential single pion v 2 's at RHIC central and semi-peripheral collisions are compatible with respective experimental measurements, as expected. More importantly, if extrapolate the same CUJET2.0 framework (with azimuthal dependence of α max ) to LHC and calculate the single particle v 2 using the same α max parameter set, theoretical results show even better agreements with ALICE, ATLAS and CMS data, in both central and semi-peripheral collisions (particularly the latter, cf. Table B.67).
The consistency, of the open heavy flavor and heavy flavor lepton's single particle azimuthal anisotropy calculated in the same CUJET2.0 framework (with azimuthal variations), with experimental measurements can shed light on the underlying physics for this azimuthal variation of path averaged coupling strength. This is a work in process.   A mini-summary of Section 3.3: by allowing the maximum coupling constant α max to vary by even less than 10% in respective reaction plane at central and semi-peripheral collisions, we can gain in CUJET2.0 the simultaneous compatibility with not only measurements of R in AA and R out AA at RHIC, but also effective R in AA and R out AA (measurements of (α in max ) and R out AA (α out max ) in the CUJET2.0 HTL scenario. Reference curves are shown in Fig. 4, Fig. 8, and Fig. 10. The p T > 8 GeV range is chosen for both RHIC R AA and LHC R AA for safer preservation of eikonal and soft approximations, p T > 8 GeV and p T > 12 GeV range is chosen for RHIC v 2 and LHC v 2 to avoid the avalanche region. The choice of α in max = 0.23, α out max = 0.26 significantly reduces the χ 2 /d.o.f. for v 2 at both RHIC and LHC, especially the latter one. By the mean time, this set of α max parameters maintains almost perfect agreement with both RHIC and LHC for azimuthally averaged R AA .

Thermalization time
At very early time of relativistic heavy-ion collisions, the matter created during the collision is characterized by extremely high energy densities. While it expands, and the gluon density decreases, the strength of interaction among particles increases, facilitating the thermalization process into quark gluon plasma. The time scale over which the thermalization process takes place is generally referred to as τ 0 and is assumed to be approximately equal to 0.5 ∼ 1 fm/c.
Once the plasma has reached the thermalized stage, the system can be approximated by a fluid and its evolution can be computed in the framework of relativistic hydrodynamics. The fundamental equations of energy-momentum and baryon number conservation are assumed to hold.
In CUJET2.0, the initial time τ 0 is set to be 0.6 fm/c to match the choice of VISH2+1, which generates 2+1D viscous hydrodynamical fields as a bulk background in the parton shower modification. Additionally, the temporal evolution of the system is parametrized in such a way that the density "seen" by the jet grows linearly until the full thermalization is reached, and decreases as 1/τ thereafter (cf. Eq. (2.15) and Appendix I).
In a recent paper by Song et al. [91], the authors found in VISHNU, where a UrQMD hadronic afterburner is coupled to VISH2+1, that to simultaneous fit RHIC and LHC particle production spectrum and all order collective flow harmonics, the initial time needs to be increased to 0.9 fm/c. We are therefore motivated to explore the effect of longer thermalization time, and to do so we modify the initial time τ 0 to be 0.9 fm/c at LHC conditions in CUJET2.0. The theoretical results compared with corresponding experimental data are shown in the bottom left and bottom right panel of Fig. 11. As illustrated in the bottom left panel of Fig. 11, with this increase in initial time, the suppression of pion at LHC slightly diminishes, which is understood by the fact that longer thermalization time in the linear scheme (cf. Appendix I) will create dilutely distributed scattering centers. Meanwhile, in the bottom right panel of Fig. 11, the best fit CUJET2.0 HTL α max for LHC Pb+Pb 2.76ATeV central collisions gains by 0.02 because of the prolonged τ 0 . Since the typical temperature range reached by LHC is higher than by RHIC, one would expect there is more reduction of density because of the τ 0 growth at LHC than at RHIC. Therefore, to maintain the same quenching magnitude with τ 0 = 0.9 GeV/c the increase of α max at RHIC will be smaller than LHC. This τ 0 effect can hence result in a lessened gap between the best fit value of CUJET2.0 HTL α max at RHIC and LHC central collisions.
However, if adding the hadronic afterburner to VISH2+1 as in VISHNU, the bulk evolution profile would be significantly different, which may or may not invalidate the above argument. Therefore a thorough test with a complete VISHNU code coupled to CUJET shall be conducted before drawing any positive conclusions.

Non-HTL scenario
Experimental R AA data favors maximum running strong coupling α max = 0.25−0.27 in the CUJET pQCD model with dynamical QCD medium and 1-loop HTL gluon propagator, i.e. a thermal gluon with mass m g = (f E ) gT 1 + n f /6 / √ 2 and f E = 1. However, lattice QCD calculation suggests this approximation breaks down in the non-perturbative region [92]. We therefore vary the HTL deformation parameters (f E , f M ) in CUJET2.0 from HTL (1, 0) to a non-HTL (2, 0) to explore this effect. The consequential pion R AA results are shown in the top left, top right, and χ 2 /d.o.f. are shown in the bottom right panel of Fig. 11.
With dynamical scattering potential and doubled thermal plasmon gluon mass in CU-JET2.0, a novel inversion of R AA for pion suppression pattern appears: when increase α max from 0.30 to 0.60, pion R AA first decreases then increases, with the strongest quenching occurs at α max =0.4, which is also the best fit to experimental data. This inversion is generated from a complex interplay between the reduced differential scattering cross section with an enlarged gluon mass and the variation of running coupling saturation scale Q min with α max (Q min = Λ QCD Exp{2π/9α max }). And at present stage there is no simple asymptotic analytical formula to track this non-HTL scenario in CUJET2.0.
On the other hand, concentrating on the χ 2 /d.o.f. fit of this non-HTL scenario whose gluon mass is doubled, in the bottom right panel of Fig. 11 the best fit non-HTL α max =0.4 still has a large χ 2 /d.o.f., implying the necessity to explore more complicated combination of HTL deformation parameters in CUJET2.0, e.g. f E = 2, f M = 0.6. In fact, a realistic non-HTL scenario would eventually have an effective running coupling and a effective scattering potential extracted from lattice QCD data on non-perturbative qq potential V (r, T ). And a comprehensive test of whether lattice QCD predicts the correct jet medium physics in the near T c non-perturbative region can be conducted in the CUJET framework. This is a work in progress.

Summary and Outlook
We presented in this paper the basic features of CUJET2.0 pQCD azimuthal jet flavor tomography model, which features DGLV opacity series with multi-scale running coupling and 2+1D viscous hydrodynamical fields, as well as the TG elastic energy loss; geometry, radiative energy loss and elastic energy loss fluctuations; and the convolutions of energy loss probablity distributions with initial production spectra and fragmentation functions.
We list our main results and conclusions with CUJET2.0, derived through focusing on CUJET2.0 itself, or through cross comparison with CUJET1.0, as follows: 3. The small value of best fit α max in the CUJET2.0 HTL scenario implies that longer jet path length in a transversely expanding medium overrides the reduction of density and contributes to overall enhanced quenching. 5. The robust crossing pattern of π, D, B, e − R AA 's is rigorously encoded in the flavor dependent energy loss structure of DGLV opacity expansion combined with TG elastic sector, and a transverse expanding medium has minor effect on this mass hierarchy.
6. Solutions to the "heavy quark puzzle" are intrinsically integrated in the framework of CUJET, through the inclusion of elastic energy loss, dynamical QCD medium effect, realistic geometry fluctuation, and energy loss fluctuations.
7. CUJET2.0 predicts a decisively less quenched B meson R AA which is well above D meson and pion at 5 GeV < p T < 15 GeV, as well as a critical alternation of R AA 's mass ordering from B > e − > D at p T < 15 GeV to e − > B > D at p T 25 GeV. 8. We explored the effect of allowing α max to vary minimally in different azimuthal directions. We find that even a less than 10% variation in the averaged coupling strength along in and out of plane paths reduces dramatically the χ 2 /d.o.f. for the v 2 fit at both RHIC and LHC, and at the same time, azimuthally averaged R AA results are consistently in agreement with experimental measurements at the level of χ 2 /d.o.f. < 1.5. The underlying physics responsible for this local effect is subject to feature explorations.
Future work of, improvement on and more thorough test with CUJET2.0 can be made in the following aspects: • Calculate the open heavy flavor and heavy flavor lepton's single particle azimuthal anisotropy in the azimuthal dependency included CUJET2.0 framework and compare with existent experimental results to explore the underlying physics responsible for azimuthal α max variations.
• Extrapolate an effective running coupling and an effective scattering potential from lattice QCD qq free energy [92] to replace the corresponding running strong coupling and differential cross section in CUJET2.0 to explore non-perturbative jet medium physics near Tc.
• Integrate Shuryak and Liao's model [93][94][95] of near T c enhancement of jet-medium interactions in CUJET2.0 to explore non-perturbative local effect. Enhancement of coupling strength originating from non-perturbative structures, created by the colorelectric jet passing a plasma of color-magnetic monopoles, dominate the near-Tc matter and could contribute to significant azimuthal variation of α max .
• Replace the presently assumed Poisson multiple gluon emission distribution with other radiative energy loss fluctuation patterns to explore the role that fluctuations play on jet quenching pattern and azimuthal anisotropy.
• Consider more carefully the effects that running scale variations have on hadron spectra and collective flow harmonics. Introduce effective running coupling for each interference term in the summation of current amplitude in the DGLV opacity expansion.
• Improve the structure and algorithm of the CUJET2.0 Monte-Carlo code to realize the calculation of jet-hardron correlation observables [96][97][98]. Constrain the hydro ambiguity and gain comprehensive information about the parton-medium interaction mechanism through azimuthal jet flavor tomography.

A Notations and conventions
We adopt the following notations and conventions throughout this paper, unless otherwise footnoted: When considering experimental observables in an A+A collision, for example hadron suppression and/or azimuthal anisotropy, z-axis is chosen along the beam direction, and azimuthal plane refers to the plane transverse to the beam axis. In particular, we always define p T as the transverse momentum perpendicular to the beam direction. Besides p T , physical quantities or concepts involving such coordinates include rapidity/pseudorapidity, bulk evolution profile, jet production distribution, etc.
When considering properties of a single jet, such as quenching in QGP, we choose the jet propagation direction as the "z-axis" and define the transverse plane accordingly. This coordinate system applies most importantly to the calculation of radiative and elastic energy loss. Take the scattering with a parton in a particular rapidity frame with a transverse momentum exchange q as an example, we denote q as four vector q µ , q as q's 3D space components, and q as q's transverse components with respect to the z-axis, i.e. the jet propagation direction. q ≡ q µ = (q 0 , q) = (q 0 , q z , q), and q ⊥ ≡ |q|. Similar notations are also applied to k µ , µ , J µ , etc.

B Review of fixed coupling DGLV
The DGLV opacity expansion [30][31][32]99] is a theory encompassing inelastic parton-medium interactions and describing gluon radiations in the pQCD framework. As in the WHDG [33] generalization of DGLV, CUJET1.0 and 2.0 supplement the radiative jet-medium interactions with TG [100] elastic collisional energy loss in the color medium (Section B.2).
The main computational task performed in CUJET via Monte-Carlo integration is to evaluate the number of radiated gluons per energy fraction dN g /dx for each initial jet production coordinates (x 0 ,n) 11 . After that the average inclusive gluon radiation distribution is calculated, fluctuations due to multiple gluon emission is computed via numerical convolution assuming uncorrelated Poisson ansatz, and the normalized radiation probability, P rad (∆E rad , E 0 ; x 0 ,n)) is evaluated via fast Fourier transform including delta function end point singularities (Section B.3.1). Normalized elastic energy loss probability, P el (∆E el , E 0 ; x 0 ,n)) is also computed with Gaussian multiple collision fluctuations (Section B.3.2). The final total energy loss probability distribution is the convolution of radiative and elastic sector, P rad ⊗ P el (Section B.4.1), it is then folded over the initial quark jet spectrum dN pp /d 2 p T dη (Section B.4.2). Finally CUJET averages over inital jet configurations via d 2 xdnT A (x + b/2)T A (x − b/2) and fragments jets into different flavor hadrons or leptons to compare with data (Section B.4.3).

B.1 Radiative energy loss in fixed coupling dynamical scattering DGLV
The GLV opacity expansion model was developed by Gyulassy, Levai and Vitev [30,31], built upon the foundations of the Gyulassy-Wang (GW) potential [29], it expresses the partonic energy loss as a series in powers of the opacity L/λ, where L indicates the size of the plasma and λ the mean free path of the parton. At nth order in opacity, one considers n scatterings between the parton and the medium. This is often referred to as a thin plasma approximation, which is valid for small values of opacity, as opposed to the thick plasma limit where multiple soft scatterings apply.
The interaction between parton and medium is modeled according to a GW [29] Debye screened potential with screening mass µ. µ is considered a fundamental property of the plasma along with the medium density ρ. Both of them are expressed as functions of the local temperature T in the system. GLV model includes the power-law tail of the scattering cross section, thus scatterings with large momentum transfer (hard) are taken into account.
There are three major kinematic assumptions made in the GLV theory: soft eikonal approximation, collinear radiation and discrete scattering centers. Details of them are listed below: • Eikonal approximation: both the parton energy E and the emitted gluon energy ω are much larger than the transverse momentum exchanged with the medium q ⊥ ≡ |q|: E q ⊥ and ω q ⊥ . Soft approximation assumes ω E.
• Collinear radiation: Gluons are emitted at small angles with respect to the parent parton: ω k ⊥ , where k ⊥ ≡ |k| represents the transverse momentum of the gluon.
Under the soft eikonal approximation, the parent parton has sufficiently high energy such that its path is approximately straight. The gluon, which is radiated at small angles, does not carry away a significant portion of the original parton energy, and consequently the jet energy is not dynamically updated during the multiple scattering process. The most remarkable features of the gluon radiation spectra in the GLV opacity expansion theory comprise the interference effects between production (vertex) radiation and induced radiation, and the interference effects among subsequent scatterings of the radiated gluon in the plasma (quantum cascade). These effects lead to an expression for the double-differential gluon multiplicity distribution in x (fractional gluon energy) and k ⊥ (gluon transverse momentum), which is later integrated to give rise to the energy loss of the parent parton assuming no further exchange of energy with the medium takes place.
An extension of the GLV model to include massive quarks kinematic effects as well as plasmon mass for the gluons was developed by Djordjevic and Gyulassy in [32] (DGLV). The full derivations can be found in the original papers [30][31][32]. Here we will only show the main results and provide their physical interpretations.
In the soft eikonal approximation used to derive DGLV, the incoming jet, gluon and exchanged four momenta read where parenthesis and square brackets denotes Minkowski spacetime and light-cone coordinates respectively. In the above expressions we have suppressed the effective gluon plasmon mass m g = µ/ √ 2 as well as the the parton mass M . In the static scattering center approximation q 0 ∼ q z |q|. The gluon fractional energy x E and fractional plus-momentum x + are related via In the pure collinear limit, they coincide. Corrections need to be made for finite emission angles, which involve variations in the upper kinematic integration limit and the Jacobian of the transformation x + → x E : A detailed discussion about this issue can be found in Appendix D. The double-differential gluon multiplicity distribution in x + and k, for DGLV opacity order n = 1, i.e. the case that the hard parton scatters one single time with weakly-coupled static quasi-particles in the deconfined thermal medium, is given by with the normalized modular squared scattering potential in a static QCD medium being |v(q)| 2 = µ 2 π(q 2 + µ 2 ) 2 . (B.5) Here C R is the quadratic Casimir of the jet (C F = 4/3 for quark jets, C A = 3 for gluon jets), α s = g 2 /4π is the strong coupling constant, and L is the jet path length. Note that the opacity is written in terms of the gluon rather than the jet mean free path, λ g , because of a simplification in the color algebra known as "color triviality" [31]. χ 2 = M 2 x 2 + +m 2 g (1−x + ) controls the "dead cone" (cf. Appendix D) and LPM destructive interference effects due to both the finite quark current, M , and the thermal gluon mass m g = µ(T )/ √ 2. ∆z 1 = z 1 −z 0 represents the distance between the scattering points z 1 and z 0 (production vertex).
The DGLV all-orders result expresses an arbitrary opacity order in a closed form, for an arbitrary collision probability along the jet path. It is applicable for both coherent and incoherent geometries (cf. Appendix C). The gluon multiplicity distribution takes the form: and the normalized modular squared potential for ith static scattering center is

(B.7)
This interacting potential has the form of the Debye screened Gyulassy-Wang potential (cf. [29]), and forward scattering unitarity correction δ 2 (q i ) is subtracted from it. ∆z k = z k − z k−1 represents the distance between adjacent scattering points z k and z k−1 . The kinematic current amplitudes in Eq. (B.6) are modified versions of the Hard, Gluon-Bertsch and Cascade terms in GLV theory (cf. [31]), for finite masses: which regulates the LPM phase of color currents. In principle, the opacity series should be calculated to sufficiently high order to generate the "exact" prediction of experimental observables. However, the numerical power required to drive the computation at such high levels of precision might prove to be insufficient to calculate more complex observables than the simple energy loss for one specific plasma setup. A limitation of this kind would indisputably hinder the capabilities of our algorithm and limit its predictive power. Thus one has to quantify the error introduced by eventually limiting the computations to lower orders in opacity, and in Appendix C we conduct such a systematic study of the convergence of the DGLV opacity series. There we demonstrate that despite another set of observables might scale differently with the opacity, for the purpose of computing the energy loss of different quark flavors, truncating the series at first order already does not add a relevant source of systematic uncertainty. Therefore, unless otherwise stated, all the calculations presented in this paper will be at n = 1 order.
After computed the differential gluon multiplicity distribution according to Eq. (B.4), we then want to integrate over k to get the gluon radiation spectrum and at higher orders, Here the elastic cross section for gluon σ el (z) can be expanded into gluon-quark and gluongluon terms, i.e.
where µ 2 (z) = g 2 T (z) 2 (1 + n f /6) = 4πα s T (z) 2 (1 + n f /6) is the squared local HTL color electric Debye screening mass in a plasma with number of flavors n f and local temperature T (z) ∝ ρ(z) 1/3 along the jet path z through the plasma. Assuming ideal gas conditions, from the boson/fermions statistics we obtain the number density of quark ρ q (z) and gluon ρ g (z) is respectively Combining Eq. (B.14) and (B.15), we have To get the DGLV gluon radiation spectrum at the first order in opacity, combine Eq. (B.4), (B.10), (B.12), and (B.16) together, we get

(B.17)
This is the DGLV n = 1 kernal in fixed coupling CUJET. Recall C R is the quadratic Casimir of the jet (C F = 4/3 for quark jets, C A = 3 for gluon jets). Local χ 2 (z) = M 2 x 2 + +m 2 g (z)(1−x + ), with M being the mass of the parton, gluon mass m g (z) = µ(z)/ √ 2, and Debye mass µ 2 (z) = g 2 T (z) 2 (1 + n f /6) = 4πα s T (z) 2 (1 + n f /6) . (B.18) For a static QCD medium, |ṽ(q)| 2 is defined via 12 The total energy ∆E carried away by the emitted gluons is obtained by integrating the radiation spectrum, Eq. (B.17). Assuming no further interaction between jet and medium, this can readily been interpreted as the energy loss that the jet suffers when propagates through a hot deconfined plasma If there are no kinematic boundaries on the dq and dk integrations, for a plasma with fixed size L and constant gluon mean free path λ g , a straightforward analytic computation for this first order in opacity leads to the asymptotic result We immediately notice the L 2 dependence of the energy loss, this is the characteristic of the LPM region, which differs from the incoherent limit whose dependence on L is linear. However, as discussed previously, kinematic boundaries exist not only for the k ⊥  Eq. (B.17) comes from the z coordinates (Eq. (B.11)), and strictly speaking, i.e. the radiated gluon spectrum for a quark/gluon jet with mass M and energy E(≡ E 0 ) created at x 0 position in the transverse plane along azimuthal angle φ has explicit dependence on the strong coupling constant α s , jet path length L and number of quarkonic flavors n f 13 . In CUJET calculation, we take into account fluctuations of the geometry by cutoff the dτ integral at τ M AX , with T (z)| τ M AX = T f . By doing so the L dependence turns into a fragmentation temperature T f dependence. We leave the discussion of the systematic uncertainties associated with varying T f and n f to Appendix E. There we show the variation of jet quenching spectra originating from both of them can be fully absorbed into a simple rescaling of α s . This makes α s 14 the only free parameter of CUJET calculation in a static QCD medium, in which case the radiative energy loss reads However, in such a deconfined medium consisting of randomly distributed static scattering centers, the collisional energy loss is exactly zero. This contradicts recent calculations which show elastic collisional contribution is important and comparable to the radiative energy loss [33]. Therefore, the natural improvement on the DGLV opacity expansion is to consider a dynamically screened QCD medium. The inclusion of these dynamical effects is achieved by computing the scattering QCD diagrams in a finite temperature field theory framework, using Hard Thermal Loop resumed propagators for all gluons. The quark gluon plasma is assumed to be thermalized at temperature T and has zero baryon density. Details of the computations can be found in original papers [45,[101][102][103][104]. 13 As well as other plasma parameters such as the formation time and thermalization scheme, which are discussed extensively in Section 2.2, Section 4.1 and Appendix I. 14 In the case of running coupling CUJET which is discussed in Section 2.1, this parameter is the maximum strong coupling constant αmax.
The dynamical QCD medium brings two major corrections to the DGLV radiated gluon number distribution Eq. (B.17): firstly, the effective dynamical mean free path for gluon λ dyn , which is defined as λ −1 dyn ≡ 3α s T , will replace its static counterpart λ g (≡ λ stat ) in Eq. (B.4). According to Eq. (B.16), they are related via with n f the number of effective quark flavors in equilibrium with the gluons in the plasma. However, in CUJET calculations we degenerate this mean free path effect on δE/E by a rescaling of the effective strong coupling constant, because the coefficient c(n f ) varies from c(0) = 0.73 to c(∞) = 1.09, and does not contribute much to the energy loss compared to the magnetically enhanced potential which will soon be discussed. Secondly and most importantly, the dynamical recoiling of color electric scattering centers induces an effective color magnetic screening mass which is smaller than the Debye mass, and the interaction potential The implications of these changes are profound: the absence of the µ 2 screening for soft momenta exchanges q makes the potential diverges and the mean free path vanishes. In the limit of q → 0, each individual Feynman diagram diverges logarithmically. These singularities however cancel out after all the contributing diagrams to the energy loss are summed over, making the gluon multiplicity finite. The combined effect of the enhanced cross section and reduced mean free path contributes to a remarkable increase for the magnitude of total energy loss and the ratio of heavy to light quark energy loss in the dynamical framework, systematic studies of this effect can be found in [103] and [105].
In the CUJET model, motivated by lattice QCD qq potential data, we introduce an effective interaction potential with deformation parameters (f E , f M ), it reads 15 .

(B.26)
Here the HTL deformation parameters (f E , f M ) are used to vary the chromo-electric and chromo-magnetic screening scales relative to HTL. In principle, HTL deformations could also change m g (T ). The default HTL plasma is (1, 0), but we also consider a deformed (2, 0) non-HTL plasma model which will be discussed in Section 4.2. (α s 16 , f E , f M ) are therefore our main model space control parameters.

B.2 Elastic energy loss
The assumption that pQCD elastic energy loss is negligible compared to radiative one is questionable. In [107,108], the authors found that radiative and elastic average energy losses for heavy quarks were in fact comparable over a very wide kinematic range accessible at the RHIC. In [33], the authors confirm these previous findings and extend them to the light quark sector, showing that elastic contributions to the total energy loss can be of the same order of magnitude of radiative ones. It is then clear that quantitative tomographic predictions cannot ignore such large contributions to jet quenching, and elastic effects need to be included in CUJET as well.
We use Thoma-Gyulassy (TG) model [100] in our calculation of the elastic energy loss. Their work was based on Bjorken's estimation of elastic energy loss in QGP (cf. appendix F). By using the hard thermal loop gluon propagators to provide a more natural infrared regulator, the TG computation leads to the following leading log result: Where x is the jet path. For ultra-relativistic particles, the velocity v can be approximated to 1 and the v-dependent factor in parenthesis becomes approximately 1. The integral over k is infrared finite due to the Debye screening mass in the denominator, but a maximal momentum k max must be set in order to screen the otherwise ultraviolet divergent logarithm. Assuming that the maximal momentum transfer comes from forward scattering against target particles with average momenta q ≈ 2T is much smaller than the projectile momentum, the value of We immediately see that this model yields a result very similar to the Bjorken computation, Eq. (F.6), i.e.
with a different Coulomb log that reflects the more natural cutoffs which are being used now: For the elastic energy loss sector in CUJET, assuming jet travels at the speed of light, combining Eq. (B.27), (B.28) and (B.29), we get the following equation to account for collisional effects: Here C R is the quadratic Casimir of the jet (C F = 4/3 for quark jets, C A = 3 for gluon jets). Note similar to radiative energy loss, the n f and T f dependence of elastic energy loss is absorbed into the α s degree of freedom. The recursive scheme is stopped when E(z) drops below M and returns maximum energy loss ∆E max = E 0 − M , even if τ does not reach τ max . If local temperature T (z)| τ =τ 0 = 0 for some τ 0 , the scheme will skip computing dE/dx and keep evolving τ . Despite its improvement over the Bjorken result, the TG model leaves the ultraviolet region unbounded, because the classical calculation has no knowledge about the particle nature of the medium and particle recoil, which becomes important when the momentum transfer q is large. The hard momentum transfer contribution is more naturally taken into account by Braaten and Thoma in [109,110], but relevant analysis shows that the differences in practical applications are almost negligible. We Numerical results for elastic energy loss, especially its effects on the ratio of light to heavy quark suppression magnitude, is discussed in Appendix F.

B.3 Fluctuations
The radiative energy loss calculated from Eq. (B.17)(B.26)(B.23) and elastic energy loss calculated from Eq. (B.30)(B.31) both perform full jet path integration 17 . The non-uniform medium's fluctuating geometry due to expanding and cooling is properly embedded in the local plasma density ρ(z) and the cutoff fragmentation temperature T f . This jet path integration provides a platform to quantify the effects of complicated heavy ion collision configurations in predicting experimental observables in the pQCD framework.
Besides the geometry, fluctuations originating from multiple gluon emissions in the radiative sector and multiple partonic collisions in the elastic sector also play important roles in jet quenching, and they may significantly influence the results of hadron multiplicity and azimuthal flow. We dedicate this section to introduce the quantification and computation of them in CUJET. The simplest assumption for multiple gluon emission is the Poisson ansatz, where the number of emitted gluons follows a Poisson distribution, with the mean number N g given by the integral of the gluon emission spectrum The gluon radiation can be thought of as a stochastic event, and it makes sense to speak of a probability distribution P rad ( ) of radiating a certain amount of energy ≡ ∆E rad /E 0 : where the maximum energy loss ratio max = 1 − M/E 0 . For simplicity, in the discussion below we suppress the n = 1 superscript of N n=1 g and E subscript of x E and write them as N g and x respectively. The probability distribution Eq. (B.34) is split into three components: The first term corresponds to the probability of zero radiation, P null r = e −N g . The second term is given by and P n+1 ( ) = 1 n + 1 Fourier transform back, we have Practically, the numerical evaluation of Eq. (B.42) uses finite discrete k i and x j series, for example, k i = −1000 + i (i = 0, 1, · · · , 2000) and x j = jσ (j = 0, 1, · · · , σ −1 ; σ = 0.0025), meaning the Fourier transform in the exp(...) of Eq. (B.42) is in fact The dNg dx (x) itself is fluctuating because of limited computing power to implement Monte-Carlo iterations. At large |k i |, this fluctuation is worsened with the highly oscillating e ik i x j , and will generate unphysical variations in P r ( ). However, if take the dk e −ik in Eq. (B.42) into account, one sees components with larger |k| will have less weight in the evaluation of P r ( ). Therefore, we smoothfy the exp(...) in Eq. (B.42) by adding a Gaussian smoother with proper width, put more weight on small k Fourier components, and modify Eq. (B.42) to with x j = jσ (j = 0, 1, · · · , σ −1 ). In CUJET, P r (0) = 0, we use Eq. (B.44) to calculate P r ( ) in the range of 0 < ≤ max , as well as in the range of max < ≤ leak = 1.75 for numerical purposes. The third and last term in Eq. (B.34) represents instead the probability of total quenching. In the soft approximation, the radiated energy ω is assumed much smaller than the initial jet energy E, and x 1. Consequently, the energy of the outgoing parton E is approximately equal to E. When the {x n } are integrated up to the kinematic limit x n = 1, a "leakage" error into the unphysical region P r ( > max ) = 0 occurs, and this error is calculated in P f ull r = ∞ max d P r ( ) 18 . For the normalization of P rad ( ) we keep the weight of the physical zero quenching probability P null unchanged, and rescale the probability distribution as follows: firstly, we calculate the norm N rad from N rad = max 0 P rad ( ). When doing this integral, the Delta functions at both boundaries are included. Secondly, we rescale the complete P rad ( ) according to P rad ( ) → 1−e −Ng N rad P rad ( ). And finally we replace the coefficient of δ( ) in P rad ( ) with zero radiation probability, i.e. 1−e −Ng N rad P null → e −Ng . Through this procedure we maintain max 0 P rad ( ) = 1, and the δ( ) at the = 0 boundary has weight e −Ng . If N g = 0, P rad ( ) = δ( ).
The effects of multiple gluon emission on the ratio of light to heavy quark energy loss is a topic of Appendix G. Note P rad ( ) inherits all the jet production coordinates, parton mass and energy, and model parameter dependencies from dNg dx . And we write down explicitly those dependencies as: (B.45)

B.3.2 Elastic energy loss fluctuation
Fluctuations of the elastic energy loss around the mean were addressed in [33] and [111]. Using a framework generally applied to diffusive processes that are characterized by a large number of soft collisions, the probability distribution to lose the collisional energy ≡ ∆E el /E 0 is represented by a Gaussian centered around the average ∆E el , with variance σ 2 = 2T /E 0 . Here ≡ ∆E el /E 0 , and the average elastic energy loss ∆E el is calculated according to Eq. (B.31), The average temperature along the jet path is The collisional energy loss probability distribution reads Here and here 0 ≤ ≤ max . For numerical purposes we also calculate P e ( ) in the range of max ≤ ≤ leak = 1.75 according to Eq. (B.51). Note integrate P el ( ) over 0 ≤ ≤ max automatically gives unity. The rearrangement of Eq. (B.49) provides great conveniences for the convolution of radiative and elastic energy loss probability distributions, which will be studied in the following section. Similar to Section B.3.1, inherited from ∆E el and N c , the elastic energy loss probability distribution has jet production coordinates, parton mass and energy, and model parameter dependency. We write down all those dependencies for P el ( ) as: (B.52)

B.4 Convolutions
In the CUJET framework, after calculated the radiative energy loss probability distribution P rad ( ) from Eq. (B.34) and elastic energy loss probability distribution P rad ( ) from Eq. (B.49), we convolute the their contributions to get the total energy loss probability distribution P tot ( ) (Section B.4.1). Then integrate P tot ( = ∆E tot /E 0 , E 0 , x 0 , φ; M ) with the pQCD p+p parton (M) spectrum and binary distribution to get the quenched A+A parton (M) spectrum (Section B.4.2). Finally, we fragment this parton spectrum to get the transverse momentum and azimuthal angle dependent production spectrum for inclusive π, D, B and e − in A+A collisions (Section B.4.3).

B.4.1 Total energy loss probability distribution
To get the total energy loss probability distribution P tot ( ), we convolute the radiative sector P rad ( )(Eq. (B.34)) and the elastic sector P el ( ) (cf. Eq. (B.49)): Technically, in CUJET, when computing the convolution for total suppression, we keep the δ function at 0 in each sector while let P r ( ) and P e ( ) spread over 0 ≤ ≤ leak = 1.75, then absorb the convoluted leak to the δ function at max . Step by step, first we rewrite P rad ( ) and P el ( ) as with P r ( ) and P e ( ) calculated over 0 ≤ ≤ leak = 1.75. Then multiply them both according to Eq. (B.53), we get here 0 ≤ ≤ max . The normalization of P tot ( ) is conducted in the normal way to ensure max 0 d P tot ( ) = 1. Strictly speaking, P tot ( ) depends on other parameters such as parton masses which are inherited from P rad ( ) and P el ( ), and explicitly, For simplicity, throughout the paper we will suppress the "tot" subscript.

B.4.2 Jet quenching spectrum
CUJET computes the quenched partonic AA spectrum by convolute the total energy loss probability distribution P ( ) calculated from Eq. (B.58) and (B.57) with partonic production cross section in p+p collisions. This is a critical improvement of CUJET over its predecessor WHDG [33], which assumes instead a simple and slowly varying power law distribution for the p+p spectra (spectral index approximation) and makes considerable simplifications in the computation of the nuclear modification factor. Given the sensitivity of the results to the details of the production cross sections, and the complex interplay between the latter and the energy loss mechanism, it is essential that no approximations are carried out in this delicate step of the computation. The partonic pp spectra for CUJET are generated from pQCD calculations. For the light sector, production is based on a leading order (LO) calculation scaled by a simple K-factor and computed from the LO pQCD CTEQ5 code of X.N. Wang [112]. For the heavy jet sector, both next-to-leading order [113] and fixed-order plus next-to-leading-log (FONLL) [114,115] computations are used. In addition to including the full NLO result [116][117][118], the FONLL calculation re-sums large perturbative terms with next-to-leading logarithmic accuracy [119]. Details about the partonic spectra used in CUJET can be found in Appendix H.
It is clear at this point, that the input and output of CUJET are: the model is given a parametrization of the plasma and a jet spectrum, and it returns a quenched spectrum after computing the energy loss of the jets in the medium. In this process, no approximations are made: each jet is evolved individually and its final momentum, or better momentum probability distribution, is stored along with the direction it came from (angular distribution). This is how CUJET performs the computation of quenched partonic spectra: 1. The algorithm starts from a jet created at x 0 in the azimuthal plane (with respect to the beam axis) with azimuthal angle φ and mass M . The distribution of jets in the transverse plane in A+A collisions is given by ρ binary (cf. Section 2.2.1). The initial transverse momentum probability distribution P 0 (p i ) 19 of the partons is proportional to the production cross section: Here dσ pp→q /dp i represents a generic p+p partonic production spectrum. A range of discrete transverse momenta [p min i , p max i ] needs to be defined for the numerical computation.

For each initial transverse momentum
Eq. (B.58)): 3. Once all the {p i } in the range specified have been computed, the {P ( ; p i )} are converted into a two-dimensional distribution map that represents the probability of a jet with initial transverse momentum p i to leave the plasma with final transverse momentum p f : The normalization is such that which is automatically ensured by max 0 d P ( ; p i ) = 1. In Eq. (B.63) we dropped the explicit dependence on the jet coordinates x 0 and φ.
4. CUJET then integrates over the production spectrum, Eq. (B.60), to obtain the "quenched" partonic p+p spectra dσ pp→q dp f dφ : 5. At last, the quenched partonic A+A spectra as a function of the observed transverse momentum p f (≡ p T ) and azimuthal angle φ are obtained by integrating over the jet transverse distribution : dσ AA→q dp f dφ

B.4.3 Fragmentation functions
Partonic spectra can provide useful information about jet quenching mechanism, nevertheless, comparison with data can only be carried out at the hadronic level. The quenched partonic spectra, Eq. (B.67), need to be convoluted with a set of fragmentation functions (FF's). The process that leads to the fragmentation of partons in the medium is not theoretically well understood, especially for heavy quarks: dissociation and recombination theories [120,121] assume that heavy D and B mesons can be formed within the plasma and lose additional energy through collisional dissociation, in a similar fashion to what has been suggested for heavy quarkonium states [122]. This, however, seems to contradict more recent lattice results [123], which indicate the complete melting of open heavy flavors occurs at temperature T 220 MeV.
Since we are dealing with high p T partons, hadronization via recombination processes is suppressed compared to fragmentation. We will assume that fragmentation takes place in vacuum, on a hypersurface parametrized by µ(x, τ f ) = Λ QCD , and our results do not show a particular sensitivity on the precise choice of fragmentation temperature T f (cf. Appendix E).
The convolution of partonic spectra over appropriate FF's takes the form (B.68) Here D i→h (y; Q) represents the probability that a parton i fragments into a hadron h which carries a fraction y of the parton energy. Q is the scale at which the FF is evaluated, here it is given by the energy of the parton. Eq. (B.68) is summed over all species i that fragment into h. For light quarks and gluons fragmenting into pions, we use leading order KKP functions [124]. For heavy quarks fragmenting into D and B mesons (c → D and b → B), we use instead the Peterson [125] function with c = 0.06 and b = 0.006, as done also in [126]. While the Peterson FF does not couple well with the FONLL production cross section [127], it was shown in [126] that similar results are produced anyway even using a more accurate fragmentation description. Finally for the decay of the heavy mesons into nonphotonic electrons (c → D → e and b → B → e), we use the same functions as in [127]. The secondary decay D → B → e is also accounted for.

C.1 Uncorrelated geometry
The DGLV opacity expansion integrand in Eq. (B.6) is a function of the distance between scattering centers ∆z k . Their distribution is connected to the mean free path λ, which in the case of a non-uniform plasma is itself a function of z, i.e. λ(z). To see the average over the target coordinates in a smooth background more clearly, we write The λ(z i ) dependence significantly complicates the DGLV integral, especially at higher order in opacity n. To study the behavior of higher order opacities more efficiently, we apply an "uncorrelated geometry" for quick DGLV evaluations. In this configuration, we neglect the interconnection between the location of the scattering centers and the mean free path, as well as the mutual dependence of the spacing of collisions.
The simplest medium one can study is an uncorrelated brick of uniform density, constant temperature T and limited length L. This configuration is realized by changing Eq. (C.1) to where the normalized distribution for scattering centers ρ(z 1 , · · · , z n ) = n! L n θ(L − z n )θ(z n − z n−1 ) · · · θ(z 2 − z 1 )θ(z 1 − z 0 ) .

(C.3)
Note z 0 = 0 is the position of the production vertex. Since the 0 and L boundaries are already contained in the integration limits, one can drop either θ(L − z n ) or θ(z 1 − z 0 ) or both in the above equation. Due to the LPM phase oscillation in Eq. (B.6), a pure brick would easily create fluctuating gluon radiation spectra. We thus make a further generalization of the brick geometry by assuming an exponential distribution of scattering centers, i.e.ρ (z 1 , · · · , z n ) = In order to fix L e (n), we require that z k − z 0 = kL/(n + 1). This constrains L e (n) = L/(n + 1).
In the discussion within this paper, unless otherwise stated, if referring to a "brick", we mean an uncorrelated brick with an exponential distribution, defined by Eq. (C.2) and (C.4).

C.2 Convergence of DGLV
The DGLV opacity series approach builds upon the Bertch-Gunion (GB) incoherent radiation and includes multiple coherent scatterings, interference with the production vertex radiation and gluon cascading. The LPM effect and the interplay between the cosine factors in the DGLV integral determine how fast the series converges to its asymptotic limit.
To understand the convergence of the DGLV opacity series quantitatively, we need to study details about the LPM phases, and the interplay between the formation time τ f and mean free path λ. Generally speaking, if denote the size of the medium as L (L > λ), the interference (coherent) effects are dominant in the region λ < τ f < L, whereas the Gunion-Bertsch incoherent limit (cf. Eq. (C.8)) and factorization limit is obtained in the region τ f < λ < L and λ < L < τ f respectively.
Formation time τ f is the time gluon spends to become on-shell, it is approximately equal to with ω = x E E the energy of the radiated gluon. (In principle, only in the strict collinear limit x E and x + coincide. For simplicity, we set x + = x E ≡ x in this section. And we will discuss this issue in Appendix D.) If there are many momentum kicks from the medium within a coherence length, then q → i q i ; however, for a qualitative estimate, we can assume k q and τ f ≈ 2ω/k 2 . In reality, the interplay between k and q makes the estimation of the real formation time difficult, and once the mass of a heavy quark is taken into account, the χ 2 = M 2 x 2 + m 2 g (1 − x) factor starts playing a relevant role by reducing the formation time and by pushing the radiation back into the incoherent regime.
Mean free path λ plays an important role in determining the effects of coherence physics. In the uncorrelated geometry assumption, the relation between λ and the distribution of scattering centers looses, while in reality they are interconnected. Coherence effects are dominant when λ < τ f , and are analytically determined by the magnitude of the LPM phases ∆z/τ f : larger phases are responsible for the oscillatory behavior characteristic of the incoherent limit, while smaller phases cause an approximate cancellation among the cosine terms, typical result of coherence physics. For instance, for n = 1 and large formation times, the LPM term cos(∆z/τ f ) approaches unity, giving rise to a neat cancellation.
In order to understand how the convergence of the opacity series is related to the coherent or incoherent radiation regime, we first compare the DGLV n = 1 result with the Gunion-Bertsch incoherent limit: .
At first order, the opacity series only includes interference effects between the creation and the induced radiation vertex. By plotting the gluon transverse momentum distribution for different brick sizes L, we demonstrate in Fig. 12 the suppression of the induced radiation due to such interference effects. Comparing solid and dashed black curve in Fig. 12, we see that this coherence effect vanishes when the size of the medium is large, as expected. To move on, we add higher order corrections to the results, shown as purple, magenta and pink curves in Fig. 12. Regardless of all these higher order modifications, the dominant contribution to the suppression of the induced radiation still comes from the the n = 1 term.
In the left panel of Fig. 12, we observe that for L = 5 fm, at n = L/λ ≈ 5 the opacity series already converges to its asymptotic value, making further corrections negligible. This can be understood by assuming the probability of hitting a given number of scattering centers follows a Poisson distribution, its average equals the opacity, and we would expect the GLV series to peak around n = L/λ.
But this convergence is only valid for short path lengths, because the interference with formation radiation is the dominant effect, on top of which the corrections due to multiple scatterings in the medium are small. As L increases, this is no longer true and the resummed result is expected to asymptotically converge to the multiple soft scattering limit, as shown in the right panel of Fig. 12.
The above analysis is restricted to a particular choice of E = 50 GeV and ω = 5 GeV. To be general, we perform a systematic study of the properties of the series, by analyzing its convergence for several coherent and incoherent regimes, varying all relevant parameters. Our goal is to understand if there is an optimal order at which the series can be truncated for most of the practical needs, and quantify the error which is eventually made.
For different sets of (E, ω and L), we compute in Fig. 13 the radiation spectrum up to ninth order in opacity.
As shown in Fig. 13, the coherent radiation is associated with faster convergence: the large formation time suppresses the magnitude of the LPM phases, leading to an approxi- mate cancellation of the cosine terms in (B.6). On the other hand, the oscillatory behavior of incoherent emission results in slower convergence of the opacity series. Interestingly, the transverse momentum distribution seems to depend significantly on the gluon energy ω, rather than the original jet energy E. Finally, the convergence is improved by the reduced medium size L, as expected from the assumption of Poisson distributed scatterings.
For completeness, in Fig. 14, we show the same simulation for a heavy quark jet in a plasma of thickness L = 5 fm: the convergence rate is almost unchanged despite the dependence of the gluon formation time on the mass of the incoming quark, manifested in the term χ 2 = M 2 x 2 + m 2 g (1 − x). The increase of M is in fact compensated by the small value of x for ω 1 GeV. However, compare with the light jet results of Fig. 13, now the suppression of the radiated gluon multiplicity depends jet energy E as well. This is because of the presence of a χ 2 term in the denominator of the DGLV radiation antenna.
Combine all the above discussions, we can conclude that except for a limited region of phase-space -when the emission mechanism is clearly incoherent -a satisfactory result can already be obtained by truncating the expansion at third order. Furthermore, when averaging over all possible path lengths in a realistic nuclear collision, 2 L 5 fm, even the first order in opacity might be regarded as a good approximation to the series (Fig.  13).

D Transverse momentum distribution of radiated gluon
In this section we concentrate on the dependence of the gluon spectrum on the transverse momentum k ⊥ , and we will see that the k ⊥ functions non-trivially on various aspects of the radiative jet energy loss.
In Appendix C we observed that ω determines how fast the series converges to its asymptotic limit. However, we did not take into account the fact that the convergence appears to be faster for larger values of the transverse momentum k ⊥ , despite the shorter formation time proportional to 1/k 2 . For instance, in Fig. 13, with E = 100 GeV, ω = 5 GeV and L = 5 fm, the first order is already a good approximation for k ≥ 4 GeV, whereas between 2 ≤ k ≤ 4 GeV the fifth order is needed; below k = 2 GeV, only n = 7 is a good  Fig. 13. DGLV opacity series calculated up to n = 1, n = 3, n = 5, n = 7 and n = 9 is plotted as black, blue, green, orange and red solid curve. The effect of the quark mass in the expression for the formation time, which intuitively would slow the convergence of the series. However, this effect is balanced by the x dependence of χ 2 = M 2 x 2 + m 2 g (1 − x): for small x, the results do not differ much from their light quark jet counterpart. On the contrary to light quark jet results, we observe a remarkable splitting between radiation distributions with same gluon energy ω but different heavy jet energy E (top-right and bottom-left figures), due to the presence of the same x dependent χ 2 in the denominator of the antenna term in Eq. (C.6), which further suppresses radiation at large x.
approximation to the series. The reason can be found in the radiation antenna term of Eq. (C.6), which determines the shape of the momentum distribution: its 1/k 3 ∼ 1/k 4 asymptotic behavior suppresses high momentum corrections and dwarfs the contribution of higher orders in opacity. This effect is very similar to what we observed for heavy quarks jets in Fig. 14, where the large contribution of χ 2 in the denominator of the antenna term offsets the increased oscillatory behavior of the integral due to shorter gluon formation times.
The DGLV opacity expansion has the ability to interpolate between single hard scattering and multiple soft scattering limit. The latter is derived assuming the radiated gluon experiences Gaussian diffusion in the transverse momentum space: for small gluon emission angles, i.e. k ⊥ qL, the momentum distribution derived from the DGLV series approaches this limit. This effect is shown in Fig. 15 for a heavy quark jet. We see that at small k ⊥ the series converges to the multiple soft scattering limit quickly. However for large k ⊥ , i.e. large angle radiation which is treated poorly in the multiple soft scattering approximation, differs from ASW, the DGLV opacity expansion includes the hard power-law Landau tails of the radiation, reproduces the gluon multiplicity more accurately.

D.1 Integration and kinematic limits
The hard 1/k 3 ∼ 1/k 4 tails of the DGLV distribution offer a relevant contribution to the total emitted radiation and become a source of concern once finite kinematic limits are taken into account. If the integrand in (C.6) were exact, the result would vanish for unphysical values of k ⊥ , and there is no need for worrying about integration limits. In reality, however, the model is derived assuming collinear approximation (k ⊥ ω), therefore kinematic limits need to be imposed to enforce physicality. The integral, for consistency, should not be sensitive to the particular choice of UV k ⊥ cutoffs, but given the hard tails of the distribution, we will see that this is not always going to be the case.
The choice of upper bounds in the k ⊥ integration depends on the particular interpretation of x in the expression for the gluon energy ω = xE: x as the fractional energy carried away by the radiated gluon (x ≡ x E , ω = x E E), or x as the fraction of plus-momentum in light-cone coordinates, in which case x ≡ x + and ω ≈ x + E + /2 20 . In the strictly collinear limit in which the DGLV integral is derived, the two definitions coincide: Equation (D.1) can be easily derived by writing explicitly the gluon four-momentum in Minkowski and light-cone coordinates, denoted respectively by parenthesis and square brackets: Depending on the interpretation of x, the upper kinematic limit on k ⊥ will vary: in the case of x + , in order to ensure forward gluon emission we need to set k M AX ⊥ = x + E + , whereas in the case of x E , to keep k ⊥ real we must set k M AX ⊥ ≈ xE sin θ, where θ is the angle between the radiated gluon and the propagating parton 21 . In Fig. 16, we plot the k ⊥ integrated gluon number distribution x dNg dx , for both interpretations of x and two different cutoff angles θ; to compare apples to apples, we add the Jacobian of the transformation x + → x E to the x + curve and integrate up to k M AX ⊥ = xE sin θ: The differences are notable, and even more prominent in the small x region, which dominates the gluon spectrum. The question of how we are going to quantify the error introduced by this systematic source of theoretical uncertainty arises immediately, and an answer will be given shortly. In the discussion above, we followed closely an in-depth analysis performed by Horowitz and Cole [128].

D.1.1 Systematic uncertainties
We approach the problem of quantifying the systematic uncertainties caused by the choice of the k ⊥ integration limits in a way which will be iterated several times throughout the construction of the CUJET model. The idea is to isolate those sources of uncertainty that have a clear impact on the observables we are going to compute from other sources whose effect is hindered by the simple rescaling of a free parameter such as strong coupling constant. In the context of the k ⊥ integration, we ask in Table 4 what is the sensitivity of the energy loss ∆E/E to the particular choice of integration limits, provided the freedom to adjust a free parameter identified as the coupling constant α s .

Curve
≡  Table 4. Fractional energy loss ∆E/E integrated from x E dN g /dx E for the curves shown in Fig. 16.
The results are indicated in the second column and range from 0.23 to 0.32. In the two rightmost columns are listed the values of the effective parameter α s needed to obtain the energy loss specified in . The free parameter α s needs to be tuned at most ±10%.
Given the interest in the ratio of light to heavy quark energy loss, we can immediately construct an error band which offers a quantitative measurement of the uncertainty generated by the choice of k ⊥ limits. In this way, α s is factored out and the results are independent of the rescaling of the free parameter. Fig. 17 shows the scaling of ∆E light /∆E heavy with the jet energy E and the plasma size L, for two distinct assumptions x ≡ x E and x ≡ x + .
The conclusion is evident: the choice of k ⊥ limits has a relevant impact at small energies E ≤ 15 GeV and long path lengths L ≥ 5 fm. Further theoretical steps to address large angle radiation and relax the collinear approximation need to be taken. Until then, in the development of the CUJET model we adhere to the collinear derivation of GLV and interpret x ≡ x + . When the gluon number distribution was needed as a differential in x E , namely x E dN g /dx E , we added the Jacobian of the transformation x + → x E to (C.6), and integrated k ⊥ up to k M AX ⊥ = xE. However, given the restricted size of such phase space, we will take our preferred assumption of k ⊥ limits and ignore the source of error coming from the x interpretation, especially when studying results in the high energy range of LHC.

D.2 Dead cone
The ability to determine the quark flavor dependence of any physical observable is not only a interesting characteristic of the DLGV integral, but also an invaluable tool used to compare predictions with data. Computing the energy loss for charm and bottom quarks within the same consistent framework, in fact allows us to put additional constraints on the model and therefore gain more insights about the nature of the quark gluon plasma. In this section we want to check the effects that the parton mass has on the transverse momentum distribution of gluon radiation.
The mass term M appears to have only a minor impact on the convergence of the series (Fig. 14), if nothing else by even improving it for certain combinations of E and ω. For very soft gluons (x 1), the heavy quark jet radiation spectrum does not differ much from its light quark counterpart, while for large values of x the radiation seems highly suppressed. The strong x dependence of the magnitude and shape of dN g /dxdk, as seen in Fig. 14, breaks the scaling with ω for typical gluon radiation from light quark jet.
Another effect is the filling of the "dead cone" characteristic of the vacuum spectrum. In vacuum, the transverse momentum distribution takes the form and the depletion of radiation takes place at angles We compare in the right panel of Fig. 18 the radiation spectrum of a heavy quark at different orders in opacity with the reference vacuum spectrum radiation, one notice immediately that the induced radiation fills in the dead cone already at first order in opacity. Since the dead cone region constitutes only a small fraction of the available phase space, the energy loss experienced by a heavy quark remains smaller than that of a light jet.  The left panel of Fig. 18 shows instead a striking feature: despite its non-vanishing mass equal to 1.2 GeV, the charm quark leads to a radiation spectrum very similar to the one of light quarks: not only the dead cone is absent and the vacuum spectrum almost divergent for k ⊥ → 0, but even the spectra have approximately the same shape and magnitude. This critical feature has vast phenomenological implications in the prediction of physical observables.

E Systematic uncertainties associated with n f , T f and dN/dy
In this section we analyze the sensitivity of CUJET to three of those parameters that govern the evolution of the medium: the number of quarkonic flavors n f , the fragmentation temperature T f , and the initial rapidity density dN/dy. 22 As usual, eventually we will hinder their effects to the rescaling fixed or running coupling constant, and consider α s or α max the only free parameter of CUJET, which will be constrained by a specific set of experimental data, typically pion R AA at a given value of transverse momentum and center of mass collision energy. Therefore, it is of great interest to show how R AA changes for different plasma assumptions, and observe if its functional form is modified once a proper rescaling of the coupling has been performed.
In Fig. 19 we change the value of n f from 0 (pure gluonic matter), to 2.5 (mix of gluonic and quarkonic degrees of freedom in chemical equilibrium). We can easily observe that a simple rescaling of α s of approximately 6% leads to a perfect agreement between the two scenarios. In [57], Zakharov reaches a similar conclusion starting from a path integral approach to the energy loss and using a running strong coupling. This simple analysis demonstrates the substantial insensitivity of CUJET to the detailed composition of the quark gluon plasma. We now focus on the late phase of plasma evolution and measure the sensitivity of R AA to the jet hadronization temperature T f . In Fig. 20, the partonic nuclear modification factor is shown for light and heavy quarks, for the default T f = 100 MeV and α s = 0.3 22 In principle, formation time and thermalization scheme also affect CUJET calculation. This issue is a topic of Appendix I. parameters, T f = 50 MeV and α s = 0.3, and finally T f = 200 MeV and α s rescaled to 0.35. We observe that jet quenching is "saturated" already at T f = 100 MeV: even if we let the jets interact until T drops to the (unphysical) value of 50 MeV, no significant changes occur in R AA . On the contrary, restricting the interaction region to T > 200 MeV alters significantly the results and a moderate ∼ 20% rescaling of the coupling constant is needed in order to reproduce the the original curve. But given the freedom to fit the coupling constant, CUJET can be regarded being insensitive to this source of theoretical uncertainty.
However, regarding T f , the study of experimental observables that are sensitive to the azimuthal anisotropy of the plasma and the angular distribution of jets would prove to be more insightful, because the fragmentation region generally resides in the outer shell (corona) of the plasma, and is more sensitive to the geometry of the collision, i.e. impact parameter, and the transverse expansion.
Finally, we study in CUJET1.0 the sensitivity of R AA to the initial rapidity density dN g /dy. This parameter is constrained by experimental observations, it fixes the initial density and temperature of the plasma according to Eq. (2.14). Intuitively, we expect the quenching to be higher for denser plasma, resulting in an increased suppression of R AA for collisions observed at the LHC. Our expectations are confirmed in Fig. 21.
Note in CUJET2.0 the medium information has been encoded in the 2+1D viscous hydrodynamic fields which presumably fit properly the hadron production spectra and bulk harmonics in the soft region, the systematic uncertainties associated with initial rapidity density are therefore irrelevant in the CUJET2.0 = rcDGLV + elastic + VISH2+1 framework.  Fig. 24) are used in this plot, as well as RHIC collision parameters. However, the initial rapidity density dN g /dy is increased from 1000 (opaque lines) to 2200 (solid lines). The increased dN g /dy is responsible for the suppression of R AA . Other parameters used in the simulation are: α s = 0.3, n f = 0, T f = 100 MeV, linear thermalization with initial time τ 0 = 1 fm/c.

F.1 Bjorken's formula
The first estimation for collisional energy loss in a quark gluon plasma was made by Bjorken [129], and his work still constitutes the benchmark against which any computation of this kind should be compared. Here we briefly outline his derivation.
In the limit E k, where k the momentum of the target particle in the medium, we can approximate the quark-quark, quark-gluon and gluon-gluon elastic cross sections as where c i,j is a numerical factor equal to 4/9, 1, 9/4 for {i, j} = {q, q}, {q, g} or {g, g} respectively. The energy loss per unit length can be written as Here E − E represents the energy lost in the collision, ρ i (k) is the quark or gluon number density, and Φ is the flux factor that accounts for the relative orientation of the target and projectile. Defining θ as the angle between the incoming parton and the target, Integrating (F.2) over dt, we obtain where B is defined by the integration limitst M AX andt M IN . If assuming B is independent of k for simplicity, we can sett M AX ≈ 2 k E ≈ 4T E andt M IN = µ 2 , with µ being the Debye screening mass of the plasma. If we further write the quark and gluon number densities as we can perform the last integration over all momenta d 3 k and finally get to the Bjorken energy loss formula In order to derive this short analytic result, several approximations were made in the way that infrared and ultraviolet divergences are regulated, i.e.t M IN andt M AX . Such divergences are physically related to the absence of collective medium effects (soft scattering) and recoil (hard scattering) in the derivation of the theory.

F.2 Numerical effects
In the left panels of Fig. 22 we observe a gain in ∆E/E after elastic collisions are taken into account. We notice three main features: (1) the energy loss is increased by up to 20%; (2) elastic losses almost do not distinguish between light and heavy quarks; (3) For sufficiently large L, ∆E/E shows signs of saturation, indicating complete quenching of the jets. The partial contribution of radiative and elastic losses to the total ∆E is given in the middle panels of Fig. 22, assuming a dynamical medium. Here we immediately appreciate the difference between light and bottom quarks: while the relative elastic contribution diminishes with L and is approximately constant with E in the case of light partons, the exact opposite behavior is observed for heavy quarks. This has a remarkable impact on the phenomenology: the ratio ∆E light /∆E heavy , shown in the right panels of Fig. 22, drops by almost 25% in the large L and small E regions.
The inclusion of dynamical effects first, and elastic collisions later, has brought the light to heavy quark energy loss ratio down from a factor of more than 2x to about 1.5x, in the range of energies E ∼ 10−30 GeV and path lengths L ∼ 4−6 fm. These improvements constitute a promising step toward closing the gap between theoretical models and experimental data, which is shown at RHIC as a surprising similarity between the quenching of light and heavy quark jets.

G Radiative energy loss probability distribution and fluctuations
We dedicate this section to discuss the flavor dependent energy loss probability distribution and fluctuation effect. The energy loss is computed by integrating d P ( ) over the range  Fig. 23 we show the radiative energy loss probability distribution P ( ) for different quark flavors that propagate in a plasma of size L. The integrated radiative energy loss ∆E/E is shown in the middle panels of Fig. 23, alongside a comparison with the same quantity obtained without the inclusion of fluctuation effects, i.e. obtained by simply integrating the gluon spectrum dx xdN/dx. We see typical energy loss probability distribution peaks at small for all flavors, and light and charm have almost negligible difference in terms of P ( ), both of them lose more energy than the bottom quark. The inclusion of fluctuation effects appears to alter only minimally the result.
The total energy loss of the jet as a function of E and L is presented in right panels of Fig. 23. The same features observed in the context of DGLV are also present in the dynamical scenario, from the coherence physics that determines the quadratic or linear L numerical NLO and FONLL initial cross sections [127], and we will estimate the error band associated with this source of systematic uncertainty. In Fig. 24, we illustrate the initial partonic production spectra for gluon, light, charm and bottom quark at RHIC and LHC energies. TeV (LHC, right). Notice how steeper the RHIC spectra are compared to LHC ones. The light spectra are computed from the LO pQCD CTEQ5 code provided in [112]. Numerical computations of the NLO and FONLL initial cross sections for the heavy sector are provided in [127].
To compare spectrum variations from RHIC to LHC, in CUJET1.0, we use separate initial rapidity density dN g /dy and production spectra for RHIC and LHC. The theoretical curves are shown in 25, superimposed on Fig. 21: The impact on R AA is large, and two separate effects can be noticed: (1) softer LHC spectra cause a vertical lift in R AA that completely counters the suppression generated by the increased density; (2) pion R AA rises faster with p T , due to the particular shape of the light quark spectra at LHC.
The uncertainties that arise from the choice of NLO or FONLL schemes for heavy quark initial spectra are shown in Fig. 26. The error bands shown in the Figure are relatively small, and we estimate that to be 5% in R AA at most. At the partonic level, in fact, any uncertainty in the normalization of the production spectra is factored out: R AA is only sensitive to changes in the slope.
Depending on what physical observables we are interested to compute, different features of the partonic spectra may or may not assume a relevant role. Since R AA is defined as a ratio of particles yields, the absolute value of the cross section matters little and the normalization drops out in the definition of the observable. What influences the computation is rather the slope of the cross section, as well as the relative normalization between different flavors.
An insightful example comes from the pion yield in p+p events at LHC, Fig. 27, which is computed by convoluting the production spectra of quarks and gluons with the appropriate fragmentation functions (more details in Section B.4.3). We can make two observations: (1) Since gluons and light quark contributions are summed together to get the pion yield,   Figure 27. left: p+p gluon and light quark production spectra, same as Fig. 24. right: p+p pion spectra from gluon only contribution (green), quark only contribution (blue), and total gluon plus quark contribution (black), assuming no "cold" nuclear effects. The pion spectra are computed using KKP fragmentation functions.
other hand, drops out once the nuclear modification factor ratio is taken. (2) Despite the high production of gluons at low p T , the gluon distribution is much steeper than the quark one. As a consequence, once fragmentation is taken into account, the gluonic contribution to the total number of pions produced sinks below the quarkonic one already at p T 25 GeV. It is then reasonable to expect R AA to depend on the light quark sector only for sufficiently high transverse momentum.
Another example where only the relative steepness between production spectra matters is given by the comparison between the (unphysical) partonic yields of light and charm quarks. In previous appendices we saw that light and charm quarks approximately lose the same amount of energy when they propagate through a deconfined medium, however the production spectrum of charm quarks is much steeper than the one of light quarks (cf. Fig. 24). The immediate consequence is that the partonic yield for charm quark is suppressed in AA collisions compared to the other, regardless of the separate normalization of the production spectra.
A similar effect applies when we compare RHIC (steeper) and LHC (flatter) spectra: the expected energy loss increase at LHC with respect to RHIC due to higher densities and temperatures, which itself would drive the particle yields down, is going to be partly compensated by the flatter production cross-sections, which in turn drive the yields up.
Finally, a comment on heavy jet quenching: the measurement of heavy flavors has often been limited to the experimental analysis of non-photonic electrons, produced mainly in the secondary decays c → D → e and b → B → e (where D and B refer to the D and B meson respectively). In such case, the relative norm between charm and bottom spectra plays a critical role, in the same way that both gluon and light quarks contribute to pion R AA . Unfortunately, the uncertainties are more significant in the heavy flavor scenario than in the pion scenario, therefore a direct measurement of the intermediate mesons would undoubtedly provide a much cleaner and insightful measurement to be compared to  Figure 28. Temperature profile of the QGP in a central (b = 0) collision at RHIC energies. The density is constrained by the observed dN/dy = 1000. The black curve represents the temperature at constant τ 0 = 1 fm/c for a radial section of plasma. The red curves represents the 1/τ 1/3 temperature probed by a quark that is created at r = 0 and propagates outward along z ≡ r (with the solid, dotted and dashed curves representing the linear, divergent and free streaming cases respectively). The dashed black T ≈ 100 MeV line corresponds to the fragmentation temperature of the jet.

I.2 Systematic uncertainties
We show in Fig. 29 how differently light and heavy quarks lose energy due to elastic and inelastic collisions are during different early stages of the plasma longitudinal expansion. For jets produced in central Au+Au events, the differential d < ∆E/E > /dz indicates the fractional energy loss during the first fm's of the jet evolution. Heavy quarks lose a larger percentage of their energy via radiative processes and its radiative energy loss rate follows the medium thermalization. The mass-dependent jet behavior observed in Fig. 29 could be used as a phenomenological indicator of the thermalization mechanism. For different parametrizations of f (τ /τ 0 ), one could expect a different relative yield between light and heavy quark jets. We can in fact expect that once the free parameters of the model (α s or α max ) are fixed by a comparison of the light sector with data, each assumptions of f (τ /τ 0 ) will yield a different result for the heavy sector. This fact is portrayed in Fig. 30, where the ratio ∆E light /∆E bottom is given as a function of L for all possible temporal envelopes.
Next, we study the hadron suppression factor's sensitivity to the thermalization phase of the plasma. The results from CUJET1.0 calculations for RHIC Au+Au 200AGeV and LHC Pb+Pb 2.76ATeV central collisions are shown in Fig. 31.
We see a great sensitivity of pion R AA to the pre-thermalization phase of the evolution in Fig. 31, which however can be counter-balanced by an adequate rescaling of the coupling constant α s , i.e. varying α s in divergent or free streaming scheme down or up by 10% recovers the linear scenario. If we constrain α s to fit a specific p T point of pion R AA at  Figure 29. Differential d < ∆E/E > /dz for light (left) and heavy (right) quarks, in a QGP defined by dN/dy = 1000, τ 0 = 1 fm/c and n f = 0. The initial energy of the quarks is 20 GeV. Blue and orange colors refer to radiative losses, whereas purple and brown to elastic ones. Notice how quickly d < ∆E/E > /dz drops for heavy quarks compared to light jets. LPM interference effects are responsible for the finite value of the energy loss at very short z in the divergent plasma scenario. Results are calculated within the framework of fixed coupling CUJET1.0 with α s =0.3.  Figure 30. Energy loss ratio ∆E light /∆E heavy as a function of L between light and bottom quarks, for the three linear (solid), divergent (dotted) and free streaming (dashed) initial conditions. The energy loss is obtained by integrating the curves in Fig. 29 up to z = L. For sufficiently long path lengths, the relative difference between the three approximations reaches approximately 10%.
RHIC initial conditions, in Fig. 31 left, we observe a complete overlap -or "degeneracy"among the linear, divergent and free streaming scenarios. The constrained fit extrapolated to LHC energies in Fig. 31 right, shows on the other hand a moderate "splitting" at high p T among the same curves. Although the difference is too small to be measured exper-  imentally, in theory this effect can be studied to discriminate among pre-thermalization phenomenological models.
The same effect is visible in Fig. 32, where pions, D and B meson R AA is plotted assuming RHIC (left) and LHC (right) initial conditions. The curves are constrained by the same RHIC fit of Fig. 31 left. We observe a moderate "splitting" of B meson R AA across all p T , which is a signature of the differences between the light and heavy quark quenching mechanism during the early evolution of the plasma (cf. Fig. 29). The splitting in the nuclear modification factor is less than 10% in Fig. 32, and it is difficult experimentally resolve this splitting in the near future. Nevertheless, what we have observed here is a clear indication of the importance of making simultaneous constrained fits to as many "orthogonal" observables as possible, and it implies that the flavor dependent quenching pattern and single particle azimuthal anisotropy are key observables of interest for refining the phase space of CUJET pQCD tomographic model.