A data-driven analysis of the heavy quark transport coefficient

Using a Bayesian model-to-data analysis, we estimate the temperature dependence of the heavy quark diffusion coefficients by calibrating to the experimental data of $D$-meson $R_{\mathrm{AA}}$ and $v_2$ in AuAu collisions ($\sqrt{s_{NN}}=200$ GeV) and PbPb collisions ($\sqrt{s_{NN}}=2.76$ TeV)~\cite{Xie:2016iwq}. The spatial diffusion coefficient $D_s2\pi T$ is found to be mostly constraint around $(1.3-1.5) T_c$ and is compatible with lattice QCD calculations. We demonstrate the capability of our improved Langevin model to simultaneously describe the $R_{\mathrm{AA}}$ and $v_2$ at both RHIC and the LHC energies, as well as the feasibility to apply a Bayesian analysis to quantitatively study the heavy flavor transport in heavy-ion collisions.


Introduction
Over the past few years, significant progress has been made to precisely describe the yields and flow of particles emitted from the QGP, and to quantitatively estimate its transport properties such as the shear viscosity to entropy density ratio η/s [2,3,4]. High energetic jet and heavy quark energy loss mechanisms and their related transport coefficients (q, D s ), however, are still not yet understood on a quantitative level. While the heavy quark transport coefficients are not directly measurable, they can be encapsulated in the parameters of theoretical models. By comparing the model calculations to the experimental data, we can estimate their values and therefore understand the interaction mechanism between heavy quarks and the medium.
In general, a comparison between the model calculations and the experimental data relies on optimizing multiple parameters -some are related to the properties of the system, some are ad hoc parameters of different models. Until now, most heavy flavor studies vary the parameters manually until the agreement with the experimental data is achieved, which leads to very limited usefulness regarding comparison with experimental data. Moreover, the simultaneous description of various observables -D-meson R AA and v 2 -is still an inevitable challenge for most models. A systematic and complete approach to optimizing the models would be to perform a random walk across the parameter space and calibrate to experimental data by applying a Bayesian analysis. In such an analysis the result is the posterior possibility distribution of each parameter, from which we can extract the optimal values of transport coefficients, the uncertainty and the correlation among different parameters.
The Bayesian statistical analysis has been applied with great success to the determination of multiple QGP properties, such as the precise estimate of the shear and bulk viscosities and the constraint of the equation of state above the parton-hardon transition [3]. In this study, we extend the analysis to quantitatively study heavy flavor evolution in heavy-ion collisions.

Modeling heavy flavor evolution in heavy-ion collisions
Our analysis utilizes the well-established framework developed by the Duke QCD group to describe the full space-time evolution of heavy quarks in heavy-ion collisions: the initial entropy density of the medium as well as the heavy quark position are generated by an effective parametric initial condition model T R ENTo, which has been shown to mimic the scaling behavior of the EKRT and IP-Glasma models [5]. The initial transverse momentum distribution is provided by FONLL [6].
The propagation of the heavy quarks in the medium is described by an improved Langevin equation [7], which takes into account of both collisional and radiative energy loss: The first two terms on the right hand side of the equation are the drag and thermal random forces inherited from the standard Langevin equation. They are responsible for the heavy quark collisional energy loss and related to the momentum diffusion coefficientq via The third term f g = −d p g /dt is the recoil force from the bremsstrahlung gluon emitted from the heavy quarks. It is added in order to take the radiative energy loss into consideration. We adopt the higher twist results for the gluon emission probability [8], which is proportional to the momentum diffusion coefficient q. In literature [7], the spatial diffusion coefficient D s = 4T 2 /q is often used to characterize the interacting strength between the heavy quarks and the medium.
The evolution of the QGP medium is described by an event-by-event (2+1)-dimensional viscous hydrodynamical model VISH(2+1), which includes both shear and bulk viscous corrections [9]. It should be noted that all the parameters related to the bulk initialization and evolution have been calibrated to experimental data of charged particles [3].
The hadronization of heavy quarks into heavy mesons is described by a hybrid fragmentation and recombination model. Currently we neglect any rescattering of the heavy mesons in the hadron gas state.
In order to estimate the heavy quark diffusion coefficient, we parametrize the spatial diffusion coefficient as linearly dependent on temperature and assume two model parameters x = (D min , D slope ): By varying the model parameters x, we are able to change the temperature dependence of D s 2πT . In this work we focus on charm quarks evolution, but the same framework would apply to bottom quarks as well.

