Production of open-charm mesons in relativistic heavy-ion collisions

We present a theoretical framework to study open charm production in relativistic heavy-ion collisions. The coupling strength between the charm quarks and the QGP constituents, quantified by the spatial diffusion coefficient $2\pi TD_{s}$, is obtained by performing a phenomenological fit analysis to the lattice QCD calculations, resulting in $2\pi TD_{s}=const.$ (\textbf{Model-A}) and $2\pi TD_{s}=1.3 + (T/T_{c})^2$ (\textbf{Model-B}). We find that the relative azimuthal distribution of the initially back-to-back generated $c\bar{c}$ pairs presents a broadening behaviour, which is more pronounced for $c\bar{c}$ pairs with small initial $p_{\rm T}$, and when the Model-B approach is adopted. The competition between the initial drag and the subsequent collective effects tends to restrict the time dependence of charm quark $R_{\rm AA}$. Concerning the theoretical uncertainty on final D-meson nuclear modification, the nuclear shadowing and pp baseline components are dominat at high and low $p_{\rm T}$ ($p_{\rm T}\lesssim3~{\rm GeV/{\it c}}$), respectively. The measured D-meson $R_{\rm AA}(p_{\rm T})$ favors Model-A assumption for the diffusion coefficient both at RHIC and LHC, while their $v_{2}(p_{\rm T})$ prefer Model-B at moderate $p_{\rm T}$. These results confirm the necessity to consider the temperature- and/or momentum-dependence of $2\pi TD_{s}$ to describe well the D-meson $R_{\rm AA}$ and $v_{\rm 2}$ simultaneously.


I. INTRODUCTION
Relativistic 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 formation of a deconfined medium, the Quark-Gluon Plasma (QGP), is expected [1,2]. Experiments with heavy-ion collisions have been carried at the Relativistic Heavy Ion Collider (RHIC) at BNL and at the Large Hadron Collider (LHC) at CERN [3][4][5] in the last two decades. Among the various probes of the QGP, heavy quarks (HQ), i.e. charm and bottom quarks, are of particular interest [6][7][8][9][10][11] since, due to their large mass, they are mainly produced in hard scattering processes at the early stages of the heavy-ion collisions. Subsequently, they interact with the QGP constituents and experience the full evolution of the QGP medium. Thermal production of cc pairs in the QGP medium is expected to be negligible at the temperatures reached in heavy-ions at the RHIC and at the LHC. Interactions with the QGP constituents do not change the flavour, which make charm quark ideal probes of the medium properties.
HQ interact with the medium constituents in two main scenarios [12]: inelastic interactions via the exchange of color charge, resulting in the gluon radiation; multiple elastic collisions with small momentum transfer. Both of them cause energy loss for the HQ, usually referred as radiative and collisional energy loss, respectively. Therefore, HQ allow one to probe the mechanisms of multiple interactions with the medium, together with the strength of the collective expansion of the fireball. Considering the fact that HQ fragmentation function is quite hard [13], its properties can be inherited well enough by the corresponding open heavy-flavour hadrons (having charm or bottom quarks among these valence quarks) such as D mesons (D 0 , D + , D * + and D + s [14,15]) and B mesons (B 0 , B + and B s [16]).
Experimentally, the mentioned energy loss effects are studied by measuring the open heavy-flavour hadron nuclear modification factor R AA , which is defined as the ratio of the binary-scaled particle production cross section in nucleus-nucleus collisions to that in nucleon-nucleon collisions at the same energy, as well as the collective effects are investigated by measuring the elliptic flow coefficient v 2 , which is the second order coefficient of the Fourier expansion of particle azimuthal distributions. A strong suppression (R AA < 1) of high p T D-meson was observed at mid-rapidity in central nucleus-nucleus collisions at the BNL-RHIC and CERN-LHC [17,18], in-dicating that charm quark energy loss effect is significant. Meanwhile, a positive v 2 was measured in semicentral collisions and intermediate p T , suggesting that charm quark participages in the collective expansion of the medium. Theoretically, models were developed [19][20][21][22][23][24] to describe the available measurements. It was realized [25][26][27][28] 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-dependence of the interaction strength, which can be quantified by the spatial diffusion coefficient 2πT D s . Also, as pointed out in Ref. [9], one should explore the propagation of theoretical uncertainties in R AA calculations, including these due to the pp baseline calculation and the (anti-)shadowing parameterization.
In this work, we try to address these questions by taking into account different models for the temperaturedependence of 2πT D s , which are phenomenologically extracted from lattice QCD calculations, and then investigate its effect on the observables (R AA and v 2 ) for both charm quarks and open charmed mesons. Additionally, based on an instantaneous approach, the typical heavylight coalescence model for charm quark is adopted with the additional feature of including the contribution from higher states of the harmonic oscillator wave function (see Sec. III D 2 for details).
The paper is organized as follows: Sec. II is dedicated to the description of the HQ transport model which is implemented via Langevin approach, as well as to the discussion of the temperature dependence of 2πT D s . In order to simulate as completely as possible the evolution of HQ in heavy-ion collisions, Sec. III presents the additional components used to build our hybrid model, including the initial conditions, the hydrodynamics expansion of the underlying thermal medium and the hadronization mechanisms for both the medium constituent and the charm quarks. Sec. IV shows the results about R AA and v 2 at parton and hadron level. A summary section can be finally found.

