A data-driven analysis for the temperature and momentum dependence of the heavy quark diffusion coefficient in relativistic heavy-ion collisions
Yingru Xu, Marlene Nahrgang, Shanshan Cao, Jonah E. Bernhard, Steffen A. Bass
AA data-driven analysis for the temperature and momentum dependence of the heavyquark diffusion coefficient in relativistic heavy-ion collisions
Yingru Xu, ∗ Jonah E. Bernhard, and Steffen A. Bass
Department of Physics, Duke University, Durham, NC 27708, USA
Marlene Nahrgang
SUBATECH, UMR 6457, IMT Atlantique, Universite de Nantes,IN2P3/CNRS, 4 rue Alfred Kastler, 44307 Nantes cedex 3, France
Shanshan Cao
Department of Physics and Astronomy, Wayne State University, Detroit, MI, 48201 (Dated: October 3, 2017)By applying a Bayesian model-to-data analysis, we estimate the temperature and momentum de-pendence of the heavy quark diffusion coefficient in an improved Langevin framework. The posteriorrange of the diffusion coefficient is obtained by performing a Markov chain Monte Carlo random walkand calibrating on the experimental data of D -meson R AA and v in three different collision systemsat RHIC and the LHC: AuAu collisions at 200 GeV, PbPb collisions at 2.76 and 5.02 TeV. Thespatial diffusion coefficient is found to be consistent with lattice QCD calculations and comparablewith other models’ estimation. We demonstrate the capability of our improved Langevin model tosimultaneously describe the R AA and v at both RHIC and the LHC energies, as well as the higherorder flow coefficient such as D -meson v . We show that by applying a Bayesian analysis, we areable to quantitatively and systematically study the heavy flavor dynamics in heavy-ion collisions. I. INTRODUCTION
The theory of the strong interaction force – QuantumChromodynamics (QCD) – predicts that at sufficientlyhigh temperature and/or baryon density, nuclear matterundergoes a phase transition from hadrons to a new stateof the deconfined quarks and gluons: the quark gluonplasma (QGP) [1–3]. Over the past two decades, ultra-relativistic heavy-ion collision experiments at the Rela-tivistic Heavy Ion Collider (RHIC) and the Large HadronCollider (LHC) have been searching and exploring thisnew state of matter under extreme conditions. Com-pelling discoveries, for instance the strong suppression ofhadrons at large transverse momenta (jet quenching), re-veal the creation of the QGP medium at RHIC and theLHC [4, 5]. The observed collective flow of low trans-verse momentum hadrons [6, 7] has provided insight intoremarkable properties of the QGP such as the strongly in-teracting, almost perfect fluid behavior with a very smallshear viscosity to entropy density ratio [8–10].Since the QGP is not directly observable, the studyof its properties relies on the measurement of final stateobservables, as well as theoretical modeling, and the com-parison between those two. For example, the relativisticviscous hydrodynamical model [11–15] – one of the mostsuccessful models in heavy-ion physics – has been uti-lized for the extraction of the temperature dependenceof the specific shear-viscosity through a model-to-datacomparison with elliptic and triangular flow data of softidentified hadrons [16]. ∗ Correspond to [email protected]
In contrast to the soft medium properties, the trans-port coefficients related to the medium interaction ofhard probes (jets and heavy quarks), such as ˆ q and ˆ e, D s ,are not yet understood on a similarly quantitative level.This is in part due to the experimental difficulty in mea-suring “rare processes”, but also due to the complex-ity of modeling the dynamics of these hard probes inter-acting with the QGP medium. Nevertheless, significantprogress has been made in recent years: a number oftransport models in the market are now able to describea selection of heavy quark observables and perform quali-tative estimates of the diffusion coefficient [17–30], whichin turn can be compared with lattice QCD calculationsof the same quantities [31–33]. In current studies of theopen heavy flavor diffusion coefficient, it is common thatthe diffusion coefficient is directly or indirectly encodedin the model and one can relate its physical propertiesto one or multiple parameters. By comparing the heavyquark observables (such as the nuclear modification fac-tor R AA and elliptic flow v ) between the theoretical cal-culation and the experimental data, these parameters canbe tuned until one finds a satisfactory fit. However, thedisadvantage of such an “eyeball” comparison is that itgets exceedingly difficult to vary multiple parameters si-multaneously or to compare with a larger selection ofexperimental measurements, as all parameters are inter-dependent and affect multiple observables at once.A more rigorous and complete approach to optimizingthe model and determining the parameters would be toperform a random walk in the parameter space and cal-ibrate on the experimental data by applying a modernBayesian statistical analysis [34–37]. In such an analysis,the computationally expensive physics model is first eval-uated for a small number of points in parameter space. a r X i v : . [ nu c l - t h ] O c t These calculations are used to train Gaussian Process em-ulators that act as fast surrogates to interpolate betweenthese points and provide model predictions for arbitraryvalues of the input parameters. The emulators thus actas substitution of the full model in order to be able toperform a Markov chain Monte Carlo exploration of thecomplete parameter space. The results of such an analy-sis are the posterior distributions of the parameters, i.e.the probability distributions of the parameter values thatoptimally describe the experimental data.This type of model-to-data analysis using Bayesianstatistics has been applied with great success in the softsector of heavy-ion physics over the last few years, e.g.for the extraction of the temperature dependence of thespecific shear and bulk viscosities of the QGP [16, 39],as well as for constraining the equation of state of QCDmatter purely from experimental measurements [36, 40].In this study, we shall extend this type of analysis to thequantitative study of heavy flavor evolution in heavy-ioncollisions. Our goal is to provide a quantitative estimateof the temperature and momentum dependence of theheavy flavor diffusion coefficient. The paper is organizedas follows: in Sec. II we will describe our physics model,the improved Langevin framework which simulates thefull space-time evolution of heavy quarks inside a QGPmedium; in Sec. III the Bayesian analysis that we utilizeis introduced; the results of the calibration and the esti-mation of the heavy quark diffusion coefficient are givenin Sec. IV; a summary and outlook can be found in thefinal section.
II. MODELING HEAVY FLAVOR EVOLUTIONIN HEAVY-ION COLLISIONSA. Full space-time evolution of heavy flavors andthe QGP medium
Our analysis utilizes the well-established frameworkdeveloped by the Duke QCD group to simulate the fullspace-time evolution of heavy quarks in heavy-ion colli-sions [22, 41]:
1. Initial production
Due to their large masses, heavy quarks are believedto be primarily produced at the beginning of the colli-sion via hard scattering. Therefore the initial momen-tum distribution can be calculated using perturbativeQCD (pQCD). In this work, we adopt the fixed-order plusnext-to-leading log formula (FONLL) [42, 43] and EPS09NLO nuclear PDFs [44] to calculate the heavy quark ini-tial momentum distribution, from which we sample themomenta of heavy quarks in our calculation (using MonteCarlo methods).The initial distribution of heavy quarks in positionspace is generated consistently with the initial condi- tion for the event-by-event hydrodynamical evolution bythe parametric initial condition model TRENTo [45]. Atthe soft medium thermalization time ( τ = 0 . s ( x, y ) | τ to the nu-cleon thickness functions T A , T B of the two projectiles byasserting a generalized ansatz: s ( x, y ) | τ = τ ∝ (cid:18) T pA + T pB (cid:19) /p . (1)where p is a free parameter in TRENTo . In this study, weutilize the calibration of the soft matter properties per-formed in [16] (Table III of the median value calibratedon charged particle yields of [16]), including TRENToinitial state parameters (except for p , which is chosen tobe 0), as well as other parameters related to soft mediumproperties, therefore the generalized function can be sim-plified as s ( x, y ) | τ = τ = (cid:112) T A T B . (2)The heavy quark initial production is based on the bi-nary collision scaling and is determined by the thicknessfunction ˆ T AB = T A T B . In this way, we are able to re-late the heavy quark initial position to the soft mediuminitial production.
2. Heavy quark evolution inside a QGP medium
After their production, heavy quarks propagatethrough the QGP medium. In the quasi-particle pictureof the QGP system, the space-time evolution of heavyquarks can generally be described by the Boltzmannequation. Since heavy quark masses are much larger thanthe typical medium temperature ( m Q (cid:29) T ), their mo-mentum change due to the scattering with light partonsin a thermally equilibrated medium is relatively small.With this assumption, the Boltzmann equation can bereduced to the Fokker-Planck equation, which is realizedstochastically by the Langevin equation (for a detailedderivation, see Refs. [46, 47]). In this study, we use animproved Langevin transport model to describe the dy-namics of heavy quarks propagating in a QGP medium,which includes not only the heavy quark collisional en-ergy loss but also the radiative energy loss due to gluonradiation [22] d(cid:126)pdt = − η D ( p ) (cid:126)p + (cid:126)ξ + (cid:126)f g . (3)The first two terms on the right hand side of the equa-tion are the drag and thermal random forces inheritedfrom the standard Langevin equation. They contributeto the collisional energy loss from quasi-elastic scatter-ing between heavy quarks and light partons, and aregeneralized to the scattering between heavy quarks andthe background medium. With the requirement that theheavy quark distribution eventually reaches equilibriumin a thermal medium, a simplified form of the Einsteinrelation η D ( p ) = κ/ (2 T E ) is adopted, where κ denotesthe heavy quark mean squared momentum change perunit time, and is usually referred as momentum-spacediffusion coefficient. Here we assume that in a “mini-mal model”, the heavy quark momentum variance in thelongitudinal direction equals that in the transverse di-rection (even though the microscopic calculation of thosetwo quantities are different): κ || = κ ⊥ = κ . The validityof such an assumption needs to be investigated in a futurecalculation. For this study we follow this assumption inorder to simplify our parameterization, and therefore theheavy quark transport coefficient is defined as ˆ q = 2 κ ⊥ ,which is the transverse momentum broadening.Assuming Gaussian shaped white noise, the ther-mal random force satisfies the relation (cid:104) ξ i ( t ) ξ j ( t (cid:48) ) (cid:105) = κδ ij δ ( t − t (cid:48) ), which indicates no correlation between ther-mal forces at different times.In order to describe the heavy quark dynamics in theintermediate and high p T region, the effective model-ing the radiative component of the heavy quark energyloss becomes necessary. A third force (cid:126)f g = − d(cid:126)p g /dt hence is introduced to account for the recoil force thatis experienced by the heavy quarks when they emitbremsstrahlung gluons, with (cid:126)p g being the emitted gluonmomentum. A higher twist calculation for the medium-induced gluon radiation is adopted from [48, 49]: dN g dxdk ⊥ dt = 2 α s P ( x )ˆ q g πk ⊥ sin (cid:18) t − t i τ f (cid:19) (cid:18) k ⊥ k ⊥ + x M (cid:19) . (4)where x is the fractional energy carried by the emit-ted gluon, k ⊥ is the gluon transverse momentum, P ( x )is the splitting function and τ f is the gluon formationtime. We relate the gluon transport coefficient ˆ q g andheavy quark transport coefficient ˆ q via the color factorsˆ q g = C A /C F ˆ q ( C F = 3 , C A = 4 / q , which characterizes the interaction strengthbetween heavy quarks and the medium. In the litera-ture, the spatial diffusion coefficient D s = 4 T / ˆ q is moreoften used in the diffusion equation, and this will be thephysical property that we are trying to extract from ouranalysis.The evolution of the QGP medium is simulated by a(2+1)-dimensional event-by-event viscous hydrodynami-cal model VISHNEW [50–52], which has been updatedto include both the shear and bulk viscosities with theshear-bulk coupling. The shear and bulk viscosities havebeen parametrized as temperature dependent. In ourcurrent study, the parameters related to properties ofthe soft medium as well as the initial condition are cali-brated through an independent Bayesian analysis of lighthadrons [16].
3. Hadronization
When the temperature of the QGP medium dropsbelow the critical temperature ( T c = 154 MeV), themedium undergoes a transition from a deconfined fluid toa confined hadron gas. The phenomenon of confinementinvolves non-perturbative processes and is not well under-stood. One often utilizes an instantaneous hadronizationmodel to convert the fluid medium into hadrons. On thetransition hypersurface, an ensemble of hadrons is gen-erated by sampling the momentum distribution from theCooper-Frye formula [53, 54]: E dN i d p = (cid:90) Σ f i ( x, y ) p µ d σ µ . (5)Heavy quarks hadronize into heavy mesons within a hy-brid model of instantaneous recombination and fragmen-tation. The momentum spectra of the meson producedby recombination is determined by the Wigner func-tion [22, 41]: dN M d p M = (cid:90) d p d p dN Q d p Q dN q d p q f WM ( (cid:126)p Q , (cid:126)p q ) δ ( (cid:126)p M − (cid:126)p Q − (cid:126)p q ) . (6)where (cid:126)p Q and (cid:126)p q are the momentum of heavy quark andlight quark that constitute the heavy meson, f WM ( (cid:126)p Q , (cid:126)p q )is the Wigner function in terms of the overlap betweenthe two initial partons and the final meson. A simplequantum harmonic oscillator is used to approximate thewave-function. For heavy quarks that do not recombinewith light partons, a fragmentation process via PYTHIAtakes place. It is found that high momentum heavyquarks tend to fragment while lower momentum heavyquarks tend to recombine with the thermal light partonsand hadronize into hadrons [41].
4. Hadronic re-scattering
After hadronization, the system continues expandingas an interacting hadron gas. Subsequent interactions be-tween heavy mesons and light hadrons (e.g scattering anddecay) after hadronization are modeled with UrQMD,which solves the Boltzmann equation for all the particlesin the system [55, 56]. UrQMD continues to evolve thesystem until the hadron gas is so dilute that all interac-tions have ceased and the system reaches its kinematicfreeze-out. The particle information is then collected tocalculate the final state observables that can be comparedwith experimental data.In our study of D mesons, the two main observablesare: the nuclear modification factor R AA , which quan-tifies the heavy quark in-medium energy loss and is ob-tained by taking the ratio of the heavy meson p T spec-tra measured in nucleus-nucleus collisions and the refer-ence spectra in proton-proton collisions, scaled by the bi-nary collision number; the harmonic flow coefficients v n ,which are the n th order coefficients in the azimuthal an-gle Fourier expansion of the emitted hadron spectra. Inour calculation, the second order harmonic elliptic flow v is calculated via both: the event-plane method [58]and the cumulant method [59], while the triangle flow iscalculated via the cumulant method [59]. B. Parameterization of the diffusion coefficient
One of the goal of the heavy-ion community for thenext few years is to quantitatively determine the heavyquark diffusion coefficient at sufficiently high precision.Since the diffusion coefficient is not a quantity that canbe directly measured, its determination requires an inter-play between both experiment and theory, meaning thevalues of the parameters which encode the heavy quarkdiffusion coefficient are obtained from a comparison be-tween experimental measurement and the correspondingtheoretical calculations.At high temperature and large momentum, the dif-fusion coefficient can be calculated using perturbativeQCD [60, 61]: The simplest possible diagram for heavyquarks interacting with light partons is given by two2 → Qq → Qq , Qg → Qg ), where the heavy quark transport coefficient equalsto [62]: ˆ q = (cid:10) ( (cid:126)p Q (cid:48) ) − (ˆ p Q · (cid:126)p Q (cid:48) ) (cid:11) . (7)where (cid:126)p Q ( (cid:126)p Q (cid:48) ) is the in-(out-)going momentum of theheavy quark. (cid:104) X (cid:105) is defined as: (cid:104) X (cid:105) = γ E Q (cid:90) d p q (2 π ) E q d p Q (cid:48) (2 π ) E Q (cid:48) d p q (cid:48) (2 π ) E q (cid:48) · X · f q ( (cid:126)p q )(2 π ) δ ( (cid:126)p Q + (cid:126)p q − (cid:126)p Q (cid:48) − (cid:126)p q (cid:48) ) (cid:88) |M| Qq → Q (cid:48) q (cid:48) . (8)where (cid:126)p q ( (cid:126)p q (cid:48) ) is the in-(out-)going momentum of the lightparton (light quark or gluon), and M → are the scat-tering matrices between heavy quarks and light partons.The leading-order pQCD calculation of diffusion coeffi-cient with respect to temperature and heavy quark mo-mentum is plotted in Fig. 1.It has been found in previous comparisons to data thatthe perturbative calculations are not sufficient to explainthe experimental findings such as the single non-photonicelectron suppression stemming from decay of the heavyflavor mesons (“single electron puzzle”) [63], or fail toto simultaneously describe both the heavy quark nu-clear modification factor R AA and elliptic flow v (“heavyquark R AA and v puzzle”) [64]. Moreover, it has beenargued that the convergence of the perturbative terms israther poor [25, 65]. In order to compensate for non-perturbative effects, one may introduce a K -factor toscale the scattering cross section by an ad-hoc param-eter, which will be adjusted until the model is able todescribe the experimental data. In this study, we use a p [GeV/c]0123456 ˆ q [ G e V / f m ] T =0.15 GeV T =0.25 GeV T =0.35 GeV T =0.45 GeV T =0.55 GeV . . . . . T [GeV] p =0 GeV/c p =10 GeV/c p =20 GeV/c p =30 GeV/c p =40 GeV/c FIG. 1. (Color online) A leading order pQCD calculationof ˆ q with respect to temperature and momentum. For thecalculation of ˆ q , a fixed coupling α s = 0 . m D = 4 πα s T . All s, u, t channel contributions for Qq → Qq, Qg → Qg are included. more generalized parametrization of the diffusion coef-ficient, which combines a linear temperature dependentcomponent and a pQCD component: D s πT ( T, p ) = 11 + ( γ p ) ( D s πT ) linear + ( γ p ) γ p ) ( D s πT ) pQCD . (9)The linear component ( D s πT ) linear = α · (1+ β ( T /T c − p = 0 GeV/c limit and can becompared to lattice QCD calculation of the spatial dif-fusion coefficient at zero momentum. The pQCD com-ponent is the contribution from perturbative processes,and is related to ˆ q pQCD , which has been calculated above,by ( D s πT ) pQCD = 8 π/ (ˆ q/T ). It should be noted thatthe standard spatial diffusion coefficient D s is defined atthe zero momentum p = 0 GeV/c limit. However, inthis work we use the notation D s to refer to the diffusioncoefficient in the full momentum range.The parameter α represents the spatial diffusion co-efficient at zero momentum near T c , parameter β is theslope of D s πT ( p = 0) above T c . The linear shape ofthe parametrization is inspired by the approximately lin-ear temperature dependence of the specific shear viscos-ity [38], as we assume an underlying relationship betweenthe transport properties of the QGP medium [66]. Theparameter γ controls the ratio between the linear com-ponent and the pQCD component. For p < /γ the lin-ear component dominates while for p > /γ the pQCDcomponent is dominant. A small value of γ indicates non-perturbative processes affect the heavy quark dynamicsinto the very high momentum region, and a large value of γ indicates a quick conversion to the pQCD dominatedregion. To better illustrate the dependence of the spatialdiffusion coefficient on γ , we plot D s πT as a functionof temperature (at fixed momentum) and momentum (atfixed temperature) for different values of γ in Fig. 2. Asshown in Fig. 2, the value of γ changes from 0 to 1 while D s π T p =10 GeV/c γ =0.0 γ =1.0 T =0.2 GeV ( D s πT ) pQCD ( D s πT ) linear . . . . T [GeV]05101520 D s π T p =50 GeV/c γ =0.0 γ =1.0 p [GeV/c] T =0.45 GeV ( D s πT ) pQCD ( D s πT ) linear FIG. 2. (Color online) An example of the spatial diffu-sion coefficient parametrization. The linear component uses( D πT ) linear = α · (1 + β ( T /T c − α, β ) = (1 . , . γ varying from 0 to1 while the color change from violet to red. the color changes from violet to red in the reverse rain-bow color scheme. For a large value of γ (red) the com-bined diffusion coefficient quickly converges to the pQCDcalculation, while for a small value of γ (violet) the dif-fusion coefficient still follows the linear contribution evenfor large momenta. III. PARAMETER CALIBRATION
In this section we summarize the workflow of theBayesian analysis that allows us to determine the highlikelihood parameter ranges of ( α, β, γ ) that govern thediffusion coefficient. More details on the Bayesian anal-ysis can be found in [16, 36, 38, 40].We first evaluate the improved Langevin model for alimited number of parameter values that are selected us-ing a Latin hypercube algorithm, and calculate the heavyquark observables as outputs from the model for thesevalues. The mapping from inputs to outputs can then beused to train a set of Gaussian process emulators, whichact as fast surrogate model of our Langevin model and areable to predict the output for any arbitrary input pointin the parameter space. A Markov chain Monte Carlo(MCMC) random walk through the parameter space isthen performed in order to calibrate the model param-eters on the experimental data. After the MCMC equi- librates, we obtain the posterior distributions of the in-put parameters, and thus the posterior estimate of theparametrized diffusion coefficient.
A. Training data preparation
Over the last few years, significant progress has beenmade related to the measurement of heavy flavor observ-ables, such as yields and/or flow cumulants of heavy fla-vor mesons, single electrons from heavy flavor hadronsemi-leptonic decays, heavy flavor tagged jets, quarko-nium yields, spectra and elliptic flow etc. However, thosedata sets differ greatly regarding the statistical and sys-tematic uncertainties and it is therefore not feasible tocombine all of them for our current analysis. For thisstudy, we only focus on the D -meson R AA and v , whichare very sensitive to the interaction mechanics betweenheavy quarks and the medium. These have been mea-sured in three different systems at RHIC and the LHC:AuAu collisions at √ s NN = 200 GeV, PbPb collisions at √ s NN = 2.76 TeV, and PbPb collisions at √ s NN = 5.02TeV. Table I summarizes the measurements and kine-matic/centrality cuts of the observables that have beenused.To illustrate the degree to which the model’s calcu-lation are affected by the particular form of the tem-perature and momentum dependence of the diffusioncoefficient, we first sample 60 design points ( ˜ X = α β FIG. 3. (Color online) 60 design points ˜ X generated by theLatin hypercube algorithm, projected in the ( α, β ) dimen-sions. The flat histograms on the edge show an uniform priordistribution of the parameters. TABLE I. D -meson variables to be compared between model calculation and experimental measurementsExperiment variables kinematic cut centrality refAuAu@200 GeV R AA ( p T ) 6 p T bins from 2 − | y | < v (EP)( p T ) 8 p T bins from 1 − | y | < v (EP)( p T ) 8 p T bins from 1 − | y | < R AA ( n part ) 6 centrality bins, 5 < p T < | y | < . R AA ( n part ) 6 centrality bins, 8 < p T <
16 GeV/c, | y | < . v (EP)( p T ) 6 p T bins from 2 −
16 GeV/c, | y | < . R AA ( p T ) 10 p T bins from 3 −
36 GeV/c, | y | < . v { } ( p T ) 11 p T bins from 1 −
40 GeV/c, | y | < v { } ( p T ) 8 p T bins from 1 −
40 GeV/c, | y | < p T [GeV]10 − − d N / d p T | | y | < . ( µ b / G e V / c ) [email protected] TeVbox syst.uncbar stat.unc D -meson, FONLL D ALICE D + ALICE D ∗ + ALICEtotal D -meson FIG. 4. (Color online) Proton-proton reference spectrum cal-culated by FONLL+PYTHIA, used as the reference spectrain order to calculate D -meson R AA . The experimental resultsdenoted total D-meson are the sum of D , D + and D ∗ + datafrom ALICE collaboration. ( x , · · · , x ) T , x i = ( α, β, γ ) T ) in the parameter space(the tilde symbol ˜ x, ˜ y represents the training datasets,whose outputs are calculated from the physical model, itis used to distinguish the predicted datasets, which arelater obtained from emulators and labeled with the starsymbol x ∗ , y ∗ ). Those design points are semi-randomlysampled in the 3-dimensional parameter space using aLatin hypercube algorithm, which aims at spreading thesamples evenly across all possible values [67], and there-fore, a small amount of samples O (10p) is sufficient totrain the Gaussian process emulators to interpolate the p -dimensional parameter space. The parameter space isdeliberately chosen to be wide enough in order to coverthe full range of likely values, and is listed in Tab. II.Figure 3 visualizes the uniform distribution the initialdesign points, projected in the ( α, β ) plane.At each of the design points, we generate 5000 mini-mum bias hydro events and for each hydro event, heavyquarks are oversampled to reduce the theoretical statis-tical uncertainty. The heavy meson observables are then TABLE II. Prior range and description for the parametersthat determines the diffusion coefficientsParameter Description Range α D s πT at T c β slope of ( D s πT ) linear above T c γ ratio between D linear s and D pQCD s calculated as the following: the events are first binnedinto different centrality classes according to the final statecharged hadron multiplicity N ch at mid-rapidity. The D -meson selection is based on the corresponding experi-mental kinematic cuts. In order to calculate the nuclearmodification factor R AA , a proton-proton collision ref-erence is needed. It is calculated using a heavy quarkFONLL distribution followed by a fragmentation processperformed by PYTHIA, and the D -meson yields in theproton-proton reference are compared with experimentaldata in Fig. 4. For the calculation of D -meson ellipticflow, we try to match the experimental methods as far aspossible, therefore for the AuAu collision and PbPb col-lision at 2.76 TeV, an event-plane (approximated by theinitial participant-plane) method is used, while for thePbPb collision at 5.02 TeV, the two-particle cumulantmethod is used, though little difference has been noticedfor the two different methods [69, 71]. In total there are69 experimental data points to calibrate against.Figure 5 compares the 60 sets of model calculation withthe experimental data (black dots with errorbars). Wecan see that the model’s outputs span a wide range inobservable space as the input parameters have been ran-domly distributed in the parameter space. These input˜ X n × p and output ˜ Y n × m matrices (where n = 60 is thenumber of input parameter points, p = 3 is the dimensionof input parameters, m = 69 is the dimension of outputat each of the input point) will then be used to train theGaussian process emulators. B. Gaussian process emulator
In order to calibrate our parameters, a random walkthroughout parameter space will be performed, where p T [GeV]0 . . . . . R AA STAR D , 0-10% p T [GeV]0 . . . . v STAR D , 0-80% p T [GeV]0 . . . . v STAR D , 10-40% A u A u @ G e V part . . . . R AA ALICE D , D ∗ , D ∗ + < p T < part . . . . R AA ALICE D , D ∗ , D ∗ + < p T <
16 GeV/c p T [GeV] − . . . . . . v ALICE D , D ∗ , D ∗ + P b P b @ . T e V p T [GeV]0 . . . . R AA ALICE R AA , 30-50% p T [GeV]0 . . . . . v CMS v ,30-50% p T [GeV]0 . . . . . v CMS v , 10-30% P b P b @ . T e V FIG. 5. (Color online) Improved Langevin model calculation of D -meson observables, compared to experimental data spanningthe full range of the explored parameter space (i.e. the “prior”). Each frame contains 60 lines, corresponding to the 60 designpoints of the analysis. From top to bottom: ( top ) AuAu collisions at 200 GeV: D -meson R AA ( p T ) in the 0-10% centrality binand v ( p T ) in the 0-80%, 10-40% centrality bins; ( middle ) PbPb collisions at 2.76 TeV: D -meson R AA as function of centralityat high momentum range 5 < p T < < p T <
16 GeV/c, v as function of p T in 30-50% centrality bin; ( bottom )PbPb collisions at 5.02 TeV: D -meson R AA ( p T ) in 30-50% centrality bin and v ( p T ) in 30-50% and 10-30% centrality bins.Experimental data are measured by STAR [68, 69], ALCIE [70–72] and CMS [73]. each step is accepted or rejected according to the rela-tive likelihood. Taking a random walk throughout the 3-dimensional parameter space requires O (1000) steps andthe number increases exponentially if we try to include more parameters. At each step one needs to generate asample of events and calculate the model’s output for theobservables in order to evaluate the likelihood and takeaction for the next step. Given the amount of compu-tational time required to evaluate the likelihood at onepoint in parameter space, such a method is not feasible.To overcome this difficulty, a set of emulators is used tofunction as a fast surrogate model that can predict thephysical model’s output at any arbitrary point in the pa-rameter space. In this study, we construct the Gaussianprocess (GP) emulator, which is an mapping from a n-dimensional input space to a normal distributed output.It not only interpolates and predicts the model outputafter being trained, but also provides the uncertaintiesof its prediction. More details of the GP emulator canbe found in Ref. [74]. Here we will briefly summarize thebasic idea of GP emulator.Consider a physical model (e.g. our improved Langevinframework), whose output of a physical process is deter-mined by a set of input parameters y = f ( x ). We sup-pose that the physical model has been evaluated at n input points in the p -dimensional parameter space, (theinput parameter matrix ˜ X = ( x , . . . , x n ) T ). At each in-put point, the model has one output y i = f ( x i ), yieldsto an n -dimensional output vector ˜ y :˜ X = x ... x p ... ...x n ... x np ⇒ ˜ y = y ...y n . (10)The output ˜ y can be viewed as a conditioned Gaussianprocess which is a collection of normal distributions:˜ y = GP ( ˜ X ) ∼ N ( µ ( ˜ X ) , K ˜ X, ˜ X ) . (11)where µ ( ˜ X ) is the mean vector of each input, and K ˜ X, ˜ X = σ ( x , x ) · · · σ ( x , x n )... . . . ... σ ( x n , x ) · · · σ ( x n , x n ) . (12)is the covariance matrix. It is constructed by the covari-ance function σ ( x , x (cid:48) ) and characterizes the correlationbetween different inputs.In order to predict the model output y ∗ at any other in-put x ∗ (the star symbol are used to represent the datasetswhose outputs are predicted from the emulators), one canwrite the joint multivariate normal distribution: (cid:18) y ∗ ˜ y (cid:19) ∼ N (cid:18)(cid:18) µ ( x ∗ ) µ ( ˜ X ) (cid:19) , (cid:18) K ∗ , ∗ K ∗ , ˜ X K ˜ X, ∗ K ˜ X, ˜ X (cid:19)(cid:19) . (13) K ∗ , ∗ , K ∗ , ˜ X and K ∗ , ˜ X have the same form as Eqn. (12)but with different x . The distribution of a predictiveoutput y ∗ can be solved by: y ∗ ∼ N ( µ, K ) (14) µ = µ ( x ∗ ) + K ∗ , ˜ X K − X, ˜ X (˜ y − µ ( ˜ X )) (15) K = K ∗ , ∗ − K ∗ , ˜ X K − X, ˜ X K ˜ X, ∗ . (16)With the knowledge of a training dataset containing a setof inputs ( ˜ X ) and outputs (˜ y ), and the covariance matrix m )0.960.981.00 v a r i a n c e e x p l a i n e d b y P C s C V ( m ) FIG. 6. (Color online) Cumulative variance explained by thefirst m (cid:48) principal components. The first few PCs are able toexplain most of the total variance ( K ), one is able to solve the equations and calculate thedistribution for any other input x ∗ .The inference of the Gaussian process is determined bythe covariance function σ ( x , x (cid:48) ). Variance choice can bemade for the covariance function, based on our knowledgeand assumption of the input parameters. In this study weuse a popular squared-exponential function with a noiseterm: σ ( x , x (cid:48) ) = σ G exp (cid:34) − m (cid:88) k =1 ( x k − x (cid:48) k ) l k (cid:35) + σ n δ x , x (cid:48) . (17)The squared-exponential covariance function implies ourintuition that the inputs in the parameter space that areclose to each other are highly correlated, whilst thosefar away are uncorrelated, the correlation strength be-tween pairs of inputs is controlled by the hyperparameter( σ G , l k ). With this covariance function, the Gaussian pro-cess is very smooth as the covariance is infinitely differen-tiable. In this study, the hyperparameters ( σ G , l k , σ n ) aredetermined in a manner of “best fit parameters” whichmaximizes the parameter likelihood function[75]. C. Principle component analysis
A Gaussian process is essentially a mapping from aninput vector to a scalar output. The output from ourLangevin mode is an m -dimensional vector ( m = 69):˜ X = x ... x p ... ...x n ... x np ⇒ ˜ Y = y ... y m ...y n ... y nm . (18)One can construct a GP emulator for each of the ob-servables. However, as the elements in the output arehighly correlated, it is useful to reduce a high dimen-sional and correlated output to a lower dimensional and . . . . v . . . . . e m u l a t o r p r e d i c t e d R AA . . . v . . . . e m u l a t o r p r e d i c t e d v P b P b @ . T e V FIG. 7. (Color online) Validation of the Gaussian emula-tors for PbPb collisions at 2.76 TeV. Each point representsthe emulator’s predicted value with respect to the model cal-culated (true) value. (left:) prediction for D -meson v at 30-50%; (right:) prediction for high momentum D -meson R AA at different centrality bins. orthogonal principal components (PCs), which are thelinear combinations of the output observables and thatprovide information on the most important componentsof the dataset.In practice, the training output ˜ Y is standardized (bysubtracting the mean and divided by the standard devi-ation for each observable), and decomposed via the sin-gular value decomposition (SVD): Y m × n = U m × n S n × n V Tn × n . (19)The columns of U ( V T ) are the left(right)-singular vec-tor of Y , which are sets of orthogonal eigenvectors of Y Y T ( Y T Y ). The output matrix Y can then be trans-formed into principal component space: Z = √ nY V. (20) S is the diagonal matrix whose diagonal elements λ i, ( i =1 , ··· ,n ) are the squared roots of the eigenvalues of Y T Y . The eigenvalues λ i are proportional to the vari-ance that contributed the i -th PC, and are sorted intodescending order. The cumulative variance explained bythe first m (cid:48) -th PCs ( m (cid:48) (cid:54) m ) then equals to: CV ( m (cid:48) ) = (cid:80) m (cid:48) i =1 λ i (cid:80) mi =1 λ i . (21)As shown in fig. 6, the first few PCs are sufficient toexplain most of the variance of the model outputs. Inthis study, we use 8 PCs and for each PC a GP emulatoris constructed, which is a significant reduction from theoriginal 69 GPs that mapped directly onto the numberof data points in the calibration dataset.Once the principal component z has been determined,we obtain the outputs in the physical observable spaceby the inverse transformation: y = 1 √ n z V. (22) In order to test the emulators’ ability to predict themodels’ output, we generate another 15 sets of test in-puts and perform the full Langevin model calculation ateach of these test inputs. Meanwhile the trained GP em-ulators also predict the output for each of these inputparameters. Figure 7 compares the prediction from theemulators to the calculation from the improved Langevinmodel for PbPb collisions at 2.76 TeV. For each test in-put point, 18 observation of R AA and v are calculatedat different centralities and p T bins. The black lines arethe y = x reference, and each dot represent the emula-tors’ prediction with respect to the models’ calculation.As visualized in Fig. 7 the GP emulators in general workvery well. We should note that the emulators provide theuncertainty for each prediction, therefore the errorbarsshown in the figure correctly capture the uncertaintiesunderlying in the emulation. D. MCMC calibration
According to Bayes theorem, the probability for thetrue parameter x given the experimental data y exp andobserved ( ˜ X, ˜ Y ) is proportional to the likelihood of theparameter x and its prior distribution: P ( x | ˜ X, ˜ Y , y exp ) ∝ P ( ˜ X, ˜ Y , y exp | x ) P ( x ) . (23)where P ( x | ˜ X, ˜ Y , y exp ) is the posterior distribution of pa-rameters x given the observation of ( ˜ X, ˜ Y , y exp ), whichis our main results from this analysis; P ( ˜ X, ˜ Y , y exp | x ) isthe likelihood of observing ( ˜ X, ˜ Y , y exp ) given the param-eter x , and P ( x ) is the prior distribution of parameter x . Our goal here is to find the posterior probability dis-tribution of the parameters P ( x | ˜ X, ˜ Y , y exp ) which wouldoptimally reproduce the experimental data using our im-proved Langevin model. In order to determine the pos-terior distribution we perform a MCMC random walk inparameter space following the Metropolis-Hasting algo-rithm [76]. During the random walk, each step is ac-cepted or rejected according to the relative likelihood.Assuming a Gaussian structure for the uncertainties, thelog likelihood function has the following form as a func-tion of the output y :log P ( ˜ X, ˜ Y , y exp | x ) = −
12 ( y − y exp ) T Σ − ( y − y exp ) −
12 log | Σ | − m π. (24)Σ is the covariance matrix which accounts for all quan-tifiable uncertainties. There are various contributions tothese uncertainties, such as the emulator prediction un-certainty, the experimental statistic and systematic un-certainties, and the physical model statistic and system-atic uncertainties. Identifying all these uncertainties can0 p T [GeV]0 . . . . . R AA STAR D , 0-10% p T [GeV]0 . . . . v STAR D , 0-80% p T [GeV]0 . . . . v STAR D , 10-40% A u A u @ G e V part . . . . R AA ALICE D , D ∗ , D ∗ + < p T < part . . . . R AA ALICE D , D ∗ , D ∗ + < p T <
16 GeV/c p T [GeV] − . . . . . . v ALICE D , D ∗ , D ∗ + P b P b @ . T e V p T [GeV]0 . . . . R AA ALICE R AA p T [GeV]0 . . . . . v CMS v , 30-50% p T [GeV]0 . . . . . v CMS v , 10-30% P b P b @ . T e V FIG. 8. (Color online) Emulator predictions for 200 random input parameters sampled from the posterior distributions. Thisfigure is similar to Fig.5 but with the input parameters chosen from the posterior distribution, and the outputs are predictionsfrom the GP emulators. be very difficult, especially for systematic uncertaintieswhich in general may have correlations among each other.For the current analysis we only consider experimentalstatistical and uncorrelated systematic uncertainties:Σ = diag( σ ) + diag( σ ) + diag( σ ) . (25)where σ stat is the experimental statistical error, σ sys isthe experimental systematic error (uncorrelated) for each observable, σ GP is the theoretical uncertainty from Gaus-sian process emulator predictions. At the current stage,all the uncertainties are assumed to be uncorrelated forthe purpose of simplicity as well as maximizing the over-all constraint.1 β . . . . α . . . . . γ β . . . . . γ yellow: AuAu@200 GeVgreen: [email protected] TeVblue: [email protected] TeV FIG. 9. (Color online) Posterior distributions of the diffu-sion coefficient parameters ( α, β, γ ) for each individual col-lision system. The diagonal plots are the histogram of theMCMC samples with other parameters integrated out. Theoff-diagonal plots display the joint distributions between pairsof parameters. The three different colors refer to the three dif-ferent analyses that calibrate on three different sets of data.
IV. RESULTSA. Posterior distributions
To evaluate the success of the calibration, 200 pointsin parameter space are randomly chosen from the equili-brated MCMC trace and evaluated by the Gaussian emu-lators. Figure 8 visualizes the corresponding observablesand compares it with the experimental data. The presen-tation is similar to Fig. 5 but with calibrated parameters.For each plot, two posterior outputs are presented, eachone corresponding to an independent Bayesian analysisbut calibrated on different experimental datasets. Thered lines correspond to the Bayesian analysis calibratedon all the output observables listed in Tab. I, whereasthe yellow, green and blue lines correspond to calibra-tions on data of a single beam energy. We find that aftercalibration, our improved Langevin approach is capableof describing the experimental data reasonably well. Thebiggest deviations are found for a few R AA points at veryperipheral centrality and low p T : peripheral collisions arenot well described by our hydrodynamical background.Also, the modeling of hadronization in the low p T regionis challenging due to significant non-perturbative effects.The main results of a Bayesian analysis are the pos-terior probability distributions of the input parameters α = 1 . +0 . − . β β = 1 . +1 . − . . . . . α . . . γ β .
15 0 .
30 0 . γ magenta: LHC energiesred: all energies γ = 0 . +0 . − . FIG. 10. (Color online) Posterior distributions of the diffu-sion coefficient parameters ( α, β, γ ) for the analysis combiningall three collision systems. The red color refers to the resultthat combines only the LHC energies: PbPb collisions at 2.76and 5.02 TeV; the magenta color refers to the analysis com-bining all three energies at RHIC and LHC. ( α, β, γ ). The posterior distributions are given in Fig. 9and 10, where the results of 5 independent Bayesian anal-yses are presented as histograms. Each analysis followsthe same procedure described in the previous few sectionsbut calibrates on different sets of experimental data, andis shown using different colors. For examples, in Fig. 9the blue histograms represent the distributions calibratedon the 5.02 TeV PbPb collision data, the green ones rep-resent the calibration on the 2.76 TeV PbPb data, andthe yellows ones represent the calibration on the 200 GeVAuAu data. In Fig. 10 the magenta histogram corre-sponds to the analysis using the data from two collisionenergies at the LHC simultaneously, while the red onesdenote the analysis using all the observable across threedifferent systems and listed in Tab. I. In each figure, thehistograms along the diagonal are the marginal posteriorprobability distributions of each parameter ( α, β, γ ), withall the other parameters integrated out. The off-diagonalcontour plots are the joint distributions which show thecorrelations among pairs of the parameters.Figure 9 and 10 indicate that parameter α is wellconstrained, peaking around (1 . ∼ .
0) for all analy-ses. This parameter determines the diffusion coefficient D s πT at 0 momentum near T c . The slope parameter β is poorly contrained, although a negative correlationbetween α and β is observed. Parameter γ controls theratio between the linear component and the pQCD com-2 . . . . . . T [GeV]05101520 D s π T p =0 GeV/c median valueprior90% C.R . . . . . . T [GeV] p =10 GeV/c ( D s πT ) pQCD median valueprior90% C.R . . . . . . T [GeV] p =50 GeV/c ( D s πT ) pQCD median valueprior90% C.R p [GeV/c]05101520 D s π T T =0.154 GeV ( D s πT ) pQCD median valueprior90% C.R p [GeV/c] T =0.35 GeV ( D s πT ) pQCD median valueprior90% C.R p [GeV/c] T =0.55 GeV ( D s πT ) pQCD median valueprior90% C.R FIG. 11. (Color online) Posterior results of the spatial diffusion coefficient from the Bayesian analysis calibrated on thecombined dataset from three different systems at RHIC and the LHC.( top ) the spatial diffusion coefficient D s πT as a functionof temperature at fixed momentum ( p = 0 GeV/c, p = 10 GeV/c and p = 50 GeV/c); ( bottom ) the spatial diffusion coefficient D s πT as a function of momentum at fixed temperature ( T = 154 MeV, T = 350 MeV, T = 550 MeV). The grey area refers tothe prior range before calibration, while the red region refers to the posterior range after calibrating on experimental data. Theblack dashed line refers to the diffusion coefficient from a leading order pQCD calculation, the red lines are the parametrizeddiffusion coefficient using the median value of the posterior parameter distributions. ponent. As shown in Fig. 9, the distribution for γ ex-tracted from AuAu collisions at 200 GeV favors a slightlysmaller value than that from the LHC energies, indicat-ing a stronger contribution from the linear componentand a slower convergent to the pQCD results in AuAucollision. For the combined analysis of the LHC ener-gies and all three energies, γ peaks around (0 . ∼ . − /γ , as γ p = 1is the momentum region where linear and pQCD compo-nent contributes equally.The width of the posterior distributions is affected bythe uncertainty we have applied in the analysis. Thesmaller the uncertainty, the stronger the constraint, andtherefore a narrower width of the posterior distributions. This may explain why among the three different colli-sion systems, the posterior distributions from 5.02 TeVPbPb collisions generally show a smaller width and thecombined calibration is mostly driven by the higher pre-cision data from PbPb 5.02 TeV collisions. Higher preci-sion experimental data and a better understanding of thetheoretical uncertainties will yield a significantly betterconstraint on the parameters. B. Heavy Quark Diffusion Coefficient
Having established the posterior distribution of the pa-rameters α, β and γ we can now utilize the parametriza-tion of the spatial diffusion coefficient Eqn. (9), to ex-tract the posterior range of this quantity. Figure 11 dis-plays the estimate of the spatial diffusion coefficient, asa function of temperature and momentum respectively.The gray area represents the prior range before the cali-3 T /T c D s π T ( p = ) MC@sHQ, elastic K=1.5MC@sHQ, ela+rad K=0.8PHSDQPM(Catania), BMQPM(Catania), LGV Duke-LGV, medianDuke-LGV, 90% C.Rc-quark T-matrix U-potc-quark lattice Ding et.alHQ lattice Banergee et.al . . . T [GeV]02468 ˆ q / T ( p = G e V / c ) c-quark, pQCDc-quark, LBTDuke-LGV, medianDuke-LGV, 90% C.RJET(light quark) FIG. 12. (Color online) Comparison of the heavy quark diffusion coefficients across multiple approaches available in theliterature. ( left ) spatial diffusion coefficient at zero momentum D s πT ( p = 0). ( right ) momentum diffusion coefficient ˆ q/T at p = 10 GeV. bration, while the red area is the posterior estimate ex-tracted from the combined analysis with 90% credibility.Here we use the result of the combined analysis that cal-ibrates on all the three systems. We note, however, theposterior ranges of D s πT do not differ much betweendifferent analysis, even though the posterior distributionsof the parameters ( α, β, γ ) show some deviation. Thedashed black lines represent the diffusion coefficient cal-culated in leading order pQCD. The solid red lines depictthe diffusion coefficient using the median value of the pa-rameters from the posterior distributions.On the upper panel of Fig. 11, the diffusion coefficientis plotted as function of temperature for 3 different val-ues of the heavy quark momentum ( p = 0 GeV/c, p = 15GeV/c and p = 50 GeV/c): for p = 0 GeV/c, the diffu-sion coefficient is solely determined by the linear compo-nent, with the D s πT ( p = 0) ∼ − T c , (whichis the range of parameter α for 5-95% percentiles). Thetemperature dependence of D s πT is not far remote froma simple linear relationship with positive slope. In ad-dition, we notice that the 90% credibility range suffersthe least uncertainties in a temperature range around T ∼ −
250 MeV, which we argue, is approximately theaverage temperature that heavy quarks experience dur-ing their propagation path. At higher temperature, theposterior range of the spatial diffusion coefficient broad-ens. A likely reason for this trend is due to the shortamount of time the bulk matter retains this high tem-perature at the beginning of the system’s evolution. Asthe system expands quickly, it rapidly cools, leaving onlya short period of time for the heavy quarks to interactwith the medium at that temperature, and therefore less information can be obtained at high temperatures.On the lower panel of Fig. 11 we explore the momen-tum dependence of the diffusion coefficient for three dif-ferent temperatures ( T = 150 MeV, T = 350 MeV and T = 550 MeV). As the heavy quark momentum increases,the uncertainties of the posterior range decrease. At highmomentum, the diffusion coefficient reflects that of thepQCD calculation, which is obtained from 2 → ( γ p ) γ p ) ,which varies only little for high mo-menta. In the low momentum region, the parameterizeddiffusion coefficient shows completely different behaviorfrom the pQCD calculation, which can be only the resultof the non-negligible contribution from non-perturbativeeffects, which are clearly needed in order to obtained arealistic description of the heavy quarks at low and in-termediate momentum.In Fig. 12 we compare our estimate of the diffusion co-efficient with the coefficient used or calculated by a num-ber of other models in the market [23, 27–30, 77] as wellas with lattice QCD calculations [31, 32]. The left frameshows the temperature dependence of the spatial diffu-sion coefficient at 0 momentum D s πT ( p = 0). Our anal-ysis is consistent with lattice QCD calculations within theuncertainties currently inherent in lattice QCD calcula-tions. Although the diffusion coefficients used in differentmodels are rather different, they all have a minimum near T c with a value range of D s πT ( p = 0) = 1 − q at p = 10 GeV/c with the resultsfrom the LBT model and the JET collaboration [30, 77].The transport coefficient ˆ q is roughly comparable with4 p T [GeV/c] − . . . . . . v AuAu200 median all-median
STAR v , 10-40% p T [GeV/c] v , 40-80% p T [GeV/c] STAR v , 0-80% A u A u @ G e V p T [GeV/c] − . . . . . v PbPb5020 median LHC-median all-median
CMS, v , 0-10% p T [GeV/c] CMS v , 10-30% p T [GeV/c] CMS v , 30-50% P b P b @ . T e V FIG. 13. (Color online) D -meson v predicted by the improved Langevin model at AuAu collisions at 200 GeV and PbPbcollisions at 5.02 TeV, compared to experimental data measured by STAR [79] and CMS [73]. The input parameters ( α, β, γ )were set to the median value of the posterior distribution for the AuAu at 200 GeV dataset (AuAu200 median), PbPb at 5.02TeV dataset (PbPb5020 median), the combined LHC energy dataset (LHC-median), and the combined data set for all threesystems (all-median). those two. A detailed investigation on the causes ofthe differences in predicted or calculated transport co-efficients by the different approaches is of great inter-est, but beyond the scope of this work. We should notethat in order to make a valid apple-to-apple compari-son among different models, especially for those trans-port approaches that depend on the surrounding matter,it would be essential for the heavy quarks to propagatein the same QGP medium evolution, which is not thecase for this particular comparison. For example, thePHSD group [23] explores a non-equilibrium microscopicdescription of the bulk evolution, the Catania group [27]describes the medium by solving a relativistic Boltzmannequation for light partons, while the others use a hydro-dynamic description. Additionally the properties of themedium and the choice for the equation of state will af-fect the extraction of the diffusion coefficients as well. C. Model validation: triangular flow
The robustness and quality of our description of in-medium heavy-quark dynamics and our extraction of theheavy quark diffusion coefficient can be tested by makingpredictions of observables that have not been part of thecalibration. Here, we focus on measurements of higher or-der flow cumulants for D mesons, in particular v , whichhas been predicted as further valuable HQ observablesin [78]. First measurements of the D meson v have re-cently become available by the CMS and STAR collabo-rations [73, 79]. We calculate the D -meson v in PbPbcollisions at 5.02 TeV using the median values of α, β and γ as determined by our calibration and compare ourresults to experimental measurements in Fig. 13. Eachframe contains three results from the improved Langevinmodel calculation stemming from our different calibra-tions: utilizing the 200 GeV AuAu data, utilizing the5.02 TeV PbPb data, utilizing both LHC data sets, andutilizing the combined datasets from all three collisionenergies and systems. As we have mentioned before, the5posterior diffusion coefficients using the median valuesdo not differ much from each other in the different analy-ses. The robustness of the analysis has been confirmed inFig. 13, where we show that the D -meson v in AuAu col-lisions and in all three measured centrality bins in PbPbcollisions as predicted by our calibration agrees very wellwith the data, irrespective of the dataset used to deter-mine the diffusion coefficient. V. CONCLUSION
In summary, we have applied state-of-the-art Bayesianmethodology to systematically extract the heavy quarkdiffusion coefficient from a model-to-data analysis of ourimproved Langevin model for the in-medium heavy quarkevolution. By calibrating to the experimental data of D -meson R AA and v measured in AuAu collisions at 200GeV, PbPb collisions at 2.76 TeV and 5.02 TeV, we areable to extract a posterior range of the diffusion coeffi-cient. Our analysis is compatible with lattice QCD cal-culations within uncertainties that are inherent in thelattice calculations. With the extracted parameters, ourimproved Langevin model has been shown to be able toreproduce the experimental data of R AA and v at bothRHIC and the LHC simultaneously, and is able to de-scribe well observables that are not included in the cali-bration, such as D -meson v .Our parametrization of the spatial diffusion coefficientscombines a linear temperature dependent component –accounting for a non-perturbative contribution – anda pQCD component – calculated from a leading orderpQCD approach with a fixed coupling of α s = 0 .
3. Itsmoothly interpolates between the linear component andthe pQCD component, and converges to the pQCD calcu- lation in the large momentum limit. The spatial diffusioncoefficient at zero momentum D s πT ( p = 0) varies be-tween 1-3 near T c and exhibits a positive slope for itstemperature dependence above T c . Even at momenta inthe range of 10 −
20 GeV/c, the non-perturbative contri-bution can not be ignored.In future work, we shall improve our treatment of thedifferent sources of uncertainty, both theoretical and ex-perimental, that can affect the outcome of our analy-sis. In addition, we plan to improve our physics modelby taking the pre-equilibrium phase of the reaction ex-plicitly into account [80] and to apply the model-to-data framework to different dynamical models of heavyquark in-medium evolution, for example a comparisonbetween Langevin and Boltzmann dynamics. Moreover,this study serves as the first application of a Bayesianmodel-to-data analysis to the heavy flavor dynamics inheavy-ion collisions and we intend to expand its applica-tion to the study of other rare probes as well.
ACKNOWLEDGMENTS
The authors thank Guang-You Qin, Jussi Auvinen,Ralf Rapp, and Weiyao Ke for valuable discussion andpositive feedback. This work has used the CPU hoursprovided by the Open Science Grid, which is supportedby the National Science Foundation and the U.S. De-partment of Energy’s Office of Science. Y.X. and S.A.B.have been supported by the U.S Department of Energyunder grant DE-FG02-05ER41367 and J.B. has been sup-ported by the National Science Foundation under grantACI-1550225. S.C is supported by the U.S DoE undergrant DE-SC0013460 and the National Science Founda-tion under grant ACI-1550300. [1] J. D. Bjorken, Phys. Rev. D , 140 (1983).[2] E. V. Shuryak, Phys. Rept. , 71 (1980).[3] L. D. McLerran, Rev. Mod. Phys. , 1021 (1986).[4] D. Teaney, J. Lauret, and E. V. Shuryak, Phys. Rev.Lett. , 4783 (2001), nucl-th/0011058.[5] S. S. Adler et al. [PHENIX Collaboration], Phys. Rev.Lett. , 052301 (2006), nucl-ex/0507004.[6] K. H. Ackermann et al. [STAR Collaboration], Phys.Rev. Lett. , 402 (2001),nucl-ex/0009011.[7] B. I. Abelev et al. [STAR Collaboration], Phys. Rev.Lett. , 192301 (2007) Erratum: [Phys. Rev. Lett. ,159902 (2011)], nucl-ex/0607012[8] H. Song and U. W. Heinz, J. Phys. G , 064033(2009),arXiv:0812.4274 [nucl-th].[9] D. T. Son and A. O. Starinets, Ann. Rev. Nucl. Part.Sci. , 95 (2007),arXiv:0704.0240 [hep-th][10] P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev.Lett. , 111601 (2005),hep-th/0405231.[11] M. Luzum and P. Romatschke, Phys. Rev. C ,034915 (2008) Erratum: [Phys. Rev. C , 039903(2009)],arXiv:0804.4015 [nucl-th] [12] U. W. Heinz, H. Song and A. K. Chaudhuri, Phys. Rev.C , 034904 (2006), nucl-th/0510014.[13] B. Schenke, S. Jeon and C. Gale, Phys. Rev. Lett. ,042301 (2011),arXiv:1009.3244 [hep-ph].[14] P. Bozek, Phys. Rev. C , 034901(2012),arXiv:1110.6742 [nucl-th].[15] I. Karpenko, P. Huovinen and M. Bleicher, Comput.Phys. Commun. , 3016 (2014),arXiv:1312.4160 [nucl-th].[16] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu,and U. Heinz, Phys. Rev. C94 , 024907 (2016),arXiv:1605.03954.[17] S. Wicks, W. Horowitz, M. Djordjevic, and M. Gyulassy,Nucl. Phys.
A783 , 493 (2007), arXiv:nucl-th/0701063.[18] R. Sharma, I. Vitev, and B.-W. Zhang, Phys. Rev.
C80 ,054902 (2009), arXiv:0904.0032.[19] M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner,Phys.Rev.
C90 , 024907 (2014), arXiv:1305.3823.[20] J. Uphoff, O. Fochler, Z. Xu, and C. Greiner, Phys. Rev.
C82 , 044906 (2010), arXiv:1003.4200.[21] M. He, R. J. Fries, and R. Rapp, Phys. Lett.
B701 , 445 (2011), arXiv:1103.6279.[22] S. Cao, G.-Y. Qin, and S. A. Bass, Phys.Rev. C88 ,044907 (2013), arXiv:1308.0617.[23] T. Song, H. Berrehrah, D. Cabrera, W. Cassing andE. Bratkovskaya, Phys. Rev. C , no. 3, 034906 (2016),arXiv:1512.00891[24] H. van Hees, M. Mannarelli, V. Greco and R. Rapp,Phys. Rev. Lett. , 192301 (2008), arXiv:0709.2884.[25] S. Caron-Huot and G. D. Moore, JHEP , 081(2008), arXiv:0801.2173 [hep-ph].[26] M. He, R. J. Fries, and R. Rapp, Phys. Rev. Lett. ,112301 (2013), arXiv:1204.4442.[27] F. Scardina, S. K. Das, V. Minissale, S. Plumari andV. Greco, (2017) arXiv:1707.05452.[28] F. Riek and R. Rapp, Phys. Rev. C82 , 035201 (2010),arXiv:1005.0769.[29] P. Gossiaux and J. Aichelin, Phys.Rev.
C78 , 014904(2008), arXiv:0802.2525.[30] S. Cao, T. Luo, G. Y. Qin and X. N. Wang, Phys. Rev.C , no. 1, 014909 (2016), arXiv:1605.06447.[31] H. T. Ding et al. , Phys. Rev. D86 , 014509 (2012),arXiv:1204.4945.[32] D. Banerjee, S. Datta, R. Gavai, and P. Majumdar, Phys.Rev.
D85 , 014510 (2012), arXiv:1109.5738.[33] A. Francis, O. Kaczmarek, M. Laine, T. Neuhaus andH. Ohno, Phys. Rev. D , no. 11, 116003 (2015),arXiv:1508.04543[34] D. Higdon, J. D. McDonnell, N. Schunck, J. Sarich,and S. M. Wild, J. Phys. G42 , 034009 (2015),arXiv:1407.3017.[35] D. Higdon, J. Gattiker, B. Williams, and M. Rightley,J.Amer.Stat.Assoc. , 570 (2008).[36] J. Novak et al. , Phys.Rev.
C89 , 034917 (2014),arXiv:1303.5769.[37] S. Habib, K. Heitmann, D. Higdon, C. Nakhleh andB. Williams, Phys. Rev. D , 083503 (2007), astro-ph/0702348 [ASTRO-PH].[38] J. E. Bernhard et al. , Phys. Rev. C91 , 054910 (2015),arXiv:1502.00339.[39] J. Auvinen, I. Karpenko, J. E. Bernhard and S. A. Bass,arXiv:1706.03666 [hep-ph].[40] S. Pratt, E. Sangaline, P. Sorensen, and H. Wang, Phys.Rev. Lett. , 202301 (2015), arXiv:1501.04042.[41] S. Cao, G. Y. Qin and S. A. Bass, Phys. Rev. C , no.2, 024907 (2015),arXiv:1505.01413[42] M. Cacciari, M. Greco, and P. Nason, JHEP , 007(1998), arXiv:hep-ph/9803400.[43] M. Cacciari, S. Frixione, and P. Nason, JHEP , 006(2001), arXiv:hep-ph/0102134.[44] K. J. Eskola, H. Paukkunen, and C. A. Salgado, JHEP , 065 (2009), arXiv:0902.4154.[45] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys.Rev. C92 , 011901 (2015), arXiv:1412.4708.[46] R. Rapp and H. van Hees, Heavy Quarks in the Quark-Gluon Plasma, in
Quark-gluon plasma 4 , pp. 111–206,2010, arXiv:0903.1096.[47] G. D. Moore and D. Teaney, Phys. Rev.
C71 , 064904(2005), hep-ph/0412346.[48] X.-f. Guo and X.-N. Wang, Phys. Rev. Lett. , 3591(2000), arXiv:hep-ph/0005044.[49] B.-W. Zhang, E. Wang, and X.-N. Wang, Phys. Rev.Lett. , 072301 (2004), arXiv:nucl-th/0309040.[50] H. Song and U. W. Heinz, Phys. Rev. C77 , 064901 (2008), arXiv:0712.3715.[51] C. Shen et al. , (2014), arXiv:1409.8164.[52] U. W. Heinz and J. Liu, Nucl. Phys.
A956 , 549 (2016),arXiv:1512.08276.[53] F. Cooper and G. Frye, Phys. Rev.
D10 , 186 (1974).[54] Z. Qiu,
Event-by-event Hydrodynamic Simulations forRelativistic Heavy-ion Collisions , PhD thesis, Ohio StateU., 2013, arXiv:1308.2182.[55] S. A. Bass et al. , Prog. Part. Nucl. Phys. , 225 (1998),nucl-th/9803035.[56] M. Bleicher et al. , J. Phys. G25 , 1859 (1999), hep-ph/9909407.[57] P. B. Gossiaux, J. Aichelin, T. Gousset and V. Guiho, J.Phys. G , 094019 (2010), arXiv:1001.4166.[58] A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C ,1671 (1998), nucl-ex/9805001.[59] A. Bilandzic, R. Snellings and S. Voloshin, Phys. Rev. C , 044913 (2011), arXiv:1010.0233 [nucl-ex].[60] R. Cutler and D. W. Sivers, Phys. Rev. D , 196 (1978).[61] B. L. Combridge, Nucl. Phys. B , 429 (1979).[62] R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigneand D. Schiff, Nucl. Phys. B , 265 (1997), hep-ph/9608322.[63] P. B. Gossiaux, J. Aichelin, T. Gousset and V. Guiho, J.Phys. G , 094019 (2010), arXiv:1001.4166.[64] S. K. Das, F. Scardina, S. Plumari and V. Greco, Phys.Lett. B , 260 (2015), arXiv:1502.03757.[65] S. Caron-Huot and G. D. Moore, Phys. Rev. Lett. ,052301 (2008),arXiv:0708.4232 [hep-ph].[66] R. Rapp and H. van Hees, arXiv:0803.0901 [hep-ph].[67] M. D. Morris and T. J. Mitchell, Journal of StatisticalPlanning and Inference , 381 (1995).[68] STAR, G. Xie, Nucl. Phys. A956 , 473 (2016),arXiv:1601.00695.[69] STAR, L. Adamczyk et al. , Phys. Rev. Lett. , 212301(2017), arXiv:1701.06060.[70] ALICE, J. Adam et al. , JHEP , 205 (2015),arXiv:1506.06604, [Addendum: JHEP06,032(2017)].[71] ALICE, B. B. Abelev et al. , Phys. Rev. C90 , 034904(2014), arXiv:1405.2001.[72] ALICE, A. M. Barbano et al. , Quark Matter 2017[https://indico.cern.ch/event/433345/contributions/2358516][73] CMS, C. Collaboration, CMS-PAS-HIN-16-007 (2016).[74] C. E. Rasmussen and C. K. I. Williams,
Gaussian Pro-cesses for Machine Learning (MIT Press, Cambridge,MA, 2006).[75] S. Ambikasaran, D. Foreman-Mackey, L. Greengard,D. W. Hogg, and M. O’Neil, ArXiv e-prints (2014),arXiv:1403.6015.[76] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Good-man, Publ. Astron. Soc. Pac. , 306 (2013),arXiv:1202.3665.[77] JET, K. M. Burke et al. , Phys.Rev.
C90 , 014909 (2014),arXiv:1312.5003.[78] M. Nahrgang, J. Aichelin, S. Bass, P. B. Gossiauxand K. Werner, Phys. Rev. C , no. 1, 014904(2015),arXiv:1410.5396[79] STAR, L. He et al.et al.