Data-driven extraction of heavy quark diffusion in quark-gluon plasma

Heavy quark production provides a unique probe of the quark-gluon plasma transport properties in heavy ion collisions. Experimental observables like the nuclear modification factor $R_{\rm AA}$ and elliptic anisotropy $v_{2}$ of heavy flavor mesons are sensitive to the heavy quark diffusion coefficient. There now exist an extensive set of such measurements, which allow a data-driven extraction of this coefficient. In this work, we make such an attempt within our recently developed heavy quark transport modeling framework (Langevin-transport with Gluon Radiation, LGR). A question of particular interest is the temperature dependence of the diffusion coefficient, for which we test a wide range of possibility and draw constraints by comparing relevant charm meson data with model results. We find that a relatively strong increase of diffusion coefficient from crossover temperature $T_c$ toward high temperature is preferred by data. We also make predictions for Bottom meson observables for further experimental tests.


I. INTRODUCTION
At extremely high temperatures such as those available at the earliest moments of cosmic evolution, normal matter turns into a new form of deconfined nuclear matter known as a quark-gluon plasma (QGP). Such a state of matter once filled the early universe when the temperature was high enough. Today the QGP is recreated in laboratories by high energy nuclear collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). A lot of measurements have been performed at RHIC and the LHC, allowing the use of empirical data to extract key properties of the QGP, which are of fundamental interests. The transport properties (such as the shear and bulk viscosity, jet transport coefficient, etc) have been found to be particularly informative for unraveling the dynamical features of QGP, leading to its identification as the nearly perfect fluid [1][2][3].
Heavy quark production provides a unique probe of the quark-gluon plasma in heavy ion collisions [4][5][6][7][8][9]. Charm and bottom quarks are very hard to be thermally produced in QGP and are dominantly produced from the initial hard scatterings. These rare objects then propagate through the QGP fireball and encode the medium information during their dynamical evolution. Experimental observables like the nuclear modification factor R AA and elliptic anisotropy v 2 of (for heavy flavor mesons as well as heavy flavor decay leptons) are sensitive to the heavy quark diffusion coefficient inside the QGP. There now exist an extensive set of such measurements, which allow a data-driven extraction/constraint of this coefficient. In this work, we make such an attempt within our recently * lish@ctgu.edu.cn † liaoji@indiana.edu developed heavy quark transport modeling framework, Langevin-transport with Gluon Radiation (LGR) [10].
A particularly interesting question about the QGP transport properties is their temperature dependence, especially how they change in the temperature region from transition temperature T c to a few times T c . This is the region accessible through RHIC and LHC collision experiments. It has been suggested that such temperature dependence could be highly nontrivial, especially close to T c . For example, it was proposed long ago that the jetmedium interaction strength (quantified by e.g. normalized jet transport coefficientq/T 3 ) may rapidly increase from high temperature down toward T c and develop a near-T c peak structure [11]. Such a scenario appears to be confirmed by many subsequent studies [12][13][14][15][16][17]. Another important transport property, shear viscosity over entropy density ratio η/s, also seems to have a visible Tdependence with a considerable increase from T c toward higher temperature []. Regarding the diffusion and drag coefficients relevant for heavy quark dynamics, there are also indications of nontrivial temperature dependence []. In this work we will focus on the diffusion coefficient and test a wide range of possibility for its temperature dependence. By comparing modeling results with experimental data of charm hadrons, we draw constraints on the behavior of this important transport property of QGP. Based on that, we further make predictions for bottom hadron observables.
The rest of this paper is structured as follows. In Section 2, we introduce the detailed setup of the LGR modeling framework and discuss the temperature dependence of diffusion constant. In Section 3, we systematically compare modeling results with data and extract optimal range of this transport coefficient based on global χ 2 analysis. Direct comparison of optimized model results with experimental observables as well as predictions for new measurements are presented in Section 4. Finally we summarize this study in Section 5.

