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

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.


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 to form a new state of nuclear matter consisting of the deconfined quarks and gluons, namely the quark-gluon plasma, QGP [1,2]. 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 [3][4][5]. The QGP was found to induce the jet quenching, as well as to exhibit the collective flow behavior among various probes [6][7][8][9][10]. The jet quenching phenomenon is known [11] 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 AA , R AA (p T ) = dσ AA /dp T dσ pp /dp T .
(1) * lish@ctgu.edu.cn † liaoji@indiana.edu 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 [12,13] 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 2 , 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 [14][15][16][17][18] since, due to their large mass, (1) m Q Λ QCD , thus, its initial production can be well described by the perturbative Quantum ChromoDynamics (pQCD) at the next-to-leading order [19][20][21], in particular at high p T ; (2) m Q 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 [22,23]. 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 calculated directly in terms of the perturbative QCD. Due to large HQ mass and moderate medium temperature, the typical momentum transfers in interactions, q ∼ gT , are assumed small, gT m Q [24], 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 derived from the lattice QCD.
Many models were developed from the Boltzmann [25][26][27][28][29] and Langevin dynamics [30][31][32][33][34] 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 [35,36]) and B mesons (B 0 , B + and B s [37]). Comparing the theoretical calculations with available data, it was realized [38][39][40][41] that the simultaneous description of R AA and v 2 of open charmed meson at low and intermediate p 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 [42,43] develops a data-based hybrid model to extract the transport coefficient by utilizing the Bayesian model-to-data analysis. See Refs. [44][45][46] 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 derived from the lattice QCD calculations. However, the condition m Q gT may not always 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 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" [44]. 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, as well as the comparison for the derived transport coefficients. 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. V contains the summary and discussion.

II. BOLTZMANN AND LANGEVIN TRANSPORT APPROACHES FOR HEAVY QUARKS
A. Linearized Boltzmann transport model The linearized Boltzmann Transport Equation (BTE) reads where, p Q , E Q and f Q are the HQ 4-momentum, energy and distribution function, respectively. C[f Q ] denotes the collision integral, including all the interaction mechanisms between heavy quarks and the medium partons, i.e. light quarks and gluons. 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 [47] is utilized to sample f Q in each cell. The collision integral is solved by using the stochastic algorithm for evaluating the collision probability [48,49].
In the local rest frame (LRF) of the cell, the heavy quark transport is performed within a given time-step ∆t. Concerning a desired scattering process l, there are n (m) incoming (outgoing) partons, and the reaction probability ∆P l is expressed as [43] (4) where, Γ l (E Q , T, t) is the relevant scattering rate; g is the spin-color degeneracy factor of the incoming medium partons; ν is the statistical factor that corrects for doublecounting when there are identical particles in the initial/final state; f i denotes the heavy quark (i = Q) and medium parton (i =q, q, g) density, while the latter one follows the Maxwell−Jüttner distribution; |M | 2 l is the initial state spin-color averaged scattering matrix element squared. In this analysis, following Ref. [43], we consider both elastic and inelastic scattering processes whose total cross section is calculated via the perturbative QCD at leading-order [25,50]. Concerning the inelastic scattering, both the 2 → 3 gluon radiation and the 3 → 2 inverse absorption processes are taken into account to guarantee the detailed balance. Meanwhile, a Debye screening mass m 2 D = 8 π (N c + N f )α s T 2 based on the Boltzmann statistics [51] is considered to regulate the t−channel gluon prpagator. dΦ(n, m) in Eq.4 is the n + m body phase space integration, where, p in (p out ) indicates the total 4-momentum of all the incoming (outgoing) partons for a given scattering process l. Within the time interval ∆t, the total reaction probability ∆P total is given by It was argued [52] that the interactions between HQs and the medium partons can be encoded into the drag and momentum diffusion coefficients: 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.