II. LANGEVIN TRANSPORT APPROACH FOR HEAVY QUARK
In this section, we summarize the kinetic transport theory, including the Langevin approach that we use, to describe the heavy quark space-time evolution in a thermal medium. We also illustrate the development to include the recoil force induced by the radiated gluon. Moreover, we introduce a phenomenological model for the temperature-dependence of the drag and diffusion coefficients. The different parton in medium energy loss mechanisms are discussed as well.
A. Heavy quark diffusion as Brownian motion with the Langevin approach While traversing the QGP, HQ experience multiple elastic scatterings with its constituents and propagate with a Brownian motion, which can be quantified by a Boltzmann Transport Equation [29]. For large quark masses and moderate medium temperatures, the typical momentum transfers in interactions between HQ and the medium are small and the Boltzmann Transport Equation reduces to the Fokker-Plank Transport Equation [30]. In the framework of Fokker-Plank Transport, the interactions between HQ and the medium constituents are conveniently encoded in the drag and diffusion coefficients, which are related to each other via the relevant dissipation-fluctuation relation (or Einstein relation). Consequently, the phase-space distribution of HQ behaves according to the Boltzmann-Jüttner approach [31] and reaches the thermodynamic equilibrium in the infinite time limit.
In ultra-relativistic heavy-ion collisions, the Fokker-Plank Transport is equivalent to the Langevin approach, which consists of a "deterministic" drag term F Drag and "stochastic" diffusion term F Diff , expressed in terms of HQ momentum and its position as [12] where, dx i and dp i are the HQ position and momentum changes in the i th time-step dt. The drag force reads where Γ(p i ) is the drag coefficient. The thermal random force which acts on the HQ is expressed as In the so-called post-point scheme, the strength of the thermal noise C ij (p i ) can be associated to the momentum diffusion coefficient κ via [32] C ik C k j = κ(p)δ ij by assuming a isotropic momentum-dependence of the diffusion coefficient. As mentioned above, Γ(p) and κ(p) are bridged via the dissipation-fluctuation relation [32] As shown in Eq. 3, C ij is weighted by a random variable ρ = (ρ 1 , ρ 2 , ρ 3 ) which folloes the Gaussian-normal distribution.
B. Temperature-dependence of the drag and diffusion coefficients The spatial diffusion coefficient [33], defined as D s = lim p→0 T mQ·Γ(p) in the non-relativistic limit, is usually em-ployed to characterize the coupling strength in HQ transport calculations. The spatial diffusion coefficient is related to the momentum diffusion coefficient κ via [33] D s = 2T 2 /κ. In addition, D s is usually scaled by the thermal wavelength λ th = 1/(2πT ), namely D s /λ th = 2πT D s . The main features concerning its temperature and momentum dependence have been developed in many models [34][35][36][37]. The drag and diffusion coefficients can be represented in terms of 2πT D s as Note that, (1) the spatial diffusion coefficient D s is defined in the zero-momentum limit, while the notation D s , as shown in Eq. 5 and 6, refers to the definition of spatial diffusion coefficient extended to larger momentum values; (2) in this work, the HQ transport coefficientq Q is related to the momentum space diffusion coefficient κ viaq Q = 2κ [32]. We discuss below two approaches to model the temperature and momentum dependence of the spatial diffusion coefficient 2πT D s (T, p).
Model-A: Following the discussion in Ref. [19,33], one may neglect both the T -and p-dependence of 2πT D s and simplify it as In this case, there is only one parameter 2πT D s , which characterizes the coupling strength of HQ with the medium. It can be adjusted according to model-to-data comparisons. A remarkable feature of this approach is that the drag coefficient behaves as Γ ∝ T 2 , which is similar to AdS/CFT and pQCD calculations [19]. Model-B: alternatively, one can neglect the momentum dependence of 2πT D s , as mentioned in Ref. [9], and parameterize its temperature dependence as where, a and b are the adjustable parameters, and T c is the critical temperature for the transition from the deconfined QGP to a hadron gas. In this approach, the drag coefficient shows a weak dependence on the temperature which is consistent with the results presented in Ref. [26,38]. Figure 1 presents the temperature dependence of the spatial diffusion coefficient 2πT D s as obtained from lattice QCD calculations, i.e. Banerjee (pink circles [34]), Kaczmarek (blue square [39]) and Ding (red triangles [40]), as well as the results from the two approaches, i.e. Model-A (dashed green curve; Eq. 7) and Model-B (solid black curve; Eq. 8) described above. The model parameters were tuned to fit the lattice QCD results and their values are summarized in Tab. I. It is interesting  1. (Color online) Spatial diffusion coefficient 2πT Ds of charm quarks (mc = 1.5 GeV) from lattice QCD calculations (pink circle [34], blue square [39] and red triangle [40] symbols) at zero momentum. The phenomenological approaches (dashed green and solid black curves) are displayed as well.
to note that the values of κ/T 3 andq Q /T 3 obtained at certain values of T /T c fall within the ranges, reported in Refs. [41,42].  Fig. 1), as well as values obtained for κ/T 3 andqQ/T 3 . The values for other predictions are shown for comparison.

