Probing the transport properties of Quark-Gluon Plasma via heavy-flavor Boltzmann and Langevin dynamics

TL;DR

This study compares Boltzmann and Langevin dynamics in modeling heavy quark transport in quark-gluon plasma, revealing differences in transport coefficients and their impact on observables like energy loss, R_AA, and v_2, with implications for experimental data interpretation.

Contribution

It provides a systematic comparison of Boltzmann and Langevin approaches for heavy quark transport, highlighting their differences and compatibility with experimental data.

Findings

Boltzmann coefficients are larger at high T and energy.

Langevin dynamics predict greater energy loss and smaller R_AA.

Boltzmann approach yields larger v_2 and azimuthal broadening.

Abstract

The heavy quark propagation behavior inside the quark-gluon plasma (QGP), is usually described in terms of the Boltzmann dynamics, which can be reduced to the Langevin approach by assuming a small momentum transfer for the scattering processes between heavy quarks and the QGP constituents. In this work, the temperature and energy dependence of the transport coefficients are calculated in the framework of both Boltzmann and Langevin dynamics. The derived transport coefficients are found to be systematically larger in the Boltzmann approach as compared with the Langevin, in particular in the high temperature and high energy region. Within each of the two theoretical frameworks, we simulate the charm quark production and the subsequent evolution processes in relativistic heavy-ion collisions. We find that the total in-medium energy loss is larger from the Langevin dynamics, resulting in a smaller (larger) $R_{\rm AA}$ at high (low) $p_{\rm T}$, for both the charm quark and heavy-flavor mesons. Meanwhile, the Boltzmann model is found to induce larger $v_{\rm 2}$, in particular at moderate $p_{\rm T}$, as well as stronger broadening behavior for the azimuthal distributions. By comparing the model calculations with available experimental measurements for D-mesons, we find that the Langevin approach is more favored by the $R_{\rm AA}$ data while the Boltzmann approach is more favored favor by the $v_{\rm 2}$ data. A simultaneous description of both observables appear challenging for both models.

Figures

Equations

$$R_{\mathrm{AA}}(p_{\mathrm{T}})=\frac{\mathrm{d}{\sigma }_{\mathrm{AA}}/\mathrm{d}{p}_{\mathrm{T}}}{\mathrm{d}{\sigma }_{\mathrm{pp}}/\mathrm{d}{p}_{\mathrm{T}}}.$$

$$v_{\rm 2}(p_{\rm T})=\biggr{\langle}\frac{p_{x}^{2}-p_{y}^{2}}{p_{x}^{2}+p_{y}^{2}}\biggr{\rangle},$$

$$\frac{p_{\mathrm{Q}}}{E_{\mathrm{Q}}}\cdot \partial f_{\mathrm{Q}}=C[f_{\mathrm{Q}}]$$

$$\frac{\Delta P_l}{\Delta t}={\Gamma }_l(E_{\mathrm{Q}},T,t)=\frac{g}{\nu }\frac{(2\pi {)}^3\delta }{\delta f_{\mathrm{Q}}}\int d\Phi (n,m)\prod _{\{in\}}{f}_i\overline{{|M|}_l^2}$$

$$d\Phi (n,m)=(2\pi {)}^4{\delta }^{(4)}(p_{\mathrm{in}}-p_{\mathrm{out}})\prod _{\{in,out\}}\frac{d^3{\vec{p}}_i}{2E_i(2\pi {)}^3}$$

$$\Delta P_{total}=\sum _l(\Delta P_l)=\sum _l({\Gamma }_l\cdot \Delta t).$$

$${\eta }_{\mathrm{D}}\equiv -\frac{d<p>}{dt}/<p>$$

$${\kappa }_{\mathrm{L}}\equiv \frac{d<(\Delta p_{\mathrm{z}}{)}^2>}{dt}$$

$${\kappa }_{\mathrm{T}}\equiv \frac{1}{2}\frac{d<(\Delta p_{\mathrm{T}}{)}^2>}{dt},$$

$$\frac{d{x}^{\mathrm{i}}}{dt}=\frac{p^{\mathrm{i}}}{E^{\mathrm{i}}}$$

$$\frac{d{p}^{\mathrm{i}}}{dt}=F_{\mathrm{Drag}}^{\mathrm{i}}+F_{\mathrm{Diff}}^{\mathrm{i}}$$

$$F_{\mathrm{Drag}}^{\mathrm{i}}=-{\eta }_{\mathrm{D}}(\vec{p},T)\cdot p^{\mathrm{i}},$$

$$F_{\mathrm{Diff}}^{\mathrm{i}}=\frac{1}{\sqrt{dt}}{C}^{\mathrm{ij}}(t,\vec{p}+\xi d\vec{p},T){\rho }^j$$

$$P(\vec{\rho })=(\frac{1}{2\pi }{)}^{3/2}exp\left\{-\frac{{\vec{\rho }}^2}{2}\right\},$$

$$C^{\mathrm{ij}}(\vec{p},T)\equiv \sqrt{{\kappa }_{\mathrm{L}}(\vec{p},T)}{\hat{p}}^i{\hat{p}}^j+\sqrt{{\kappa }_{\mathrm{T}}(\vec{p},T)}({\delta }^{ij}-{\hat{p}}^i{\hat{p}}^j),$$

$${\eta }_{\mathrm{D}}=$$

$$\frac{d-1}{2p^2}\left[\xi (\sqrt{{\kappa }_{\mathrm{T}}}+\sqrt{{\kappa }_{\mathrm{L}}}{)}^2-(3\xi -1){\kappa }_{\mathrm{T}}-(\xi +1){\kappa }_{\mathrm{L}}\right],$$

$${\eta }_{\mathrm{D}}=\frac{{\kappa }_{\mathrm{L}}}{2TE}-\frac{d-1}{2p^2}(\sqrt{{\kappa }_{\mathrm{T}}}-\sqrt{{\kappa }_{\mathrm{L}}}{)}^2$$

$${\eta }_{\mathrm{D}}=\frac{\kappa }{2TE},$$

$$2\pi T{D}_s=\underset{p\to 0}{\lim }\frac{2\pi T^2}{m_{\mathrm{Q}}\cdot {\eta }_{\mathrm{D}}(p)}=\underset{E\to m_{\mathrm{Q}}}{\lim }\frac{2\pi T^2}{m_{\mathrm{Q}}\cdot {\eta }_{\mathrm{D}}(E)},$$

$${\eta }_{\mathrm{D}}=\frac{1}{2\pi T{D}_s}\cdot \frac{2\pi T^2}{E}$$

$$\kappa =\frac{1}{2\pi T{D}_s}\cdot 4\pi T^3.$$

$${\tau }_{charm}\equiv \underset{p\to 0}{\lim }{{\eta }_{\mathrm{D}}(p)}^{-1}=\frac{m_{charm}}{2\pi T_c^2}\cdot \frac{(2\pi T{D}_s)}{(T/T_c{)}^2},$$

$$x(t+\Delta t)-x(t)=\frac{p(t)}{E_p(t)}\Delta t$$

$$\frac{d{N}_{\mathrm{M}}}{d^3{\vec{p}}_{\mathrm{M}}}=$$

$$\times {\overline{W}}_{\mathrm{M}}^{(\mathrm{n})}({\vec{y}}_{\mathrm{M}},{\vec{k}}_{\mathrm{M}}){\delta }^{(3)}({\vec{p}}_{\mathrm{M}}-{\vec{p}}_{\mathrm{Q}}-{\vec{p}}_{\bar{\mathrm{q}}})$$

$${\overline{W}}_{\mathrm{M}}^{(\mathrm{n})}({\vec{y}}_{\mathrm{M}},{\vec{k}}_{\mathrm{M}})$$

