A data-drive analysis for heavy quark diffusion coefficient

We apply a Bayesian model-to-data analysis on an improved Langevin framework to estimate the temperature and momentum dependence of the heavy quark diffusion coefficient in the quark-gluon plasma (QGP). The spatial diffusion coefficient is found to have a minimum around 1-3 near Tc in the zero momentum limit, and has a non-trivial momentum dependence. With the estimated diffusion coefficient, our improved Langevin model is able to simultaneously describe the D-meson RAA and v2 in three different systems at RHIC and the LHC.


Introduction
Over the past few decades, facilities at RHIC and the LHC have been exploring the strongly coupled QGP medium by conducting heavy-ion collision experiments. Since the QGP is not directly observable, the study of its properties relies on the comparison between experimental measurements and theoretical modeling. Open heavy flavor, as one of the hard probes of the QGP medium, is valuable yet difficult to study, since it involves measuring 'rare process', as well as modeling the complex evolution of heavy quarks during different stages of the collision. Nevertheless, significant progress has been made in recent years: a number of theoretical models are able to describe some of the heavy quark observables, and provide qualitative agreement on the heavy quark diffusion coefficient [1]. However, the simultaneous description of different heavy quark observables in different collision systems still poses challenges for some of the models. Moreover, the current studies of the heavy quark often rely on 'eye-ball fit' or 'χ 2 -goodness fit' comparisons between the model calculations and the experimental measurements, leading to very limited distinguishability between the performance of different models, and possible bias on the estimate of the heavy quark diffusion coefficient.
A rigorous and systematic approach to a model-to-data comparison requires optimizing the model and simultaneously describing different experimental observables (which are often highly interdependent and have correlated uncertainties). A Markov chain Monte Carlo (MCMC) algorithm and Bayesian statistics can be used to optimize the process, and perform the comparison and constrain the parameters in a data-driven manner. A Bayesian analysis has already been applied to the study of soft medium properties with great success, for example, the determination of the temperature dependence of shear and bulk viscosities η/s(T ), ζ/s(T ) [2], or the constraint of the QCD matter equation of state purely from experimental measurements [3]. In this work, we shall extend this type of analysis to the heavy quark sector.

