Determining the jet transport coefficient \hat{q} from inclusive hadron suppression measurements using Bayesian parameter estimation
S. Cao, Y. Chen, J. Coleman, J. Mulligan, P. M. Jacobs, R. A. Soltz, A. Angerami, R. Arora, S. A. Bass, L. Cunqueiro, T. Dai, L. Du, R. Ehlers, H. Elfner, D. Everett, W. Fan, R. J. Fries, C. Gale, F. Garza, Y. He, M. Heffernan, U. Heinz, B. V. Jacak, S. Jeon, W. Ke, B. Kim, M. Kordell, A. Kumar, A. Majumder, S. Mak, M. McNelis, C. Nattrass, D. Oliinychenko, C. Park, J.-F. Paquet, J. H. Putschke, G. Roland, A. Silva, B. Schenke, L. Schwiebert, C. Shen, C. Sirimanna, Y. Tachibana, G. Vujanovic, X.-N. Wang, R. L. Wolpert, Y. Xu
DDetermining the jet transport coefficient ˆ q from inclusive hadron suppressionmeasurements using Bayesian parameter estimation S. Cao,
1, 2
Y. Chen, J. Coleman, J. Mulligan,
5, 6
P. M. Jacobs,
5, 6
R. A. Soltz,
1, 7
A. Angerami, R. Arora, S. A. Bass, L. Cunqueiro,
10, 11
T. Dai, L. Du, R. Ehlers,
10, 11
H. Elfner,
13, 14, 15
D. Everett, W. Fan, R. J. Fries,
16, 17
C. Gale, F. Garza,
16, 17
Y. He, M. Heffernan, U. Heinz, B. V. Jacak,
5, 6
S. Jeon, W. Ke,
5, 6
B. Kim,
16, 17
M. Kordell II, A. Kumar, A. Majumder, S. Mak, M. McNelis, C. Nattrass, D. Oliinychenko, C. Park,
18, 1
J.-F. Paquet, J. H. Putschke, G. Roland, A. Silva, B. Schenke, L. Schwiebert, C. Shen,
1, 23
C. Sirimanna, Y. Tachibana,
1, 24
G. Vujanovic, X.-N. Wang,
19, 5, 6
R. L. Wolpert, and Y. Xu (The JETSCAPE Collaboration) Department of Physics and Astronomy, Wayne State University, Detroit MI 48201. Institute of Frontier and Interdisciplinary Science,Shandong University, Qingdao, Shandong 266237, China Laboratory for Nuclear Science, Massachusetts Institute of Technology, Cambridge MA 02139. Department of Statistical Science, Duke University, Durham NC 27708. Department of Physics, University of California, Berkeley CA 94270. Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley CA 94270. Lawrence Livermore National Laboratory, Livermore CA 94550. University of Texas at San Antonio, 1 UTSA Circle, San Antonio, TX 78249. Department of Physics, Duke University, Durham NC 27708. Department of Physics and Astronomy, University of Tennessee, Knoxville TN 37996. Physics Division, Oak Ridge National Laboratory, Oak Ridge TN 37830. Department of Physics, The Ohio State University, Columbus OH 43210. GSI Helmholtzzentrum f¨ur Schwerionenforschung, 64291 Darmstadt, Germany. Institute for Theoretical Physics, Goethe University, 60438 Frankfurt am Main, Germany. Frankfurt Institute for Advanced Studies, 60438 Frankfurt am Main, Germany. Cyclotron Institute, Texas A&M University, College Station TX 77843. Department of Physics and Astronomy, Texas A&M University, College Station TX 77843. Department of Physics, McGill University, Montr´eal QC H3A-2T8. Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China. Department of Physics, Massachusetts Institute of Technology, Cambridge MA 02139. Department of Physics, Brookhaven National Laboratory, Upton NY 11973. Department of Computer Science, Wayne State University, Detroit MI 48202. RIKEN BNL Research Center, Brookhaven National Laboratory, Upton NY 11973. Akita International University, Yuwa, Akita-city 010-1292, Japan (Dated: February 24, 2021)We report a new determination of ˆ q , the jet transport coefficient of the Quark-Gluon Plasma.We use the JETSCAPE framework, which incorporates a novel multi-stage theoretical approach toin-medium jet evolution and Bayesian inference for parameter extraction. The calculations, basedon the Matter and
Lbt jet quenching models, are compared to experimental measurements ofinclusive hadron suppression in Au + Au collisions at RHIC and Pb + Pb collisions at the LHC. Thecorrelation of experimental systematic uncertainties is accounted for in the parameter extraction.The functional dependence of ˆ q on jet energy or virtuality and medium temperature is based on aperturbative picture of in-medium scattering, with components reflecting the different regimes ofapplicability of Matter and
Lbt . In the multi-stage approach, the switch between
Matter and
Lbt is governed by a virtuality scale Q . Comparison of the posterior model predictions to theRHIC and LHC hadron suppression data shows reasonable agreement, with moderate tension inlimited regions of phase space. The distribution of ˆ q/T extracted from the posterior distributionsexhibits weak dependence on jet momentum and medium temperature T , with 90% Credible Region(CR) depending on the specific choice of model configuration. The choice of Matter + Lbt , withswitching at virtuality Q , has 90% CR of 2 < ˆ q/T < p T , jet >
40 GeV /c . The value of Q ,determined here for the first time, is in the range 2.0-2.7 GeV. I. INTRODUCTION
The Quark-Gluon Plasma (QGP) is the state of mat-ter in conditions of extreme temperature and density,similar to those of the universe a few microseconds af- ter the Big Bang, with structure and dynamics governedby interactions of sub-hadronic quanta [1]. The QGP isgenerated in the laboratory by collisions of heavy nuclei( A + A ) at the Relativistic Heavy Ion Collider (RHIC)and the Large Hadron Collider (LHC), and its properties a r X i v : . [ nu c l - t h ] F e b have been measured extensively by large experiments atthose facilities ([2] and references therein). Comparisonof these measurements with theoretical calculations showthat the QGP exhibits collective behavior with very lowspecific viscosity [3]. The QGP is likewise found to beopaque to penetrating probes carrying color charge, aphenomenon known as “jet quenching” [4].Jets in high-energy hadronic collisions are generatedby hard (high Q ) interactions of partons (quarks andgluons) from the incoming projectiles. The scattered par-tons are initially highly virtual and evolve through QCDbremsstahlung and pair production, forming a partonshower that hadronizes into a collimated spray of stableparticles that is observable experimentally (a “jet”). Jetproduction rates and jet structure have been measuredextensively in p (¯ p ) + p collisions [5–12], and theoreticalcalculations based on perturbative QCD (pQCD) are inexcellent agreement with such measurements [13–17].Jets are likewise produced in high-energy A + A col-lisions, in parallel with formation of the QGP. Jets gen-erated in A + A collisions propagate during their showerevolution through the QGP and interact with it, therebymodifying final-state jet structure and jet distributionsrelative to production in vacuum. Measurements of thesemodifications provide unique and sensitive probes of theQGP.A key experimental observable of jet quenching is thesuppression of inclusive hadron production at high trans-verse momentum (high p T ) [4, 18]. Hadron suppressionis measured using the ratio R AA of the inclusive hadronyield in A + A collisions to that in a reference system, usu-ally pp collisions at the same collision energy, whose yieldis scaled to account for nuclear geometry [19]; R AA = 1corresponds to the absence of nuclear effects in high- p T hadron production. Such effects arise from both nuclearParton Distributions Functions and from jet quenchingin the QGP. Inclusive hadron suppression at high p T has been measured extensively at RHIC [20, 21] and theLHC [22–25] .Jet quenching is understood theoretically to arise fromelastic and inelastic interactions of the partons in the jetshower as it traverses the QGP, with coherence effectsplaying an important role [4, 18, 26–30]. Several differ-ent formalisms have been developed to describe this pro-cess, as reviewed in [31]. Calculations based on these for-malisms have been carried out for inclusive hadron sup-pression [32–35], di-hadron production [36–38], γ -hadroncorrelations [39–43], reconstructed jets [41, 44–49], andjet substructure [50–54].These formalisms are based on various approximationsthat are applicable in limited ranges of shower energyand virtuality scales. Several of these formalisms have re-cently been implemented in a unified analysis framework,called JETSCAPE [55], providing a multi-stage model ofjet evolution in which each jet quenching formalism is ap-plied only in its appropriate range of validity in showerenergy and virtuality.Models of in-medium jet-thermal parton interactions have parameters, known as jet transport coefficients, thatcan be determined by comparison of their calculations tojet quenching measurements [27–29, 56–65]. Phenomeno-logically, the most significant transport coefficient is ˆ q ,which denotes the mean square of the momentum trans-fer between the propagating hard jet and the soft mediumper unit length, ˆ q ≡ d (cid:104) k ⊥ (cid:105) /dL , where (cid:104) . . . (cid:105) indicates theaverage over all jet propagation paths for the event pop-ulation.Quantitative determination of ˆ q by comparison of the-ory models with experimental data has previously beencarried out by the JET Collaboration, which comparedseveral jet quenching model calculations to inclusivehadron R AA measurements at RHIC and the LHC [66].This was the first extraction of ˆ q using data from bothcolliders, which are expected to generate a QGP withdifferent initial temperature. Non-perturbative contribu-tions to the value of ˆ q can also be calculated using first-principles lattice QCD [67] though challenges remain toinclude quark dynamics into a full QGP calculation. Arecent 2+1 flavor calculation with physical quark masseson N τ =8 lattices [68, 69] yields a value for ˆ q/T in therange of 2.5–3.5 for the highest temperatures generatedin heavy-ion collisions at RHIC and the LHC, consistentwith the results reported by the JET Collaboration.Each model in the JET calculation of ˆ q has a single freeparameter, either ˆ q or the effective strong coupling pa-rameter α s , which is determined separately for RHIC andLHC data. In order to go beyond separate extractions ofˆ q at RHIC and LHC, and instead obtain a distributionfor ˆ q that is a smooth function of the medium temper-ature and jet energy, parameter extraction incorporat-ing data from both colliders is required. However, thatapproach is beyond the scope of the least-squares mini-mization approach used by the JET analysis, and a morecomprehensive approach is needed, based on Bayesianinference [70–72].Bayesian inference has been utilized previously to ex-tract parameters of the QGP from heavy-ion collisiondata, in particular the specific shear viscosity η/s [73–76] and the Equation of State (EoS) [77], with the latterresult agreeing well with Lattice QCD calculations. Seealso [78]. This approach has likewise been used to studythe heavy quark diffusion coefficient of the QGP [79].The theoretical description of jet modification includesloss of energy-momentum by hard partons in the shower,generation of softer partons by medium excitation, andexcitation of the medium due to energy and momentumexchanges at a scale below that describable by pertur-bative techniques. A complete description of the entirejet, in terms of all available jet observables, requires themodeling of several transport coefficients. Focus on high- p T hadrons constrains discussion to the hardest partonsin the shower, whose distribution is modified by (mainlytransverse) momentum exchanges with the medium thatdepends primarily on ˆ q .This manuscript presents a quantitative extractionof the temperature and momentum dependence of ˆ q in the QGP, using Bayesian inference methods in theJETSCAPE framework. The analysis extends that of theJET collaboration [66] by determining the functional de-pendence of ˆ q on the jet energy and virtuality, and the lo-cal temperature. Comparison of model calculations withdata from both RHIC and the LHC provides a broad scanin jet energy and virtuality. Two collision centralities areutilized, providing variation in the medium temperatureprofile.As the jet shower propagates through the QGP, it ex-changes energy and momentum with the dense medium.Momentum exchanges above a certain scale are describedusing perturbative QCD (pQCD), whereas softer mo-mentum exchanges are modeled by partial thermalizationof the exchanged four-momentum with a hydrodynamicbackground. This calculation incorporates multi-stagejet evolution, using the Matter model [80, 81] at highvirtuality scale, and the
Lbt model [35, 42, 43, 55] atlow virtuality scale, with the switching between modelsgoverned by a free parameter Q . Two different analyticparametrizations of ˆ q are explored which are functions ofthe medium temperature and either jet energy or jet vir-tuality, and which are based on perturbative treatmentof jet-medium interactions.The experimental data used to for Bayesian parame-ter extraction are measurements of inclusive hadron sup-pression in central and semi-central Au + Au collisions at √ s NN =200 GeV for p T , h > /c [21], and Pb + Pbcollisions at √ s NN =2.76 and 5.02 TeV for p T , h > /c [23, 24].The manuscript is organized as follows: Sect. IIpresents the jet evolution models; Sect. III presents theˆ q parametrizations; Sect. IV presents the experimen-tal data and treatment of their uncertainties; Sect. Vpresents training of the Gaussian process emulator;Sect. VI presents implementation of Bayesian inference;Sect. VII presents closure tests; Sect. VIII presents theresults in terms of posterior distributions for ˆ q from theBayesian parameter extraction; and Sect. IX gives a sum-mary and outlook. II. MODELING JET-MEDIUM INTERACTIONS
JETSCAPE provides a general numerical frameworkfor simulating jet-medium interactions, with several al-ternative models to simulate the QGP and jet evolution.QGP evolution is modeled using relativistic hydrodynam-ics; in this study, jet evolution is calculated using the
Matter model to describe interactions with the QGPat high virtuality and the
Lbt model to describe inter-actions with the QGP at low virtuality.
Matter and
Lbt are also combined in a multi-stage approach. Thissection discusses each model in turn.
A. QGP evolution
The evolution of the QGP is simulated with second-order relativistic hydrodynamics as implemented inVISH2+1 [82–86]. An initial entropy density profilefrom the Monte-Carlo Glauber model [19, 87] is evolvedwith dissipative fluid dynamics, starting at longitudi-nal proper time τ = 0 . c . For Au-Au colli-sions at √ s NN = 200 GeV and Pb-Pb collisions at √ s NN = 2 .
76 TeV, a partial chemical equilibrium equa-tion of state was used with T chem = 165 MeV (s95p-v0-PCE165 [88, 89]), along with a constant specific shearviscosity η/s = 0 .
08 and no bulk viscosity; these choiceswere found to provide a good description of data atRHIC and the LHC in Ref. [90]. For Pb-Pb collisions at √ s NN = 5 .
02 TeV, the hydrodynamic profiles were latertuned in Ref. [86] with a different equation of state andshear viscosity: partial chemical equilibrium equation ofstate with T chem = 150 MeV (s95p-v1-PCE150 [88, 89])and a temperature dependent η/s [86]. The hydrody-namic model provides the space-time evolution profileof the temperature T and flow velocity u µ ; the viscouspart of the energy-momentum tensor is not used. Thein-medium jet evolution equations are solved using theseprofiles, as discussed in the following subsections.In this work, we apply event-averaged hydrodynamicprofiles for jet evolution and only consider jet energyloss inside the QGP phase; interactions in both the briefpre-hydrodynamic stage ( τ < . c ) and the dilutehadronic stage ( T < T stop ) are neglected. Jet-mediuminteractions are stopped at T stop = 165 MeV, close thethe pseudo-critical chiral temperature T c . If a lower valueof T stop is used, the jet energy loss will increase and asmaller value of ˆ q will be extracted from the jet quench-ing data, though the effect is expected to be small dueto the minor enhancement of jet energy loss at such lowtemperature. B. Jet production
Energetic jets are generated in hard scatterings,with rate based on a leading-order perturbative QCD(LO pQCD) calculation in momentum space using theCTEQ5 parametrization of parton distribution func-tions [91]. In nucleus-nucleus collisions, the EPS09parametrizations of nuclear PDF modification [92] aretaken into account for the momentum space distributionof hard partons. Their position space distribution is sam-pled according to the Monte-Carlo (MC) Glauber model.
C. Jet-medium interaction at high virtuality:MATTER
The M odular A ll T wist T ransverse-scattering E lastic-drag and R adiation ( Matter ) model [80, 81, 93] simu-lates the splitting of highly virtual partons, i.e. jet par-tons whose virtuality is much larger than the multiple-scattering scale of the medium it probes ( ∼ √ ˆ qE ), where E is the jet energy. At high virtuality, the number ofsplittings dominates over the number of scatterings insidethe medium, and the parton splitting process is describedby a medium-modified virtuality-ordered shower [64, 94–96], where the scattering in the medium provides an ad-ditional contribution to the splitting functions.A jet shower is initiated by a single hard parton pro-duced at a point r with a forward light-cone momentum p + = ( p + ˆ n · (cid:126)p ) / √
2, in which ˆ n = (cid:126)p/ | (cid:126)p | specifies thedirection of the jet. The virtuality ( Q ) of a particularparton is sampled based on the Sudakov form factor thatdetermines the virtuality distribution [80, 97],∆( Q , Q ) = (cid:89) a ∆ a ( Q , Q )= (cid:89) a exp − Q (cid:90) Q dQ Q α s ( Q )2 π − z c (cid:90) z c dzP a ( z, Q ) . (1)Here, a represents the channels through which the jetparton can split, Q varies from the maximum possiblevalue Q max that initiates at the parton energy downto the minimum allowed value of Q below which thevirtuality-order parton shower breaks. In the equationabove, z c is taken as Q /Q , and the splitting func-tion contains both vacuum and medium-induced contri-butions, P a ( z, Q ) = P vac a ( z ) + P med a ( z, Q ) . (2)Here, the medium-induced part is adopted from thehigher-twist formalism [60, 64, 98, 99] and treated as aperturbation to the vacuum part: P med a ( z, Q ) = P vac a ( z ) z (1 − z ) Q ζ +max (cid:90) dζ + ˆ q g ( r + ζ ) × (cid:34) − (cid:32) ζ + τ + f (cid:33) − ζ + τ + f sin (cid:32) ζ + τ + f (cid:33) + 2 (cid:32) ζ + τ + f (cid:33) cos (cid:32) ζ + τ + f (cid:33) (cid:35) . (3)Here ˆ q g denotes the gluon transport coefficient; it is eval-uated locally at (cid:126)r + ˆ nζ + and is related to the quark trans-port coefficient ˆ q q by color factors. The maximum lengthsampled ζ +MAX is taken as 1 . τ + f , where τ + f = 2 p + /Q isthe mean light-cone formation time.After Q of the parent parton is determined, z is chosenby sampling the splitting function P ( z ). The maximumpossible virtualities of the two daughters are thus zQ and(1 − z ) Q , from which the virtualities of the two daugh-ters Q and Q are similarly assigned by sampling the form factor in Eq. (1). The transverse (to ˆ n ) momen-tum of the produced pair is then calculated according tothe difference in invariant mass between the parent anddaughters: k ⊥ = z (1 − z ) Q − (1 − z ) Q − zQ . (4)The actual time/length for each splitting process issampled using a Gaussian distribution with a mean valueof τ + f [81]. This process is iterated until virtualities ofall partons within the jet shower reaches the predeter-mined value of Q . This virtuality-ordered parton showermethod is similar to the time-like shower implemented in Pythia , except that here the medium effect is includedin a consistent way.Hard partons evolve through multiple splittings in
Matter starting with maximum possible virtualities( Q = E ) until their virtualities drop below Q . When the Matter model is applied alone, Q is fixed at 1 GeV.For proton-proton collisions, only the vacuum contribu-tion to the splitting function Eq. (2) is taken into account.As shown in Ref. [81], this approach provides a good de-scription of the single inclusive hadron and jet spectraat high p T in proton-proton collisions, serving as a reli-able baseline for studying their nuclear modification inheavy-ion collisions. For nucleus-nucleus collisions, boththe vacuum and medium-induced parts are implemented.At Q = 1 GeV, all partons are converted into hadronsusing Pythia fragmentation.Partons are fragmented independently using the py1ent function of PYTHIA [100]. We note that there issizable uncertainty in parametrized fragmentation func-tions at LHC energies [101]. However, since the com-bined calculation of initial production and hadronizationin JETSCAPE accurately describes jet spectra in proton-proton collisions [102], we assume that is can also be usedto calculate in-medium modification in heavy-ion colli-sions.For the medium-induced part of the splitting func-tion in Eq. (2), the local fluid velocity of the dynamicalmedium is taken into account by rescaling ˆ q in Eq. (3) viaˆ q = ˆ q local · p µ u µ /p [103], where p µ is the four-momentumof the jet. The value of ˆ q is zero before jets enter the ther-mal medium ( τ < . q local asˆ q . It should be understood that a boost is invoked withinthe calculation for the case of a moving frame.The jet transport coefficient of the QGP medium isthe sole parameter of the Matter model. As discussedin [81], the minimal assumption that ˆ q is proportional tothe entropy density s in the local rest frame, ˆ q/s = ˆ q /s ,is able to describe single inclusive hadron and jet R AA ,but distinct values of ˆ q at a given reference point s are required at RHIC and the LHC. The present workexplores a more general form of ˆ q as a function of mediumtemperature, jet energy, and virtuality scale, which canbe uniformly applied to data from RHIC and the LHC. D. Jet-medium interactions at low virtuality: LBT
The Linear Boltzmann Transport (
Lbt ) model calcu-lates the time evolution of jets in a thermal QGP gen-erated in relativistic heavy-ion collisions, by using a ki-netic approach that includes elastic and inelastic colli-sions [35, 41–43, 47, 104]. The evolution of the phasespace distribution of a jet parton a with p µa = ( E a , (cid:126)p a ) isdescribed using the Boltzmann equation p a · ∂f a ( x a , p a ) = E a ( C el a + C inel a ) , (5)where C el a and C inel a are the collision integrals for elasticand inelastic scatterings.For elastic scattering of a with a thermal parton b fromthe medium background, the collision term C el a is evalu-ated with the leading-order vacuum matrix elements forall possible ab → cd channels. The collinear ( u, t → S ( s, t, u ) = θ ( s ≥ µ ) θ ( − s + µ ≤ t ≤ − µ ), inwhich µ = g T ( N c + N f / / a reads:Γ el a = (cid:88) b,c,d γ b E a (cid:90) (cid:89) i = b,c,d d [ p i ] f b ( (cid:126)p b ) S ( s, t, u ) × (2 π ) δ (4) ( p a + p b − p c − p d ) |M ab → cd | , (6)where d [ p i ] = d p i / [2 E i (2 π ) ], and γ b and f b representthe spin-color degeneracy and thermal distribution ofparton b respectively. The probability of elastic scatter-ing of parton a in each time step ∆ t is thus P el a = Γ el a ∆ t .During the propagation of a jet through the QGPmedium, each elastic scattering changes its transversemomentum (perpendicular to its initial direction). Thisresults in an increase of the average transverse momen-tum square (or transverse momentum broadening) overtime. This transverse momentum broadening per unittime (or length) is known as the jet transport coefficientˆ q , characterizing both the local thermal parton densityand the strength of jet-medium interaction. In LBT thevalue of ˆ q may be inferred by evaluating Eq. (6) weightedby the transverse momentum broadening of the jet par-ton:ˆ q a = (cid:88) b,c,d γ b E a (cid:90) (cid:89) i = b,c,d d [ p i ] [ (cid:126)p c − ( (cid:126)p c · ˆ p a )ˆ p a ] f b ( (cid:126)p b ) × S ( s, t, u )(2 π ) δ (4) ( p a + p b − p c − p d ) × |M ab → cd | (7) In addition to elastic scattering, inelastic scatteringwhich generates medium-induced gluon radiation is alsoincluded in the Lbt model. The inelastic scattering rateat a given time t is defined as the average number of emit-ted gluons from parton a per unit time, and is evaluatedas [104–106]Γ inel a ( E a , T, t ) = 11 + δ ag (cid:90) dzdk ⊥ dN ag dzdk ⊥ dt , (8)in which the δ ag term is imposed to avoid double countingfor the g → gg process. The medium-induced gluon spec-trum is taken from the higher-twist energy loss formalism[60, 64, 107], dN ag dzdk ⊥ dt = 2 α s ( k ⊥ ) P vac a ( z ) k ⊥ π ( k ⊥ + x m a ) ˆ q g sin (cid:18) t − t i τ f (cid:19) . (9)Here, z and k ⊥ are the fractional energy and trans-verse momentum of the emitted gluon with respect toits parent parton a , and P vac a ( z ) is the vacuum split-ting function, and ˆ q g is the gluon transport coefficient.The initial time t i denotes the production time of theparent parton a from which the gluon is emitted, and τ f = 2 E a z (1 − z ) / ( k ⊥ + z m a ) is the formation time ofthe radiated gluon with m a being the mass of the parton.In this work, we assume zero mass for light quarks andgluons. Note that the gluon spectrum is proportional toˆ q g , which is related to the medium parameters in Lbt through Eq. (7). To avoid possible divergence as z → z min = 2 πT /E is imple-mented for the energy of the emitted gluon [105]. Notethat Eq. (9) is consistent with the medium-induced split-ting function Eq. (3) used in Matter , except that the( ζ/τ f ) and ( ζ/τ f ) terms are ignored here in Lbt . Thecontribution from these two terms has been discussedin [80] and shown to be small when ζ (cid:46) τ f . Multiplegluon emissions are allowed in each time step ∆ t . Differ-ent medium-induced gluon emissions are assumed to beindependent of each other; their number n is therefore aPoisson distribution with mean as (cid:104) N ag (cid:105) = Γ inel a ∆ t , P ( n ) = (cid:104) N ag (cid:105) n n ! e −(cid:104) N ag (cid:105) n . (10)Thus, the probability of an inelastic scattering processoccurring is P inel a = 1 − e −(cid:104) N ag (cid:105) . Interference betweenmultiple gluon emissions have not been taken into ac-count yet, which should be improved in our future work.To combine the elastic and inelastic processes, the to-tal scattering probability is divided into two regions: pureelastic scattering with probability P el a (1 − P inel a ) and in-elastic scattering with probability P inel a . The total scat-tering probability is thus P tot a = P el a + P inel a − P el a · P inel a .Based on these probabilities, the Monte Carlo approachis used to determine whether a given jet parton a scat-ters in the thermal medium, and whether the scattering ispurely elastic or inelastic. For a selected scattering chan-nel, the energy and momentum of the outgoing partonsare sampled using the corresponding differential spectragiven by Eq. (6) and (9).For realistic nuclear collisions, we initialize jet partonsfrom hard scatterings in the same way as for the Matter model (Sect. II C). To account for the effect of mediumflow effect on jet transport during the QGP phase, in eachtime step we first boost each jet parton into the local restframe of the fluid cell in which its energy and momentumare updated based on the
Lbt model, and then boost itback to the global collision frame where it propagatesto the spacetime of the next time step. Note that theprevious rescaling of ˆ q in Matter ( p ˆ q = ( p · u )ˆ q local )has the same effect as the boost method here when themedium-induced gluon spectrum is written in a boost-invariant form. On the freeze-out hypersurface of theQGP ( T stop = 165 MeV), high p T jet partons are passedto Pythia 6 for conversion into hadrons.In the original work of
Lbt the strong coupling con-stant α s was the sole parameter determining both theelastic [Eq. (6)] and inelastic [Eq. (8)] scattering pro-cesses. In this work we take an alternative approach,parametrizing ˆ q directly (Sect. III). E. Multi-stage evolution with MATTER+LBT
Both the medium-modified virtuality shower with
Matter alone, and the vacuum virtuality shower with
Lbt low-virtuality parton transport, can be used to de-scribe the modification of inclusive hadron and jet distri-butions in heavy-ion collisions. However, the applicationof either of these models alone to the entire jet evolutionis not theoretically complete; the
Matter formalism isnot applicable for parton virtuality below the mediumscale, and
Lbt ignores the in-medium modification of jetsin the highly virtual stage. JETSCAPE has therefore de-veloped a multi-stage approach to calculating in-mediumjet evolution, in which
Matter is applied for partonswith high virtuality and
Lbt is applied for partons withlow virtuality [55].For an energetic parton generated by a hard scatter-ing,
Matter is used to simulate its virtuality-orderedsplitting process (Sect. II C). In each splitting, the vir-tuality of each of the daughter partons is smaller thanthat of the parent. When the virtuality of a parton inthe shower falls below a specified scale Q , it is passed to Lbt for the subsequent time-ordered in-medium evolu-tion. In this combined approach,
Matter largely deter-mines the spectrum of final state partons for high-energyjets or short in-medium path length, while
Lbt predom-inantly governs low energy parton scattering, especiallywhen the in-medium path length is large [55]. A key parameter of this multi-stage approach is theseparation (or switching) scale Q between Matter and
Lbt evolution, whose value is expected to be similar tothat of the medium scale, Q ∼ ˆ qτ f . A similar separa-tion scale was explored in [108]. This analysis providesthe first phenomenological determination of Q from ex-perimental data, using Bayesian inference.Note that in the Lbt stage we assume zero virtualityfor thermal QGP partons, and recoiling partons are scat-tered out of the background medium by the jet. In con-trast, the virtuality of jet partons is fully tracked since itis fed from
Matter . When a parton splits in
Lbt , we as-sume that its two daughters share its virtuality in propor-tion to their z fraction. At temperature below T stop ∼ T c , Lbt partons with
Q <
Pythia 6 , while those with virtualities still above1 GeV are passed back to
Matter for subsequent vac-uum showering until all partons satisfy
Q < τ < . T Matter vacuum shower is appliedwith the value of Q = 1 GeV. This same vacuum showeris used for p + p collisions. It is only in the hydrodynamicstage ( τ > . T > T stop ) of A + A collisions thatthe Matter ( Q > Q ) + Lbt ( Q < Q ) model is usedwith a Q value greater than 1 GeV. III. ˆ q PARAMETRIZATION As discussed in Sec. II, the jet transport coefficientˆ q is the sole quantity constrained by fitting Matter and Lbt calculations to experimental data. We em-ploy three different parametrizations of ˆ q . The form ofthe parametrization is derived from Eq. (7), which isbased on perturbative scattering of a jet parton insidea medium. Assuming a thermal distribution for f b andtaking the small-angle approximation for elastic scatter-ing gives [109–112]ˆ q ≈ C R ζ (3) π α T ln (cid:18) CET µ (cid:19) , (11)where C R is the color factor of the jet parton (4/3 forquark and 3 for gluon), T is the medium temperature,and C is a constant depending on the kinematic cutsimplemented in Eq. (11). For the Lbt model using aconstant α s = 0 . C is approximately 5.6 for gluons and5.8 for quarks [112] . In this work ˆ q refers to the light-quark jet transport coefficient; the gluon jet transportcoefficient is obtained by scaling with the relative colorfactor.The definition of ˆ q need not be limited to perturbativescattering of a jet parton with a thermal medium at thescale T . Thus, we extend Eq. (11) to a more general formas follows:ˆ q ( E, T ) | A,B,C,D T = 42 C R ζ (3) π (cid:18) π (cid:19) (cid:40) A (cid:2) ln (cid:0) E Λ (cid:1) − ln( B ) (cid:3)(cid:2) ln (cid:0) E Λ (cid:1)(cid:3) + C (cid:2) ln (cid:0) ET (cid:1) − ln( D ) (cid:3)(cid:2) ln (cid:0) ET Λ (cid:1)(cid:3) (cid:41) , (12)where ( A, B, C, D ) are parameters that will be deter-mined from the experimental data using Bayesian pa-rameter extraction. If the first part in the braces { ... } (orparameter A ) is set to zero, the second part reduces toEq. (11) if the coupling constant α s = 4 π/ / ln( ET / Λ )is assumed to run with both jet energy and medium tem-perature scales at leading order. For the parameter Λ weuse Λ = 0 . ∼ T ), the value of ˆ q is controlledsolely by the scale of the jet parton itself, and not bythe medium temperature. This first part in { ... } , withparameters A and B , represents the physics assumed bythe Matter model. The second part in { ... } is expectedto represent the physics of an on-shell jet parton scat-tering with quasi-particles inside a thermal medium, asassumed by the Lbt model. The arguments of the log-arithms in Eq. (12) involve additional constant factorsthat depend on the particular cut-off value implementedin the t -channel scattering. We treat these as parameters, called B and D .We consider Eq. (12) to be a a sufficiently generalansatz of the energy and momentum dependence of ˆ q within the perturbative picture of jet-medium interac-tion. We use this parametrization consistently in both Matter and Lbt when they are applied separately todescribe experimental data. We expect that the physicsof the high virtuality stage in Matter is described pre-dominantly by the first term (with A and B ), while thephysics of the thermal stage in Lbt is described by thesecond term (with C and D ).For the multi-stage calculation combining Mat-ter + Lbt we utilize two different parametrizations of ˆ q .The first parametrization uses Eq. (12) to calculate ˆ q in both Matter and Lbt stages, while introducing anadditional parameter Q that represents the virtualityboundary between the two stages. This five-parameterformulation is denoted “ Matter + Lbt q , it can becompared directly to the parametrization in which Mat-ter and Lbt are applied separately.To reduce the number of parameters and capture thejet physics of virtuality evolution in Matter more pre-cisely, we introduce a second ˆ q parametrization for themulti-stage Matter + Lbt model, as follows:ˆ q ( Q, E, T ) | Q ,A,C,D T = 42 C R ζ (3) π (cid:18) π (cid:19) A (cid:104) ln (cid:16) Q Λ (cid:17) − ln (cid:16) Q Λ (cid:17)(cid:105)(cid:104) ln (cid:16) Q Λ (cid:17)(cid:105) θ ( Q − Q ) + C (cid:2) ln (cid:0) ET (cid:1) − ln( D ) (cid:3)(cid:2) ln (cid:0) ET Λ (cid:1)(cid:3) . (13)This parametrization is denoted “ Matter + Lbt Q asthe scale in the first term instead of the jet energy E .The motivation behind this parametrization is that the Matter model better characterizes the parton showeras a function of virtuality. However, the value of ˆ q de-termined using this parametrization cannot be directlycompared to that from Lbt .The parameter B is replaced by the switching virtu-ality Q , so that this formulation likewise has four pa-rameters. The θ function ensures that, during the Mat-ter stage ( Q > Q ), ˆ q receives contributions from bothterms, while during the Lbt stage ( Q < Q ) only thesecond term contributes. In this parametrization, thedistribution of ˆ q is continuous at Q = Q . IV. EXPERIMENTAL DATA This analysis carries out Bayesian parameter ex-traction using experimental measurements of inclusivehadron production in A + A collisions at RHIC and LHC( R AA ). Selection of experimental data for this process re-quires consideration of the p T range suitable for compari-son to theoretical calculations of jet quenching, in partic-ular the possible role of medium-modified hadronizationat low p T .The energy loss formalism in this manuscript involvesthe convolution of initial state and hard scattering dis-tributions with energy loss calculations applied to hardpartons as they propagate through the medium. The fi-nal parton distributions are then convoluted with vacuumfragmentation functions to calculate the hadron distribu-tions to be compared to data. The calculations are there-fore based on the assumption that the hadronization ofleading hadrons takes place outside the dense medium.The space-time distribution of jet hadronization in thepresence of the QGP is currently an open issue, to beresolved using both experimental data and theoreticalmodeling. Relevant experimental data to address thisquestion include high- p T di-hadron correlations in centralAu + Au collisions at RHIC [113], whose jet-like angulardistributions indicate that charged hadrons with p T (cid:38) /c arise predominantly from vacuum fragmentation;and particle-identified relative yields in reconstructed jetsin central Pb + Pb collisions at the LHC, which are simi-lar to those for jets in vacuum for p T > /c , in con-trast to a striking enhancement in the baryon/meson ra-tio for bulk (non-jet) production [114]. These experimen-tal observations suggest that hadrons with p T (cid:38) /c are generated in central A + A collisions predominantlyby jet fragmentation in vacuum; in other words, that theprocesses of jet-medium interactions and hadron forma-tion largely factorize for hadrons with p T > /c . Onthe other hand, parametric theoretical arguments suggestthat vacuum hadronization occurs only for hadrons with p T > 10 GeV /c , at both RHIC and the LHC.In order to simplify the analysis presented in thismanuscript, we therefore restrict the p T range of the in-clusive hadron R AA measurements considered for com-parison to the theory calculations to p T > /c atRHIC and p T > 10 GeV /c at the LHC. We note thatthis cut limits significantly the statistical weight of theRHIC data relative to that at the LHC, due to the muchnarrower kinematic range accessible at the lower √ s NN ofRHIC (see Sect. VIII). Lowering of this limit, to enablegreater statistical weight of RHIC data, will be exploredin future work.The experimental datasets used in this analysis, whichcover a wide range in hadron p T and medium tempera-ture, are as follows: • Au-Au collisions at √ s NN = 200 GeV, 0-10% and40-50% centrality [21]; • Pb-Pb collisions at √ s NN = 2 . 76 TeV, 0-5% and30-40% centrality [23]; • Pb-Pb collisions at √ s NN = 5 . 02 TeV, 0-10% and30-50% centrality [24]. A. Experimental uncertainties Bayesian parameter extraction requires specification ofexperimental uncertainties, optimally the full covariancematrix Σ E . However, the full covariance matrix of mea-surement uncertainty is difficult to determine, and it isusually not reported in experimental publications. We fo-cus here on measurements of inclusive hadron R AA , anddiscuss how the covariance matrix is determined for thereported measurements used in this analysis.For these measurements, the experimental uncertain-ties are specified as a function of hadron p T . The CMS [24] and ATLAS [23] publications report the fol-lowing four types of uncertainty:1. uncorrelated statistical error and systematic uncer-tainty on each data point;2. luminosity uncertainty, fully correlated in all cen-trality bins for a given collision system;3. Glauber scaling ( (cid:104) T AA (cid:105) ) uncertainty, fully corre-lated in p T for a given collision system and cen-trality bin; and4. other correlated errors of unspecified origin, withonly qualitative dependence on hadron p T specified.Because the luminosity and (cid:104) T AA (cid:105) uncertainties are in-dependent of p T , it is straightforward to calculate theircontribution to the off-diagonal terms of the covariancematrix. However, the other correlated uncertainties arisefrom sources such as track selection, momentum resolu-tion, and efficiency correlations, which vary in differentways with p T . To account for this complexity we intro-duce a correlation-length parameter (cid:96) (defined below) torepresent the range in p T over which these uncertaintiescontribute.For the R AA distribution of a specific collision systemand centrality from a specific experimental publication(indexed by k ), let Σ Ek be the corresponding covariancematrix block constructed from “uncorrelated”, “fully-correlated”, and “length-correlated” uncertainty vectorsrespectively, { σ uncorr k } , { σ fcorr k } , { σ lcorr k } , as reported bythe experiments. Then the uncertainty covariance blockΣ Ek is given byΣ Ek = Σ uncorr k + Σ fcorr k + Σ lcorr k Σ uncorr k,ij = σ uncorr k,i σ uncorr k,j δ ij Σ fcorr k,ij = σ fcorr k,i σ fcorr k,j Σ lcorr k,ij = σ lcorr k,i σ lcorr k,j exp (cid:20) − (cid:12)(cid:12)(cid:12)(cid:12) p k,i − p k,j (cid:96) k (cid:12)(cid:12)(cid:12)(cid:12) α (cid:21) . (14)Here p k,i is the i th p T value in block k , and δ ij = 1 if i = j and 0 otherwise. Thus, Σ uncorr k is a diagonal ma-trix, representing the combined, uncorrelated statisticaland systematic experimental uncertainties. Σ fcorr k corre-sponds to the fully correlated, p T -independent luminosityand (cid:104) T AA (cid:105) uncertainties, and Σ lcorr k is constructed fromthe correlated experimental uncertainties using a powerexponential covariance function. We set the exponent α = 1 . 9, similar to the common choice α = 2 but com-putationally more stable [115]. The p k,i transverse mo-mentum values and correlation length (cid:96) k in Eq. (14) arelinearly rescaled so that all values lie within [0,1]. Therescaled correlation length (cid:96) k is nominally set to a valueof (cid:96) k = 0 . 2. Other values were used to study the sensi-tivity of our results to this parameter choice.The PHENIX publication [21] reports uncertainties ina similar fashion but with different labels. Uncorrelatederrors are denoted as Type A, fully correlated and p T -independent are reported as Type C, and Type B refersto correlated systematic errors with an unspecified p T -dependence. Therefore Type A, B, and C errors aretreated as σ uncorr , σ lcorr , σ fcorr , respectively, accordingto Eq. (14). V. GAUSSIAN PROCESS EMULATORS Because JETSCAPE model calculations are computa-tionally expensive, we use Gaussian Process Emulators(GPEs) to interpolate the model-parameter space [70,116]. GPEs offer a non-parametric method of regression,providing a statistical surrogate for the computationallyexpensive model by using a limited set of training pointsto predict with a defined uncertainty any untried value ofinput parameters. This allows us to make rigorous sta-tistical comparisons to the experimental data efficiently,and perform inference on the input parameters.Our implementation of the GPE is identical to thatused in [117], except for the choice of the covariancefunction which controls the correlation between pairs ofpoints. To improve emulator stability, we replace thesquared-exponential function used in [117] with a Mat´ern5/2 covariance function, c ( x i , x j ) = (cid:32) √ d i,j (cid:96) + 5 d i,j (cid:96) (cid:33) exp (cid:32) − √ d i,j (cid:96) (cid:33) (15)where d i,j = | x i − x j | denotes the difference betweenpairs of points. The correlation length parameters { (cid:96) } are found through hyperparameter optimization in thescikit-learn package [118]. A. Design points Performance of the GPE depends critically on thechoice of design points. We base our initial choice of thedesign points on the method of a space-filling Latin Hy-percube Design (LHD) [119, 120], which ensures marginaluniformity and optimizes the distance between points.This was implemented with the function optimumLHS inthe R package lhs . In order to reduce the emulator in-terpolation uncertainty in the most relevant regions ofphase space, the choice is then revisited and improvedby adding more design points. For instance, for the Lbt model calibration described later, we start with 60 train-ing points uniformly sampled within A × C ∈ [0 , × [0 , , . × [0 , . , × [0 , 1] and then 20 points in [0 , . × [0 , . 75] – andrepeat the calibration procedure several times to ensurethat sufficient training points have been sampled within A C Design Points of Inputs A,C FIG. 1. Latin Hypercube Design for input parameters A and C for the Lbt model. the region where the peaks of the posterior distributionsof our model parameters locate. Note that the designpoint parameter space has four or five dimensions, whichcannot be directly visualized. To illustrate the final setof design points, the distribution of inputs A and C forthe Lbt model is shown in Fig. 1. B. Multivariate output For each design point, we run the computer model forthree collision systems and two centralities (Sect. IV) todetermine the inclusive hadron R AA at various p T val-ues. The set of R AA values at each p T provides a 66-dimensional output for each design point. Instead ofpassing the 66-dimensional output directly to a high-dimensional GPE, we first employ a Principal Compo-nent Analysis (PCA). The PCA both reduces the outputdimensionality and provides a linearly independent de-scription, making the output data more tractable andallowing the application of independent GPEs.PCA rotates the data onto an orthogonal space, utiliz-ing the Singular Value Decomposition (SVD) of the data,as follows. Let Y be the centered and scaled output with n rows and p columns; then for a diagonal matrix S andorthogonal matrices U and V , the SVD of Y is Y = USV (cid:48) , (16)where V (cid:48) denotes the transpose of V . The rotation of Y by V gives the matrix US , which has uncorrelatedcolumns. If we assume normality of the data, then thecolumns are independent as well. Thus, we apply thetransformation Z = YV (17)0 . . . . Fraction of Variance Explained Number of PC r F R FIG. 2. Percent of the output data variance as a function ofthe number of components determined in the PCA. 99% ofthe variance is described with two components, while 99.9%is described with three components. and train independent GPs on the columns of Z . Topredict a new point, we take the GP predictions androtate them by V (cid:48) . PCA can also be used for dimension reduction. Thevalues { s r } of S are the square roots of the eigenval-ues { λ r } of the scaled sample covariance matrix of Y ,in non-increasing order, with associated eigenvectors in V . The fraction of variance corresponding to the first R eigenvectors, for R ≤ p , is F R = (cid:80) Rr =1 λ r (cid:80) r λ r . (18)If we use only the first R columns of V in our transforma-tion (call this matrix V R ), we capture F R of the varianceexplained by Y . Thus the goal is to choose a value of R that is large enough to explain a suitable amount of vari-ance, and small enough for a tractable number of GPEs.Figure 2 depicts the variance corresponding to R for the Lbt data. For these data we choose R = 3, which cap-tures 99.9% of the variation. C. Emulator validation To validate our GPEs, we perform “hold out” tests inwhich we remove a design point from the emulator train-ing procedure, and compare the emulator predictions atthat design point to the corresponding model calculation.We repeat this procedure for all design points, in orderto validate the emulator performance broadly across thedesign space. Figure 3 depicts examples of these compar-isons for the LBT model, along with a comparison to the R trueAA R e m u l a t o r AA Au Au200GeV,40 50%emulator =7.8% Pb Pb2.76TeV,30 40%emulator =10.3% Pb Pb5.02TeV,30 50%emulator =8.5% LBT Au Au 200 GeV,40 50%Pb Pb 2.76 TeV,30 40%Pb Pb 5.02 TeV,30 50% ( R t r u e AA R e m u l a t o r AA ) / e m u l a t o r FIG. 3. Validation of the Gaussian process emulator predic-tions. For each design point, the emulators are re-trainedwithout constraint from that holdout point, and the emula-tor predictions are compared with the model calculations atthe design point. emulator uncertainties in the right panel. We find thatthe GPEs generally predict the model well, and that theemulator uncertainties capture the deviations reasonablywell; the emulator uncertainties in fact slightly overesti-mate the observed deviations. There exist a small num-ber of design points which are poorly predicted, as shownby the off-diagonal scatter points. These points originateat the boundaries of the parameter space, where interpo-lation is not possible; we verified that they do not impactour results. We validated the GPE performance for allmodels, with the average emulator uncertainties in therange (cid:104) σ emulator (cid:105) ≈ − 22% ( Lbt ), (cid:104) σ emulator (cid:105) ≈ − Matter ), (cid:104) σ emulator (cid:105) ≈ − 14% ( Matter + Lbt (cid:104) σ emulator (cid:105) ≈ − 24% ( Matter + Lbt VI. BAYESIAN CALIBRATION With a validated emulator, we can proceed to calibra-tion - using the experimental data to perform inferenceon the input parameters. We use Bayesian inference,treating model parameters as random variables charac-terized by probability distributions, and use Bayes’ Ruleto update the prior distribution of input parameters θ (e.g. for LBT, these are θ = { A, B, C, D } ) to the poste-rior distribution θ conditional on the experimental values Y E [121, 122]. Let m ( θ ) denote the computationally ex-pensive computer model; then, the posterior distribution f ( θ | Y E ) is f ( θ | Y E ) ∝ f ( Y E | m ( θ )) f ( θ ) . (19)Because the posterior f ( θ | Y E ) is not analyticallytractable, we employ an affine invariant Markov ChainMonte Carlo algorithm [123] to draw samples from f ( θ | Y E ).1Recall that we train R independent GPEs on the first R columns of Z = YV . Let m ∗ r ( θ ) be the GPE interpo-lation for given inputs θ . Then m ∗ r ( θ ) has Normal distri-bution with mean µ ∗ r ( θ ) and variance σ ∗ r ( θ ).Note that each predictive mean µ ∗ r ( θ ) and variance σ ∗ k ( θ ) are implicitly conditioned on column r of Z (i.e.transformed design output) and design input. Becausethe GPEs are independent, we can easily write down thejoint distribution of m ∗ ( θ ) = [ m ∗ ( θ ) , . . . , m ∗ R ( θ )] (cid:48) : m ∗ ( θ ) ∼ N ( µ ∗ ( θ ) , Σ ∗ ( θ )) µ ∗ ( θ ) = [ µ ∗ ( θ ) , . . . , µ ∗ R ( θ )] (cid:48) Σ ∗ ( θ ) = diag (cid:18)(cid:104) σ ∗ ( θ ) , . . . , σ ∗ R ( θ ) (cid:105) (cid:48) (cid:19) (20)Since we emulate in PCA space, we must rotate ourpredictive interpolations back into the observable space,i.e. multiply m ( θ ) by V (cid:48) R . However, even though wecapture over 99% of the variance with our choice of R ,we have found calibration to be more stable if we addback the extra variation lost when transforming back tothe physical space. From the SVD decomposition Y = USV (cid:48) , we see that Y (cid:48) Y = VS V (cid:48) . Additionally, if we let V b denote the matrix comprised of the columns of V from R +1 onward (and similarly to S ) then we can decomposethe VS V (cid:48) into the sum VS V (cid:48) = V R S R V R (cid:48) + V b S b V b (cid:48) . (21)Noting that the sample covariance matrix is n Y (cid:48) Y , wedenote Σ extra = n V b S b V b (cid:48) the covariance matrix of ex-tra variation lost when transforming back and forth fromthe PCA space.Initially, we model Y E (centered and scaled to match Y ) as multivariate Normal, centered at m ∗ ( θ ) V (cid:48) with co-variance matrix Σ E + Σ extra . However, because m ∗ ( θ ) V (cid:48) is also multivariate Normal, we can analytically integrateover m ∗ ( θ ). Our final calibration model is thus f ( Y E | θ ) ∼ N ( µ ∗ ( θ ) V (cid:48) R , V R Σ ∗ ( θ ) V (cid:48) R + Σ E + Σ extra ) f ( θ ) ∼ unif ( θ ) , (22)where we assign a uniform prior on the design space for θ . To sample the posterior distribution, we discard thefirst 30,000 samples (“burn-in”) of the Markov ChainMonte Carlo algorithm, for which the sampler has notyet reached equilibration, and then save the next 100,000as draws from the posterior distribution f ( θ | Y E ). VII. CLOSURE TESTS In order to validate the end-to-end analysis procedure,we perform a set of closure tests. We “hold out” a design / T q TruthExtracted 90% CR p = 100 GeV/cLBT (a) / T q Prior / T q TruthExtracted 90% CR T = 300 MeVLBT (b) / T q Prior FIG. 4. Example closure test for a single design point. Thetruth (solid line) is compared to the inferred 90% range ofˆ q (colored band), as a function of (a) temperature and (b) jetmomentum. point from the emulator training, as described in SectionV C, and instead use the model predictions at that de-sign point to generate “pseudo-data” equivalent to theexperimentally measured datasets. We then perform theBayesian calibration procedure using this pseudo-data inplace of the experimental measurements, and comparethe inferred parameters to the original parameters of thedesign point.Figure 4 shows an example of such a closure test for asingle design point, in which the inferred credible regionfor ˆ q is compared to the true value from the design point.We repeat these closure tests for each design point, andstatistically evaluate their consistency using a p -value.The p -value is defined as the percentage of posterior sam-ples that are more compatible with the pseudo-data thanthe truth, using a χ taking into account the correlationin the uncertainties. We generally find consistent per-formance, as can be seen as an example in Fig. 5. Thedistribution deviates from a flat distribution with a shifttoward high p -values, indicating that the uncertainty ob-2 N u m be r o f ho l dou t t e s t s LBT FIG. 5. Distribution of p -values from the closure tests per-formed using all the design points. tained is conservative.Furthermore, we examine the closure differentially inˆ q and θ . As described in Sect. V C, we observe consis-tency except for occasional failure at the boundary of theparameter space, which are found not to be near the ex-tracted solutions. The boundary points consist of thosewith p -value very close to 0. One caveat is that for the Matter + Lbt q/T (cid:38) 5. Thisbecomes relevant in the low momentum region, and ac-cordingly the results in that region should be interpretedcautiously. VIII. RESULTS In this section we discuss the posterior distribu-tion from the Bayesian parameter extraction for theparametrizations in Sec. III, and the corresponding val-ues of ˆ q . We first discuss the analysis using the Mat-ter and Lbt models separately, and then the analysis ofthe combined model of Matter + Lbt with two differentchoices of ˆ q parametrization. A. Parameter extraction using MATTER and LBTseparately We first carry out Bayesian parameter extraction for Matter and Lbt , using Eq. (12). Figure 6 shows thedistribution of inclusive hadron R AA for the three mea-sured datasets, compared to calculations based on Lbt at the initial design points prior to parameter extraction.The prior distribution of the parameter space covers allexperimental data and serves as the training data for the Parameter A B C D Q Matter Lbt Matter + Lbt Matter + Lbt GPE. The analogous distributions for the other calcula-tions discussed below look similar to Fig. 6 and will notbe shown.Figures 7 and 8 show the same data and the posteriordistributions for Lbt and Matter . The dashed linesindicate median values, corresponding to the median pa-rameters values given in Tab. I. The models describe thedata moderately well compared to the experimental un-certainties, but exhibit systematic deviations at high p T in central Pb + Pb collisions at √ s NN =5.02 TeV, and atall p T for semi-central collisions in both √ s NN =2.76 and5.02 TeV.We compare Matter results with the LHC data(Fig. 8) only above p T ∼ 30 GeV /c because, when Mat-ter is applied alone, after the first few splittings theparton virtuality may drop below Q while its energy isstill high. Such partons, when modeled by Matter , donot interact further with the medium and therefore yieldan inaccurate description of the p T dependence of R AA . Matter is therefore expected to work well at high p T and to fail at low p T . However, experimental data uncer-tainties are smaller at low p T than high p T , and the low p T data therefore contribute with higher weights to thecalibration. This generates the deviations between Mat-ter and data at high p T for the LHC data in (Fig. 8).The implementation in the Bayesian inference analysis ofa confidence measure for a model calculation in differentregions of phase space will be explored in a future studyto address this issue.Figure 9 shows the posterior distribution of the param-eter space from this procedure, with the median value ofeach parameter given in Tab. I. The diagonal panels show1-D projections onto each parameter; a clear differencecan be seen between the Matter and the Lbt models.The off-diagonal panels show 2-D projections for Mat-ter (upper right) and Lbt (lower left).For Matter , the extracted value of A is significantwhile that of C peaks close to zero, indicating that theextracted value of ˆ q is due primarily to the first term inthe braces in Eq. (12). In contrast, Lbt results in a moresignificant contribution from the second term in Eq. (12).This is consistent with the respective domains of applica-bility of the two models: parton splitting inside Matter is driven by high virtuality and is insensitive the thermalscale of the medium, while Lbt describes the scatter-ings between jet partons with a thermal medium with3 10 20 30 40 50 60 70 80 9000.20.40.60.811.21.41.6 10 20 30 40 50 60 70 80 90 10000.20.40.60.811.21.41.6 AA R (GeV/c) T p (GeV/c) T p (GeV/c) T p FIG. 6. (Color online) Inclusive hadron R AA for the three measured datasets [21, 23, 24], together with prior calculations basedon Lbt using design points of the parameter space. Inner error bars on experimental data points are statistical errors; outererror bars are the quadrature sum of statistical error and systematic uncertainty. 10 20 30 40 50 60 70 80 9000.20.40.60.811.21.41.6 10 20 30 40 50 60 70 80 90 10000.20.40.60.811.21.41.6 AA R (GeV/c) T p (GeV/c) T p (GeV/c) T p FIG. 7. (Color online) Posterior predictive distributions of inclusive hadron R AA using Lbt compared to the same data asFig. 6. an on-shell approximation. Note that the domain of ex-perimental data was restricted according to the expectedregime of validity of each model.Figure 9 also shows the correlation between pairs ofparameters in the off-diagonal panels. A marked anti-correlation between parameters A and C is observed inthe first column of the third row, because both A and C contribute positively to the overall normalization of ˆ q inthis parametrization. On the other hand, a weaker cor-relation is seen between B and D , which is also expectedin this parametrization.Figure 10 shows the 90% credible region (C.R.) forˆ q , determined from the posterior distributions in Fig. 9.The dotted and solid lines show the median values forfixed quark momentum and medium temperature inthe upper and lower panels, respectively, illustratingmore differential information than in Tab. I. This new constraint on ˆ q is consistent within uncertainties withthe value determined previously by the JET Collabora-tion [66], although the median value is smaller. Thisis expected, since the semi-analytical calculations usedin the JET Collaboration analysis did not include elas-tic scattering processes, and some calculations consideredonly single gluon emission for the inelastic process. Theinclusion of multiple gluon emission channels and elasticscattering in the both Matter and Lbt reduces the ex-tracted ˆ q value relative to these simpler approximations.The extracted value of ˆ q/T has only weak T -dependence for both the Matter - and Lbt -based analy-ses. Figure 10(b) shows a slight decrease in ˆ q/T at highjet p T for both Matter and Lbt . The uncertainty islarger at low p T due to the p T -range of R AA data con-sidered in this work.4 10 20 30 40 50 60 70 80 9000.20.40.60.811.21.41.6 10 20 30 40 50 60 70 80 90 10000.20.40.60.811.21.41.6 AA R (GeV/c) T p (GeV/c) T p (GeV/c) T p FIG. 8. (Color online) Posterior predictive distributions of inclusive hadron R AA using Matter compared to the same dataas Fig. 6. Data points at lower p T values are excluded from this comparison due to the applicability of the model. AA BB CC DD LBTMATTER FIG. 9. (Color online) Posterior distribution of the 4-D spacefor ˆ q when Matter and Lbt are applied separately. Off-diagonal panels show correlations of posterior distributionsfor Lbt (lower left, red) and Matter (upper right, blue). Tocompare the distributions for the two models, parameters Aand C are on similar scales, as are B and D. B. Parameter extraction using MATTER+LBTcombined: five parameter Figure 11 shows inclusive hadron R AA with poste-rior parameter distributions for the combined Mat-ter + Lbt q in Eq. 12 and the switching virtuality Q ,for a total of five parameters. The bands show the pos- / T q MATTER 90% CRLBT 90% CRJET Collaboration p = 100 GeV/c(a) / T q Prior MATTER/LBT / T q MATTER 90% CRLBT 90% CR T = 300 MeV(b) / T q Prior MATTER/LBT FIG. 10. (Color online) The (quark) jet transport coefficient ˆ q from Bayesian parameter extraction using Matter and Lbt separately: (a) as function of the medium temperature, and(b) as function of quark momentum. The solid and dashedlines indicate the median value for Matter and Lbt , respec-tively. Matter or Lbt alone,there is no significant improvement when fitting with thecombined model. The level of agreement of the posteriordistributions is similar to that seen in Figs. 7 and 8, in-dicating that the simple model of a virtuality scale Q for switching between Matter and Lbt may not fullycapture the virtuality dependence of jet quenching.Figure 12 shows the correlation of posterior parameterdistributions for the Matter + Lbt Q ,the virtuality scale at which the calculation switches be-tween Matter and Lbt . The median value is 2.02 GeV,with 90% credible region [1 . , . 72] GeV. It is evidentthat the RHIC data have significantly less impact thanthe LHC data on the posterior distributions in this anal-ysis, because of the low relative statistical weight of theselected RHIC data (Sect. IV). Line 3 of Table I givesthe median values of the five parameters.The value of Q in this analysis is taken to be the samefor RHIC and LHC data, though in practice the value of Q may be smaller at RHIC than at the LHC because ofthe lower average QGP temperature. Exploration of thisdegree of freedom requires consideration of additional ex-perimental data, however, and is beyond the scope of thecurrent analysis.Figure 13 shows the 90% C.R. of the extracted quarkjet ˆ q from the Matter + Lbt Matter and Lbt takenfrom Fig. 10. The value of ˆ q determined using the multi-stage approach is lower than those determined from fits to Matter or Lbt separately. This is because Matter iseffective for jet energy loss at high virtuality but is blindto parton evolution at low virtuality, while the oppositeis the case for Lbt . Combining the models for a multi-stage evolution approach leads to larger jet energy lossthan found by applying only one of them across the entirephase space. Matter + Lbt therefore requires a smallerˆ q value than Matter or Lbt does to describe the samejet quenching data. C. Parameter extraction using MATTER+LBTcombined: four parameter Finally, we discuss results using the Matter + Lbt R AA distribution and the model calculationutilizing the median values of the parameters. Quali-tatively, the model calculations describe the overall p T -dependence of the data well although, as also seen inFigs. 7, 8 and 11, the posterior distributions fall outsideof the systematic uncertainty limits of the data. Thisagain indicates that introduction of a virtuality scale Q for switching between Matter and Lbt may not be suffi-cient to describe the virtuality dependence of jet quench-ing.Figure 15 shows the correlation of posterior param-eter distributions of Matter + Lbt Mat-ter + Lbt Q from com-bined RHIC and LHC data is 2.70 GeV, with 90% CR[1 . , . 41] GeV. While the median value is larger thanthat determined using Matter + Lbt 1, the differenceis not significant, as seen from the mutually compati-ble 90% confidence regions. This comparison providesan estimate of the systematic uncertainty due to differ-ent model parametrization. Line 4 of Table I gives themedian values of the four parameters.Figure 16 shows ˆ q as a function of medium temperatureand quark momentum using Matter + Lbt 2. Referringto Eq. (13), ˆ q has a larger value in the Matter stagethan in the Lbt stage, since both terms in { . . . } con-tribute to Matter but only the second term contributesto Lbt . The contribution from the first term dependson the virtuality of the parton Q , which varies for dif-ferent partons with the same momentum. In order toestimate the ˆ q value in Fig. 16 we use the average valueof Q , obtained from Eq. (1). Comparing with Fig. 13,we also observe that the value of ˆ q extracted from Mat-ter + Lbt Matter + Lbt IX. SUMMARY AND OUTLOOK We have reported the application of state-of-the-artBayesian inference methodology to determine the QGPjet transport coefficient ˆ q from inclusive hadron suppres-sion data measured at RHIC and the LHC. Two jet en-ergy loss models were utilized, Matter and Lbt . Mat-ter is applicable to modeling the medium-modified split-ting of highly virtual partons, while Lbt is applicablefor the in-medium transport of nearly on-shell partons.The models are first applied separately, and then com-bined to form a multi-stage evolution approach. Twodifferent parametrizations were used for the functionaldependence of ˆ q on jet momentum or virtuality scale andthe medium temperature, based on the picture of per-turbative scattering between jets and a thermal medium.A novel treatment of experimental uncertainties is em-ployed, taking into account their covariance for the firsttime in the determination of ˆ q .Such model calculations are computationally expen-sive. Gaussian process emulators are therefore employedto render this process computationally efficient, trainedat design points selected starting with a Latin Hypercubein parameter space. The resulting gain in computationalefficiency enabled calibration of the multi-dimensionalparameter space.The Bayesian inference process generates posterior pa-6 10 20 30 40 50 60 70 80 9000.20.40.60.811.21.41.6 10 20 30 40 50 60 70 80 90 10000.20.40.60.811.21.41.6 AA R (GeV/c) T p (GeV/c) T p (GeV/c) T p FIG. 11. (Color online) Posterior predictive distributions of Matter + Lbt AA CC BB DD QQ RHIC+LHCLHCRHIC MATTER+LBT1 FIG. 12. (Color online) Posterior distribution of the 4-D spacefor for Matter + Lbt 1. Off-diagonal panels show correlationsof posterior distributions for RHIC+LHC (lower left, red) andLHC only (upper right, blue). rameter distributions for each model configuration. Toconstrain the model parameters, we used 66 inclusivehadron R AA datapoints at two centralities for Au + Aucollisions at √ s NN =200 GeV and Pb + Pb collisions at √ s NN =2.76 and 5.02 TeV. For both the Matter -onlyand Lbt -only configurations, the extracted value of ˆ q/T has only weak dependence on the medium temperature T . The value of ˆ q determined using these approachesis consistent with a previous determination by the JETCollaboration [66].A multi-stage jet evolution approach, combining Mat- / T q MATTER 90% CRLBT 90% CRMATTER+LBT1 90% CRJET Collaboration p = 100 GeV/c(a) / T q Prior MATTER+LBT1 / T q MATTER 90% CRLBT 90% CRMATTER+LBT1 90% CR T = 300 MeV(b) / T q Prior MATTER+LBT1 FIG. 13. (Color online) 90% CR regions for the quark-jettransport coefficient ˆ q using Matter and Lbt (Fig. 10) and Matter + Lbt 10 20 30 40 50 60 70 80 9000.20.40.60.811.21.41.6 10 20 30 40 50 60 70 80 90 10000.20.40.60.811.21.41.6 AA R (GeV/c) T p (GeV/c) T p (GeV/c) T p FIG. 14. (Color online) Posterior predictive distributions of R AA using Matter + Lbt 2, compared to the same data as Fig. 6.Dashed lines show model calculation using median values of parameters. AA CC DD QQ RHIC+LHCLHCRHIC MATTER+LBT2 FIG. 15. (Color online) Posterior distribution of the 4-D spacefor for Matter + Lbt 2. Off-diagonal panels show correlationsof posterior distributions for RHIC+LHC (lower left, red) andLHC only (upper right, blue). ter and Lbt , is applied here for the first time. Thetransition between Matter and Lbt , based on par-ton virtuality, is controlled by the virtuality parame-ter Q which separates the virtuality-ordered-splittingdominating region and the time-ordered-transport dom-inating region for jet quenching inside a medium. Theposterior distribution of ˆ q/T from the combined model( Matter + Lbt ) is systematically lower than that deter-mined using Matter or Lbt alone, since the combinedmodel more accurately describes energy loss over the full / T q MATTER in MATTER+LBT 2LBT in MATTER+LBT 2JET Collaboration p = 100 GeV/c(a) / T q Prior MATTER+LBT2 / T q MATTER in MATTER+LBT 2LBT in MATTER+LBT 2 T = 300 MeV(b) / T q Prior MATTER+LBT2 FIG. 16. (Color online) 90% CR regions for the quark-jettransport coefficient ˆ q using Matter + Lbt p T -dependence. The twodifferent ˆ q parametrizations give consistent results, al-though with differences in the median extracted param-eter values; the median value of Q is 2.0 and 2.7 GeV,respectively.This analysis represents a significant advance in quan-titative understanding of jet-medium interactions in theQuark-Gluon Plasma. However, the posterior distri-butions from this analysis do not fully describe themagnitude and p T -dependence of inclusive hadron R AA measurements in the datasets considered. This tensionindicates that additional components in the modelingof jet quenching are needed, for instance a more de-tailed parametrization of the virtuality dependence of jetquenching than the single switching scale Q used here.Future work will also provide more detailed account-ing of experimental and theoretical uncertainties andtheir covariance, and incorporate additional measure-ment channels, in particular those involving recon-structed jets. While the ˆ q parametrization employed inthis analysis is derived from the perturbative approachto jet-medium scattering, non-perturbative effects mayrequire additional dependence of ˆ q on jet energy andmedium temperature, and additional transport param-eters may be required to fully describe the jet measure-ments. Other models of the plasma, incorporating quasi-particle degrees of freedom [124, 125], will also be con-sidered. ACKNOWLEDGMENTS This work was supported in part by the National Sci-ence Foundation (NSF) within the framework of theJETSCAPE collaboration, under grant numbers ACI-1550172 (Y.C. and G.R.), ACI-1550221 (R.J.F., F.G.,M.K. and B.K.), ACI-1550223 (D.E., M.M., U.H., andL.D.), ACI-1550225 (S.A.B., J.C., T.D., W.F., R.W.,S.M., and Y.X.), ACI-1550228 (J.M., B.J., P.J., W.K.,X.-N.W.), and ACI-1550300 (S.C., L.C., A.K., A.M.,C.N., C.P., A.S., J.P., L.S., C.Si., R.A.S. and G.V.).It was also supported in part by the NSF under grantnumbers OAC–2039142 (R.A.), PHY-1516590, PHY-1812431 and PHY-2012922 (R.J.F., B.K., F.G., M.K., and C.S.), and by the US Department of Energy, Officeof Science, Office of Nuclear Physics under grant num-bers DE-AC02-05CH11231 (D.O., X.-N.W.), DE-AC52-07NA27344 (A.A., R.A.S.), DE-SC0013460 (S.C., A.K.,A.M., C.S. and C.Si.), DE-SC0004286 (L.D., M.M., D.E.and U.H.), DE-SC0012704 (B.S. and C.S.), DE-FG02-92ER40713 (J.P.) and DE-FG02-05ER41367 (T.D., J.-F.P., S.A.B. and Y.X.). The work was also supportedin part by the National Science Foundation of China(NSFC) under grant numbers 11935007, 11861131009and 11890714 (Y.H. and X.-N.W.), by the NaturalSciences and Engineering Research Council of Canada(C.G., M.H., S.J., C.P. and G.V.), by the Fonds derecherche du Qu´ebec – Nature et technologies (FRQNT)(G.V.), by the Office of the Vice President for Research(OVPR) at Wayne State University (C.P. and Y.T.), bythe S˜ao Paulo Research Foundation (FAPESP) underprojects 2016/24029-6, 2017/05685-2 and 2018/24720-6(M.L.), and by the University of California, Berkeley -Central China Normal University Collaboration Grant(W.K.). U.H. would like to acknowledge support by theAlexander von Humboldt Foundation through a Hum-boldt Research Award. Allocation of supercomputingresources (Project: PHY180035) were obtained in partthrough the Extreme Science and Engineering Discov-ery Environment (XSEDE), which is supported by Na-tional Science Foundation grant number ACI-1548562.Calculation were performed in part on Stampede2 com-pute nodes, generously funded by the National ScienceFoundation (NSF) through award ACI-1134872, withinthe Texas Advanced Computing Center (TACC) at theUniversity of Texas at Austin [126], and in part on theOhio Supercomputer [127] (Project PAS0254). Compu-tations were also carried out on the Wayne State Gridfunded by the Wayne State OVPR, and on the super-computer Guillimin from McGill University, managed byCalcul Qu´ebec and Compute Canada. The operation ofthe supercomputer Guillimin is funded by the CanadaFoundation for Innovation (CFI), NanoQu´ebec, R´eseaude M´edicine G´en´etique Appliqu´ee (RMGA) and FRQ-NT. Data storage was provided in part by the OSIRISproject supported by the National Science Foundationunder grant number OAC-1541335. [1] E. Shuryak, Rev. Mod. Phys. , 035001 (2017),arXiv:1412.8393 [hep-ph].[2] W. Busza, K. Rajagopal, and W. van der Schee, Ann.Rev. Nucl. Part. Sci. , 339 (2018), arXiv:1802.04801[hep-ph].[3] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. , 123 (2013), arXiv:1301.2826 [nucl-th].[4] G.-Y. Qin and X.-N. Wang, Int. J. Mod. Phys. E24 ,1530014 (2015), arXiv:1511.00790 [hep-ph].[5] G. Arnison et al. (UA1 Collaboration), Phys. Lett. B , 461 (1986). [6] J. Appel et al. (UA2 Collaboration), Phys. Lett. B ,349 (1985).[7] T. Aaltonen et al. (CDF Collaboration), Phys. Rev. D , 052006 (2008).[8] V. M. Abazov et al. (D0 Collaboration), Phys. Rev.Lett. , 062001 (2008).[9] B. Abelev et al. (STAR Collaboration), Phys. Rev. Lett. , 252001 (2006).[10] G. Aad et al. (ATLAS Collaboration), J. High EnergyPhys. , 153 (2015).[11] Phys. Rev. C , 034911 (2020). [12] A. M. Sirunyan et al. (CMS), JHEP , 082 (2020),arXiv:2005.05159 [hep-ex].[13] M. Dasgupta, F. A. Dreyer, G. P. Salam, and G. Soyez,J. High Energy Phys. , 57 (2016).[14] J. Currie, E. W. N. Glover, and J. Pires, Phys. Rev.Lett. , 072002 (2017).[15] X. Liu, S.-O. Moch, and F. Ringer, Phys. Rev. D ,056026 (2018).[16] M. Czakon, A. van Hameren, A. Mitov, and R. Pon-celet, J. High Energy Phys. , 262 (2019).[17] A. J. Larkoski, I. Moult, and B. Nachman, PhysicsReports , 1 (2020).[18] X.-N. Wang and M. Gyulassy, Phys. Rev. Lett. , 1480(1992).[19] M. L. Miller, K. Reygers, S. J. Sanders, and P. Stein-berg, Ann. Rev. Nucl. Part. Sci. , 205 (2007),arXiv:nucl-ex/0701025.[20] J. Adams et al. (STAR), Phys. Rev. Lett. , 172302(2003), arXiv:nucl-ex/0305015 [nucl-ex].[21] A. Adare et al. (PHENIX), Phys. Rev. C87 , 034911(2013), arXiv:1208.2254 [nucl-ex].[22] S. Chatrchyan et al. (CMS), Eur. Phys. J. C72 , 1945(2012), arXiv:1202.2554 [nucl-ex].[23] G. Aad et al. (ATLAS), JHEP , 050 (2015),arXiv:1504.04337 [hep-ex].[24] V. Khachatryan et al. (CMS), JHEP , 039 (2017),arXiv:1611.01664 [nucl-ex].[25] S. Acharya et al. (ALICE), JHEP , 013 (2018),arXiv:1802.09145 [nucl-ex].[26] E. Braaten and M. H. Thoma, Phys. Rev. D44 , 2625(1991).[27] M. Gyulassy and X.-N. Wang, Nucl. Phys. B420 , 583(1994), nucl-th/9306003.[28] R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne,and D. Schiff, Nucl. Phys. B483 , 291 (1997), hep-ph/9607355.[29] B. G. Zakharov, JETP Lett. , 952 (1996), arXiv:hep-ph/9607440.[30] G.-Y. Qin et al. , Phys. Rev. Lett. , 072301 (2008),arXiv:0710.0605 [hep-ph].[31] N. Armesto et al. , Phys. Rev. C , 064904 (2012),arXiv:1106.1106 [hep-ph].[32] S. A. Bass et al. , Phys. Rev. C79 , 024901 (2009),arXiv:0808.0908 [nucl-th].[33] N. Armesto, M. Cacciari, T. Hirano, J. L. Nagle,and C. A. Salgado, J. Phys. G37 , 025104 (2010),arXiv:0907.0667 [hep-ph].[34] X.-F. Chen, C. Greiner, E. Wang, X.-N. Wang, andZ. Xu, Phys. Rev. C81 , 064908 (2010), arXiv:1002.1165[nucl-th].[35] S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys.Lett. B777 , 255 (2018), arXiv:1703.00822 [nucl-th].[36] A. Majumder, E. Wang, and X.-N. Wang, Phys. Rev.Lett. , 152301 (2007), arXiv:nucl-th/0412061.[37] H. Zhang, J. F. Owens, E. Wang, and X.-N.Wang, Phys. Rev. Lett. , 212301 (2007), arXiv:nucl-th/0701045.[38] T. Renk, Phys. Rev. C78 , 034904 (2008),arXiv:0803.0218 [hep-ph].[39] H. Zhang, J. F. Owens, E. Wang, and X.-N. Wang,Phys. Rev. Lett. , 032302 (2009), arXiv:0902.4000[nucl-th].[40] G.-Y. Qin, J. Ruppert, C. Gale, S. Jeon, and G. D.Moore, Phys. Rev. C80 , 054909 (2009), arXiv:0906.3280 [hep-ph].[41] X.-N. Wang and Y. Zhu, Phys. Rev. Lett. , 062301(2013), arXiv:1302.5874 [hep-ph].[42] W. Chen, S. Cao, T. Luo, L.-G. Pang, and X.-N. Wang,Phys. Lett. B777 , 86 (2018), arXiv:1704.03648 [nucl-th].[43] T. Luo, S. Cao, Y. He, and X.-N. Wang, Phys. Lett. B782 , 707 (2018), arXiv:1803.06785 [hep-ph].[44] G.-Y. Qin and B. Muller, Phys. Rev. Lett. , 162302(2011), [Erratum: Phys. Rev. Lett.108,189904(2012)],arXiv:1012.5280 [hep-ph].[45] W. Dai, I. Vitev, and B.-W. Zhang, Phys. Rev. Lett. , 142001 (2013), arXiv:1207.5177 [hep-ph].[46] N.-B. Chang and G.-Y. Qin, Phys. Rev. C , 024902(2016), arXiv:1603.01920 [hep-ph].[47] Y. He, S. Cao, W. Chen, T. Luo, L.-G. Pang,and X.-N. Wang, Phys. Rev. C , 054911 (2019),arXiv:1809.02525 [nucl-th].[48] J.-W. Qiu, F. Ringer, N. Sato, and P. Zurita, Phys.Rev. Lett. , 252301 (2019), arXiv:1903.01993 [hep-ph].[49] Z.-B. Kang, F. Ringer, and I. Vitev, Phys. Lett. B (2017), 10.1016/0370-2693(86)90290-X.[50] Y.-T. Chien and I. Vitev, Phys. Rev. Lett. , 112301(2017), arXiv:1608.07283 [hep-ph].[51] Y. Mehtar-Tani and K. Tywoniuk, JHEP , 125(2017), arXiv:1610.08930 [hep-ph].[52] N.-B. Chang, S. Cao, and G.-Y. Qin, Phys. Lett. B781 ,423 (2018), arXiv:1707.03767 [hep-ph].[53] J. Casalderrey-Solana, Y. Mehtar-Tani, C. A. Salgado,and K. Tywoniuk, Phys. Lett. B , 357 (2013).[54] P. Caucal, E. Iancu, and G. Soyez, JHEP , 273(2019).[55] S. Cao et al. (JETSCAPE), Phys. Rev. C96 , 024909(2017), arXiv:1705.00050 [nucl-th].[56] R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne,and D. Schiff, Nucl. Phys. B484 , 265 (1997), hep-ph/9608322.[57] R. Baier, Y. L. Dokshitzer, A. H. Mueller, and D. Schiff,Nucl. Phys. B531 , 403 (1998), arXiv:hep-ph/9804212[hep-ph].[58] M. Gyulassy, P. Levai, and I. Vitev, Nucl. Phys. B571 ,197 (2000), arXiv:hep-ph/9907461.[59] U. A. Wiedemann, Nucl. Phys. B588 , 303 (2000),arXiv:hep-ph/0005129.[60] X.-F. Guo and X.-N. Wang, Phys. Rev. Lett. , 3591(2000), arXiv:hep-ph/0005044.[61] P. Arnold, G. D. Moore, and L. G. Yaffe, JHEP ,030 (2002), hep-ph/0204343.[62] N. Armesto, C. A. Salgado, and U. A. Wiedemann,Phys. Rev. D69 , 114003 (2004), arXiv:hep-ph/0312106[hep-ph].[63] M. Djordjevic and U. W. Heinz, Phys.Rev.Lett. ,022302 (2008), arXiv:0802.1230 [nucl-th].[64] A. Majumder, Phys. Rev. D85 , 014023 (2012),arXiv:0912.2987 [nucl-th].[65] S. Caron-Huot and C. Gale, Phys. Rev. C , 064902(2010), arXiv:1006.2379 [hep-ph].[66] K. M. Burke et al. (JET), Phys. Rev. C90 , 014909(2014), arXiv:1312.5003 [nucl-th].[67] A. Majumder, Phys. Rev. C , 034905 (2013),arXiv:1202.5295 [nucl-th].[68] A. Kumar, A. Majumder, and C. Nonaka, PoS Hard-Probes2018 , 051 (2019), arXiv:1901.03878 [nucl-th]. [69] A. Kumar, A. Majumder, and J. H. Weber, (2020),arXiv:2010.14463 [hep-lat].[70] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn,Statistical Science , 409 (1989).[71] C. Currin, T. Mitchell, M. Morris, and D. Ylvisaker,Journal of the American Statistical Association , 953(1991).[72] M. C. Kennedy and A. O’Hagan, Journal of the RoyalStatistical Society: Series B (Statistical Methodology) , 425 (2001).[73] J. Novak, K. Novak, S. Pratt, J. Vredevoogd,C. Coleman-Smith, and R. Wolpert, Phys. Rev. C89 ,034917 (2014), arXiv:1303.5769 [nucl-th].[74] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu,and U. Heinz, Phys. Rev. C94 , 024907 (2016),arXiv:1605.03954 [nucl-th].[75] D. Everett et al. (JETSCAPE), (2020),arXiv:2010.03928 [hep-ph].[76] D. Everett et al. (JETSCAPE), (2020),arXiv:2011.01430 [hep-ph].[77] S. Pratt, E. Sangaline, P. Sorensen, and H. Wang, Phys.Rev. Lett. , 202301 (2015), arXiv:1501.04042 [nucl-th].[78] G. Nijs, W. Van Der Schee, U. G¨ursoy, and R. Snellings,(2020), arXiv:2010.15134 [nucl-th].[79] Y. Xu, J. E. Bernhard, S. A. Bass, M. Nahrgang,and S. Cao, Phys. Rev. C97 , 014907 (2018),arXiv:1710.00807 [nucl-th].[80] A. Majumder, Phys. Rev. C88 , 014909 (2013),arXiv:1301.5323 [nucl-th].[81] S. Cao and A. Majumder, Phys. Rev. C , 024903(2020), arXiv:1712.10055 [nucl-th].[82] H. Song and U. W. Heinz, Phys. Lett. B658 , 279 (2008),arXiv:0709.0742 [nucl-th].[83] H. Song and U. W. Heinz, Phys. Rev. C77 , 064901(2008), arXiv:0712.3715 [nucl-th].[84] H. Song, S. A. Bass, U. Heinz, T. Hirano, and C. Shen,Phys. Rev. Lett. , 192301 (2011), arXiv:1011.2783[nucl-th].[85] C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, et al. ,(2014), arXiv:1409.8164 [nucl-th].[86] U. G¨ursoy, D. Kharzeev, E. Marcus, K. Rajagopal,and C. Shen, Phys. Rev. C , 055201 (2018),arXiv:1806.05288 [hep-ph].[87] https://github.com/chunshen1987/superMC .[88] C. Shen, U. Heinz, P. Huovinen, and H. Song, Phys.Rev. C82 , 054904 (2010), arXiv:1010.1856 [nucl-th].[89] P. Huovinen and P. Petreczky, Nucl. Phys. A837 , 26(2010), arXiv:0912.2541 [hep-ph].[90] Z. Qiu, C. Shen, and U. Heinz, Phys. Lett. B707 , 151(2012), arXiv:1110.3033 [nucl-th].[91] H. L. Lai et al. (CTEQ), Eur. Phys. J. C12 , 375 (2000),arXiv:hep-ph/9903282.[92] K. J. Eskola, H. Paukkunen, and C. A. Salgado, JHEP , 065 (2009), arXiv:0902.4154 [hep-ph].[93] M. Kordell and A. Majumder, in (2017)arXiv:1702.05862 [nucl-th].[94] A. Majumder, (2009), arXiv:0901.4516 [nucl-th].[95] X.-N. Wang and X.-F. Guo, Nucl. Phys. A696 , 788(2001), arXiv:hep-ph/0102230.[96] A. Majumder and C. Shen, Phys. Rev. Lett. ,202301 (2012), arXiv:1103.0809 [hep-ph]. [97] A. Majumder and J. Putschke, Phys. Rev. C93 , 054909(2016), arXiv:1408.3403 [nucl-th].[98] P. Aurenche, B. G. Zakharov, and H. Zaraket, JETPLett. , 605 (2008), arXiv:0804.4282 [hep-ph].[99] P. Aurenche, B. Zakharov, and H. Zaraket, (2008),arXiv:0806.0160 [hep-ph].[100] T. Sjostrand, S. Mrenna, and P. Z. Skands, JHEP ,026 (2006), arXiv:hep-ph/0603175 [hep-ph].[101] D. d’Enterria, K. J. Eskola, I. Helenius, andH. Paukkunen, Nucl. Phys. B , 615 (2014),arXiv:1311.1415 [hep-ph].[102] A. Kumar et al. (JETSCAPE), Phys. Rev. C ,054906 (2020), arXiv:1910.05481 [nucl-th].[103] R. Baier, A. H. Mueller, and D. Schiff, Phys. Lett. B649 , 147 (2007), arXiv:nucl-th/0612068 [nucl-th].[104] S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys. Rev. C94 , 014909 (2016), arXiv:1605.06447 [nucl-th].[105] S. Cao, G.-Y. Qin, and S. A. Bass, Phys. Rev. C88 ,044907 (2013), arXiv:1308.0617 [nucl-th].[106] S. Cao, G.-Y. Qin, and S. A. Bass, Phys. Rev. C92 ,024907 (2015), arXiv:1505.01413 [nucl-th].[107] B.-W. Zhang, E. Wang, and X.-N. Wang, Phys. Rev.Lett. , 072301 (2004), arXiv:nucl-th/0309040 [nucl-th].[108] P. Caucal, E. Iancu, A. H. Mueller, and G. Soyez, Phys.Rev. Lett. , 232001 (2018), arXiv:1801.09703 [hep-ph].[109] X.-N. Wang, Phys. Rept. , 287 (1997), arXiv:hep-ph/9605214.[110] G.-Y. Qin and A. Majumder, Phys. Rev. Lett. ,262301 (2010), arXiv:0910.3016 [hep-ph].[111] J. Auvinen, K. J. Eskola, and T. Renk, Phys. Rev. C , 024906 (2010), arXiv:0912.2265 [hep-ph].[112] Y. He, T. Luo, X.-N. Wang, and Y. Zhu, Phys. Rev. C91 , 054908 (2015), arXiv:1503.03313 [nucl-th].[113] J. Adams et al. (STAR), Phys. Rev. Lett. , 162301(2006), arXiv:nucl-ex/0604018 [nucl-ex].[114] V. Kuˇcera (ALICE), Nucl. Part. Phys. Proc. ,181 (2016), arXiv:1511.02766 [hep-ex].[115] M. Gu, X. Wang, and J. Berger, Ann. Statist. , 3038(2018).[116] C. E. Rasmussen and C. K. I. Williams, Gaussian Pro-cesses for Machine Learning (The MIT Press, 2006).[117] J. E. Bernhard, P. W. Marcy, C. E. Coleman-Smith,S. Huzurbazar, R. L. Wolpert, and S. A. Bass, Phys.Rev. C (2015), 10.1103/physrevc.91.054910.[118] F. Pedregosa et al. , Journal of Machine Learning Re-search , 2825 (2011).[119] M. D. McKay, R. J. Beckman, and W. J. Conover,Technometrics , 239 (1979).[120] R. Stocki, Computer Assisted Mechanics and Engineer-ing Sciences , 393 (2005).[121] M. C. Kennedy and A. O’Hagan, Journal of the RoyalStatistical Society: Series B (Statistical Methodology) , 425 (2001).[122] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson,A. Vehtari, and D. B. Rubin, Bayesian Data Analysis3rd Edition (Chapman and Hall/CRC, 2013).[123] D. Foreman-Mackey, D. W. Hogg, D. Lang, andJ. Goodman, Publications of the Astronomical Societyof the Pacific , 306–312 (2013).[124] S. K. Das, F. Scardina, S. Plumari, and V. Greco, Phys.Lett. B747 , 260 (2015), arXiv:1502.03757 [nucl-th]. [125] J. Xu, J. Liao, and M. Gyulassy, Chin. Phys. Lett. ,9 (2015), arXiv:1411.3673 [hep-ph].[126] . [127] “Ohio Supercomputer Center,” http://osc.edu/ark:/19495/f5s1ph73http://osc.edu/ark:/19495/f5s1ph73