$$=\int \frac{d^3{\vec{x}}_{\mathrm{Q}}^{\prime }{d}^3{\vec{p}}_{\mathrm{Q}}^{\prime }}{(2\pi {)}^3}\frac{d^3{\vec{x}}_{\bar{\mathrm{q}}}^{\prime }{d}^3{\vec{p}}_{\bar{\mathrm{q}}}^{\prime }}{(2\pi {)}^3}{W}_{\mathrm{Q}}({\vec{x}}_{\mathrm{Q}}^{\prime },{\vec{p}}_{\mathrm{Q}}^{\prime })W_{\bar{\mathrm{q}}}({\vec{x}}_{\bar{\mathrm{q}}}^{\prime },{\vec{p}}_{\bar{\mathrm{q}}}^{\prime })$$

$$\times W_{\mathrm{M}}^{(\mathrm{n})}({\vec{y}}_{\mathrm{M}}^{\prime },{\vec{k}}_{\mathrm{M}}^{\prime })$$

$$\displaystyle=\biggr{[}\frac{1}{2}\biggr{(}\frac{\vec{y}^{\;2}_{\rm M}}{\sigma_{\rm M}^{2}}+\sigma_{\rm M}^{2}\vec{k}^{\;2}_{\rm M}\biggr{)}\biggr{]}^{n}exp\biggr{\{}-\frac{1}{2}\biggr{(}\frac{\vec{y}^{\;2}_{\rm M}}{\sigma_{\rm M}^{2}}+\sigma_{\rm M}^{2}\vec{k}^{\;2}_{\rm M}\biggr{)}\biggr{\}}\biggr{/}n!$$

$${\sigma }_{\mathrm{M}}^2$$

$$\displaystyle\frac{d^{3}\vec{k}}{2k}\rightarrow\frac{d^{3}\vec{k}}{2k}\cdot\biggr{\{}2[1-cos\biggr{(}\frac{t-t_{0}}{\tau_{f}}\biggr{)}]\biggr{\}},$$

$${\tau }_f=\frac{2x(1-x)E}{k_{\perp }^2+(x{m}_{\mathrm{Q}}{)}^2+(1-x)m_{\mathrm{D}}^2/2},$$

$$x=\frac{k+k_z}{E+p_z}$$

$$\frac{d{p}^{\mathrm{i}}}{dt}=F_{\mathrm{Drag}}^{\mathrm{i}}+F_{\mathrm{Diff}}^{\mathrm{i}}+F_{\mathrm{Gluon}}^{\mathrm{i}}.$$

Full text

I INTRODUCTION

In ultrarelativistic collisions of heavy nuclei such as Au or Pb, an extreme high temperature and energy density environment can be produced around the collision point, where allows form a new state of nuclear matter consisting of the deconfined quarks and gluons, namely the quark-gluon plasma, QGP Gyulassy and L. McLerran (2005); E. Shuryak (2005). To investigate its properties, the experiments using the Au and Pb as the colliding beams have been carried at the Relativistic Heavy Ion Collider (RHIC) at BNL and at the Large Hadron Collider (LHC) at CERN, respectively, in the past two decades B. Muller, J. Schukraft, and B. Wyslouch (2012); J. Schukraft (2017); E. Shuryak (2017). The QGP was found to induce the jet quenching, as well as to exhibit the collective flow behavior among various probes PHENIX Collaboration (2004, 2005); STAR Collaboration (2005); ALICE Collaboration (2013a, b). The jet quenching phenomenon is known X. N. Wang and M. Gyulassy (1992) as the energy loss of the fast partons traversing the QGP medium, and it can be investigated by measuring the suppression behavior of the cross-section of the desired particles produced in nucleus-nucleus collisions to that in binary-scaled nucleon-nucleon collisions at the same energy, which is the so-called nuclear modification factor, $R_{\rm AA}$,

[TABLE]

The collective effect can be interpolated as the strong collective expansion of QGP when its (local) thermal equilibrium state is achieved, and it can be studied by a Fourier expansion S. Voloshin and Y. Zhang (1996); A. M. Poskanzer and S. A. Voloshin (1998) of the particle azimuthal distributions with respect to the reaction plane, which is defined as the plane including impact parameter and beam axis. Normally, the second coefficient, $v_{\rm 2}$,

[TABLE]

is called elliptic flow coefficient, which allows to describe the anisotropy of the transverse momentum.

Heavy quark (HQ), including charm and bottom, are of particular interest E. Braaten and M. H. Thoma (1991); Rapp and van Hees (2010); A. Andronic et al (2016); G. Aarts et al (2017); V. Greco (2017) since, due to their large mass, (1) $m_{\rm Q}\gg\Lambda_{\rm QCD}$, thus, its initial production can be well described by the perturbative Quantum ChromoDynamics (pQCD) at the next-to-leading order M. Cacciari, M. Greco, and P. Nason (1998); M. Cacciari, S. Frixione, and P. Nason (2001); M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason, and G. Ridolfi (2012), in particular at high $p_{\rm T}$; (2) $m_{\rm Q}\gg T$, resulting in the negligible thermal production of HQ pairs in QGP medium with the temperature reached at RHIC and LHC energies. In addition, HQ flavor is conserved throughout the interactions with the surrounding QGP constituents, i.e. gluons and (anti-)light quarks. Therefore, the initially produced HQ pairs will experience the full evolution of QGP, and serve as its ideal probes.

During the propagating through the QGP medium, the HQ dynamics is usually described by the Boltzmann or Langevin model R. Rapp and H. van. Hess (2008); F. Prino and R. Rapp (2016). For the Boltzmann approach, the evolution of the HQ distribution function behaves the Boltzmann Transport Equation (BTE), where the elastic and inelastic scattering processes between HQs and the quasi-particles of QGP are quantified by the relevant scattering matrix. Consequently, it can be given with the help of the perturbative QCD. Due to large HQ mass and moderate medium temperature, the typical momentum transfers in interactions, $q\sim gT$, are assumed small, $gT\ll m_{\rm Q}$ W. M. Alberico, A. Beraudo, A. De Pace, A. Molinari, M. Monteno, M. Nardi, and F. Prino (2011), therefore, the HQ trajectory will be changed significantly only after receiving lots of soft momentum kicks from the surrounding QGP constituents, resulting in the Brownian motion. Based on this assumption, BTE is reduced to the Fokker-Plank Transport Equation (FPTE), which can be realized stochastically by a Langevin Transport Equation (LTE). In the framework of LTE, all the interactions are conveniently encoded into three transport coefficients, satisfying the dissipation-fluctuation relation. Therefore, with LTE, all the problem reduced to the evaluation of three transport coefficients, which can be extracted from the lattice QCD at zero momentum limit.