Model-to-data comparison
To calibrate the model to the data, i.e. to determine the optimal values of all the parameters, requires a random walk through the parameter space using the Markov chain Monte Carlo (MCMC) algorithm [10]. However, evaluating the full Langevin model in a fine grid in parameter space is computationally highly expensive, therefore directly performing the MCMC random walk with the Langevin model is not feasible. In this situation Gaussian process emulators can be used as an alternate surrogate model to fast interpolate the output of the Langevin model at any arbitrary point in the parameter space. In practice, a small number of parameter sets (( x 1 , x 2 , ..., x n ), i = 1, ..., 80) are uniformly sampled via a Latin hypercube algorithm across the parameter space [11] and evaluated in the Langevin model. The model outputs y = (R AA , v 2 ) are then calculated at each of the parameter set. Those prior results (( x i , y i ), i = 1, .., 80) are used to train the Gaussian process emulators such that we can calibrate our parameters to experimental data through the MCMC random walk. The posterior possibility distribution P( x * |X, Y, y exp ) for the true parameter x * is calculated according to Bayes' theorem [4]: where P(X, Y, y exp | x * ) is the likelihood of observing (X, Y, y exp ) with x * , and is given by [12]:

Posterior Results
To verify our analysis, 200 input parameter values are randomly selected from the posterior distribution. In Fig. 1 we compare our model outputs predicted from the Gaussian process emulator with the experimental data. Three individual analyses are performed in this work, and each is labeled with different color: the blue one calibrating our parameters to the experimental data in AuAu collisions, the green one calibrating to the experimental data in PbPb collisions, and the red one calibrating to both. The spread of the posterior outputs visualizes the remaining uncertainty in our analysis. In addition, the median value of each parameter x = (D min , D slope ) is applied in the full Langevin simulation and the corresponding results are shown as the dashed lines. Overall the results demonstrate good agreement between the calibration and experimental values, as well as the validity of the Gaussian process emulator to predict the output from a physical model.
The main result of the analysis is the posterior possibility distribution of each parameter, which is shown in Fig. 2. The diagonal is the marginal distribution of the parameter x = (D min , D slope ) with the other one integrated out, and the off-diagonal shows the correlation between them. We find a narrow distribution for the parameter D min for all the three analysis. From the off-diagonal histograms, we observe a nontrivial negative correlation between D min and D slope . However, D slope is relatively unconstrained, indicating a possible missing piece in our parametrization or too large experimental data uncertainties. Figure 3 presents our estimate of the charm quark diffusion coefficient D s 2πT as a function of temperature, compared to other model calculations [13]. The posterior estimate for D s 2πT is significantly constrained compared to the prior range. In addition, the uncertainty of the diffusion coefficient is smallest  between (1.3 − 1.5)T c -a temperature range where charm quarks spend most of the time. This feature of D s 2πT is consistent with the Bayesian analysis on the temperature dependence of η/s, indicating that this is the temperature range where the medium spends most of the time. Though in general charm quarks tend to have a larger D s 2πT in Pb+Pb collisions than in Au+Au collisions, the combined analysis still gives an estimate within the overlapped region of those two. As expected, the combined analysis is the most constrained. Our result is compatible with the lattice QCD calculations.

Conclusion
In this study, we have performed a Bayesian model-to-data analysis to extract the temperature dependence of the charm quark diffusion coefficient in the QGP medium. The spatial diffusion coefficient D s 2πT is found to be compatible to lattice QCD calculation, with the most constrained region between 1.3T c to 1.5T c . This study has demonstrated the applicability of Bayesian analysis to the heavy flavor analysis in heavy-ion physics. Future work will improve the model description of the heavy quark evolution, such as applying Boltzmann transport, including hardonic scattering, and extend to more collision systems.