II. METHODOLOGY
In this Section we present the details of our modeling framework. The heavy quark evolution is described by the following Langevin transport equation that incorporates gluon radiation: [18] where the deterministic drag force reads with η D ( p, T ) being the drag coefficient . The two-point temporal correlation of the stochastic thermal force F T is given by [19] indicating the uncorrelated random momentum kicks from the medium partons. P ij = p i p j /p 2 and P ij ⊥ = δ ij − p i p j /p 2 are the projection operators for momentum components parallel and perpendicular to the direction of the HQ motion, respectively. The relation between the drag coefficient (η D ), the longitudinal (κ ) and transverse momentum diffusion coefficients (κ ⊥ ) follows (4) The third term on the right hand side of Eq. 1, denotes the total recoil force induced by the emitted gluons. The emission rate of gluons is estimated with the following Higher-Twist model formula [20]: In the above, z denotes the fraction of energy carried away by the emitted gluon, and P (z) represents the quark splitting function; α s (k ⊥ ) is the strong coupling constant of QCD at leading order approximation; is the gluon formation time;q q is the quark jet transport coefficient. The recoil force becomes important when the heavy quark is in the high energy regime E m Q , where HQ velocity v Q = 1 − (m Q /E) 2 ∼ 1 and radiative energy loss becomes significant. In this regime, theq q in the above can be approximated bŷ As one can see at this point, the key parameters controlling all the forces in Eq. 1 are the momentum diffusion coefficients (κ and κ ⊥ ). It is customary in heavy quark phenomenological modelings [18,[21][22][23][24]] to further connect these parameters to the spatial diffusion constant under reasonable approximations. In the low momentum regime where diffusion dynamics is most important, one could assume approximate isotropy for the momentum diffusion coefficients, i.e. κ = κ ⊥ ≡ κ. Thus the Eq. 4 is further reduced to i.e. the so-called dissipation-fluctuation relation in the non-relativistic approximation. The connection to the (scaled) spacial diffusion constant, strictly speaking, is valid at zero-momentum limit [25], 2πT . Such relation has been phenomenologically generalized to finite momentum and widely used in heavy quark modelings [18,23,24], allowing the expressions of both the drag and the momentum diffusion coefficients (Eq. 8) in terms of spacial diffusion constant: We can see that now there is only one key transport parameter, the spatial diffusion constant (2πT D s ) that quantifies all relevant components: the drag force (Eq. 2), thermal random force (Eq. 3) and the recoil force (Eq. 5 and 6) in the Langevin approach (Eq. 1). Thus, the dynamical interactions between the heavy quarks and the QGP medium are conveniently encoded into 2πT D s . We note that, in a naive way, a small/large spatial diffusion corresponds to a short/long mean-free path and thus strong/weak HQ-medium coupling strength. Indeed, many past studies have demonstrated sensitivity of experimental observables (such as R AA and v 2 of heavy flavor mesons) to this key parameter. It appears that, very similar to the situation of jet energy loss, the R AA is mainly controlled by the HQ-medium interaction on average while the v 2 is strongly influenced by the temperature dependence of the diffusion constant [22,[26][27][28]. In this study, we aim to investigate such temperature dependence. Let us focus on the temperature range (1 ∼ 3)T c and frame the question in a model-independent way. Consider D s (T /T c ) as an arbitrary function in this range, it can always be expressed via a series of polynomials (as long as one can include enough terms): .. without the need of assuming any theoretically-motivated temperature dependence. In principle, with sufficient experimental data and adequate computing power, one could exploit data-driven extraction of all these coefficients term by term. As a first step, we take only the first two terms, i.e. a constant plus a linear dependence, with the following ansatz: The two dimensionless parameters in Eq. 10, the slope α and the intercept β, will be explored in a very wide range without presuming any reasonable value. Our approach is to compute observables for any given (α, β) and let the large set of experimental data decide what would be preferred via χ 2 analysis. We will then compare so-extracted spatial diffusion constant with various other results [29]. Finally we describe a few detailed aspects of the numerical implementation. When solving the Langevin transport equation (Eq. 1), we need the space-time evolution of the medium temperature and the fireball velocity field. It is simulated in terms of a 3+1 dimensional relativistic viscous hydrodynamics based on the HLLE algorithm ( -see details in [30]). Concerning the hadronization of HQ, a "dual" approach, including both fragmentation and heavy-light coalescence mechanisms, is utilized when the local temperature is below T c = 165 GeV. Following our previous work [24,31], the Braaten fragmentation functions [32] is employed with the parameter r = 0.1 [33]. Within the instantaneous coalescence approach, the momentum distributions of heavy-flavor mesons (M ) composed of a heavy quark (Q) and a light anti-quark (q) reads M represents the coalescence probability for Qq combination to form the heavy-flavor meson in the n th excited state, and it is defined as the overlap integral of the Wigner functions for the meson and Qq pair [34] are the relative coordinate and the relative momentum, respectively, in the center-of-mass frame of Qq pair. Note that the parton Wigner functions are defined through the Gaussian wave-function, while for heavyflavor meson, it is quantified by a harmonic oscillator one [35]. The width parameter σ M is expressed as [24] where, K = 2/3 (K = 2/5) for the ground state n = 0 (1 st excited state n = 1); r 2 M is the mean-square charge radius of a given species of D-meson, which is predicted by the light-front quark model [36]; e i and m i are the charge and mass of a given parton, respectively. See Ref. [24] for more details.

III. CONSTRAINING DIFFUSION CONSTANT
In this Section we focus on constraining diffusion constant using experimental data. With any given set of parameters (α, β) in Eq. 10, we can calculate the corresponding final observable y for the desired species of D-meson. Then, a χ 2 analysis can be performed by comparing the model predictions with experimental data In the above σ i is the total uncertainty in data points, including the statistic and systematic components which are added in quadrature. n = N −1 denotes the degree of freedom (d.o.f ) when there are N data points used in the comparison. In this study, we use an extensive set of LHC data: D 0 , D + , D * + and D + s collected at mid-rapidity (|y| < 0.5) in the most central (0−10%) and semi-central (30 − 50%) Pb-Pb collisions at √ s NN = 2.76 TeV [37,38] and √ s NN = 5.02 TeV [39], as well as the v 2 data in semicentral (30 − 50%) collisions [40,41]. We scan a wide range of values for (α, β) in Eq. 10: 0 ≤ α ≤ 9 and −8.5 ≤ β ≤ 4. We note this covers a significantly broader span than existing studies and than commonly conceived reasonable values of (2πT )D s . It would be highly unlikely, if not impossible, for the actual QGP diffusion constant to fall outside this range. A total of 25 different combinations were computed and compared with experimental data, and we summarize these in Tab. I. The χ 2 values were computed separately for R AA and v 2 as well as for all data combined. To better visualize the results, we also show them in Fig. 1, with panel (a) for R AA analysis and panel (b) for v 2 analysis. In both panels, the y-axis labels the slope α and x-axis labels the 2πT D s (T = 2T c ): basically y-axis quantifies a model's temperature dependence while x-axis calibrates the average diffusion in that model. The different points (filled circles) represent the different combinations of parameters (α, β) in Tab. I, with the number near each point to display the relevant χ 2 /d.o.f for that model. A number of observations can be drawn from the comprehensive model-data comparison. For the R AA , several models achieve χ 2 /d.o.f ∼ 1 with widespread values of slope parameter. This suggests that R AA appears to be more sensitive to the average diffusion constant while insensitive to the temperature dependence. For the v 2 , it clearly shows a stronger sensitivity to the temperature dependence. There also exist several models with χ 2 /d.o.f ∼ 1 and it appears that a small value of (2πT )D s near T c is crucial for a better description of v 2 data. Taken all together, we are able to identify two particular models that outperform others in describing both R AA and v 2 data simultaneously with χ 2 /d.o.f ∼ 1. These two will be the parameter-optimized models from the LGR framework: the LGR (Model-A) with (α, β) = (3, −1) and the LGR (Model-B) with (α, β) = (6.5, −5.5). While both models give similarly nice good description of R AA , the Model-B has a much stronger temperature dependence and gives a better description of v 2 . Let us make a comparison with various existing modeling frameworks, e.g. TAMU [22], PHSD [26], LTB [42], POWLANG [43], BAMPS (eastic) [44], BAMPS (elas-tic+radiative) [44] and CUJET3 [45,46]. The published results from these models for relevant observables are taken from pertinent references and used to evaluate the corresponding χ 2 for each model. The analysis results χ 2 /d.o.f (R AA ) and χ 2 /d.o.f (v 2 ) are then shown and compared in the panel (c) of Fig. 1. One can see that to describe simultaneously both R AA and v 2 data is challenging in general. The LGR Model-A and Model-B, featuring a moderate to strong temperature dependence and a small diffusion constant very close to T c , demonstrate a successful description of current measurements.  [47], red triangle [48] and blue square [49]), CUJET3(red region [50]), a Bayesian analysis in 95% CL from LIDO (shadowed gray band [51]) and a LO calculation with a Boltzmann dynamics (long dashed blue curve [52]). The result for D-meson (dot dashed green curve [53]) in the hadronic phase is also shown for comparison. Figure 2 we present the spacial diffusion constant 2πT D s of charm quark from various phenomenological extractions and theoretical calculations. Very close to T c , the results from both LGR Model-A (solid black curve) and Model-B (dashed orange curve) are compatible with the LQCD calculations within their significant uncertainties (pink circle [47], red triangle [48] and blue square [49]) as well as consistent with other models [52,54,55]. Toward the higher temperature end, the region spanned by our Model-A and Model-B compare well with the band from the Bayesian analysis based on the Duke model (gray region [51]). Combing various information together, we observe that: (1) a small value 2πT D s (2 ∼ 4) appears to be much preferred in the vicinity of T c ; (2) a relatively strong increase of its value toward higher temperature is favored, albeit still with large uncertainty for T 2T c .

IV. LGR RESULTS FOR OBSERVABLES
In this Section, we present the results for various observables to be compared with experimental data. Specifically we use the optimized LGR Model-B based on analysis from the previous Section.  are done with FONLL initial charm quark spectra and EPS09 NLO parametrization for the nPDF in Pb [24], and the green band reflects the theoretical uncertainties coming from these inputs. It can be seen that the model calculations provide a very good description of the measured p T -dependent R AA data for various charm mesons. The same conclusion can be drawn for the comparison in Pb-Pb collisions at √ s NN = 5.02 TeV, as shown in   p T ∼ 3 − 5 GeV, suggests that charm quarks actively participate in the collective expansion of the fireball. Given that our model has provided a very good description of charm meson data, it is tempting to further test it with bottom meson measurements. Here we present LGR Model-B results for the strange and non-strange bottom mesons. To do that, one would need the relevant transport coefficient for bottom quarks. It has been suggested [43,52,57] that the ratio of bottom quark spacial diffusion constant to that for charm quark exhibits a weak T -dependence and varies within ∼ 0.8 − 0.9 in the range T c < T < 4T c . We therefore use a constant factor 0.85 to give a temperature dependent spatial diffusion constant 2πT D s (bottom) = 0.85 × 2πT D s (charm) for calculating the nuclear modification factor of openbottom hadrons. Figure 6 shows the obtained results in the 30 − 50% centrality Pb-Pb collisions at √ s NN = 5.02 TeV. It is found that R AA (B 0 s ) is significant larger than R AA (B + ), in particular at p T ∼ 4 − 6 GeV. This difference decreases toward high p T . Similar to previous results for the open-charm systems [31], the enhancement behavior is mainly induced by the heavy-light coalescence effect, which is more pronounced for the B 0 s (bs) than for the B + (bu). The observation is consistent with the Bmeson measurements (0 − 100%) reported by the CMS Collaboration [58].

V. SUMMARY
In summary we have used a recently developed heavy quark transport modeling framework (Langevintransport with Gluon Radiation, LGR) to study the heavy flavor spatial diffusion constant in the quark-gluon plasma in a data-driven approach. In particular we've examined the temperature dependence of this transport co-efficient by systematically scanning a wide range of possibilities. Our global χ 2 analysis using extensive set of LHC data on charm meson R AA and v 2 has allowed us to constrain the preferred range of this parameter. It is found that R AA is more sensitive to the average value in the relevant temperature region while v 2 is more sensitive to the temperature dependence. Taken together, our analysis suggests that a small value 2πT D s (2 ∼ 4) appears to be much preferred in the vicinity of T c while a relatively strong increase of its value toward higher temperature is favored. The extracted temperaturedependent 2πT D s curve is shown in Fig. 2 and consistent with other phenomenological analyses as well as lattice calculations. With the optimized LGR model calculations we've demonstrated a simultaneous description of charm meson R AA and v 2 observables. We've further made predictions for bottom meson observables in the same model, for which an enhancement of the ratio R AA (B 0 s )/R AA (B + ) is found in the low to intermediate p T region with its maximum around p T ∼ 4 − 6 GeV.
We end with a discussion on further extending the present analysis for future studies. An important step further is to go beyond the linear ansatz Eq. 10 for the spatial diffusion constant. We plan to further employ a multi-term nonlinear ansatz and to use Bayesian inference for efficiently extracting the full temperature dependence. Another improvement would be the inclusion of RHIC data in the analysis. Compared with LHC, the RHIC fireball shall be more sensitive to the near-Tc region and would help further constrain the transport coefficient there. It has also been noticed that the major uncertainty in nailing down the temperature dependence lies in the high temperature end, which has not been well constrained by current data. One possibility is to explore the extremely central collisions (e.g. top 1% events of multiplicity) at the highest LHC energies, which should produce a fireball that has more fraction of its spacetime evolution in the high temperature region and thus becomes more sensitive to the behavior of QGP in that region. We also plan to explore this idea in the future.
versity Pervasive Technology Institute, and in part by the Indiana METACyt Initiative. The Indiana META-Cyt Initiative at IU was also supported in part by Lilly Endowment, Inc.