Many models were developed from the Boltzmann O. Fochler, Z. Xu, and C. Greiner (2010); S. S. Cao, T. Luo, G. Y. Qin,and X. N. Wang (2016); M. Djordjevic and M. Djordjevic (2016); J. C. Xu, J. F. Liao and M. Gyulassy (2016); T. Song, H. Berrehrah, D. Cabrera, W. Cassing, and E. Bratkovskaya (2016) and Langevin dynamics Y. Akamatsu, T. Hastuda and T. Hirano (2009); S. K. Das, J. Alam and P. Mohanty (2010); M. He, R. J. Fries, and R. Rapp (2013); S. S. Cao, G. Y. Qin, and S. A. Bass (2015); T. Lang, H. van Hees, G. Inghirami, J. Steinheimer, and M. Bleicher (2016) to study the suppression and collective effect of the final heavy-flavor productions (having the charm or bottom quarks among these valence quarks) such as D mesons ($D^{0}$, $D^{+}$, $D^{*+}$ and $D^{+}_{s}$ S. Li (2014a) (for the ALICE Collaboration); S. Li (2014b) (for the ALICE Collaboration)) and B mesons ($B^{0}$, $B^{+}$ and $B_{\rm s}$ CMS Collaboration (2017)). Comparing the theoretical calculations with available data, it was realized J. F. Liao and E. Shuryak (2009a); S. K. Das, F. Scardina, S. Plumari, and V. Greco (2015); J. C. Xu, J. F. Liao and M. Gyulassy (2015); S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018) that the simultaneous description of $R_{\rm AA}$ and $v_{\rm 2}$ of open charmed meson at low and intermediate $p_{\rm T}$ is sensitive to the temperature and energy dependence of the transport coefficients. It is necessary to mention that, in order to improve the description of the measurements, the Duke group Y. R. Xu, J. E. Bernhard, S. A. Bass, M. Nahrgang, and S. S. Cao (2018) develops a data-based hybrid model to extract the transport coefficient by utilizing the Bayesian model-to-data analysis. See Refs. R. Rapp $et.~{}al.$ (2018); S. S. Cao $et.~{}al.$ (2018); Y. R. Xu $et.~{}al.$ (2018) for the recent review.

As mentioned, the Langevin approach is a very convenient and widely used model, and it allows to establish, directly, a link between the observables and transport coefficients, which can be extracted from the lattice QCD calculations. However, the condition $m_{\rm Q}\gg gT$ may not always be justified, in particular for charm quark with the medium temperature close to its initial value, resulting in the possible modification of the heavy meson $R_{\rm AA}$. So, in this work, we focus on the discussion related to the “benefits and limitations for Boltzmann vs. Langevin implementations of the heavy-flavor transport in an evolving medium” R. Rapp $et.~{}al.$ (2018). Both the BTE and LTE will be employed to investigate the temperature and energy dependence of the various transport coefficients, as well as to study the charm quark transport behaviors in the QGP medium.

The paper is organized as follows. In Sec. II we summarize the employed Boltzmann and Langevin dynamics, together with the comparison for the extracted transport coefficients including only the elastic interactions. Sec. III is dedicated to the description of the hybrid model, including the initial state configuration, the hydrodynamic expansion of the underlying medium, heavy quark propagation and hadronization via fragmentation and “heavy-light” coalescence mechanisms. Sec. IV shows the results obtained at parton and hadron levels with only the elastic processes, while Sec. V with both the elastic and inelastic contributions. Sec. VI contains the summary and discussion.

II Boltzmann and Langevin Dynamics with only elastic processes

II.1 Linearized Boltzmann transport model

The Boltzmann Transport Equation (BTE) reads

[TABLE]

where, $p_{\rm Q}$, $E_{\rm Q}$ and $f_{\rm Q}$ are the HQ 4-momentum, energy and distribution function, respectively. $C[f_{\rm Q}]$ denotes the collision integral, including all the interaction mechanisms between heavy quarks and the medium partons. Equation 3 can be linearized by ignoring the change of thermal parton distribution in the medium due to the heavy quark propagation, and thus, $C[f_{\rm Q}]$ becomes a linear function of $f_{\rm Q}$. Based on the Monte Carlo techniques, Eq. 3 can be solved numerically by slicing the coordinate space into a 3-dimension grid, and then the test particle method C. Y. Wong (1982) is used to sample $f_{\rm Q}$ in each cell. The collision integral is solved by using the stochastic algorithm for evaluating the collision probability Z. Xu and C. Greiner (2005); F Scardina, M Colonna,S Plumari,and V Greco (2013). In this work, we utilize only the linearized Boltzmann module in the Lido hybrid model W. Y. Ke, Y. R. Xu, and S. A. Bass (2018) with all the default parameters, except the charm quark mass $m_{\rm c}=1.5~{}{\rm GeV}$.

In the local rest frame (LRF) of the cell, the heavy quark transport is performed within a given time-step $\Delta t$. Concerning a desired scattering process $l$, there are $n$ ($m$) incoming (outgoing) partons, and the reaction probability $\Delta P_{l}$ is expressed as W. Y. Ke, Y. R. Xu, and S. A. Bass (2018)

[TABLE]

where, $\Gamma_{l}(E_{\rm Q},T,t)$ is the relevant scattering rate; $g$ is the spin-color degeneracy factor of the incoming medium partons; $\nu$ is the statistical factor that corrects for double-counting when there are identical particles in the initial/final state; $f_{i}$ denotes the heavy quark ($i=Q$) and medium parton ($i=\bar{q},q,g$) density, while the latter one follows the $\rm Maxwell-J{\ddot{u}}ttner$ distribution; $\overline{{|M|}_{l}^{2}}$ is the initial state spin-color averaged scattering matrix element squared for two-body interactions, which can be calculated via the perturbative QCD at leading-order B. L. Combridge (1978). $d\Phi(n,m)$ in Eq.4 is the $n+m$ body phase space integration,

[TABLE]

where, $p_{\rm in}$ ($p_{\rm out}$) indicates the total 4-momentum of all the incoming (outgoing) partons for a given $2\rightarrow 2$ scattering process $l$. Within the time interval $\Delta t$, the total reaction probability $\Delta P_{total}$ is given by

[TABLE]

It was argued Guy D. Moore and D. Teaney (2005) that the interactions between HQs and the medium partons can be encoded into the drag and momentum diffusion coefficients:

[TABLE]

which describes the average momentum/energy loss, momentum fluctuations in the direction that parallel (i.e. longitudinal) and perpendicular (i.e. transverse) to the propagation, respectively.

II.2 Langevin transport model

While traversing the Quark-Gluon Plasma (QGP), HQ suffers frequent but soft momentum kicks from the medium partons, therefore, HQ behaves the Brownian motion, which can be described by the Langevin Transport Equation (LTE) M. He, H. van Hees, P. B. Gossiaux, R. J. Fries, and R. Rapp (2013)

[TABLE]

The deterministic drag force reads

[TABLE]

where $\eta_{\rm D}(\vec{p},T)$ is the drag coefficient.

The stochastic force which acts on the HQ is expressed as

[TABLE]

with the Gaussian noise $\rho^{j}$ follows a normal distribution

[TABLE]

resulting in $<\rho^{i}>_{\rho}=0$ and $<\rho^{i}\rho^{j}>_{\rho}=\delta^{ij}$. Therefore, there is no correlation for the random force between two different time scales $<F^{\rm i}_{\rm Diff}(t)F^{\rm j}_{\rm Diff}(t^{\prime})>_{\rho}\equiv C^{\rm ik}C^{\rm kj}\delta(t-t^{\prime})$, indicating the uncorrelated random momentum kicks from the medium partons. During the numerical implementation, as shown in Eq. 10, the stochastic process depends on the specific choice of the momentum argument of the covariance matrix, $C^{\rm ij}(t,\vec{p}+\xi d\vec{p},T)$, via a parameter $\xi\in[0,1]$. Typically, $\xi=0$ for pre-point Ito, $\xi=1/2$ for mid-point and $\xi=1$ for post-point discretization scheme of the stochastic integral. Finally, $C^{\rm ij}$ can be represented in terms of the longitudinal ($\kappa_{\rm L}$) and transverse momentum diffusion coefficients ($\kappa_{\rm T}$) A. Beraudo, A. De Pace, W. M. Alberico, and A. Molinari (2009), i.e.

[TABLE]

therefore, the relation between $\eta_{\rm D}$, $\kappa_{\rm L}$ and $\kappa_{\rm T}$ is given by

[TABLE]

where $d=3$ denotes the spatial dimension. As pointed, HQ diffusions are conveniently encoded in the three coefficients $\eta_{\rm D}$, $\kappa_{\rm L}$ and $\kappa_{\rm T}$. Note that Eq. 13 can be reduced to

