Cone size dependence of jet suppression in heavy-ion collisions

The strong suppression of high-$p_T$ jets in heavy ion collisions is a result of elastic and inelastic energy loss suffered by the jet multi-prong collection of color charges that are resolved by medium interactions. Hence, quenching effects depend on the fluctuations of the jet substructure that are probed by the cone size dependence of the spectrum. In this letter, we present the first complete, analytic calculation of the inclusive $R$-dependent jet spectrum in PbPb collisions at LHC energies, including resummation of energy loss effects from hard, vacuum-like emissions occurring in the medium and modeling of soft energy flow and recovery at the jet cone. Both the geometry of the collision and the local medium properties, such as the temperature and fluid velocity, are given by a hydrodynamic evolution of the medium, leaving only the coupling constant in the medium as a free parameter. The calculation yields a good description of the centrality and $p_T$ dependence of jet suppression for $R=0.4$ together with a mild cone size dependence, which is in agreement with recent experimental results. Gauging the theoretical uncertainties, we find that the largest sensitivity resides in the leading logarithmic approximation of the phase space of resolved splittings, which can be improved systematically, while non-perturbative modeling of the soft-gluon sector is of relatively minor importance up to large cone sizes.

Introduction. Jets are collimated sprays of energetic particles produced in collider experiments that act as proxies of accelerated quark and gluon degrees of freedom originating from elementary large momentum-transfer processes or decays of massive bosons. In this context, precision computations of QCD jet events play a crucial role in a wide range of fundamental measurements at colliders [1,2], including measurements of the Higgs boson properties [3] and searches beyond the Standard Model.
In contrast, jet physics in heavy ion collisions probes the discovery frontier to potentially reveal and detail new emergent QCD phenomena in dense partonic systems. The creation of a short-lived, hot and dense state of deconfined matter, also known as the quark-gluon plasma (QGP), leaves a strong imprint on high-p T probes [4,5]. This phenomenon, commonly referred to as "jet quenching", was observed for the first time at RHIC [6][7][8] and later at the LHC [9][10][11][12][13][14]. Currently, the exact mechanisms responsible for jet modifications, including details of the energy transport from high-energy to low-energy modes and color/quantum decoherence of multi-partonic states, are under intense investigation.
The basic mechanisms of parton energy loss were understood and formalized in the 90's and implemented for RHIC phenomenology [15][16][17][18][19], where a number of approximations, in particular for the medium induced ra-diative spectrum, were then necessary to allow for analytic computations. This introduced a theoretical bias on model calculations, absent from full numerical approaches [20][21][22], that could be alleviated by incorporating the two main scattering regimes: the Rutherford scattering regime, dominated by a single hard momentum transfer, and the low momentum regime where multiple scatterings contribute with order one probability. Furthermore, with the measurements of fully reconstructed jets at the LHC and RHIC, it was soon recognized that higher order corrections accounted for by parton cascades are not negligible. The need to address these effects spurred the rapid development of Monte Carlo (MC) event generators [23][24][25][26][27][28][29][30][31][32][33] which, to some extent, rely on modeling of the quantum nature of jet evolution. In parallel with this computational effort, tremendous conceptual progress has been made in addressing these questions by analyzing the interference structure of two successive splittings within the medium [34][35][36][37][38]. To leading logarithmic accuracy, it has been shown that the in-medium jet evolution is characterized by an early vacuum parton cascade whose constituents either get resolved by the medium due to color decoherence, whereas unresolved splittings factorize from the in-medium evolution, losing energy coherently as a single color charge [39,40]. This work aims to address these two challenges within a firstprinciple analytic framework.
In addition to the modification of the hard components of the jet and their interactions with the plasma constituents, there are non-universal contributions to jet observables pertaining to how soft jet constituents thermalize in the plasma. In analogy with hadronization ef-fects, these non-perturbative contributions are bound to be modeled. This leads us to one of the most important questions in jet quenching physics: what is the relative magnitude of the uncertainties associated with describing the hard, perturbative structures-that are systematically improvable-and the soft, infrared features of medium-modified jets as a function of their kinematics? Providing a quantitative answer to this question is crucial if one aims to establish the predictive power of weak coupling techniques in jet quenching phenomenology and probe the transport properties of the QGP. With this work, we aim to provide an answer to this fundamental question.
In this letter, we report a first-principle calculation of the single-inclusive jet spectrum and its cone size dependence in heavy ion collisions where high density effects are resummed to all orders. Even though jets with a larger cone do retain a larger fraction of the lost energy, a priori reducing jet suppression compared to a smaller one, we show that resumming the additional energy loss induced by the cone-size dependent jet substructure fluctuations yields a final jet suppression that is very mildly dependent on R. The well established connection between energy loss dynamics and coherence effects, which determine the actual resolved phase space of the jet in the medium, allows us to confront our results with highstatistics experimental data merely by constraining the strength of the QCD coupling in the medium. Additional fluctuations on the path and medium density explored by the jet, which vary event by event, are taken care of by embedding our framework into a realistic heavy-ion environment in which the medium is described by the explosion of a liquid droplet of deconfined QCD matter.
Theoretical formalism. The spectrum of jets with cone size R in proton-proton collisions is given by the convolution of the initial hard parton spectra with the corresponding fragmentation function. The latter describes the energy remaining within the jet at different angular resolutions R, starting from a large value R 0 ∼ 1. For a steeply falling initial spectrum, it can be written as [41] where n ≡ n k (p T , R 0 ) is the power-index of the crosssection of the hard parton with flavor k. This is calculated at leading order (LO) at the factorization scale Q 2 fac , such thatσ k = f i/A ⊗ f j/A ⊗σ ij→k(l) , and involves a convolution of parton distribution functions (PDFs) f i/A (x, Q 2 fac ) with the 2 → 2 QCD scattering cross section σ ij→kl . The moment of the fragmentation function of an initial hard parton with flavor k, i.e. f i/k , due to flavor conversion during the DGLAP evolution [41][42][43][44].
Correspondingly, the cross section in nucleus-nucleus collisions (AA) are convolved with a probability distribution P ( ) describing medium-induced energy loss out of the jet cone, and reads whereσ pp i corresponds to the quark/gluon contribution to the total cross section in Eq. (1) (the tilde serves as a reminder that the proton PDFs are replaced by nuclear PDFs). Finally, the flavor dependent resummed quenching factors (QF) account for the energy loss by a jet with momentum p T and size R during the passage of a background medium [45]. In the limit of large power index n, we we invoke the asymptotic expansionσ pp ∝ (p T + ) −n ∼ p −n T e −n /p T 1 + O(n 2 /p 2 T ) , which to leading order allows us to identify the QF with the Laplace transform (LT) of the energy loss probability, i.e., Q(p T ) = ∞ 0 d P ( )e −ν | ν=n/p T , where we have omited the flavor subscript for clarity.
The factorization (2) reduces trivially to the jet production cross section in the absence of final-state interactions, Eq. (1), by setting the quenching factors to unity, Q i → 1, and replacing the nuclear PDFs by standard proton PDFs. It is justified by the fact that out of cone vacuum evolution takes place at much shorter times than energy loss and was used as a basis for the extraction of the quenching weights from the data [46,47].
A novel ingredient of our setup are the quenching factors Q i (p T , R) that resum contributions to the total energy loss of a jet consisting of many color charges that interact with the medium. Every splitting that occurs at short time scales within the medium, gives rise to an additional color current that can scatter with the plasma constituents and source further medium-induced energy loss. The magnitude of this effect can be gauged by comparing the formation time of a splitting, t f = 2/[z(1 − z)p T θ 2 ], to the characteristic time scale the medium needs to resolve the product of the splitting, namely t d = [q 0 θ 2 /12] −1/3 [34][35][36][37]. Here,q 0 ≡ d k 2 ⊥ /dt is the transport coefficient that encodes medium properties, the so-called jet quenching parameter [16]. Hence, jet splittings occurring at time scales much shorter than the related medium time scale, that is if t f t d L, are unaffected by the medium and obey the same properties as vacuum splittings [39,40]. The latter inequality implies that a splitting with θ < θ c , where the critical angle is θ c = (q 0 L 3 /12) −1/2 , will not be resolved by the medium. With these considerations in mind, one can show that in the large n limit, owing to the fact that the convolution of energy loss probability distributions reduces to a direct product of quenching factors (in Laplace space), the evolution equation for the resummed quenching weight is [39] where ji (z) are the un-regularized Altarelli-Parisi splitting functions and the phase space constraint is given by Θ res (p, R) = Θ(t f < t d < L). Above, it is understood that p is evaluated at p ≡ p T . This distinction is necessary when solving Eq. (3) since the initial condition also depends on p T . The non-linear evolution equations account for the energy loss of the multi-prong jet substructures that are generated by early collinear splittings.
The initial conditions for the resummed quenching factors Q i (p, R) at R = 0 are the bare quenching factors for single partons. In this work, we have where the two bare quenching factors are the LT of the corresponding probability distributions that describe radiative and elastic energy loss [45,48], For their precise definitions, see Eqs. (8) and (9) below. The radiative and elastic energy loss are driven by the transport coefficientsq andê [49], respectively, which are related by the fluctuation-dissipation relationê =q/(4T ) in a weakly-coupled plasma [50] (whereê g =ê for gluons, andê q = C F Ncê g for quarks). The quenching factor due to radiative energy loss off a single parton is simply the exponential of the LT for a single inclusive gluon radiative spectrum [45,48]. For our purposes, we should rather consider how single partons contribute to the energy loss of the jet by accounting for the energy that is transported outside of the jet reconstruction cone. To this aim, we exploit the wide parametric angular separation between the regime of soft emissions that undergo a rapid turbulent cascade responsible for transporting energy from the jet scale to the medium temperature where dissipation forces take over, and the regime of collimated semi-hard emissions, which experience broadening through collisions with the medium constituents [51][52][53][54][55].
The medium-induced gluon radiation spectrum has been computed up to next-to-leading order (NLO) within the improved opacity expansion (IOE) in the soft limit [56][57][58] and unifies both the BDMPS approach with the GLV/higher-twist formalism [59,60], which has proven to be an important ingredient for phenomenological studies [61]. The IOE was also shown to be very accurate when compared to exact numerical solutions [62]. It is expressed as where Ω = (1 − i) q/(4ω), k 2 (s) = i ωΩ 2 [cot Ωs − tan Ω(L−s)], and the strong coupling constant runs with the typical transverse momentum of the emission, i.e. α s = α s (qω) 1/4 [63]. In this expansion, the effective transport coefficientq differs from the bareq 0 by a factor that reflects the full leading logarithmic contribution, i.e.q ATeV, compared to ATLAS data [65], for different centralities.
whereq 0 = g 2 med N c m 2 D T /(4π) for a thermal medium in the Hard Thermal Loop (HTL) theory and the lower cutoff scale is µ 2 [57,58]. The Debye mass m D computed at LO in a thermal medium reads m 2 D = 3g 2 med T 2 /2 (for three active quark flavors). The effective medium scale Q 2 depends itself on the amount of rescattering in the medium and can be found by solving the transcendental equation Q 4 =q 0 ω ln Q 2 /µ 2 * [64]. In our framework, the medium coupling g med , the only free parameter that determines energy loss, is to be extracted from the comparison to experimental data.
We first consider semi-hard gluons that are emitted within the range ω s < ω ω c , where ω c ≡ q 0 ln(q 0 L/µ 2 * )L 2 /2 corresponds to the maximum accumulated energy through multiple soft scatterings, and ω s ≡ (g 2 med N c /(2π) 2 ) 2 πq 0 L 2 is the energy scale at which emission probability is of order one, determining the onset of turbulent energy loss [66]. Their broadening distribution reflects the typical transverse momentum kicks received in the plasma. The fact that the two terms entering the full NLO spectrum are dominated by different kinds of processes has to be reflected in the typical behavior of the respective broadening distribution. In this way, the softer gluons from Eq. (4), with ω ω c and small transverse kicks k 2 ⊥ ∼qL, will experience Gaussian broadening, while the harder emissions from Eq. (5), with ω ω c and typically large transverse momenta k 2 ⊥ > ω/L qL, where the first inequality arises from demanding that t f = ω/k 2 ⊥ < L, are governed by a power-law behavior, ∼q 0 L/k 4 ⊥ . We assume that the effect of broadening appears as a multiplicative factor B ωR; , representing the probability for the emitted gluon to be transported out to an angle larger than the jet cone, θ > R, where Q 2 broad denotes a characteristic broadening scale. This is concretely realized by integrating the broadening probability distribution, derived in [67], for angles larger than the jet cone. This distribution and the proposed factorized form correctly interpolate between the multiple-scattering and higher-twist regimes. The full out-of-cone spectrum consists of two terms, from the IOE expansion up to NLO with their corresponding broadening factors [52], such that Since emissions can take place anywhere along the inmedium path, one also has to average over the radiation time. This is approximated by simply setting Q 2 broad =qL/2 [51,52] in the first term. The choice of scale in the second term reproduces the correct behavior of the full GLV spectrum in the high ω and k ⊥ regime up to a logarithmic factor that we neglect. In brief, dI/(dωdk 2 ⊥ ) 2ᾱq0L (1) dω , where the ellipses represent the logarithmic contributions and dI (1) /dω ∼ πᾱq 0 L 2 /ω is the limiting behavior of Eq. (5) for ω ω c and k ⊥ qL. The B distribution was used in the second term to insure proper normalization. For more details see the supplemental material.
Soft gluons, with T < ω < ω s , cascade quasiinstantaneously to the thermal scale [66] and should effectively be treated within hydrodynamics. Their emission rate is therefore not affected by transverse momentum broadening. Assuming that their distribution becomes approximately uniform in the solid angle around the jet, we account for the possibility that a fraction of this energy ends up back in the jet cone by modifying ω → ω(1 − (R/R rec ) 2 ), where the recovery angle R rec is a free parameter [68]. An analogous modification is applied to the elastic quenching factor. Emissions at ω < T belong to the Bethe-Heitler regime, and are not relevant for our present phenomenological application [62,69].
Putting all the pieces together, the final expression for the radiative bare quenching factor reads where the parton flavor index is implicit and ν ≡ n/p T . We have approximated dI NLO /dω dI (0) /dω in the soft regime. The bare quenching factor for elastic energy loss is also with implicit parton flavor dependence. It results from taking the LT of δ( −êL(1−R 2 /R 2 rec )). We shall see that our results at small cone sizes are not very sensitive to the above modeling of energy recovery.
Numerical results. Using the bare quenching factors for the radiative (8) and elastic (9) contributions to energy loss, we numerically solve the coupled evolution equations in Eq. (3). The cone-size dependence of the bare quenching factors, through the broadening effects encoded in dI > /dω and resulting in more energy loss for smaller R, are to a large extent washed away by the evolution. This is because wider jets have a larger resolved phase space and hence comprise more radiating charges than the narrower jets, effectively hampering energy recovery.
In our numerical computations, we fix the values of the two free parameters of our setup, g med and R rec . The energy recovery parameter R rec has been varied between R rec = π/2 and R rec = (5/6) π/2, which was estimated from a linearized approach to model the QGP wake for [70,71]. To constrain g med we have compared our results for the widely used nuclear suppression factor for jet production, also known as R AA , for jets with R = 0.4 around p T ∼ 100 GeV against high-statistics experimental data from ATLAS for the 0-10% centrality class of PbPb collisions at √ s = 5.02 ATeV [65]. In order to compute R AA , we have taken the ratio between the nuclear and the vacuum spectra, both defined through Eqs. (1) and (2), which comprise a weighted sum of the quark and gluon jet contributions to the full spectrum [72]. Eventby-event in-medium path fluctuations of a jet through the QGP have been taken into account by embedding our framework into a realistic heavy ion environment as simulated in the VISHNU hydrodynamical model [73], see the supplemental material for further details. The value of g med is thus constrained by the experimental data to be within the range g med ∈ {2.2, 2.3}. We emphasize that the magnitude of quenching is predominantly driven by the emission of copious soft gluons [45], see also Tab. I below. The extracted parameters yield an average value q 0 0.41 GeV 2 /fm in 0-10% central PbPb collisions that is well within the perturbative regime, see the supplemental material for more information on the centrality dependence of key parameters. However, the logarithmic corrections to the bare medium parameters, resulting in Q 2 = 14.2 GeV 2 for the factorization scale andq = 2.46 GeV 2 /fm, produce a relatively large maximal medium energy scale ω c ≈ 65 GeV. These effects only become apparent when carefully treating the dominant scattering regimes of the full spectrum, as achieved with Eqs. (4) and (5).  We show results for R AA as a function of jet p T and centrality in Fig. 1, confronted against data from AT-LAS [65] for R = 0.4 jets. We can see that our theoretical results, constrained to describe R AA around p T ∼ 100 GeV for 0-10% only, give an excellent description of both the p T dependence and centrality evolution of jet suppression. We note that the downturn of R AA at the highest jet p T is due to the nuclear PDF modifications imprinted in the initial hard parton spectra of PbPb collisions. Finally, we quantify the R (in)dependence of jet suppression by taking double ratios of the full results for R AA as in Fig. 2, with the largest size R = 1 in the denominator. Such notably mild dependence of jet suppression with R is in agreement with ALICE results at low-p T [74] and with recent experimental preliminary data from CMS at high-p T [75].
Summary and discussion. We have provided a first, analytical description of the cone-size dependent jet spectrum in heavy ion collisions at the LHC implemented in a realistic event-by-event setup including nuclear geometry and hydrodynamic expansion of the quark-gluon plasma and accounting for multiple scattering effects. By adequately introducing the notion of single-parton energy loss within the context of a multi-parton object such as a jet, in which the phase space is determined not only by the jet p T and size R but also on whether splittings are resolved by the medium as determined by color coher-ence effects, we have shown that the resummation of the quenching factors yield results for jet suppression that are very mildly dependent on R. Our procedure to embed our formalism into a realistic heavy ion environment has proven successful, given in particular our good description of the centrality evolution of high-statistics experimental data. The various tools and formalism introduced in this work can be systematically improved and will be applied to other jet substructure observables in heavy ion collisions in the future.
While the error bands provided in the plots above stem mainly from the constraining power of the experimental data at R = 0.4 on the medium coupling g med , we would currently like to discuss the sensitivity of our results on the various assumptions made in the setup in order to identify the main sources of uncertainty. Our main findings are summarized in Tab. I for moderate cone sizes 0.2 ≤ R ≤ 0.6, see also the supplemental material for a full scan of parameters and their centrality dependence. First, the inclusion of the higher-twist radiative spectrum in the IOE is of mild importance for this observable, since such emissions typically occur at small angles, but it improves the description at high-p T . Furthermore, as expected, notable bias effects can be identified through the strong sensitivity to the power of the steeply falling spectrum n, which point to the importance of higher order terms in the large n expansion that can be calculated systematically. More importantly, comparing the effect of changing the hard phase space (through θ c ) and the parameters governing the behavior and recovery of soft gluons (through ω s and R rec ), we note that an increased precision in the perturbative sector is still needed before the sensitivity to non-perturbative effects start to dominate. For instance, going beyond leading logarithmic accuracy to compute θ c will be important to rigorously study the interesting marked centrality dependence of this critical angle. The importance of the recovery parameter R rec has been gauged between two limiting scenarios of R rec = 1 corresponding to almost complete energy recovery for large-R jets and R rec = ∞ corresponding to no energy recovery. Surprisingly, we find relatively little sensitivity to this parameter at these moderate cone sizes. Conversely, the sensitivity becomes the dominant source of uncertainty only at large-R, i.e. R ≈ 1, jets.
In conclusion, we have demonstrated how the conesize dependent jet spectrum and R AA is largely governed by the energy loss off hard, resolved jet splittings in the medium through copious, soft gluon radiation and their subsequent broadening out of the jet cone and elastic drag. The results obtained by analyzing the various components of our setup emphasize the importance of analytical tools to guide more sophisticated numerical models, such as MC parton showers.