C. Heavy quark in-medium energy loss
As introduced in the previous sub-sections, the multiple scattering of heavy quarks (HQ) off the thermal partons inside a hot and dense QCD medium results in a Brownian motion, which can be described by the Langevin Transport Equation in the small momentum transfer limit. This accounts to the so-called collisional energy loss of the HQ. However, after traversing the QGP, heavy quarks can interact with the medium constituents via inelastic scattering, resulting in gluon radiation [43,44]. This medium-induced gluon radiation leads to the so-called radiative energy loss. In this analysis, we follow the strategy proposed in Ref. [21,45] to incorporate in the Langevin Equation both the collisional and the radiative energy loss of HQ propagating through the QGP medium. Equation 1 is therefore modified as with where, F Gluon i is the recoil force which acts on the HQ, and p ij indicates the momentum of the radiated gluon. The transverse momentum and radiation time dependence of the radiated gluon is quantified by pQCD Higher-Twist calculations [46].
It should be noticed that the Langevin Equation (Eq. 1) is modified to include the recoil force induced by the emitting gluon (Eq. 9 and 10), resulting in the possible violation of the fluctuation-dissipation relation (Eq. 4) by a certain amount: moreover, in this approach, the collisional and radiative energy loss effects are treated as independent, while, as pointed in Ref. [47,48], they are not entirely independent since the transport coefficients for collisional and radiative processes are correlated, which is not taken into account in this work. Finally, note that a lower cut-off is imposed on the gluon energy (ω πT [32]) to balance the gluon radiation and the inverse absorption, and to constrain the evolution of low-energy heavy quarks to follow the soft multiple scattering scenario, where the detailed balance is well defined. We follow Ref. [32] by assuming that the fluctuationdissipation relation (Eq. 4) is still valid between the drag and the diffusion terms of the modified Langevin Transport Equation (Eq. 9).

III. HYBRID MODELING OF HEAVY QUARK EVOLUTION
In order to simulate open charm hadron production in heavy-ion collisions, one needs to employ a hybrid model including the initial conditions, the hydrodynamics expansion of the underlying medium and the hadronization mechanisms for both the medium and the charm quarks.
A. Initial distribution of heavy quarks

Spatial-space initialization via Glauber model
The initial spatial distributions of heavy quark pairs are determined by simulating the initial entropy density distributions in heavy-ion collisions ¶ . The relevent transverse profile, i.e. perpendicular to the beam direction, is modeled by the MC-Glauber model (SuperM C [49]) which allows one to sample randomly the position of each ¶ Exactly, the spatial distributions of heavy quark pairs are sampled according to a event-averaged smooth initial transverse profile, which will be discussed in detail in Fig. 4 (Sec. III B 1). nucleon inside the projectile and the target nuclei according to their Woods-Saxon distributions, while the longitudinal profile, i.e. parallel to the beam direction, is described by a data-inspired phenomenological function.
At the initial time of the collision, τ = √ t 2 + z 2 ≡ 0, the entropy density, s(τ = 0, r ⊥ , η s ), can be factorized as where, s(τ = 0, r ⊥ ) is the initial entropy density deposited in the transverse plane [49]. The function ρ(η s ) (Eq. 11) allows one to quantify the longitudinal profile of initial entropy density as a function of the spatial pseudorapidity η s = 0.5ln(t + z)/(t − z). Experimentally, charged particle pseudorapidity distributions exhibit a plateau behaviour in the central region (η ∼ 0), followed by a rapid drop-off toward forward/backward regions (i.e. at large η) [50]. It was argued [51] that this observation can be reproduced by composing the initial entropy density into two regions: the initial entropy density is flat near η s ∼ 0 and smoothly fall-off as a half part of a Gaussian approach in the forward/backward space-time rapidity regions. Therefore, we parameterize the longitudinal distribution ρ(η s ) as where, y beam is the beam rapidity; ∆η and σ η describe the plateau and Gaussian fall-off behaviour, respectively; H is the Heaviside step function. Using the typical parameters such as the initial time scale τ 0 = 0.6 fm/c, the shear viscosity η/s = 1/(4π) corresponding to the predicted low-limit and the critical temperature T c = 165 MeV at both RHIC and LHC energies, we can compare the calculated charge particle multiplicity with the available measurements, and fix the parameters of Eq. 12. For instance, we take ∆η = 0.5 and σ η = 0.7 in central (0 − 10%) Au-Au collisions at √ s NN = 200 GeV, and the resulting mean multiplicity per participant pair is < 2N ch /N part > model = 3.745, which is consistent with the available measurements < 2N ch /N part > data = 3.64 ∼ 3.82 [52].

Momentum-space initialization via pQCD-based calculation
The initial momentum distributions of heavy quarks are determined according to FONLL (Fixed Order Nextto-Leading Logarithms [53][54][55]) calculations in the desired rapidity intervals, considering also the related systematic uncertainties on the calculations. The differential production cross section of charm quarks calculated in the range 0 < y < 0.5 for pp collisions at √ s = 5.02 TeV is shown in the panel-a of Fig. 2. The corresponding central values (solid red curve) of FONLL calculations are obtained with [53]  where µ R and µ F are the renormalization and factorization scales, respectively; m Q is the heavy quark mass, and its central value is m c = 1.5 GeV and m b = 4.75 GeV for charm and bottom, respectively. The upper (dashed black) and lower (dotted blue) curves represent the systematic uncertainties which are estimated by varying the renormalization and factorization scales and the quark mass in a conservative approach. The common variations are [56,57] The ratios to the central values are shown in the panel-b of Fig. 2. The uncertainty on parton distribution functions (PDFs) is given by different sets of inputs from CTEQ6 [58]. One can see that the uncertainty on QCD scales (solid red curve) dominates in the considered p T region, while the one on PDFs (long dashed green curve) is negligible for 2 < p T < 15 GeV/c. Note that the different sources mentioned above are considered in this analysis.
The charm quantum numbers are conserved in strong interactions, therefore, the charm quark c is always created together with its anti-quarkc. Then we assume the back-to-back azimuthal correlations, where, i = x, y, z. Consequently, the p T -and ydependence of the cc pair yields are sampled according to the FONLL calculations (e.g. panel-a in Fig. 2) via Monte-Carlo, and then, they are restricted to satisfy the above conditions.