[TABLE]

with the post-point scheme, i.e. $\xi=1$. Following our previous analysis S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018); S. Li and C. W. Wang (2018), a “minimum model” by assuming a isotropic momentum dependence of the diffusion coefficient, $\kappa_{\rm L}=\kappa_{\rm T}\equiv\kappa$, is adopted in this work, although it is just validated at $p=0$ and, they not exactly the same at $p\neq 0$ region from the analytical calculations A. Beraduo, A. De Pace, W. M. Alberico, and A. Molinari (2009). Eq. 14 is therefore further reduced to

[TABLE]

which is the so-called dissipation-fluctuation relation (or Einstein relation) in the non-relativistic approximation.

II.3 Boltzmann vs. Langevin

In this sub-section, we mainly focus on the comparison of the transport coefficients obtained via the Boltzmann and Langevin approaches with considering only the elastic scattering ($2\rightarrow 2$) off the QGP constituents. We show before that the scattering rate (Eq. 4) for $c+q\rightarrow c+q$ process in Fig. 1, which is presented as a function of charm quark energy and the medium temperature.

It is found that the energy dependence is weak, while the temperature dependence is stronger.

II.3.1 Boltzmann vs. Langevin: spatial diffusion coefficient

The spatial diffusion coefficient $D_{s}$ Guy D. Moore and D. Teaney (2005) scaled by the thermal wavelength $1/(2\pi T)$,

[TABLE]

is defined at $p\rightarrow 0$ limit, which can be obtained directly by substituting Eq. LABEL:eq:BTECoef with the Boltzmann approach. $2\pi TD_{s}$ is available from the Lattice QCD calculation, moreover, it is found S. Li and C. W. Wang (2018) that, according to a phenomenological fitting analysis with the Langevin approach, model predictions based on $2\pi TD_{s}=7$ allow to reproduce all the measured $p_{\rm T}$ dependence of the nuclear modification factor at both RHIC and LHC energies. Therefore, in the Langevin approach, the drag and the momentum diffusion coefficients (Eq. 15) can be obtained via Eq. 16 by setting $2\pi TD_{s}=7$. Note that, in this case, (1) the definition of spatial diffusion coefficient is extended to larger momentum values. Similar strategy is adopted in Ref. Y. Akamatsu, T. Hatsuda, and T. Hirano (2009); S. S. Cao, G. Y. Qin, and S. A. Bass (2015); M. Greif, F. Senzel, H. Kremer, K. Zhou, C. Greiner, and Z. Xu (2017); and (2) the drag and momentum diffusion coefficients in Eq. 15 can be represented in terms of $2\pi TD_{s}$ as

[TABLE]

The temperature dependence of the spatial diffusion coefficient $2\pi TD_{s}$ is presented in Fig. 2. The results from the Boltzmann (only $c+q\rightarrow c+q$ and $c+g\rightarrow c+g$) and Langevin (only collisional) approaches are displayed as the dashed red and solid black curves, respectively. Lattice QCD and Ads/CFT calculations, i.e. Banerjee (pink circles D. Banerjee, S. Datta, R. Gavai, and P. Majumdar (2012)), Kaczmarek (blue square O. Kaczmarek (2014)), Ding (red triangles H. T. Ding, A. Francis, O. Kaczmarek, F. Karsch, H. Satz, and W. Soeldner (2012)) and Oleg (dotted blue curve O. Andreev (2018)) are shown as well as for comparison. Within the significant systematic uncertainties, both Boltzmann (dashed red curve) and Langevin predictions (solid black curve) are consistent with the Banerjee and Oleg calculations. Similar behavior can be found by comparing with the other model predictions, such as the CUJET3 (red region) J. C. Xu, J. F. Liao and M. Gyulassy (2016) and Lido (green region) W. Y. Ke, Y. R. Xu, and S. A. Bass (2018).

With Eq. 16, the thermalization time of charm quark, defined in $p\rightarrow 0$ limit Guy D. Moore and D. Teaney (2005), can be expressed as

[TABLE]

which is about 2.27 and 3.07 ${\rm fm/{\it c}}$ for Boltzmann and Langevin approach, respectively, with $T=2T_{c}\approx 330~{}{\rm MeV}$ and $m_{charm}=1.5~{}{\rm GeV}$.

II.3.2 Boltzmann vs. Langevin: transport diffusion coefficients

In Fig. 3, the drag coefficient (left), longitudinal (middle) and transverse momentum diffusion coefficients (right), are presented as a function of the medium temperature (upper) at a given energy $E\approx 10~{}{\rm GeV}$, and as a function of the charm quark energy (bottom) at fixed temperature $T=0.3~{}{\rm GeV}$. The results obtained with the Boltzmann (Eq. LABEL:eq:BTECoef) and Langevin model (Eq. LABEL:eq:LTECoef) are shown as dashed red and solid black curves, respectively, in each panel.

Concerning the drag coefficient $\eta_{\rm D}$ [(a), (d)], the two models show an increasing temperature dependence and a decreasing behavior for the energy. The results with the Boltzmann approach (dashed red curve) is systematically larger than that with the Langevin approach (solid black curve). Both the longitudinal $\kappa_{\rm L}$ [(b), (e)] and transverse momentum diffusion coefficient $\kappa_{\rm T}$ [(c), (f)] increases strongly with increasing the energy and temperature via the Boltzmann approach, while they change slowly, as expected (Eq. LABEL:eq:LTECoef), via the Langevin approach.

III Methodology

In the previous analysis S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018), we construct a theoretical framework to study the charm quark propagation in ultrarelativistic heavy-ion collisions. The general modules of the hybrid model are discussed in the following.

III.1 Initial state configuration

The initialization of the heavy quark pairs is performed in the spatial and momentum space, respectively. In the transverse direction, the initial spatial distribution is sampled according to the initial binary collision density that is modeled by a Glauber-based approach C. Shen, Z. Qiu, H. C. Song, J. Bernhard, S. A. Bass, and U. Heinz (2016), while in the longitudinal direction, it is described by a data-inspired phenomenological function S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018). The initial momentum distribution of $c/\bar{c}$ is predicted by the FONLL calculations M. Cacciari, M. Greco, and P. Nason (1998); M. Cacciari, S. Frixione, and P. Nason (2001); M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason, and G. Ridolfi (2012), assuming a back-to-back azimuthal correlation between $c$ and $\bar{c}$ ($|\Delta\phi^{c\bar{c}}|=\pi$). Fo nucleus-nucleus collisions, e.g. Pb–Pb, the nuclear modification of the parton distribution functions (nPDFs) is taken into account by utilizing the EPS09 NLO parametrization approach K. J. Eskola, H. Paukkunen, and C. A. Salgado (2009).

The above initial state configuration allows providing the relevant entropy density distribution, which will be taken as the input of the subsequent hydrodynamical evolution. All the parameters in this procedure are tuned by the model-to-data comparison S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018).

III.2 Hydrodynamic description

The underlying medium evolution is modeled by a 3+1D relativistic viscous hydrodynamics, vHLLE I. Karpenko, P. Huoviven and M. Bleicher (2014), with the initial time scale $\tau_{0}=0.6~{}{\rm fm}/{\it c}$, shear viscosity ${\eta/s=1/(4\pi)}$ and critical temperature $T_{c}=165~{}{\rm MeV}$ in both Au–Au and Pb–Pb collisions. Note that the hydrodynamic simulation provides the space-time evolution of the temperature and the flow velocity field, which will be used in the HQ Boltzmann and Langevin dynamics.