Implementing broadening in the medium-induced spectrum
The broadening distribution is related to the dipole scattering amplitude via a Fourier transform. In a static medium it reads within the leading logarithmic approximation. Here, µ * represents a model-dependent, infrared scale, see below Eq. (6) for its precise definition in the HTL model of medium interactions. Its asymptotic behavior are as follows In the latter case, we have neglected additional logarithmic terms. In the regime of multiple scattering, parton splitting take place as a two-stage process where, first, the splitting products separate and decohere and, second, they broaden independently [37,82]. Therefore, we can easily propose a similar ansatz for the integrated, out-of-cone spectrum where B ωR; Q 2 broad = Q 2 broad 4π ∞ y dx P(x) with y = (ωR) 2 /Q 2 broad and Q broad is the characteristic broadening scale.
The gluon emission regime where ω < ω c is dominated by multiple scatterings. For ω > ω c , the dominant contribution to the spectrum is provided by a single, hard scattering with the medium. In this case, the unintegrated GLV (N = 1) spectrum reads where k = (k x , k y ) is the transverse momentum vector and k ⊥ ≡ |k|. At large energy and transverse momentum, i.e. k 2 ⊥ ω/L and ω µ 2 L, the spectrum behaves as Given that the integrated spectrum in the limit ω ω c becomes ω dI dω =ᾱ π 4q 0L 2 2ω , we can thus express the above unintegrated spectrum as a function of the integrated one where the ellipses correspond to the logarithmic corrections in Eq. (14). For the purpose of our applications, in the following we will neglect such corrections. The two discussed limits can be accounted for by expanding the full potential around the harmonic oscillator [67], where the x-dependent logarithm is neglected in Eq. (10). The formula for broadening up to next-toleading order in the improved opacity expansion reads where we have defined x ≡ k 2 /Q 2 s , and the expansion parameter is λ = (ln Q 2 found implicitly. For our purposes, λ is restricted to be within 0 and 1. After integrating out the angle, we get B(y) = −λ + e −y 1 + λ 1 + y Ei(y) − log 4y (17) with y ≡ (ωR) 2 /Q 2 s . In addition, we should also take an average of the production point. A good estimate, at least for the leading order part, is to consider Q 2 s → Q 2 s /2. The behavior of the N = 1 spectrum at high ω is well described by considering instead max Q 2 s , 16ω/(π 2 L) .

Bare and resummed quenching factors
We plot the bare and resummed quenching factors for quarks and gluons in Fig. 3 as a function of p T for values of the medium parameters corresponding to 0-5% centrality PbPb collisions at √ s = 5.02 ATeV. As discussed above, these include both the effect of radiative and elastic energy loss, incorporated in the bare quenching factors as where i = q, g. The p T dependence appears mainly through the ratio n i (p T )/p T , where n i (p T ) is the powerlaw indices of the spectra of quark and gluon initiated jets. These weight factors also constitute the initial conditions for the non-linear evolution in Eq. (3) for the fully resummed weights Q i (p T , R). The error bands correspond to varying the medium coupling g med ∈ [2.2, 2.3] and recovery angle R rec ∈ [5π/12, π/2], as described in the main text. In Fig. 3 we see that the clear ordering with R of the bare quenching factors (dashed), with more energy loss for smaller R, is no longer present in the resummed case (solid). The energy loss is even enhanced when comparing R = 0.2 to R = 0.4. Both for quark initiated jets (top) and gluon initiated jets (bottom), resummed quenching factors show only slight differences among the different R.

Embedding into a heavy ion environment
In order to account for the fluctuations in jet energy loss due to the different paths that the jet can explore through the expanding QGP, we need to embed our theoretical framework into a realistic heavy ion background. The necessary steps are the following: • Sample the production point in the transverse plane (x, y) using the overlap of the thickness functions of the two nuclei, T AB (x, y; b) = T A (x − b/2, y)T B (x + b/2, y), where b is the impact parameter of the nuclear collision and the thickness function is the transverse density of nucleons of the Lorentz contracted nuclei, distributed according to the Woods-Saxon density function (see [83] for more details on Glauber modeling).
• Assign a random orientation in the transverse plane and a random value of rapidity within the range −2 ≤ y ≤ 2.
• Following the path of the jet, compute the integrated values of the necessary physical variables, which in general depend on the local temperature T and fluid velocity u, until the jet exits the QGP phase at T c = 145 MeV (possible energy loss effects during the hadron gas phase, which could be more relevant for the softer particles [84][85][86][87], have been ignored for the moment and will be addressed in the future). The values of T and u are read from event averaged hydrodynamic profiles for the evolution of an expanding droplet of liquid QGP [73] in PbPb collisions at √ s = 5.02 ATeV for different centrality classes.
Given that quantities like the fluid temperature T are given in the local fluid rest frame, we need to consider the actual distance traversed by the jet at each time step within such reference frame [88]: where v ≡ p/E is the jet axis and u and γ F are the local fluid velocity and Lorentz factor, respectively. Up to numerical factors, the set of physical variables that we need are found by integrating along the path Γ(t) of a jet such as, for example: and analogously for the rest of variables: T , m D , θ c ,ê and w c . Note that, due to the presence of a flowing medium, transport coefficients get a dilution factor (p·u)/p 0 , where p is the four-momentum of the jet and u the fluid fourvelocity [89].
In this way, we obtain a set of representative histories for the jet in-medium path so we can use the event-byevent values of the path integrated physical variables to compute jet energy loss jet-by-jet. Our final results are obtained by averaging over all configurations.
To provide the reader some guidance on the typical magnitude of these quantities, we show in Table II the average values for several of the relevant physical variables that enter our calculation, as a function of centrality, where we set g med = 2.25. We point out the interesting correlation between the presence of decoherence effects, for θ c < R, and complete coherence, θ c ∼ R, as a function of centrality.

Sensitivity of results to systematic uncertainties
The calculation presented in this work is based on a perturbative description of hard, vacuum-like and semihard, medium-induced emissions. In addition, the regime of soft gluon emissions where thermalization plays an important role are modelled using parametric estimates. In order to systematically quantify the theoretical uncertainties imbued in our modelling, we have performed a scan varying most of the parameters that go into the computation. This is presented as a relative change in R AA at p T 110 GeV for several jet radius R for the 0-10% centrality class in PbPb collisions at √ s = 5.02 ATeV, when varying the parameters around the central values g med = 2.25 and R rec = π/2, as shown in Table III. Clearly, R AA is sensitive to the index n of the hard, steeply falling spectrum. This is a reflection of the bias effect related to jet selection. Next, a large sensitivity is reported in the variation of θ c in the phase space of the resummed quenching factors. This points to the importance of pushing beyond the leading logarithmic precision presently implemented. We also report little sensitivity to the precise value of the transition between the semi-hard and soft gluon regimes, encoded in the parameter ω s . Finally, we observe that even though we vary R rec between two extreme values such as R rec = 1 (full recovery of energy at R = 1) and R rec = ∞ (no recovery at all, for any R), these non-perturbative effects remain moderate ( 10%) up to relatively large cones of around R = 0.6.