Shadowing effect in nucleus-nucleus collisions
The nuclear modification of the parton distribution functions (nPDFs) should be taken into account in nucleus-nucleus collisions since the nucleons are bound in a nucleus. The most relevant effect at RHIC and LHC energy is a depletion at small Bjorken-x, usually called shadowing [59], which reduces the production cross section of charm quarks at low p T . At large Bjorken-x, shadowing is replaced by an enhancement of the PDF, usually called anti-shadowing in the literature [59]. In this work, we employ EPS09 NLO parameterization [60] for the gold (Au) and lead (Pb) nucleus PDFs. Figure 3 shows the effect on HQ production as the ratio of the production cross sections with and without EPS09 modification on lead (Pb) nuclear PDFs at √ s = 2.76 TeV. It is found that the effect of shadowing is more pronounced for charm (filled black circle) than that for bottom quarks (open red circle), due to the smaller Bjorken-x (∝ m T = m 2 Q + p T 2 at a given rapidity) for charm in this region probed by charm. The effect becomes similar towards high p T , induced by the similar Bjorken-x values probed when p T m Q . It is found that the ratio is slightly larger than unity (∼ 15% at maximum) in the range 10 p T 40 GeV, which will enhance the heavy-flavor production at high p T .
The corresponding systematic uncertainties on the EPS09 NLO parameterization are defined baed on various nPDFs sets which are obtained by tuning fit param-eters ¶ .

B. Underlying QGP medium
The hot and dense strongly-interacting medium produced in heavy-ion collisions is in pre-equilibrium state until it reaches local thermalization. We assume that the QCD medium undergoes a rapid thermalization and forms a QGP in equilibrium at τ 0 = 0.6 fm/c, at which the hydrodynamical evolution commences. The thermalization time scale is much shorter than the total life time of the QGP. Therefore, we neglect the pre-equilibrium evolution and thermalization in this analysis, assuming s(τ = 0) ≈ s(τ 0 = 0.6) hereafter. As discussed below, we utilize the initial conditions described above to model the initial entropy density distribution at the starting time scale of the hydrodynamical evolution, as well as an Equation of State (EoS) obtained via the lattice QCD calculations to describe the phase transition from the deconfined partons to the hadronic state.

Hydrodynamic description
The description of the QGP medium evolution is implemented by means of a 3+1 dimensional relativistic viscous hydrodynamics based on the HLLE algorithm [61], with τ 0 = 0.6 fm/c, shear viscosity η/s = 1/(4π) and critical temperature T c = 165 MeV in Au-Au and Pb-Pb collisions. It provides the space-time evolution of the temperature and the flow velocity field.
x (fm) Concerning the initial state simulation for the hydrodynamic medium evolution, we rely on the Glauber-based model introduced in Sec. III A 1 (see Eq. 11). However, by considering that the full event-by-event hydrodynamic ¶ see Eq. 2.12 and 2.13 in Ref. [60] for details. We use the nPDFs sets up to k = 7 in this analysis. simulation requires a huge computational time and disk space, we overcome these issues by utilizing a weighting approach to have an event-averaged smooth initial transverse profile of the entropy density distribution. Figure 4 shows the results obtained for the centrality interval 30 − 50% for Pb-Pb collisions at √ s NN = 2.76 TeV. The initial entropy density distribution for a single event is presented in the panel-a (left), while the result after averageing all the events belonging to 30 − 50% is displayed in the panel-b (right). As expected, event-by-event fluctuations are largely suppressed.

Isothermal freeze-out
In this work, we neglect the chemical freeze-out procedure and consider, only, the kinetic freeze-out (or freezeout since now) occurs at T c = 165 MeV.
To model the freeze-out of the QGP medium, we utilize an instantaneous approach across a hypersurface of constant temperature, namely isothermal freeze-out [62]. We employ a widely used model, cornelius [63], to reconstruct the isothermal particlization hypersurface, on which the momentum distributions of the different hadron species are evaluated using the Cooper-Frye formalism [64].

C. Simulation of heavy quark Brownian motion
In this sub-section, we describe the numerical framework utilized for the Langevin evolution of heavy quarks (HQ) coupled with the expanding underlying hydrodynamic medium. Generally, in the local rest frame (LRF) of the fluid cell, the HQ motion follows the modified Langevin Transport Equation, and the local temperature and the local flow velocity at the considered cell position are provided by the relativistic hydrodynamics simulations. The steps of this numerical procedure are (1) Sample the HQ pairs at the position x µ and momentum p µ , in the laboratory frame (LAB), according to the initial phase space configurations (τ ∼ 0); (2) Move all the HQ from τ ∼ 0 to τ 0 = 0.6 fm/c as free streaming particles, and modify the positions x µ correspondingly; (3) Search the fluid cell at x µ , and extract its temperature T and velocity u µ from the hydrodynamic simulations; then, boost the HQ to the 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 µ where the three terms in the right hand side are driven by • the drag force term F Drag i : drag coefficient Γ, which is determined by substituting Eq. 7 or Eq. 8 into Eq. 5, with the fluid cell temperature T obtained in the previous step; • the thermal force term F Diff i : the relevant time correlation profile behaves as: where, the momentum diffusion coefficient κ is given by Eq. 6. The above correlation is implemented by applying a momentum deflection sampled according to a Gaussian distribution with the width κ/∆t; • the recoil force term F Gluon i : during each time step ∆t, the Higher-Twist model gives the average number of radiated gluons, which is assumed to follow the Poisson distributions. The resulting total probability to radiate at least one gluon is used to determine whether or not the radiation process is triggered; (5) Update the HQ position after the time-step ∆t with the four-momentum p i obtained in the previous step, and boost back the HQ from the LRF to the LAB reference frame; As discussed above, the QGP medium hadronizes in our model when the local temperature reaches the critical one T c = 165 GeV. When the temperature T reaches T c , the heavy quark (HQ) will hadronize into the relevant heavy-flavor hadrons. It is known that the hadronization is an intrinsically non-perturbative process, which is treated as phenomenological models. Two approaches are usually employed to describe the HQ hadronization processes, namely "fragmentation" [65] and "heavy-light coalescence" [66]. In this work, we adopt a "dual" approach [67,68], including both fragmentation and coalescence, to model the HQ hadronization in heavy-ion collisions.