The QGP medium expands and cools down, and the (local) temperature drops below the critical one $T_{c}$, resulting in the transition from the QGP phase to hadrons gas, namely hadronization. After the transition, the hadron gas can in principle continue to interact inelastically until the chemical freeze-out, subsequently, the hadronic system continues to expand and interact elastically until the kinetic freeze-out. In this work, we neglect the chemical freeze-out procedure and consider, only, the kinetic freeze-out (or freeze-out since now) occurs at $T_{c}=165~{}{\rm MeV}$. An instantaneous approach across a hypersurface of constant temperature, namely isothermal freeze-out, is utilized and modeled by a widely used approach, Cornelius P. Huovine and H. Petersen (2012).

III.3 Heavy quark propagation in medium

We refer to Ref. S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018) for the detailed discussion about the numerical framework of charm quark Langevin evolution, which is coupled with the expanding hydrodynamic medium. For the Boltzmann case, it is quite similar except the procedure to update the charm quark momentum in a discrete time-step. In the following, we show the general strategy for both cases:

1. (1)

sample a given number of HQ pairs at the position and momentum ($x^{\mu},p^{\mu}$), in the laboratory frame (LAB), according to the previous initial phase space configurations ($\tau\sim 0$); 2. (2)

move all the HQs from $\tau\sim 0$ to $\tau_{0}=0.6~{}{\rm fm/{\it c}}$ as free streaming particles, and modify the positions $x^{\mu}$ correspondingly; 3. (3)

search the fluid cell at the same position as HQ, $x^{\mu}$, and extract its temperature $T$ and velocity $u^{\mu}$ from the hydrodynamic simulations; then, boost the current HQ to the local rest frame (LRF) of the fluid cell and get the HQ momentum in this frame; 4. (4)

make a discrete time-step $\Delta t=0.01~{}{\rm fm/{\it c}}$ for the HQ in order to update its momentum $p^{\mu}$

• •

Boltzmann dynamics: for the current HQ with $p^{\mu}_{old}$, calculate its reaction probability $\Delta P_{l}$ for each possible scattering channel $l$ (Eq. 4); the target channel is selected according to the relative reaction probabilities $\Delta P_{l}/\Delta P_{total}$ (Eq. 6), meanwhile, the 4-momentum $p^{\mu}_{new}$ of the heavy quark after the scattering can be obtained according to the relevant scattering kinematics;

• •

Langevin dynamics: fix the drag and momentum diffusion coefficient with the fluid cell temperature $T$ (Eq. LABEL:eq:LTECoef), as well as the drag (Eq. 9) and thermal force (Eq. 10); and then, modify the HQ momentum $p^{\mu}$ according to the Langevin transport equation (Eq. LABEL:eq:LTE_Col); 5. (5)

update the HQ position after the time step $\Delta t$

[TABLE]

with the $p^{\mu}$ obtained in the previous step, and then boost back the HQ to the LAB frame; 6. (6)

repeat the above steps (3)-(5) when the local temperature $T\geqslant T_{c}$.

III.4 Heavy quark hadronization via fragmentation and coalescence

The heavy quark will suffer the instantaneous hadronization procedure via a “dual” approach, including fragmentation and heavy-light coalescence mechanisms, when the local temperature drops below the critical one $T_{c}=165~{}{\rm MeV}$. In this work, we follow the previous analysis S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018) and use this “dual” model for the final heavy-flavor meson productions.

Concerning the universal fragmentation function, the Braaten approach E. Braaten, K. Cheung, and T. C. Yuan (1993) is employed in this work. Due to the limitation of the measurements, it is difficult separate the open charmed hadrons produced in decays of each excited charmed hadrons. Practically, the relevant contributions can be treadted together with the fragmentation, by including their contributions in the fragmentation fraction of a particular open charmed hadron. Finally, the fragmentation fractions for the various hadron species are given by S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018) $f(c\rightarrow D^{0})=0.566$, $f(c\rightarrow D^{+})=0.227$, $f(c\rightarrow D^{\ast+})=0.230$, $f(c\rightarrow D_{s}^{+})=0.081$ and $f(c\rightarrow\Lambda_{c})=0.080$. The open charmed hadrons listed above are all the species included in the fragmentation model, and the higher state contributions are considered, which is consistent with the heavy-light coalescence model (to-be discussed below; Eq. 24).

The momentum distributions of heavy-flavor mesons ($Q\bar{q}$) reads

[TABLE]

where, $g_{\rm M}$ is the spin-color degeneracy factor; $f_{\rm Q}(\vec{x}_{\rm Q},\vec{p}_{\rm Q})$ is the phase-space distributions of heavy quark, which can be obtained after the HQ propagate through the underlying QGP medium; $f_{\rm\bar{q}}(\vec{x}_{\rm\bar{q}},\vec{p}_{\rm\bar{q}})$ is the one for light anti-quark, which follows the Boltzmann-J$\ddot{\rm u}$ttner distribution in the momentum space and it is spatially distributed on the freeze-out hypersurface. The coalescence probability for $Q\bar{q}$ combination to form the heavy-flavor meson in the $n^{th}$ excited state, is quantified by the overlap integral of the Wigner function of the meson and the $Q\bar{q}$ pair K. C. Han, R. J. Fries, and C. M. Ko (2016),

[TABLE]

where, $\vec{y}_{\rm M}=(\vec{x}_{\rm Q}-\vec{x}_{\rm\bar{q}})$ and $\vec{k}_{\rm M}=(m_{\rm\bar{q}}\vec{p}_{\rm Q}-m_{\rm Q}\vec{p}_{\rm\bar{q}})/(m_{\rm Q}+m_{\rm\bar{q}})$ are the relative coordinate and the relative momentum, respectively, in the center-of-mass frame of $Q\bar{q}$ pair; $W_{\rm Q}(\vec{x}^{\;\prime}_{\rm Q},\vec{p}^{\;\prime}_{\rm Q})$ and $W_{\rm\bar{q}}(\vec{x}^{\;\prime}_{\rm\bar{q}},\vec{p}^{\;\prime}_{\rm\bar{q}})$ are, respectively, the Wigner functions of heavy quark and light anti-quark with their centroids at $(\vec{x}_{\rm Q},\vec{p}_{\rm Q})$ and $(\vec{x}_{\rm\bar{q}},\vec{p}_{\rm\bar{q}})$, and they are both defined by taking the relevant wave function to be a Gaussian wave packet C. B. Dover, U. Heinz, E. Schnedermann, and J. Zimányi (1991). $W_{\rm M}^{\rm(n)}(\vec{y}^{\;\prime}_{\rm M},\vec{k}^{\;\prime}_{\rm M})$ denotes the Wigner function of heavy-flavor meson, which is based on the well-known harmonic oscillator C. B. Dover, U. Heinz, E. Schnedermann, and J. Zimányi (1991). The width parameter $\sigma_{\rm M}$ is expressed as S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018)

[TABLE]

where, $\langle r_{\rm M}^{2}\rangle\approx(0.9~{}{\rm fm})^{2}$ is the mean-square charge radius of D-meson; $e_{\rm Q}$ and $e_{\rm\bar{q}}$ are the absolute values of the charge of heavy quark and light anti-quark, respectively; the light (anti-)quark mass takes $m_{\rm u/\bar{u}}=m_{\rm d/\bar{d}}=300~{}{\rm MeV}$ and $m_{\rm s/\bar{s}}=475~{}{\rm MeV}$. We consider the various heavy-flavor meson species up to their first excited states ($n\leqslant 1$), see the Tab. II as shown in Ref. S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018) for details.