B. Gluon radiation incorporated 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) [33,41,53] The deterministic drag force reads where η D ( p, T ) is the drag coefficient.
The stochastic force which acts on the HQ is expressed as with the Gaussian noise ρ j follows a normal distribution resulting in < ρ i > ρ = 0 and < ρ i ρ j > ρ = δ ij . Therefore, there is no correlation for the random force between two different time scales , 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 ij (t, p+ξd p, T ), via a parameter ξ ∈ [0, 1]. Typically, ξ = 0 for pre-point Ito, ξ = 1/2 for mid-point and ξ = 1 for post-point discretization scheme of the stochastic integral [54]. Finally, C ij can be represented in terms of the longitudinal (κ L ) and transverse momentum diffusion coefficients (κ T ) [55], i.e. (12) therefore, the relation between η D , κ L and κ T is given by where d = 3 denotes the spatial dimension. As pointed, HQ diffusions are conveniently encoded in the three coefficients η D , κ L and κ T . Note that Eq. 13 can be reduced to with the post-point scheme, i.e. ξ = 1. Following our previous analysis [41,53], a "minimum model" by assuming a isotropic momentum dependence of the diffusion coefficient, κ L = κ T ≡ κ, is adopted in this work, although it is just validated at p = 0 and, they not exactly the same at p = 0 region from the analytical calculations [56]. Eq. 14 is therefore further reduced to which is the so-called dissipation-fluctuation relation (or Einstein relation) in the non-relativistic approximation. The recoil force induced by the emitted gluons where, p 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 calculations [57]: where, z = ω/E Q denotes the energy fraction of the emitted gluon with respect to its parent HQ, and k ⊥ indicates the transverse momentum of gluon; α s (k ⊥ ) = 4π 11Nc/3−2N f /3 (ln k 2 ⊥ Λ ) −1 is the strong coupling constant of QCD at leading order approximation; P (z) = (z 2 − 2z + 2)/z is the splitting function for process "Q → Q + g";q q is the jet transport coefficient for quarks ¶ ; t 0 is the initial time for gluon radiation; is the gluon formation time, with E Q and m Q are the HQ energy and mass, respectively. It was argued [33] that an additional lower cutoff was imposed on the emitted gluon energy, ω π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.

C. Boltzmann vs. Langevin
In this sub-section, we mainly focus on the comparison of the transport coefficients obtained via the Boltzmann and Langevin approaches. We show before that the scattering rate (Eq. 4) for c + q → 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 E (G eV ) 5 10 15 20 T ( G e V ) dependence is weak, while the temperature dependence is stronger.

Boltzmann vs. Langevin: spatial diffusion coefficient
The spatial diffusion coefficient D s [52] scaled by the thermal wavelength 1/(2πT ), is defined at p → 0 limit, which can be obtained directly by substituting Eq. 7 with the Boltzmann approach. 2πT D s is available from the Lattice QCD calculation, moreover, it is found [53] that, according to a phenomenological fitting analysis with the Langevin approach, model predictions based on 2πT D s = 7 allow ¶ According to its definition,qq = 2κ to reproduce all the measured p 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. 18 by setting 2πT D 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. [33,51,58]; and (2) the drag and momentum diffusion coefficients in Eq. 15 can be represented in terms of 2πT D s as The temperature dependence of the spatial diffusion coefficient 2πT D s is presented in Fig. 2. The results from the Boltzmann (only c + q → c + q and c + g → c + g) and Langevin (only collisional) approaches are displayed as the dashed red and solid black curves, respectively. Lattice QCD calculations, i.e. Banerjee (pink circles [59]), Kaczmarek (blue square [60]) and Ding (red triangles [61]) 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 calculations. Similar behavior can be found by comparing with the CUJET3 predictions (red region) [28], as well as the results based on the Bayesian model-to-data analysis (green region) [43].  [59], blue square [60] and red triangle [61] symbols) at p = 0. The CUJET3 (shadowing red region [28]) and Lido model predictions (shadowing green region [43]), together with the results obtained from the Boltzmann (dashed red curve) and Langevin dynamics (solid black curve) are displayed for comparison.
With Eq. 18, the thermalization time of charm quark, defined in p → 0 limit [52], can be expressed as which is about 2.27 and 3.07 fm/c for Boltzmann and Langevin approach, respectively, with T = 2T c ≈ 330 MeV and m charm = 1.5 GeV.

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 ≈ 10 GeV, and as a function of the charm quark energy (bottom) at fixed temperature T = 0.3 GeV. The results obtained with the Boltzmann (Eq. 7) and Langevin model (Eq. 19) are shown as dashed red and solid black curves, respectively, in each panel.
Concerning the drag coefficient η 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 κ L [(b), (e)] and transverse momentum diffusion coefficient κ T [(c), (f)] increases strongly with increasing the energy and temperature via the Boltzmann approach, while they change slowly, as expected (Eq. 19), via the Langevin approach.