Fragmentation model
The HQ fragmentation can be implemented by using the "Lund symmetric fragmentation function" (PYTHIA 6.4 [13]) with all the defaults parameters. Alternatively, in this analysis, we utlize two other phenomenological models: • Peterson fragmentation function [69]: with the parameter fixed to c = 0.06 and b = 0.006 for charm and bottom, respectively; • Braaten fragmentation functions [70]: with the parameter r = 0.1 for charm quarks with m c = 1.5 GeV [71]. Figure 5 shows the normalized fragmentation functions as a function of the fragmentation fraction z, which is defined as the momentum fraction taken away by the fragmented heavy-flavor hadron with respect to that of its mother HQ. The average value of z is larger for Braaten model (dashed blue curve) as compared to Peterson model (solid black curve), resulting in a harder transverse momentum distribution at hadron level for the former one. Note that the fragmentation functions are assumed to be universal, i.e. to be the same, for different colliding systems and at different colliding energies.
To model the hadronization, we need also to provide the fragmentation fractions for the various hadron species, i.e. the fraction of charm quarks hadronizing in the different hadron species, except for the Lund-PYTHIA approach. Reference [72]

Heavy-light coalescence model
Within the instantaneous hadronization approach [73,74], the heavy-light quark coalescence is commonly modeled in terms of the overlap among the Wigner functions, which are based on Gaussian wave packets for the heavy quark and the light anti-quark, and the harmonic oscillator wave function the for charm hadron. However, in the calculations, some groups [23,75,76] consider only the harmonic oscillator wave functions of the ground state. The heavy-light coalescence probability for the excited charm hadron, such as c → D 1 (2420) 0 , is then obtained with some artificial assumptions. This was recently updated, for light quarks [77], by including the contribution from higher states of the harmonic oscillator wave functions. We follow this strategy and further extend it to charm quarks in this analysis.
According to the heavy-light coalescence model, the momentum distributions of heavy-flavor mesons (M ) composed of a heavy quark (Q) and a light anti-quark (q) are given as where, g M is the degeneracy factor accounting for the spin-color degrees of freedom; f Q ( x Q , p Q ) and fq( xq, pq) are the phase-space distributions of heavy quark and light anti-quark, respectively. For the heavy quark, f Q ( x Q , p Q ) can be obtained after the HQ propagate through the underlying QGP medium. For the thermal light anti-quark, fq( xq, pq) follows the Boltzmann-Jüttner distribution in the momentum space and it is spatially distributed on the freeze-out hypersurface ¶ . The coalescence probability is quantified by W are the relative coordinate and the relative momentum, respectively, in the center-of-mass (CMS) frame of the Qq pair; W Q ( x Q , p Q ) and Wq( x q , p q ) are, respetively, 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 [66]. W Finally, the overlap integral function for a heavy-flavor meson (Eq. 18) in the n th excited state in the CMS of the Qq pair is re-written as [77] W (n) M ( y, k) = υ n n! e −υ , υ = 1 2 Note that, in this work, we just consider the open charmed mesons up to their first excited states (n 1) according to the PDG data [56]. The width parameter σ M in the harmonic oscillator wave function is determined by the radius of the formed heavy-flavor meson. The charge radius of the Qq system is given by [32] where, e Q and eq are the absolute values of the heavy quark and light anti-quark charges, respectively. r 2 denotes the average squared distance, and it can be calculated from the Wigner function By substituting Eq. 20 into Eq. 23, we can relate σ 2 M to r 2 M via