In Fig. 4 the coalescence probabilities obtained in central ($0-10\%$) Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02~{}{\rm TeV}$, are presented as a function of the charm quark transverse momentum. The contributions of the ground states and the first excited states are shown separately as the long dashed blue and short dashed black curves, respectively. It is found that the coalescence into a ground state has maximum probability at $p_{\rm T}^{\rm HQ}\sim 0$, and it decreases towards high $p_{\rm T}$, due to the difficulty to find a coalescence partner in this region. The coalescence probability into the first excited states shows similar behavior. The total coalescence probability is shown as a solid red curve, which decreases from $\sim 0.7$ at $p_{\rm T}^{\rm HQ}\sim 0$ to $0.2$ at $p_{\rm T}^{\rm HQ}=10~{}{\rm GeV}$. Moreover, the total coalescence probability is larger than $0.5$ in the range $p_{\rm T}^{\rm HQ}\lesssim 4~{}{\rm GeV}$, reflecting its dominance in this region.

As displayed in Fig. 4, the hadronization of charm quark is divided into three channels: fragmentation, coalescence to form D mesons at ground state and at first excited state. During the implementation, we generate a random number, using the Monte Carlo techniques, with flat distribution between zero and one, and then compare it to the above three probabilities. Finally, the target channel can be selected, and the momentum of the relevant heavy-flavor meson will be obtained by assuming the energy-momentum conservation in the $Q$ and $\bar{q}$ combination procedure. See Ref. S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018) for details.

IV Results with considering only elastic processes

IV.1 Momentum distribution inside a static medium

In order to study the difference between the Boltzmann and Langevin dynamics, in this sub-section, we focus on the time evolution of the charm quark momentum distribution, which is obtained inside a static medium with temperature fixed at $T=0.3~{}{\rm GeV}$, as well as the momentum initialized at $p=10~{}{\rm GeV}$.

In Fig. 5, the charm quark momentum distribution $dN/dp$ based on the Boltzmann model [(a)], is calculated at various times during the hydrodynamic evolution of the medium, showing as different styles. At the starting time $\tau_{0}=0.6~{}{\rm fm/\it{c}}$ (solid black curve), as expected, the initial $dN/dp$ behaves a delta distribution at $p=10~{}{\rm GeV}$. During the evolution up to $\tau=12~{}{\rm fm/\it{c}}$, $dN/dp$ is broadened comparing with the initial distribution, meanwhile, the average momentum is shifted toward low $p$, which is mainly induced by the drag force. This is caused by the fact that the initial momentum spectrum of charm quark is much harder than that of medium parton, and the multiple elastic scatterings are therefore dominated by the drag rather than the diffusion term S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018). The results based on the Langevin approach [(b)] present a different broadening behavior, which follows a Gaussian-like shape, as expected in the construction (Eq. 11). Similar results can be found in Ref. S. K. Das, F. Scardina, S. Plumari, and V. Greco (2014). Comparing Boltzmann with Langevin calculations, it is observed that the momentum broadening profile is stronger with the Boltzmann model, since the scatterings with large momentum transfer are allowed in this approach, which are discarded with the Langevin approach. Consequently, the relevant azimuthal angle distribution with the Boltzmann model, is expected to show a stronger broadening behavior as compared to Langevin. Note that, for both Boltzmann and Langevin dynamics, $dN/dp(\tau=3~{}{\rm fm/\it{c}})$ (dotted red curve) is followed by a tail in the range $p>10~{}{\rm GeV}$, where the interaction processes allow the charm quark to gain more energy respect to the lost term.

IV.2 Elastic energy loss inside a realistic medium

Figure 6 shows the average in-medium energy loss of charm quark, due to elastic scatterings, as a function of its initial energy. The results with the Boltzmann and Langevin dynamics are presented as the thick and thin curve, respectively. When comparing Boltzmann with Langevin results, they are similar at low energy ($E\lesssim 10~{}{\rm GeV}$) where the interactions with small momentum transfer are dominated, while the former one is systematically larger than the latter one at higher energy, resulting in a softer charm quark spectrum in this region.

IV.3 $R_{\rm AA}$ and $v_{\rm 2}$ for charm quarks

Figure 7 shows the nuclear modification factor $R_{\rm AA}$ ($v_{\rm 2}$) of charm quark obtained with the Boltzmann (solid red curves) and Langevin approach (dashed black curve) in central (semicentral) Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02~{}{\rm TeV}$. It is observed that $R_{\rm AA}$, as displayed in the panel (a), is suppressed at high $p_{\rm T}$ with the Boltzmann approach as compared to the Langevin approach. Therefore, charm quark loses more its initial energy while traversing the medium in the Boltzmann dynamics, which is consistent with results shown in Fig. 6. The elliptic flow coefficient $v_{\rm 2}$, as presented in the panel (b), with the Boltzmann approach is systematically larger as compared to the Langevin approach, which means that the Boltzmann dynamics is more efficient in producing $v_{\rm 2}$. Similar behavior is observed for different centrality classes and at different energies. It is most probably due to fact that the drag coefficient is larger in Boltzmann, resulting in a larger drag force acted on the charm quarks, which is able to introduce more significant interactions with the QGP partons, as well as to transfer more $v_{\rm 2}$ from the medium partons to the charm quarks.

Concerning the relative azimuthal angle distribution, the yields of the initially back-to-back generated $c\bar{c}$ pairs can be described by a delta distribution at $|\Delta\phi|=\pi$. After propagating through the medium, it is found that the above $|\Delta\phi|=\pi$ distribution is broadened within different initial transverse momentum interval $p_{\rm T}^{\rm c/\bar{c}}$, as shown in different curves in Fig. 8. With the Boltzmann approach (thick curves), It is clear to see that there is an almost flag behavior with the lower initial transverse momentum $p_{\rm T}^{\rm c/\bar{c}}<1.5~{}{\rm GeV}$ (dotted black curve), indicating the corresponding initially back-to-back properties are largely washed out throughout the interactions with the surrounding medium constituents S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018). Meanwhile, the broadening behavior tends to decrease with increasing $p_{\rm T}^{\rm c/\bar{c}}$ (dashed pink curve). Similar results can be found with the Langevin approach (thin curves). Note that the nuclear (anti-)shadowing effect is not included. Comparing Boltzmann with Langevin approach, they are close within small $p_{\rm T}^{\rm c/\bar{c}}$ region, while the former one shows stronger broadening behavior at larger $p_{\rm T}^{\rm c/\bar{c}}$. This is because, as explained above, with the larger initial drag force, the interactions in the Bolatzmann model are stronger and more powerful to pull the $c\bar{c}$ pairs from high momentum to low momentum.

IV.4 $R_{\rm AA}$ and $v_{\rm 2}$ for heavy-flavor mesons

Figure 9 presents the average $R_{\rm AA}$ [(a)] and $v_{\rm 2}$ [(b)] of the nonstrange D-meson ($D^{0}$, $D^{+}$ and $D^{*+}$) in central ($0-10\%$) and semicentral ($30-50\%$) Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02~{}{\rm TeV}$, respectively, with the Boltzmann (solid red curve) and Langevin approach (dashed black curve). It is observed that $R_{\rm AA}$ is suppressed at high $p_{\rm T}$ for the Boltzmann dynamics as compared to the Langevin, while $v_{\rm 2}$ is systematically higher in the while $p_{\rm T}$ region. This behavior is consistent with the results found at parton level (see Fig. 7). The available measurements for $R_{\rm AA}$ and $v_{\rm 2}$ (boxes) are shown for comparison. The calculations with both the Boltzmann and Langevin approaches fail to reproduce both the $R_{\rm AA}$ and $v_{\rm 2}$ data, which may due to the missing effects such as radiative energy loss mechanisms.

Based on the Bayesian model-to-data analysis, the original Lido hybrid model W. Y. Ke, Y. R. Xu, and S. A. Bass (2018) is developed to study the fundamental interaction mechanisms between heavy quark and the QGP constituents. However, in this work, we utilize only its Boltzmann module to describe the charm quark propagation inside the underlying thermal medium. Therefore, one cannot expect same $R_{\rm AA}$ and $v_{\rm 2}$ results between us, since the other used modules are different such as the initial charm quark momentum spectra, hydrodynamic modeling and the heavy-light coalescence in the subsequent hadronization procedure.

