Langevin dynamics of heavy quarks in a soft-hard factorized approach

Ultrarelativistic heavy-ion collisions provide a unique opportunity to create and investigate the properties of strongly interacting matter in extreme conditions of temperature and energy density, where the normal matter turns into a new form of nuclear matter, consisting of deconfined quarks and gluons, namely quark-gluon plasma (QGP Bzdak et al. (2020)). Such collisions allow us to study the properties of the produced hot and dense partonic medium, which are important for our understanding of the properties of the universe in the first few milliseconds and the composition of the inner core of neutron stars E. Shuryak (1980); Braun-Munzinger and Wambach (2009); Braun-Munzinger et al. (2016). Over the past two decades, the measurements with heavy-ion collisions have been carried at the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron Collider (LHC) at CERN Gyulassy and L. McLerran (2005); E. Shuryak (2005); B. Muller, J. Schukraft, and B. Wyslouch (2012), to search and explore the fundamental properties of QGP, notably its transport coefficients related to the medium interaction of hard probes.

Heavy quarks (HQs), including charm and bottom, provide a unique insight into the microscopic properties of QGP Zhou et al. (2014); Tang et al. (2014); Andronic et al. (2016); X. Dong and V. Greco (2019); Dong et al. (2019); Zhao et al. (2020). Due to large mass, they are mainly produced in initial hard scatterings and then traverse the QGP and experience elastic and inelastic scatterings with its thermalized constituents Rapp and van Hees (2010); F. Prino and R. Rapp (2016). The transport properties of HQ inside QGP are encoded in the HQ transport coefficients, which are expected to affect the distributions of the corresponding open heavy-flavor hadrons. The resulting experimental observables like the nuclear modification factor $R_{\rm AA}$ and elliptic anisotropy $v_{\rm 2}$ of various $D$ and $B$-mesons are therefore sensitive to the HQ transport coefficients, in particular their energy and temperature dependence. There now exist an extensive set of such measurements, which allow a data-based extraction of these coefficients. In this work, we make such an attempt by using a soft-hard factorized model (see Sec. III) to calculate the diffusion and drag coefficients relevant for heavy quark Langevin dynamics.

A particularly important feature of the QGP transport coefficients is their momentum and temperature dependence, especially how they change within the temperature region accessed by the RHIC and LHC experiments. For instance the normalized jet transport coefficient $\hat{q}/T^{3}$ was predicted to present a rapidly increasing behavior with decreasing temperature and develop a near-$T_{c}$ peak structure J. F. Liao and E. Shuryak (2009). The subsequent studies  J. C. Xu, J. F. Liao and M. Gyulassy (2015, 2016); JET Collaboration (2014); S. Z. Shi, J. F. Liao, and M. Gyulassy (2018, 2019); Ramamurti and Shuryak (2018) seem to confirm this scenario. Another important transport property, shear viscosity over entropy density ratio $\eta/s$, also presents a visible $T$-dependence with a considerable increase above $T_{c}$ J. E. Bernhard, J. S. Moreland and S. A. Bass (2019). Concerning the HQ diffusion and drag coefficients, there are also indications of nontrivial temperature dependence from both the phenomenological extractions S. S. Cao, G. Y. Qin, and S. A. Bass (2015); Cao et al. (2016); M. Greif, F. Senzel, H. Kremer, K. Zhou, C. Greiner, and Z. Xu (2017); S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018); S. Li and C. W. Wang (2018); Prado et al. (2020); Wang et al. (2020) and theoretical calculations M. He, R. J. Fries, and R. Rapp (2013); T. Song, H. Berrehrah, D. Cabrera, W. Cassing, and E. Bratkovskaya (2016); W. M. Alberico, A. Beraudo, A. De Pace, A. Molinari, M. Monteno, M. Nardi, and F. Prino (2011); O. Fochler, Z. Xu, and C. Greiner (2010); S. Z. Shi, J. F. Liao, and M. Gyulassy (2018, 2019). The difference among the derived hybrid models is mainly induced by the treatment of the scale of QCD strong coupling, hadronization and non-perturbative effects Gossiaux (2019); Cao (2021). See Ref. Beraudo et al. (2018); Xu et al. (2019); Cao et al. (2019) for the recent comparisons.