Implementation of the "dual" hadronization model
Since we focus on the open charmed meson production, composed of a c (c) and its partnerq (q), in this work, a "dual" hadronization approach is implemented as described in the following. We will take as example the cq combination.
(1) Extract the three-vectors r i,c and p i,c for the i th charm quark position and momentum at T c , after the propagating through the QCD medium; (2) Sample a number of associatedq candidates, N i,partners , for the considered charm quark, and initialize their positions and momentum according to: • position: set the three-vector for the j th partner, r ij,q , according to the coordinate of the hypersurface cells which the considered current charm quark is located; • momentum initialization: sample the partner momentum, p ij,q , in the LRF of the fluid cell, according to the Boltzmann-Jüttner distribution; boost the partner to the LAB frame; (3) Calculate coalescence probabilities up to the first excited state, i.e. W M ( y ij , k ij ), via Eq. 21, in the CMS of each cq pair, with y ij ( r i,c , r ij,q ) and k ij ( p i,c , p ij,q ) defined in Eq. 19; get the relevant total coalescence probability: : the fragmentation process will be triggered for the considered charm quark; • rdm < P Max i : the coalescence approach will be implemented * rdm < W In the following, we will show the coalescence probability for c quarks into D-meson, including D 0 , D * 0 , D + , D * + , D + s and their first excited states listed in Tab. II, for different centrality classes and at different energies.
In the panel-a (upper) of Fig. 6 the coalescence probabilities obtained in central (0 − 10%) Pb-Pb collisions at √ s NN = 2.76 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 separatedly 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. On the other case the coalescence probability into the first excited states shows a slightly increasing behaviour in the range p HQ T 3 GeV, followed by a decreasing trend at higher p HQ T . This behaviour may be induced by the fact that energetic charm (anti-)quarks are needed to form D mesons in the highly excited states. The total coalescence probability is shown as a solid red curve, which decreases from ∼ 0.75 at p HQ T ∼ 0 to 0.15 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. Similar behaviour was found for Pb-Pb and Au-Au collisions in different centrality classes.
The panel-b (bottom) of Fig. 6 shows the results calculated for Pb-Pb collisions at √ s NN = 2.76 TeV in the 0 − 10% (solid curve) and 30 − 50% (dashed curve). The coalescence probability is systematically larger for more central collisions. This is because the parton density is larger in the 0 − 10% than in the 30 − 50%, resulting in a larger probability to form heavy-light combinations.

E. Experimental observables
We investigate the nuclear modification factor R AA which is defined as the ratio of the binary-scaled particle production cross section in nucleus-nucleus collisions to that in nucleon-nucleon collisions at the same energy, R AA (p T , y) = d 2 σ AA /dp T dy d 2 σ pp /dp T dy , where, d 2 σ AA /dp T dy is the p T and y double-differential production cross section in nucleus-nucleus collisions, scaled by the number of binary nucleon-nucleon collisions; d 2 σ pp /dp T dy is the double-differential result in nucleon-nucleon collisions. The deviation of R AA from unity is sensitive to the effects such as initial (anti-)shadowing and the in-medium energy loss, consequently, it can theorefore be used to quantify the nuclear effects in heavy-ion collisions. The elliptic flow coefficient v 2 is defined as the second harmonic when representing the particle azimuthal distributions via a Fourier expansion: Therefore, v 2 allows to describe the anisotropy of the transverse momentum distribution of the produced particles. It is sensitive to the EoS and to the initial conditions in the low p T region, while at hight p T it originates for the path-length-dependence of in-medium energy loss. As discussed in Sec. III A 2 (Eq. 15 and 16), the cc pairs are initially back-to-back generated before including the nuclear matter effects such as (anti-)shadowing and in-medium energy loss. Therefore, the initial relative azimuthal distribution dN cc pair /d|∆φ| can be described by a delta function at |∆φ| = π, with the relative azimuthal angle |∆φ| defined as, where, φ c (φc) denotes the azimuthal angle of the c (c) quark. However, the relative azimuthal distribution will be broadened by a certain amount after the propagation of the c quarks through the medium and this behaviour can be inherited by the corresponding heavyflavour hadrons dN DD /d|∆φ|. Note that, in this work, we neglect the hadronic rescatterings in the late stages, which can slightly reduce the open charmed meson R AA at high p T and enhance its v 2 at moderate p T [75], but, it is not expected to significantly affect the azimuthal correlation of DD pairs [78].

IV. NUMERICAL RESULTS
In this section, we summarize the results obtained at parton and hadron level, respectively. The comparisons with available measurements are discussed as well.
A. Results for heavy quarks

Profile of heavy quark energy loss
The average in-medium energy loss of charm quarks is shown in Fig. 7 as a function of the initial energy, displaying separately the contributions of collisional (long dashed blue curve) and radiative (dashed black curve) mechanisms. The results based on Model-A (Eq. 7) are shown in the panel-a (upper). It can be seen that collisional energy loss is significant at low energy, while radiative energy loss is the dominant mechanism at high energy. The crossing point between collisional and radiative contributions is around E = 7 ∼ 8 GeV. In the panel-b (bottom) of Fig. 7, the results based on Model-B (Eq. 8) are shown. The energy loss with Model-B is slightly larger than that with Model-A; the crossing point between collisional and radiative contributions is around E = 8 ∼ 9 GeV. This is caused by the fact that, (1) the initial transverse momentum spectrum of HQ is much more harder that of medium constituent, hence, the multiple elastic scatterings are dominated by the drag term; (2) the drag coefficient with Model-B is larger than Model-A around T c (see Fig. 1 together with Eq. 5), resulting in a stronger interaction strength between the HQ and the incident medium constituents. Consequently, the HQ lose more energy with Model-B approach.  Figure 8 shows the (raw) yields of the initially backto-back generated cc pairs, after propagating through the medium, with Model-A approach (Eq. 7) in central (0 − 10%) Pb-Pb collisions at √ s NN = 2.76 TeV.