V Boltzmann and Langevin Dynamics with both elastic and inelastic processes

Concerning the scattering inelastically with the light (anti-)quarks and gluons of QGP in Boltzmann, both the $2\rightarrow 3$ gluon radiation and the $3\rightarrow 2$ inverse absorption processes are taken into account to guarantee the detailed balance. The relevant scattering matrixes are derived in an improved Gunion-Bertsch model in the Soft-Eikonal limit O. Fochler, Z. Xu, and C. Greiner (2010). Meanwhile, a Debye screening mass $m_{\rm D}^{2}=\frac{8}{\pi}(N_{c}+N_{f})\alpha_{s}T^{2}$ based on the Boltzmann statistics M. Greif, F. Senzel, H. Kremer, K. Zhou, C. Greiner, and Z. Xu (2017) is considered to regulate the $t-$channel gluon prpagator. The LPM effect is included by restricting the momentum space integration of the emission-absorption gluon with a coherence factor W. Y. Ke, Y. R. Xu, and S. A. Bass (2018),

[TABLE]

where, $t_{\rm 0}$ is the initial time for gluon radiation/absorption, $k_{\perp}$ is the transverse momentum of gluon, $\tau_{f}$ the gluon formation time

[TABLE]

with $E$ and $m_{\rm Q}$ are the HQ energy and mass, respectively, and

[TABLE]

characterize the light-cone momentum fraction of the radiated/sbsorbed gluon. Note that the coherence factor shown in Eq. 25, is obtained by requiring the radiation rate reduces to the Higher-Twist prediction in the limits W. Y. Ke, Y. R. Xu, and S. A. Bass (2018): soft-emission ($x\ll 1$); large gluon transverse momentum comparing with the momentum transfer ($k^{2}_{\perp}\gg q^{2}_{\perp}$).

The gluon radiation incorporated Langevin transport model is expressed as S. S. Cao, G. Y. Qin, and S. A. Bass (2015); S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018); S. Li and C. W. Wang (2018)

[TABLE]

Comparing with Eq. LABEL:eq:LTE_Col, the additional term $F^{\rm i}_{\rm Gluon}$ is the recoil force induced by the emitted gluons

[TABLE]

where, $p^{\rm ij}$ indicates the momentum of the radiated gluon. The transverse momentum together with the radiation time dependence of the radiated gluon is quantified by pQCD Higher-Twist model B. W. Zhang, E. K. Wang, and X. N. Wang (2004):

[TABLE]

$x=k/E$ denotes the fraction of energy carried away by the emitted gluon, which is equivalent to Eq. 27 in the high-energy ($E\sim p_{z}$) and collinear-radiation ($k\sim k_{z}$) limt; $\alpha_{s}(k_{\perp})=\frac{4\pi}{11N_{c}/3-2N_{f}/3}({\rm ln}\frac{k_{\perp}^{2}}{\Lambda^{2}})^{-1}$ is the strong coupling constant of QCD at leading order approximation; $P(x)=(x^{2}-2x+2)/x$ is the splitting function for process “$Q\rightarrow Q+g$”; $\hat{q}_{\rm q}$ is the jet transport coefficient for quarks555According to its definition, $\hat{q}_{\rm q}=2\kappa_{\perp}/v_{\rm Q}\approx 2\kappa_{\perp}$ at high energy $E\gg m_{\rm Q}$, where HQ velocity $v_{Q}=\sqrt{1-(m_{\rm Q}/E)^{2}}\sim 1.$; $\tau_{\rm f}=2x(1-x)E/[k_{\perp}^{2}+(xm_{\rm Q})^{2}]$ is the gluon formation time without considering the contribution of the gluon thermal mass ($m^{2}_{g}=m_{\rm D}^{2}/2$; see Eq. 26). It was argued S. S. Cao, G. Y. Qin, and S. A. Bass (2015) that an additional lower cutoff was imposed on the emitted gluon energy, $k\geqslant\pi T$, to balance the gluon radiation and the inverse absorption, so as to ensure that HQ equilibrium state can be reached after sufficiently long evolution time.

We can see that the implementations of radiative energy loss are different in the Boltzmann and Langevin approaches, which will apparently introduce the source of uncertainty when comparing these two models. However, it is still necessary to check further the modifications for each dynamics.

Figure 10 displays the elastic (or collisional) and inelastic (or radiative) energy loss as the dashed black and long dashed blue curves, respectively, while the combined results are shown as the solid red curves. The results with the Boltzmann and Langevin dynamics are presented as the thick and thin curves, respectively. We can see that the inelastic contribution (thick long-dashed blue curve) with the Boltzmann approach, is dominated at high energy, while the elastic component (thick dashed black curve) is significant at low energy. Similar behavior is observed with the Langevin approach (thin curves). The energy loss due to elastic scattering is larger with Boltzmann (similar with Fig. 6), while the total in-medium energy loss is larger with Langevin, in particular at high energy region.

Figure 11 shows the $R_{\rm AA}$ ($v_{\rm 2}$) of charm quark with Boltzmann (solid red curves) and Langevin approach (dashed black curve) but including both the elastic and inelastic scattering processes. When comparing with the results including only the elastic component (see Fig. 7), $R_{\rm AA}$ is suppressed (enhanced) at high (low) $p_{\rm T}$ for both these two models, while $v_{\rm 2}$ is enhanced in the range $2\lesssim p_{\rm T}\lesssim 7~{}{\rm GeV}$, in particular with the Boltzmann approach. This is mainly due to the fact that, as discussed in Fig. 10, inelastic component dominates at high $p_{\rm T}$, meanwhile, it introduces more interactions between charm quarks and QGP partons, transfering more $v_{\rm 2}$ from QGP partons to charm quarks. Similar behavior was observed in Ref. S. S. Cao, G. Y. Qin, S. A. Bass, and B. Muller (2013).

Figure 12 presents the average $R_{\rm AA}$ [(a)] and $v_{\rm 2}$ [(d)] of the nonstrange D-meson ($D^{0}$, $D^{+}$ and $D^{*+}$) in central ($0-10\%$) and semicentral ($30-50\%$) Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02~{}{\rm TeV}$, respectively, with the Boltzmann (solid red curves) and Langevin approach (dashed black curves). The central values are obtained in terms of the central predictions of the initial heavy quark spectra and the nPDFs. The bands are determined according to the total theoretical uncertainties, which are contributed by the FONLL model predictions on the initial charm momentum spectra, as well as the EPS09 NLO parameterization for the nPDF in Pb S. Li, C. W. Wang, X. B. Yuan, and S. Q. Feng (2018). We take the maximum derivation with respect to the central values in each $p_{\rm T}$ bin, and add in quadrature the above two components to get a conservative range. It is observed that $R_{\rm AA}$ is suppressed at high $p_{\rm T}$ for the Langevin dynamics as compared to the Boltzmann, while $v_{\rm 2}$ is systematically higher at moderate $p_{\rm T}$ ($2\lesssim p_{\rm T}\lesssim 4~{}{\rm GeV}$). This behavior is consistent with the results found at parton level (see Fig. 11). The available measurements for $R_{\rm AA}$ and $v_{\rm 2}$ (boxes) are shown for comparison. We find that, (1) when comparing with Fig. 9, both the elastic (or collisional) and inelastic (or radiative) energy loss mechanisms are needed to reduce the discrepancy between model and data, as concluded in Ref. S. Wicks, W. Horowitz, M. Djordjevic, and M. Gyulassy (2007); (2) the calculations with the Langevin approach seem to give a better description of the measured $R_{\rm AA}$ as compared to those with Boltzmann approach, in particular in the range $p_{\rm T}\gtrsim 10~{}{\rm GeV}$. Meanwhile, nonstrange D-meson $v_{\rm 2}$ calculated with the Boltzmann approach is closer to the available data at $p_{\rm T}\lesssim 4~{}{\rm GeV}$. The comparison of $R_{\rm AA}$ and $v_{\rm 2}$ gives the opposite indications about the two models, confirming that it is challenging to describe well $R_{\rm AA}$ and $v_{\rm 2}$ simultaneously, as observed in Ref. S. K. Das, F. Scardina, S. Plumari, and V. Greco (2015). A similar behavior can be observed in Pb–Pb collisions at $\sqrt{s_{\rm NN}}=2.76~{}{\rm TeV}$ (panel-b and e in Fig. 12) and Au–Au collisions at $\sqrt{s_{\rm NN}}=200~{}{\rm GeV}$ (panel-c and f in Fig. 12).