The paper is organized as follows. In Sec. II we introduce the general setup of the employed Langevin dynamics. Section III is dedicated to the detailed calculation of the heavy quark transport coefficients with the factorization model. In Sec. IV we show systematic comparisons between modeling results and data and optimize the model parameters based on global $\chi^{2}$ analysis. With the parameter-optimized model, the energy and temperature dependence of the heavy quark transport coefficients are presented in Sec. V, as well as the comparisons with data and other theoretical calculations. Section VI contains the summary and discussion.

The classicial Langevin Transport Equation (LTE) of a single HQ reads Rapp and van Hees (2010)

with $i,k=1,2,3$. The first term on the right hand side of Eq. 1b represents the deterministic drag force, $F^{i}_{drag}=-\eta_{D}p^{i}$, which is given by the drag coefficient $\eta_{D}(E,T)$ with the HQ energy $E=\sqrt{\vec{p}^{\;2}+m^{2}_{Q}}$ and the underlying medium temperature $T$. The second term denotes the stochastic thermal force, $F^{i}_{thermal}=C^{ik}\rho^{k}/\sqrt{dt}$, which is described by the momentum argument of the covariance matrix $C^{ik}$,

together with a Gaussian-normal distributed random variable $\rho^{k}$, resulting in the uncorrelated random momentum kicks between two different time scales

$\kappa_{T}$ and $\kappa_{L}$ are the transverse and longitudinal momentum diffusion coefficients, respectively, which describe the momentum fluctuations in the direction that perpendicular (i.e. transverse) and parallel (i.e. longitudinal) to the propagation. Considering the Einstein relationship, which enforces the drag coefficient starting from the momentum diffusion coefficients as S. Li and J. F. Liao (2020)

The parameter $\xi$ denotes the discretization scheme of the stochastic integral, which typically takes the values $\xi=0,0.5,1$, representing the pre-point Ito, the mid-point Stratonovic, and the post-point discretization schemes, respectively; $d=3$ indicates the spatial dimension. In the framework of LTE the HQ-medium interactions are conveniently encoded into three transport coefficients, i.e. $\eta_{D}$, $\kappa_{T}$ and $\kappa_{L}$ (Eq. 4). All the problems are therefore reduced to the evaluation of $\kappa_{T/L}(E,T)$, which will be mainly discussed in this work.

Finally, we introduce a few detailed setups of the numerical implementation. The space-time evolution of the temperature field and the velocity field are needed to solve LTE (Eq. 1a and 1b). Following our previous analysis S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018); S. Li, C. W. Wang, R. Z. Wan, and J. F. Liao (2019), they are obtained in a 3+1 dimensional viscous hydrodynamic calculation I. Karpenko, P. Huoviven and M. Bleicher (2014), with the local thermalization started at $\tau_{0}=0.6~{}{\rm fm}/{\it c}$, the shear viscosity ${\eta/s=1/(4\pi)}$ and the critical temperature $T_{c}=165~{}{\rm MeV}$ (see details in Ref. I. Karpenko, P. Huoviven and M. Bleicher (2014)). When the medium temperature drops below $T_{c}$, heavy quark will hadronize into the heavy-flavor hadrons via a fragmentation-coalescence approach. The Braaten-like fragmentation functions are employed for both charm and bottom quarks E. Braaten, K. Cheung, and T. C. Yuan (1993); M. Cacciari, P. Nason and R. Vogt (2005). An instantaneous approach is utilized to characterize the coalescence process for the formation of heavy-flavor mesons from the heavy and (anti-)light quark pairs. The relevant coalescence probability is quantified by the overlap integral of the Winger functions for the meson and partons, which are defined through a harmonic oscillator and the Gaussian wave-function C. B. Dover, U. Heinz, E. Schnedermann, and J. Zimányi (1991), respectively. See Ref. S. Li and J. F. Liao (2020) for more details.

When propagating throughout the QGP, the HQ scattering off the gluons and (anti-)partons of the thermal deconfined medium, can be characterized as the two-body elementary processes,