III. METHODOLOGY
In the previous analysis [41], 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.

A. 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 [62], while in the longitudinal direction, it is described by a data-inspired phenomenological function [41]. The initial momentum distribution of cc pairs is predicted by the FONLL calculations [19][20][21], assuming a back-toback azimuthal correlation between c andc (|∆φ cc | = π). 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 [63].
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 [41].

B. Hydrodynamic description
The underlying medium evolution is modeled by a 3+1D relativistic viscous hydrodynamics, vHLLE [64], with the initial time scale τ 0 = 0.6 fm/c, shear viscosity η/s = 1/(4π) and critical temperature T c = 165 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 MeV. An instantaneous approach across a hypersurface of constant temperature, namely isothermal freeze-out, is utilized and modeled by a widely used approach, Cornelius [65].

C. Heavy quark propagation in medium
We refer to Ref. [41] 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) sample a given number of HQ pairs at the position and momentum (x µ , p µ ), in the laboratory frame (LAB), according to the previous initial phase space configurations (τ ∼ 0); (2) move all the HQs from τ ∼ 0 to τ 0 = 0.6 fm/c as free streaming particles, and modify the positions x µ correspondingly; (3) search the fluid cell at the same position as HQ, x µ , and extract its temperature T and velocity u µ 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) make a discrete time-step ∆t = 0.01 fm/c for the HQ in order to update its momentum p µ • Boltzmann dynamics: for the current HQ with p µ old , calculate its reaction probability ∆P l for each possible scattering channel l (Eq. 4); the target channel is selected according to the relative reaction probabilities ∆P l /∆P total (Eq. 6), meanwhile, the 4- momentum p µ 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. 19), as well as the drag (Eq. 9) and thermal force (Eq. 10); the momentum of the radiated gluon is given by the higher-twist model (Eq. 17), together with the recoil force (Eq. 16); and then, modify the HQ momentum p µ according to the Langevin transport equation (Eq. 8); (5) update the HQ position after the time step ∆t with the p µ obtained in the previous step, and then boost back the HQ to the LAB frame; (6) repeat the above steps (3)-(5) when the local temperature T T c .

D. 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 MeV. In this work, we follow the previous analysis [41] and use this "dual" model for the final heavy-flavor meson productions.
where, g M is the spin-color degeneracy factor; f Q ( x Q , p Q ) is the phase-space distributions of heavy quark, which can be obtained after the HQ propagate through the underlying QGP medium; fq( xq, pq) is the one for light antiquark, which follows the Boltzmann-Jüttner distribution in the momentum space and it is spatially distributed on the freeze-out hypersurface. The coalescence probability for Qq 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 Qq pair [67], W Q ( x Q , p Q ) and Wq( x q , p q ) are, respectively, the Wigner functions of heavy quark and light anti-quark with their centroids at ( x Q , p Q ) and ( xq, pq), and they are both defined by taking the relevant wave function to be a Gaussian wave packet [68]. W (n) M ( y M , k M ) denotes the Wigner function of heavy-flavor meson, which is based on the well-known harmonic oscillator [68]. The width parameter σ M is expressed as [41] where, r 2 M ≈ (0.9 fm) 2 is the mean-square charge radius of D-meson; e Q and eq are the absolute values of the charge of heavy quark and light anti-quark, respectively; the light (anti-)quark mass takes m u/ū = m d/d = 300 MeV and m s/s = 475 MeV. We consider the various heavy-flavor meson species up to their first excited states (n 1), see the Tab. II shown in Ref. [41] for details. In Fig. 4 the coalescence probabilities obtained in central (0 − 10%) Pb-Pb collisions at √ s NN = 5.02 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 HQ T ∼ 0, and it decreases towards high p 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 ∼ 0.7 at p HQ T ∼ 0 to 0.2 at p HQ T = 10 GeV. Moreover, the total coalescence probability is larger than 0.5 in the range p HQ T 4 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 andq combination procedure. See Ref. [41] for details.