Correlation in relative azimuthal angle
The different dashed curves refer to different intervals of c quark p T . Note that the (anti-)shadowing effects are not included in this plot. It is clearly observed that an almost flat |∆φ| distribution with the lower initial transverse momentum interval p c/c T < 1.5 GeV (dotted black curve), indicating the initially back-to-back properties are largely washed out throughout the interactions with the surrounding medium constituents. The broadening of the distributions tends to decrease with increasing p c/c T , reflecting a larger survival probability, for the initial back-to-back correlation, towards high p c/c T . The results in the whole momentum range p c/c T < 80 GeV are shown as solid red curve. Similar conclusions are obtained with Model-B (Eq. 8), which the broadening is more pronounced as compared to Model-A. As explained above, the larger initial drag term cases, the stronger interactions in Model-B, which are more powerful to pull the cc pairs from high momentum to low momentum.  The results including both collisional and radiative contributions (solid red curve) are shown as well.
3. Nuclear modification factor and elliptic flow Figure 9 shows the nuclear modification factor R AA of charm quarks obtained by considering only the collisional (dashed black curve) and radiative (long dashed blue curve) energy loss mechanisms, with Model-A approach in semi-central (30 − 50%) Pb-Pb collisions at √ s NN = 2.76 TeV. The result including the two effects (solid red curve) is closer to the one with only collisional energy loss at low p T , while it is closer to the radiative curve at high p T . This is consistent with what we discussed above in Fig. 7, which indicated that collisional energy loss is the dominant contribution at low momentum/energy of the charm quark, while at high p T radiative processes dominate. Similar behaviour can be found for the results with Model-B.
In Fig. 10 the charm quark R AA (panel-a, c and e) and v 2 (panel-b, d and f) are calculated, including both the collisional and radiative energy loss mechanisms, at various times during the hydrodynamic evolution of the medium. At the starting time τ 0 = 0.6 fm/c, R AA (panela; upper left) is determined by the (anti-)shadowing effect (see also Fig. 3), and the corresponding v 2 (panel-b; bottom left) is close to zero (even though the statistics is limited) for both Model-A (panel-a and b; left two panels) and Model-B (panel-c and d; middle two panels). During the evolution up to τ = 7 fm/c, R AA rises in the low p T region, while it drops at high p T , because the initial drag term is dominant with respect to the diffusion term, as discussed in Fig. 7. The variation between neighboring time-windows exhibits a decreasing trend, and the modification of R AA is less pronounced after τ = 7 fm/c. This may be induced by the late stage collective flow, which allows to transport the HQ from low momentum to high momentum, as mentioned in Ref. [79]. This means that the competition between the initial drag and the subsequent collective flow tends to restrict the time dependence of R AA . This can be confirmed by studying the time evolution of v 2 , as displayed in the panel-b (bottom left). It clearly shows that v 2 develops mostly at late times, reaching the maximum at τ ∼ 7 fm/c. The results of Model-A and Model-B are qualitatively similar.

Production cross section
The p T -differential production cross section of D * + mesons, in the range |y| < 1 and p T 7 GeV for pp collisions at √ s = 200 GeV is shown in the panel-a (upper) of Fig. 11. The curves in different styles are the model calculations, for the central values (Eq. 13), relying on various fragmentation functions: Lund-Pythia (dashed blue curve), Peterson (solid red curve) and Braaten (long dashed purple curve). As discussed in Sec. III D, the spectrum with the Braaten fragmentation function is found to be harder than that with with Peterson function. The experimental data (black boxes) are shown as well for comparison, which is obtained by scaling the cc production cross section reported by the STAR experiment [80] by the factor f (c → D * + ) = 0.230. Within the experimential uncertainties, the measured p T differential cross section is better described by the central prediction with the Braaten fragmentation function.
However, one should consider simultaneously the theoretical uncertainties on the FONLL calculation due to the perturbative QCD scales and the heavy quark mass (Eq. 14). The resulting D * + cross section is displayed in the panel-b (bottom) of Fig. 11. The curves in solid, dotted and long dashed curves denote the lower, central and upper bands of the model calculations, respectively. Within both the theoretical and the experimental uncertainties, the results based on the different fragmentation functions can provide a good description of the measured D meson corss section in the whole p T region [80]. Same conclusions can be drawn for pp collisions at √ s = 2.76 and 7 TeV. Hereafter, all the results are based on the Braaten fragmentation function. pT-differential production cross section for D * + mesons with Braaten fragmentation function, including the theoretical uncertainties. Experimental data arederived from Ref. [80].
In Fig. 12 the results of the calculations of the p Tdifferential production cross section, of D 0 (panel-a; upper), D + (panel-b; middle) and D * + (panel-c; bottom) mesons at mid-rapidity (|y| < 0.5) in pp collisions at √ s = 5.02 TeV are shown. They can be compared with the upcoming measurements at the LHC. Figure 13 shows the (raw) yields of DD pairs, produced from the initially back-to-back generated cc pairs, as function of the relative azimuthal angle |∆φ|, obtained with the Model-A approach for central (0 − 10%) Pb-Pb collisions at √ s NN = 2.76 TeV. Note that:

Correlation in relative azimuthal angle
(1) both the in-medium energy loss and nuclear (anti-)shadowing effects are included; (2) the heavy-light coalescence (Sec. III D 2) is not considered and all D mesons are produced via fragmentation. The broadening observed at hadron level is similar to the one observed at quark level (Fig. 8). A significant broadening is observed at low transverse momentum and it decreases with increasing transverse momentum. A qualitatively similar trend is found with Model-B, but with a more pronounced broadening.   It is found that the relative ratio (="with/without" nPDFs) between them is about 0.7 (1.1) at p T =1 (10) GeV. This behaviour is consistent with the observation made for charm quarks when discussing the results shown in Fig. 3. When comparing the model calculations with the corresponding measurements, we can see that, within the experimental uncertainties, the data can be described by the results including shadowing in the whole p T region, even without taking into account the heavy-light coalescence effect. The measurements with higher precision are needed to quantify this effect, in particular in the intermediate p T region. After propagating the theoretical uncertainties on the pp reference and on the nuclear (anti-)shadowing to the nuclear modification factor, we find that the uncertainty in the pp reference is significant at low p T (p T 3 GeV/c), while the one one the nuclear PDFs is dominated at higher p T . Figure 15  at low (high) p T as compared to Model-A, while v 2 is significantly higher at moderate p T . This is due to the different temperature-dependence of the spatial diffusion coefficient, as is was also pointed out when discussing the results at parton level (panel-b of Fig. 10). The calculations with Model-A seem to give a better description of the measured R AA [81] as compared to those with Model-B, in particular in the range p T 4 GeV even though the theoretical uncertainties are large. On the other hands, D-meson v 2 calculated with Model-B approach is closer to the available data [82] at p T 5 GeV. The comparison of R AA and v 2 gives the opposite indications about Model-A and Model-B, confirming that it is challenging to describe well R AA and v 2 simultaneously. A similar behaviour was observed in Ref. [26].
The results at RHIC energy are displayed in Fig. 16. Within the experimental uncertainties, the measured D 0 R AA but the data samples collected in 2010/2011 [83] and 2014 [84], for Au-Au collisions at   The contributions of the different hadronization mechanism are displayed separately: fragmentation (Braaten) as long dashed black curve and coalescence as dotted green curve. For the coalescence, the contributions of the ground states (dashed blue curve) and first exciting states (dot-dashed purple curve) are displayed separately. Experimental data are taken from Ref. [81] .
curves) and Model-B approaches (dashed blue curves).
The results with Model-A are closer to the measurements in the range 2 < p T < 4 GeV. The same conclusion can be drawn for the D + meson R AA [84]. The data-to-model comparison can be improved with the future high precision measurements. As observed at LHC energy (panel-b of Fig. 15), the temperature dependent coupling strength of Model-B allows to improve the description of the measured D-meson v 2 . Figure 17 shows the D 0 meson R AA obtained with Model-A approach at mid-rapidity (|y| < 0.5) in central (0 − 10%) Pb-Pb collisions at √ s NN = 2.76 TeV. Both the fragmentation and the coalescence mechanisms are considered in the hadronization model. The fragmentation (long dashed black curve) is dominant in the range p T 6 GeV, while the coalescence contribution (dotted green curve) is significant in 1 p T 5 GeV. Concerning the different components of the coalescence mechanism, the contributions due to the ground states (dashed blue curve) cannot be neglected in 2 p T 4 GeV, while the excited states contribution (dot-dashed purple curve) is dominant in the range p T 2 GeV. This is due to the associated coalescence probabilities (Fig. 6), as well as the related degeneracy factors (g M in Eq. 17).  Fig. 18, respectively. The charm quark hadronization is implemented by using the "dual" model, including both the fragmentation and the coalescence. The available data points for v 2 (boxes [86]) are shown for comparison. As already observed at the other collisions energies ( Fig. 15 and 16), the v 2 results with Model-B (dashed blue curve) give a better description of the data.

V. SUMMARY AND CONCLUSIONS
In this analysis, we investigated the charm quark evolution through the QGP together with the relevant open charmed meson observables in relativistic heavy-ion collisions. We tried to study the temperature dependence of the coupling strength of the charm quark and the medium, as well as its effects on the nuclear modification factor R AA and the elliptic flow v 2 of open charmed mesons in Au-Au collisions at √ s NN = 200 GeV and Pb-Pb collisions at √ s NN = 2.76 and 5.02 TeV. We built a theoretical framework to achieve this goal, and all the adjustable parameters were tuned according to the comprehensive sets of available data at RHIC and LHC energies. The coupling strength for charm quark 2πT D s , is determined according to lattice QCD calculations using two different assumptions: 2πT D s = const.
(Model-A) and 2πT D s = 1.3 + (T /T c ) 2 (Model-B). It is found that: (1) charm quark in-medium energy loss due to gluon radiation is dominant at high p T , while the quasielastic scattering is significant at low p T ; the energy loss is stronger with Model-B approach because the relevant 2πT D s is smaller, and therefore the initial drag term is larger, resulting in stronger interactions. Hence, charm quarks will lose more energy while traversing the QGP; (2) the azimuthal angle (|∆φ|) distribution of the initially back-to-back generated cc pairs presents a broadening behaviour, which is mainly due to quark pairs with small initial p T ; this broadening effect is more pronounced with Model-B due to the larger drag force, which is more powerful to pull cc pairs from high momentum to low momentum; (3) the charm quark R AA is mostly determined by interactions occuring in the time window 0.6 τ 7 fm/c, due to the competition between initial drag and subsequent collective effect; (4) hadronization due to fragmentation is dominant at high p T , while the coalescence is significant at moderate p T ; the coalescence probability induced by the higher state component is relevant at moderate-low p T , resulting in an enhancement of D-meson yield in this region; (5) the theoretical uncertainty on the D-meson R AA is dominated by the pp reference uncertainty at p T 3 GeV/c, and by the nuclear (anti-)shadowing parameterization at higher p T ; (6) model-to-data comparisons for D-meson R AA favor Model-A assumption for the dependence of 2πT D s on temperature, while the measured v 2 prefer Model-B. This conclusion holds true for all the available measurements at RHIC and LHC energy, suggesting the need for a temperature dependent 2πT D s , as well as a possible momentum dependent 2πT D s to describe simultaneously R AA and v 2 .
Some effects such as pre-equilibrium interactions and hadronic rescatterings are missing in this model. More detailed checks and results on this developments will be discussed in forthcoming publications.