with $p_{1}=(E_{1},\vec{p}_{1})$ and $p_{2}$ are the four-momentum of the injected HQ ($Q$) and the incident medium partons $i=q,g$, respectively, while $p_{3}$ and $p_{4}$ are for the ones after scattering. Note that the medium partons are massless ($m_{2}=m_{4}\sim 0$) in particular comparing with the massive HQ ($m_{1}=m_{3}=m_{Q}$ in a few times $\rm GeV$). The corresponding four-momentum transfer is $(\omega,~{}\vec{q}^{\;})=(\omega,~{}\vec{q}_{T},~{}q_{L})$. The Mandelstam invariants read

with $q\equiv|\vec{q}\;|\ll E_{1}$ for small momentum exchange. The transverse and longitudinal momentum diffusion coefficients can be determined by weighting the differential interaction rate with the squared transverse and longitudinal momentum transfer, respectively. It yields

and

As the momentum transfer vanishes ($|t|\rightarrow 0$), the gluon propagator in the $t$-channel of the elastic process causes an infrared divergence in the squared amplitude $\overline{|\mathcal{M}^{2}}|_{t-channel}\propto 1/t^{2}$, which is usually regulated by a Debye screening mass, i.e. $t\rightarrow t-\lambda m^{2}_{D}$ with an adjustable parameter $\lambda$ B. Svetitsky (1988); P. B. Gossiaux and J. Aichelin (2008). Alternatively, it can be overcome by utilizing a soft-hard factorized approach S. Peign$\acute{\rm e}$ and A. Peshier (2008a, b), which starts with the assumptions that the medium is thermal and weakly coupled, and then the interactions between the heavy quarks and the medium can be computed in thermal perturbation theory. Finally, this approach allows to decompose the soft HQ-medium interactions with $t>t^{\ast}$, from the hard ones with $t<t^{\ast}$. For soft collisions the gluon propagator should be replaced by the hard-thermal loop (HTL) propagator E. Braaten and T. C. Yuan (1991); J. P. Blaizot and E. Iancu (2002), while for hard collisions the hard gluon exchange is considered and the Born approximation is appropriate. Therefore, the final results of $\kappa_{T/L}$ include the contributions from both soft and hard components.

As discussed in the QED case S. Peign$\acute{\rm e}$ and A. Peshier (2008a), $\mu+\gamma\rightarrow\mu+\gamma$, the complete calculation for the energy loss of the energic incident heavy-fermion is independent of the intermediate scale $t^{\ast}$ in high energy limit. However, in the QCD case, there is the complication that the challenge of the validity of the HTL scenario, $m_{D}^{2}\ll T^{2}$ S. Peign$\acute{\rm e}$ and A. Peshier (2008b), due to the temperatures reached at RHIC and LHC energies. Consequently, in the QCD case, the soft-hard approach is in fact not independent of the intermediate scale $t^{\ast}$ P. B. Gossiaux and J. Aichelin (2008); W. M. Alberico, A. Beraudo, A. De Pace, A. Molinari, M. Monteno, M. Nardi, and F. Prino (2011). In this analysis we have checked that the calculations for $\kappa_{T/L}$ are not sensitive to the choice of the artificial cutoff $t^{\ast}\sim m_{D}^{2}$.

In the next parts of this section, we will focus on the energy and temperature dependence of the interaction rate ${\Gamma}$ at leading order in $g$ for the elastic process, as well as the momentum diffusion coefficients (Eq. 7 and 8) in soft and hard collisions, respectively.

Following the strategies utilized in our previous work S. Li and J. F. Liao (2020); S. Li, W. Xiong, and R. Z. Wan (2020), the two key parameters in this study, the intermediate cutoff $t^{\ast}$ and the scale $\mu$ of running coupling (Eq. 31a), are tested within a wide range of possibility and drawn constrains by comparing the relevant charm meson data with model results. We calculate the corresponding final observable y for the desired species of D-meson. Then, a $\chi^{2}$ analysis can be performed by comparing the model predictions with experimental data