A. Momentum distribution inside a static hydrodynamic 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 hydrodynamic medium with temperature fixed at T = 0.3 GeV, as well as the momentum initialized at p = 10 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 τ 0 = 0.6 fm/c (solid black curve), as expected, the initial dN/dp behaves a delta distribution at p = 10 GeV. During the evolution up to τ = 8 fm/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 that of medium constituent, and the multiple elastic scatterings are therefore dominated by the drag term. 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. [69]. 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 angular distribution with the Boltzmann model, is expected to show a stronger broadening behavior as compared to the Langevin model. Note that, for both Boltzmann and Langevin dynamics, dN/dp(τ = 2 fm/c) (dotted red curve) is followed by a tail in the range p > 10 GeV, where the interaction processes allow the charm quark to gain more energy respect to the lost term.  Figure 6 shows the average in-medium energy loss of charm quark as a function of initial energy, displaying separately the contributions of elastic (or collisional) and inelastic (or radiative) mechanisms 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 energy loss (think blue curve) with the Boltzmann approach, is dominated at high energy, while the elastic component (thick black curve) is significant at low energy. Similar behavior is observed with the Langevin approach (thin curves). When comparing the Boltzmann with Langevin results, the total energy loss for the latter one is slightly larger at high energy, resulting in a softer charm quark spectrum in this region.   , with the Boltzmann approach is systematically larger in the range 2 p T 5 GeV, as compared to the Langevin approach, thus, the Boltzmann dynamics is more efficient in producing v 2 . Note that this is expected according to the discussion in Fig. 5.

B. Energy loss inside a realistic medium
Concerning the relative azimuthal angle distribution, the yields of the initially back-to-back generated cc pairs can be described by a delta distribution at |∆φ| = π. After propagating through the medium, it is argued [41] that the above |∆φ| = π distribution is broadened within different initial transverse momentum interval p    √ s NN = 5.02 TeV, respectively, with the Boltzmann (solid red curves) and Langevin approach (dashed black curves). The bands denote the theoretical uncertainties, which are mainly contributed by the QCD scales and the EPS09 NLO parameterization [41]. It is observed that R AA is enhanced (suppressed) at low (high) p T for the Langevin dynamics as compared to the Boltzmann, while v 2 is systematically higher at moderate p T (2 p T 3 GeV). This behavior is consistent with the results found at parton level (see Fig. 7). The available measurements for R AA and v 2 (boxes) are shown for comparison. The calculations with the Langevin approach seem to give a better description of the measured R AA [70] as compared to those with Boltzmann approach, in particular in the range p T 10 GeV. Meanwhile, nonstrange D-meson v 2 calculated with the Boltzmann approach is closer to the available data [71] at p T 4 GeV. The comparison of R AA and v 2 gives the opposite indications about the two models, confirming that it is challenging to describe well R AA and v 2 simultaneously, as observed in Ref. [39]. Experimental data taken from Ref. [70,71].
A similar behavior can be observed in Pb-Pb collisions at √ s NN = 2.76 TeV (panel-a and d in Fig. 10) and Au-Au collisions at √ s NN = 200 GeV (panel-b and e in  Fig. 10).

V. CONCLUSION AND DISCUSSION
In this work, we investigated the charm quark evolution via the Boltzmann and Langevin dynamics in relativistic heavy-ion collisions. The derived drag coefficient (η D ), momentum diffusion coefficients (κ L and κ T ) and spatial diffusion coefficient (2πT D 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 AA ) and elliptic flow coefficient (v 2 ) at parton and hadron level, are discussed and compared with the available measurements at RHIC and LHC energies.
It is found that η D , κ L and κ T calculated from the Boltzmann dynamics (2πT D s 7 in the range 1 < T /T c < 3), are systematically larger than the ones obtained with the Langevin approach (2πT D s = 7). The total in-medium energy loss is larger with the Langevin approach, resulting in a smaller (larger) charm quark R AA at high (low) p T , as compared to the Boltzmann. Meanwhile, Boltzmann dynamics is more efficient in producing v 2 , in particular at moderate p T , as well as in developing the broadening effect for the relative azimuthal angle distributions. The above R AA and v 2 behaviors observed at parton level are well inherited by the corresponding heavy-flavor hadrons. Comparing with the available data at RHIC and LHC energies, the model calculations for D-meson R AA favor the Langevin approach, while v 2 prefer the Boltzmann approach. A simultaneous description of both R AA and v 2 remains a challenge for both models.
The resolution of such 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 chromoelectric and chromo-magnetic degrees of freedom [76,77] and that leads to a strong temperature dependence of transport coefficients [28,78,79]. Whether a similar strategy may help address the R AA and v 2 challenge in the heavy flavor sector would be an interesting problem for future investigation. More detailed studies will be reported in forthcoming publications. . Experimental data taken from Refs. [72][73][74][75].