Heavy quark in-medium evolution
Our study of the full space-time evolution of heavy quarks utilizes the well-established improved Langevin framework developed by the Duke QCD group [4]: The initial entropy density and the heavy quark positions are generated by an effective parametric initial condition model T R ENTomimicking the scaling behavior of IP-Glasma [5]. The heavy quark initial transverse momentum distribution is provided by FONLL [6]; After being produced, heavy quarks propagate in the QGP medium according to an improved Langevin equation [4], which takes into account both collisional and radiative energy loss; at the critical temperature T c = 154 MeV, heavy quarks hadronize into heavy mesons via a hybrid model of recombination and fragmentation. After hadronization, the heavy mesons continuously experience hadronic re-scattering with light hadrons -all hadronic interactions are modeled by UrQMD [8].
The improved Langevin equation that describes the dynamics of heavy quarks inside the QGP medium is: where −η D (p)⃗ p and ⃗ ξ are the drag and thermal random forces. They are responsible for the heavy quark collisional energy loss and related to the transport coefficientq via η The third term ⃗ f g = −d⃗ p g /dt is the recoil force when the heavy quark emits bremsstrahlung gluons, with ⃗ p g being the gluon momentum. We adopt the higher twist formalism for the medium-induced gluon radiation [7] and the recoil force is dependent onq. With sufficient large gluon energy limit, the system will equilibrate when the Einstein relationship is satisfied [4]. Under this construction, the drag force, the thermal random force and the recoil force are dependent on the heavy quark transport coefficientq. In literature, the spatial diffusion coefficient D s = 4T 2 /q is often used to characterize the interacting strength between heavy quarks and the medium.
At high temperature or large momentum, the heavy quark transport coefficient can be calculated using the perturbative QCD approach (pQCD). However, it has been found that the pQCD calculation is not sufficient to explain the suppression experienced by heavy quarks when propagating in the medium, nor able to simultaneously describe both the heavy flavor R AA and v 2 . In order to introduce a non-perturbative effect, we parameterize the spatial diffusion coefficient as a combination of linear temperature dependent component and a pQCD component: The soft component (D s 2πT ) soft = α * (1 + β * (T/T c − 1)) is the heavy quark diffusion coefficient in the p = 0 GeV/c limit. It accounts for the non-perturbative effects. The (D s 2πT ) pQCD = 8πT 3 /q pQCD is calculated by summing over the leading-order 2 → 2 scattering between heavy quarks (M c = 1.3 GeV) and massless light partons at fixed coupling α s = 0.3. The value of heavy quark diffusion coefficient is determined by three parameters x = (α, β, γ), which can be interpreted as following: Parameter α is D s 2πT (p = 0) around T c , parameter β is its slope above T c . Parameter γ controls the ratio between the linear component and the pQCD component. A small value of γ indicates the non-perturbative effects are non-negligible towards the high momentum region, and a large value of γ indicates a quick conversion to the pQCD dominated region. The evolution of heavy quarks depends on the surrounding matter. In this study, the QGP medium is simulated by an event-by-event (2+1)-dimensional viscous hydrodynamical model VISHNEW [9], which also includes shear and bulk viscous corrections. All parameters related to the medium evolution (such as η/s, ζ/s) as well as the initial entropy deposition have been calibrated to experimental data of charged particles by an independent Bayesian analysis [2].
According to Bayes' theorem, the posterior probability distribution of the model parameters (P(x|y exp )) is proportional to the product of the parameters' prior distribution (P(x)), and the likelihood of observing experimental results given parameter (L(y exp |x)). In this case, we assume the prior parameter have a uniform distribution , which implies that we don't have any initial constraint/knowledge on the parameters. Thus the posterior distribution of the parameters is proportional to the likelihood: where y exp is the experimental data (D-meson R AA , v 2 at different centralities at AuAu collisions (200 GeV) and PbPb collisions (2.76 and 5.02 TeV)) [10], y is the model calculations at parameter x = (α, β, γ), Σ is the covariance matrix which includes quantified uncertainties (experimental and theoretical).
In practice, we explore the parameter space and sample the Bayesian posterior distribution by performing a MCMC random walk, which requires evaluating the model at thousands or even millions of parameter values. However, the simulation of heavy quark evolution in heavy-ion collision is computationally highly expensive. Therefore, efficient techniques such as Gaussian process emulators are applied to solve the problem -these act as a fast surrogate model of the full Langevin simulation. A detailed description of the construction of the emulators and performance of the MCMC algorithm can be found in Ref. [2]. To verify our analysis, 200 sets of parameters are randomly selected from the posterior distributions and evaluated by the GP emulators. Figure 1 compares a selection of observables with experimental data. Two different colors correspond to two independent Bayesian analysis: the blue lines correspond to the results obtained by calibrating on a single collision system -PbPb collision at 5.02 TeV, the red lines correspond to the results obtained by calibrating on the three collision systems at RHIC and the LHC. After the calibration, our improved Langevin approach is capable of simultaneously describing the experimental data. The spreads of the posterior R AA and v 2 visualize the uncertainties kept in our analysis, which come from the experimental statistic and uncorrelated systematic uncertainties, as well as the uncertainty arising from the prediction made by the GP emulators.

Posterior results and discussion
The posterior distributions of parameters (α, β, γ) are plotted in Fig. 2, where the diagonal histograms are the marginal distribution of each parameter with the others integrated out. The off- diagonal plots show the correlation between each pair. We observe a narrow distribution for parameters α and γ, implying those two are constraint with current experimental data and uncertainties. Although parameter β is poorly constraint, the strong negative correlation between α and β indicates some collinear interdependency and should not be ignored.
In Fig. 2 we also estimate the diffusion coefficient by taking the posterior samples from the combined analysis (with 90% credibility). Our estimate of D s 2πT (p = 0) shows a minimum between (1-3) around T c , and has a positive temperature dependent. Comparing our estimate with the coefficient used/calculated by a number of other models in the market [1,11], as well as the lattice QCD calculations [12], our estimate is consistent with lattice QCD calculations within the large uncertainties current inherent in lattice QCD calculations. Although the diffusion coefficient used in different models are rather different, they all have a minimum near T c with value range D s 2πT (p = 0) ∼ (1−7). Nevertheless, we should note that in order to make an apple-to-apple comparison among different models, it is essentially important for heavy quarks to propagate in the same QGP medium, especially for those transport models, which is not the case for the current comparison. A detailed investigation is needed, but beyond the scope of this work.