In the above $\sigma_{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 in the range $p_{\rm T}\leq 8~{}\rm{GeV}$: $D^{0}$, $D^{+}$, $D^{\ast+}$ and $D_{s}^{+}$ $R_{\rm AA}$ data collected at mid-rapidity ($|y|<0.5$) in the most central ($0-10\%$) and semi-central ($30-50\%$) Pb–Pb collisions at $\sqrt{s_{\rm NN}}=2.76~{}{\rm TeV}$ ALICE Collaboration (2016a, b) and $\sqrt{s_{\rm NN}}=5.02~{}{\rm TeV}$ ALICE Collaboration (2018a), as well as the $v_{\rm 2}$ data in semi-central ($30-50\%$) collisions ALICE Collaboration (2014); ALICE Collaboration (2018b); CMS Collaboration (2018a).

We scan a wide range of valus for ($t^{\ast}$, $\mu$): $1\leq\frac{|t^{\ast}|}{m_{D}^{2}}\leq 3$ and $1\leq\frac{\mu}{\pi T}\leq 3$. A total of 20 different combinations were computed and compared with the experimental data. The obtained results are summarized in Tab. 1. The $\chi^{2}$ values are computed separately for $R_{\rm AA}$ and $v_{\rm 2}$ as well as for all data combined. To better visualize the results, we also show them in Fig. 1, with left panels for $R_{\rm AA}$ analysis and right panels for $v_{2}$ analysis. In both panels, the y-axis labels the desired parameters, $\frac{|t^{\ast}|}{m_{D}^{2}}$ (upper) and $\frac{\mu}{\pi T}$ (lower), and x-axis labels the $\chi^{2}/d.o.f$ within the selected ranges⁵⁵5$\chi^{2}/d.o.f$ is shown in the range $0.5<\chi^{2}/d.o.f<2.5$ for better visualization.. The different points (filled gry circles) represent the different combinations of parameters ($\frac{|t^{\ast}|}{m_{D}^{2}}$, $\frac{\mu}{\pi T}$) in Tab. 1, with the number on top of each point to display the relevant “$Model~{}ID$” for that model. A number of observations can be drawn from the comprehensive model-data comparison. For the $R_{\rm AA}$, several models achieve $\chi^{2}/d.o.f\simeq 1$ with $\frac{|t^{\ast}|}{m_{D}^{2}}\lesssim 1.5$ and widespread values of $\mu$: $1\leq\frac{\mu}{\pi T}\leq 3$. This suggests that $R_{\rm AA}$ appears to be more sensitive to the intermediate cutoff while insensitive to the scale of coupling constant. For the $v_{\rm 2}$, it clearly shows a stronger sensitivity to $\mu$, which seems to give a better description ($\chi^{2}/d.o.f\simeq 1.5-2.0$) of the data with $\frac{\mu}{\pi T}=1$. It is interesting to see that $R_{\rm AA}$ data is more powerful to constrain $\frac{|t^{\ast}|}{m_{D}^{2}}$, while $v_{\rm 2}$ data is more efficient to nail down $\frac{\mu}{\pi T}$. Taken all together, we can identify a particular model that outperforms others in describing both $R_{\rm AA}$ and $v_{\rm 2}$ data simultaneously with $\chi^{2}/d.o.f=1.3$. This one will be the parameter-optimized model in this work: $|t^{\ast}|=1.5m_{D}^{2}$ and $\mu=\pi T$.

In this section we will first examine the $t^{\ast}$ dependence of $\kappa_{T/L}$ (Eq. 15) for both charm and bottom quark. Then, the relevant energy and temperature dependence of $\kappa_{T/L}(E,T)$ will be discussed with the optimized parameters. For the desired observables we will perform the comparisons with the results from lattice QCD at zero momentum limit, as well as the ones from experimental data in the low to intermediate $p_{\rm T}$ region.

In this work we have used a soft-hard factorized model to investigate the heavy quark momentum diffusion coefficients $\kappa_{T/L}$ in the quark-gluon plasma in a data-driven approach. In particular we’ve examined the validity of this scenario by systematically scanning a wide range of possibilities. The global $\chi^{2}$ analysis using an extensive set of LHC data on charm meson $R_{\rm AA}$ and $v_{\rm 2}$ has allowed us to constrain the preferred range of the two parameters: soft-hard intermediate cutoff $|t^{\ast}|=1.5m_{D}^{2}$ and the scale of QCD coupling constant $\mu=\pi T$. It is found that $\kappa_{T/L}$ have a mild sensitivity to $t^{\ast}$, supporting the validity of the soft-hard approach when the coupling is not small. With this factorization model, we have calculated the transport coefficient $\hat{q}$, drag coefficient $\eta_{D}$, spatial diffusion coefficient $2\pi TD_{s}$, and then compared with other theoretical calculations and phenomenological extractions. Our analysis suggests that a small value $2\pi TD_{s}\simeq 6$ appears to be much preferred near $T_{c}$. Finally we’ve demonstrated a simultaneous description of charm meson $R_{\rm AA}$ and $v_{\rm 2}$ observables in the range $p_{\rm T}\leq 8$ GeV. We’ve further made predictions for bottom meson observables in the same model.

We end with discussions on a few important caveats in the present study that call for future studies:

• $\bullet$

  In the framework of Langevin approach, the heavy flavor dynamics is encoded into three coefficient, $\kappa_{T}$, $\kappa_{L}$ and $\eta_{D}$, satisfying Einstein’s relationship. It means that two of them are independent while the third one can be obtained accordingly. The final results therefore depend on the arbitrary choice of which of the two coefficients are calculated with the employed model F. Prino and R. Rapp (2016); Xu et al. (2019). For consistency, we calculate independently the momentum diffusion coefficients $\kappa_{T/L}$ with the factorization approach, and obtain the drag coefficient via $\eta_{D}=\eta_{D}(\kappa_{T},\kappa_{L})$ (see Eq. 4). A systematic study among the different options may help remedy this situation.

• $\bullet$

  According to the present modeling, the resulting spatial diffusion coefficient exhibits a relatively strong increase of temperature comparing with lattice QCD calculations, in particular at large $T$, as shown in Tab. 2. It corresponds to a larger relaxation time and thus weaker HQ-medium coupling strength. This comparison can be improved by (1) determining the key parameters based on the upcoming measurements on high precision observables (such as $R_{\rm AA}$, $v_{\rm 2}$ and $v_{\rm 3}$) of both $D$ and $B$-mesons at low $p_{\rm T}$; (2) replacing the current $\chi^{2}$ analysis with a state-of-the-art deep learning technology.

  	$T=T_{c}$	$T=2T_{c}$	$T=3T_{c}$
  This work	$\sim$6	$\sim$22	$\sim$31
  LQCD (median value) D. Banerjee, S. Datta, R. Gavai, and P. Majumdar (2012); N. Brambilla, V. Leino, P. Petreczky, and A. Vairo (2020)	$\sim$5	$\sim$10	$\sim$13

• $\bullet$

  It is realized Cao et al. (2020) that the heavy quark hadrochemistry, the abundance of various heavy flavor hadrons, provides special sensitivity to the heavy-light coalescence mechanism and thus plays an important role to understand the observables like the baryon production and the baryon-to-meson ratio. A systematic comparison including the charmed baryons over a broad momentum region is therefore crucial for a better constraining of the model parameters, as well as a better extraction of the heavy quark transport coefficients in a model-to-data approach.

• $\bullet$

  The elastic scattering ($2\rightarrow 2$) processes between heavy quark and QGP constituents are dominated for heavy quark with low to moderate transverse momentum S. Li, C. W. Wang, R. Z. Wan, and J. F. Liao (2019). Thus, here we consider only the elastic energy loss mechanisms to study the observables at $p_{\rm T}\leq 8$ GeV. The missing radiative ($2\leftrightarrow 3$) effects may help to reduce the discrepancy with $R_{\rm AA}$ data in the vicinity of $p_{\rm T}=8$ GeV, as mentioned above (see Fig. 8). With the soft-hard factorized approach, it would be interesting to include both elastic and radiative contributions in a simultaneous best fit to data in the whole $p_{\rm T}$ region. We also plan to explore this idea in the future.