VI Conclusion and Discussion

In this work, we investigated the charm quark evolution via the Boltzmann and Langevin dynamics in relativistic heavy-ion collisions. By including only the elastic scattering contributions, the extracted drag coefficient ($\eta_{\rm D}$), momentum diffusion coefficients ($\kappa_{\rm L}$ and $\kappa_{\rm T}$) and spatial diffusion coefficient ($2\pi TD_{s}$) are calculated as a function of charm quark energy and the medium temperature, and further compared between the two approaches. The relevant in-medium energy loss together with its effect on the nuclear modification factor ($R_{\rm AA}$) and elliptic flow coefficient ($v_{\rm 2}$) at parton and hadron level, are discussed and compared with the available measurements at RHIC and LHC energies.

It is found that $\eta_{\rm D}$, $\kappa_{\rm L}$ and $\kappa_{\rm T}$ calculated from the Boltzmann dynamics ($2\pi TD_{s}\lesssim 7$ in the range $1<T/T_{c}<3$), are systematically larger than the ones obtained with the Langevin approach ($2\pi TD_{s}=7$). The collisional energy loss is larger with the Boltzmann approach, resulting in a smaller charm quark $R_{\rm AA}$ at $p_{\rm T}\gtrsim 10~{}{\rm GeV}$, as compared to the Langevin. Meanwhile, due to the larger drag force and stronger interactions in Boltzmann, it is more efficient in producing larger $v_{\rm 2}$, as well as in developing the broadening effect for the azimuthal angle distributions. The above $R_{\rm AA}$ and $v_{\rm 2}$ behaviors observed at parton level are well inherited by the corresponding heavy-flavor hadrons. When comparing the model with available data, it is realized that the calculations including only the contributions from the elastic processes, are unable to describe both the $R_{\rm AA}$ and $v_{\rm 2}$ measured at RHIC and LHC energies. This discrepancy can be reduced by including the inelastic contributions in both the Boltzmann and Langevin dynamics, even though the relevant implementations are different between these two models. Finally, we find that the model calculations for non-strange D-meson $R_{\rm AA}$ favor the Langevin approach, while $v_{\rm 2}$ prefer the Boltzmann approach. A simultaneous description of both $R_{\rm AA}$ and $v_{\rm 2}$ remains a challenge for both models.

It is necessary to mention that Ref. S. K. Das, F. Scardina, S. Plumari, and V. Greco (2014) is also a systematical study of Boltzmann versus Langevin, by considering only the elastic scattering processes. We obtain the similar conclusions for charm quarks, for instance, (1) drag coefficients show a decreasing momentum/energy dependence, while the momentum diffusion coefficients present an increasing behavior from the Boltzmann transport equation; (2) after the in-medium evolution, charm quark spectra is harder with the Boltzmann approach in the range $p_{\rm T}\lesssim 7-10~{}{\rm GeV}$, resulting in a larger (smaller) $R_{\rm AA}$ at $2\lesssim p_{\rm T}\lesssim 7~{}{\rm GeV}$ ($p_{\rm T}\lesssim 2~{}{\rm GeV}$); (3) as explained above, Boltzmann model gives larger $v_{\rm 2}$ at both parton and hadron levels. On the other hand, few differences are observed between us: (1) the calculations for $R_{\rm AA}$ and $v_{\rm 2}$ with both the Boltzmann and Langevin dynamics, including only the elastic processes, are failed to describe the available data in this analysis, while it is not true in Ref. S. K. Das, F. Scardina, S. Plumari, and V. Greco (2014), in particular for the Boltzmann approach, which reproduce well the measured $p_{\rm T}$ dependence of both $R_{\rm AA}$ and $v_{\rm 2}$ (see references therein); (2) in this analysis, the additional inelastic (or radiative) contributions are powerful to reduce the discrepancy with data, in particular at $p_{\rm T}$, however, this effect is not discussed in Ref. S. K. Das, F. Scardina, S. Plumari, and V. Greco (2014); (3) the theoretical uncertainty such as the one on the initial charm quark production, is taken into account in this work, which is missing in Ref. S. K. Das, F. Scardina, S. Plumari, and V. Greco (2014). These differences could be induced by the following sources: (1) comparing with the hybrid model utilized in this analysis (Sec. III), Ref. S. K. Das, F. Scardina, S. Plumari, and V. Greco (2014) takes different approaches in the relevant modules, such as the initial charm quark spectra is given by a parameterized power-law function, which works better only at high momentum region; nuclear (anti-)shadowing and heavy-light coalescence effects are missing; (2) with the Boltzmann approach, only quark-gluon scattering ($Q+g\rightarrow Q+g$) is considered in Ref. S. K. Das, F. Scardina, S. Plumari, and V. Greco (2014), while quark-quark ($Q+q\rightarrow Q+q$) is also included in the two-body interactions in this work; constant running coupling and Debye mass are used in Ref S. K. Das, F. Scardina, S. Plumari, and V. Greco (2014), but a momentum and temperature dependent scenario is adopted for us.

Finally, it is interesting to note that the resolution of the above model-to-data challenge may require the inclusion of nonperturbative dynamics in the medium. It may be noted that a similar challenge was previously investigated for light flavor jet energy loss and a viable solution was previously proposed by introducing a nontrivial medium color structure that includes both chromo-electric and chromo-magnetic degrees of freedom J. F. Liao and E. Shuryak (2008); J. F. Liao and E. Shuryak (2009b) and that leads to a strong temperature dependence of transport coefficients S. Z. Shi, J. F. Liao, and M. Gyulassy (2018a); J. C. Xu, J. F. Liao and M. Gyulassy (2016); S. Z. Shi, J. F. Liao, and M. Gyulassy (2018b). Whether a similar strategy may help address the $R_{\rm AA}$ and $v_{\rm 2}$ challenge in the heavy flavor sector would be an interesting problem for future investigation. More detailed studies will be reported in forthcoming publications.

Acknowledgements.

The authors are grateful to Weiyao Ke for providing the Boltzmann module in Duke model, as well as the data shown in Fig. 2. S. Li is supported by National Science Foundation of China (NSFC) under Grant Nos.11847014 and 11875178, China Three Gorges University (CTGU) Contracts No.1910103, Hubei Province Contracts No.B2018023, China Scholarship Council (CSC) Contract No.201807620007, and the Key Laboratory of Quark and Lepton Physics Contracts No.QLPL2018P01. C. W. Wang acknowledges the support from the NSFHB No.2012FFA085. R. Z. Wan acknowledges support from NSFC under Project No.11505130. J. F. Liao is supported by the National Science Foundation under Grant No.PHY-1352368. The computation of this research was performed on IU’s Big Red II cluster, which was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute, and in part by the Indiana METACyt Initiative. The Indiana METACyt Initiative at IU was also supported in part by Lilly Endowment, Inc.