Directed flow of D mesons at RHIC and LHC: non-perturbative dynamics, longitudinal bulk matter asymmetry and electromagnetic fields

We present a study of the directed flow $v_1$ for $D$ mesons discussing both the impact of initial vorticity and electromagnetic field. Recent studies predicted that $v_1$ for $D$ mesons is expected to be surprisingly much larger than that of light charged hadrons; we clarify that this is due to a different mechanism leading to the formation of a directed flow with respect to the one of the bulk matter at both relativistic and non-relativistic energies. We point out that the very large $v_1$ for $D$ mesons can be generated only if there is a longitudinal asymmetry between the bulk matter and the charm quarks and if the latter have a large non-perturbative interaction in the QGP medium. A quite good agreement with the data of STAR and ALICE is obtained if the diffusion coefficient able to correctly predict the $R_{AA}(p_T)$, $v_{2}(p_T)$ and $v_{3}(p_T)$ of $D$ meson is employed. Furthermore, the mechanism for the build-up of the $v_1(y)$ is associated to a quite small formation time that can be expected to be more sensitive to the initial high-temperature dependence of the charm diffusion coefficient. We discuss also the splitting of $v_1$ for $D^0$ and $\bar D^0$ due to the electromagnetic field that is again much larger than the one observed for charged particles and in agreement with the data by STAR that have however still error bars comparable with the splitting itself, while at LHC standard electromagnetic profile assuming a constant conductivity is not able to account for the huge splitting observed.


I. INTRODUCTION
The theoretical and experimental studies on the formation and evolution of the Quark-Gluon Plasma (QGP) in ultra-Relativistic Heavy-Ion Collisions (uRHICs) have launched a new investigation stage, in which the main discoveries of the last decades on the hot QCD matter as a nearly inviscid fluid with strongly non-perturbative behaviour and the development of collective flows lay the foundations for dealing with new interesting features, such as the intense vorticity [1][2][3][4][5][6] and electromagnetic fields [7][8][9][10][11][12][13] generated in non-central collisions. Indeed, the fraction of the orbital angular momentum of the colliding system transferred to the hot plasma manifests itself with swirls of the order of 2 − 3 fm −1 whose observable effects could manifest in the chiral vortical effect [14][15][16][17], the spin polarization of baryons and vector mesons [18][19][20][21][22][23][24] and the "wiggle" slope in the directed flow of hadrons [25][26][27][28]. Furthermore, electromagnetic fields of the order of few to tens of m 2 π are produced in the overlap area mainly due to the motion of spectator charges and can be probed by observations of the chiral magnetic effect and related quantum phenomena [7,[29][30][31], the splitting in the spin polarization [32][33][34] and the splitting of the directed flow (a dipole asymmetry) [35][36][37][38][39][40][41][42][43][44]; more recently also a possible impact on v 2 , v 3 and the average transverse momentum < p T > has been pointed out [39]. In the last 15 years there has been a large effort to study the interaction of heavy quarks (HQs) with the bulk medium [45][46][47][48][49][50][51][52][53][54][55][56][57] and more recently a first determination of its interaction in terms of a drag coefficient γ(T ) or space-diffusion coefficient has been achieved, as expounded also in several reviews of the last two years [58][59][60]. There is a general consensus that the interac-tion is non-perturbative and around the pseudo-critical temperature T c can even get close to the value predicted in the infinite coupling limit by AdS/CFT [57,58]. A critical analysis of the several sources of indetermination can be found in recent joint activities comparing the features of the different approaches [61][62][63]. Heavy quarks and antiquarks have a very short formation time with respect to that of light quarks, therefore they are probably the earliest charged particles appearing in the uRHICs created matter. So one may expect they keep traces of both the initial stage and the subsequent evolution into a thermalized QGP. Furthermore, due to the large masses, their thermalization time is comparable or larger than the QGP lifetime.
A surprising novel theoretical prediction has been that the heavy quarks can manifest a directed flow v 1 = p x /p T that is more than one order of magnitude larger than the one of light hadrons despite their large mass [28,38]. This has been predicted for both the average charm v 1 (y) [28] and for the splitting ∆v 1 (y) between particle and anti-particles (D 0 andD 0 ) [38]. In Ref. [38] a Langevin approach coupled to the Maxwell equations has been used to describe the effect of the intense initial electromagnetic fields on the dynamics of heavy quarks in the QGP medium. It was predicted for the first time an electromagnetically-induced splitting in the directed flow of D 0 and D 0 .
A similar approach has been used in Ref. [28,40] to study the directed flow of neutral D mesons considering the tilt of the fireball in the reaction plane with respect to the beam axis. The main merit is to predict a dv 1 /dy 0.04 which is about half the one measured by STAR [64], but a successful prediction of a directed flow of D mesons can be more than one order of magnitude larger than the one of light hadrons.
Recent experimental efforts on measuring the directed flow of neutral D mesons at both RHIC [64] and LHC energies [65] have measured flow signals much larger than the one of light hadrons and even larger than early theoretical prediction for v 1 (y) or comparable with them for ∆v 1 (y) = v 1 (D 0 ) − v 1 (D 0 ) that could be associated to the presence of early electromagnetic fields. Among the D mesons, the neutral D 0 and D 0 are excellent particles to study the influence of the electromagnetic fields in the QGP phase, since it is natural to expect that the influence on them of the electromagnetic fields could have origin only in the deconfined phase.
We aim at investigating the effect of both the initial vorticity and electromagnetic fields on the directed flow of D 0 and D 0 in comparison to the recent experimental measurements by STAR [64] and ALICE [65] collaborations, once the longitudinal distribution of the bulk is fixed to reproduce correctly the light hadron v 1 and the charm-bulk interaction is the one able to describe with good accuracy the D mesons observables. The study we present here is performed by the relativistic Boltzmann transport equation for heavy quarks developed to describe the elliptic flow v 2 and the nuclear modification factor R AA of D mesons at both RHIC and LHC energies; hence it will be used as standard bulk-charm interaction the one corresponding to the determination of the nonperturbative 2πT D s coefficient [58,66]. Moreover, we have modified our standard boost-invariant initial conditions breaking the forward-backward symmetry by means of a tilted distribution of the fireball which leads to a finite directed flow of charged particles in agreement with the experimental data. Furthermore, with respect to the first prediction on the splitting ∆v 1 (y) [38] at LHC energy, we include the bulk vorticity and present prediction also and mainly for RHIC energy. We find a satisfying agreement with experimental data for both v 1 (D 0 ) and v 1 (D 0 ) predicting a dv 1 /dy 0.07 and for their difference ∆v 1 finding d∆v 1 /dy 0.01 at RHIC energy. A main part is dedicated to clarify the origin of the large v 1 of D mesons pointing out that its origin is different from the one of the bulk QGP and that its large value manifest only due to their non-pertubative interaction with the bulk QGP matter.
The article is organized as follows. In Sec. II we introduce our transport approach, focusing on our description of the heavy quarks dynamics and the electromagnetic fields. In Sec. III we explain in detail our modified initialization of the tilted fireball, while in Sec. IV we discuss the space-time evolution of the QGP vorticity profiles. Our results on the directed flow of neutral D mesons are presented in Sec. V and Sec. VI for RHIC and LHC energies respectively. Finally, we draw our conclusions.

II. CHARM QUARKS TRANSPORT EQUATION IN THE ELECTROMAGNETIC FIELD
We use a relativistic transport code developed to perform studies of the dynamics of heavy-ion collisions at both RHIC and LHC energies and different collision systems [66][67][68][69][70][71][72][73][74][75]. The dynamical evolution of gluons and light quarks as well as of charm quarks in the QGP bulk medium is described by the Relativistic Boltzmann Transport equations given by is the phase-space one-body distribution function of the parton j, k = g, q, q (gluon, quark or antiquark) and f Q (x, p) is the phase-space one-body distribution function of the heavy quark Q = c, c (charm or anticharm); on the left-hand sides F µν is the electromagnetic strength tensor and on the right-hand sides C = C 22 represents the relativistic collision integral accounting for 2 → 2 scattering processes.
For the bulk composed by quarks and gluons the evolution is described by Eqs. (1) and in this paper we assume that it is independent of f Q . In the collision integral C[f j , f k ] the total cross section is determined in order to keep fixed the ratio η/s = 1/(4π) during the evolution of the QGP, see Refs [70,73,74,76] for more details. However, this is equivalent to simulate the dynamical evolution of a fluid with specified η/s by means of the Boltzmann equation. In the present paper we have employed a bulk with massive quarks and gluons in the framework of a Quasi-Particle Model (QPM) [77] where quarks and gluons are dressed with thermal masses m g,q (T ) ∝ g(T ) T and the T -dependence of the coupling g(T ) is tuned to lattice QCD thermodynamics [78]. In this approach the dynamical evolution of the system has approximately the lattice QCD equation of state, with a softening of the equation of state consisting in a decreasing speed of sound approaching the cross-over region [78].
In the evolution equations for heavy quarks, Eqs.
(2), we treat the phase-space distribution functions of the bulk medium f j,k (x, p) as external quantities through C[f j , f k , f Q ] and we adopt the established approximation of neglecting collisions between heavy quarks. The HQs interact with the medium by means of only 2 → 2 elastic processes using scattering matrices calculated at Leading-Order in pQCD; nevertheless, a successful way to treat non-perturbative effects in heavy-quark scattering is given by the QPM [77]. In our approach the effective coupling g(T ) leads to effective vertices and a dressed massive gluon propagator for g+HQ → g+HQ and massive quark propagator for q + HQ → q + HQ scatterings. The detail of the calculations can be found in Ref. [66].
The hadronization process plays a crucial role in determining the final spectra, R AA (p T ) and v 2 (p T ). When the temperature of the fireball goes below the quark-hadron transition temperature, T c = 155 MeV, we hadronize the charm quark to D-meson by means of a fragmentation model, see Ref. [79] for more details. The multiplicity is chosen by matching it with the experimental value of the dN/dy for the considered centrality class.
We perform simulations with our default initialization and simulations with modified initial conditions that allow to implement the onset of vorticity in the fireball due to the angular momentum of the system. The latter modelling is explained in detail in Section III, whereas our standard initialization is as follows.
The initial conditions for the bulk in the coordinate space are given by the standard Glauber model assuming boost invariance along the longitudinal direction. In momentum space they are given by a Boltzmann-Jüttner distribution function up to transverse momentum p mj T while at larger momenta we adopt mini-jet distributions as calculated by pQCD at NLO order in [80]: we take p mj T = 2.0 GeV at RHIC and p mj T = 3.5 GeV at LHC energy. The charm and anticharm quark distributions are initialized in coordinate space according to the number of binary nucleon-nucleon collisions N coll . In momentum space we use charm quark production according to the Fixed Order + Next-to-Leading Log (FONLL) calculation [81] which describes the D-meson spectra in protonproton collisions after fragmentation [66].
In the last decade a lot of studies have been performed to determine the electromagnetic field generated in uRHICs [7][8][9][10][11][12][13]. In the present work we aim to study the impact of the electromagnetic field on heavy quark dynamics. The electromagnetic fields appearing in the kinetic equations (1) and (2) are computed as done in Ref. [38], following Ref. [35]. In our convention for the reference frame the nucleus A whose centre is at x ≥ 0 moves towards the positive z -axis and the nucleus B whose centre is at x ≤ 0 moves towards the negative z -axis. Therefore, the magnetic field B generated in noncentral collisions points on average along the negative y direction. The time evolution of B y is described assuming a constant electrical conductivity σ el of the QGP in order to obtain analytic solutions of the Maxwell equations. The values of σ el employed in our calculations are typical values that can be reached in the range of temperatures explored in a heavy-ion collision, in agreement with the lattice QCD calculations [82][83][84]. Due to Faraday induction an electric field E along the x axis is generated.
The total magnetic field is the sum of B +,− y generated by Z point-like charges as [7] In the above equation B + y (τ, η s , x ⊥ , φ) and B − y (τ, η s , x ⊥ , φ) are the elementary magnetic fields generated by a single charge e located in the transverse plane at x ⊥ = (x ⊥ , φ) and moving towards +z and −z respectively with speed β related to the longitudinal space-time rapidity η s ≡ arctan(z/t); x in and x out are the endpoints of the x ⊥ integration regions given by where R is the radius of the nucleus and b is the impact parameter of the collision. The elementary electromagnetic fields are obtained by solving the Maxwell equations and can be written as where α = e 2 /(4π) is the electromagnetic coupling and y ≡ arctan(β) is the beam rapidity of the + mover. In the above equation A and ξ are given by . In a similar way an electric field produced by the charges moving along the z direction can be obtained as this field has to be convoluted with the transverse charge distribution ρ ± (x ⊥ ) as for the magnetic field. In Fig. 1 [83,85,86]. From the left panel we see that at η s = 0 the electric field vanishes due to symmetry. However, at forward and backward rapidity E x become huge and comparable to B y , as we can see from the middle and right panels in which the results at positive rapidities are shown; for symmetry reasons, at the corresponding points in the negative η s plane B y is the same and E x has opposite sign: B y (η s ) = B y (η s ) and E x (−η s ) = E x (η s ). At η s = 0 the magnetic field |eB y | reduces initially by about one order of magnitude in 1 fm/c and the decrease is higher for smaller electric conductivity. At larger η s the decrease of |eB y | becomes smaller, especially at η s = 1.5.
It is important to notice that the electromagnetic field at RHIC and LHC is quite similar and has a maximum value that is about 0.1 GeV/fm, however in vacuum the maximum at RHIC is some factor larger while at LHC is more than 50 times larger.
Aside from E x and B y , the other components of the electromagnetic field averaged over many initial conditions vanish or are very small. However, in event-byevent collisions the large fluctuations in the positions of the protons inside the two colliding nuclei can generate non-zero values of the other components of the electromagnetic field with magnitude comparable but generally smaller than B y and E x [11,12]. Furthermore, for a more complete calculation we should also include the electromagnetic field generated by the participant protons; however, it has been shown in [9,35] that it is subdominant in magnitude with respect to that produced by spectators, especially in the early stage that is crucial for the formation of a directed flow of heavy mesons that is the focus of this work.
Within our transport approach we are able to simultaneously describe the R AA (p T ) and the v 2 (p T ) of D mesons both at RHIC and LHC energies [66]. The inclusion of the electromagnetic field does not affect visibly those quantities, while it produces a sizeable effect on the rapidity dependence of the directed flow v 1 of D 0 and D 0 as we will see in the following sections.

III. SET UP OF THE INITIAL CONDITION
In this subsection we discuss how we set up the initial vortical structure of the fireball produced in uRHICs.
In a non-frontal collision the system constituted by the two incoming nuclei possesses an angular momentum which depends mainly by the energy, the size and the impact parameter of the collision. Only a fraction of this angular momentum is transferred to the plasma created after the collision, since most of it is carried away by the spectator nucleons. The angular momentum retained by the QGP manifest itself as a nonzero vorticity of the system [1][2][3][4][5]. Even if the total angular momentum J of the plasma, which quantifies its global rotation, could be precisely determined, several phase-space configurations of the plasma would correspond to that value of J. Moreover, it is not well known so far how this vortical structure is produced in the plasma.
In hydrodynamical simulations J has been introduced as an asymmetric initial energy density distribution respect to the reflection η s → −η s [3,27]. In those initial conditions the initial flow velocity Bjorken components are zero, but the longitudinal asymmetry of the initial energy density profile generates fluid vorticity. At variance with the common use of hydrodynamic models, many simulations based on relativistic kinetic theory generate automatically a vorticity in the QGP since positions and momenta of partons are determined by the primary collision of the two incoming nuclei [4,5,87,88]. Hence, the initial longitudinal velocity distribution of the fireball is asymmetric with respect to the reflection η s → −η s and the induced shear flow produces a nonzero local vorticity.
In this paper we extend our model in order to take into account the nonzero angular momentum that is transferred to the QGP by the colliding nuclear system. Comparing it to our default initialization that does not include any effect of the angular momentum of the collision we are able to disentangle the effect of vorticity on dynamical evolution and final observables of charged particles and heavy mesons. The default type of initial conditions, introduced in the previous section, is given by a full longitudinally boost invariant distribution, i.e. a uniform pseudorapidity distribution in the range η s ∈ [−2, 2] for RHIC and η s ∈ [−3, 3] for LHC and a standard transverse distribution based on the Glauber model at the initial time t 0 of the simulation. The second initialization is an extension of the previous one in order to include the angular momentum generated in noncentral collisions, inspired by the modelling adopted within the hydrodynamic framework in Ref. [3,27].
Within our convention, the nucleus A whose centre is at x ≥ 0 moves towards the positive z -axis respectively and the nucleus B whose centre is at x ≤ 0 moves towards the negative z -axis. We implement a vortical structure in the QGP modifying the longitudinally boost invariant Bjorken picture with an initial density profile asymmetric with respect to the reflection η s → −η s : where x ⊥ is the transverse coordinate, ρ 0 = ρ(0, 0) is the density at the centre of the fireball and W is the wounded nucleon weight function given by N A and N B are the number of participant nucleons of the nuclei A and B respectively; they are determined from the Glauber model in the following way: where σ in N N is the inelastic nucleon-nucleon cross section taken to be 4.
determines in the initial density distribution a central plateau of length 2η s0 and gaussian tails at larger rapidity of width σ η . These two parameters are fixed to better reproduce the experimental pseudorapidity density of charged particles dN ch /dη. A complete description of the experimental tails in the fragmentation regions cannot be achieved within our model which involves only the partonic phase, however our initial condition allow us to reproduce fairly well the dN ch /dη in the central rapidity region |y| 1.5 which is of interest for our work.
The functions f ± (η s ) in Eq. (8), representing the emission contributions respectively from forward-going and backward-going participant nucleons, are given by The parameter η m in the functions f ± (η s ) determines an asymmetry in the contribution to the local participant density from the forward and backward-going nuclei and this leads to a tilted fireball in the reaction plane [27] which owns a huge angular momentum reflecting that transferred by the two initial nuclei during the collision. This initial spatial asymmetry of the fireball along the longitudinal direction is translated, by means of the collective evolution and expansion, to a final momentum distribution of the particles with a nonzero directed flow v 1 that is odd in pseudorapidity η. The angular momentum J of the plasma is not a directly measurable quantity and, moreover, the same value of J can be obtained by means of different realizations of the initial conditions. It has been shown within the hydrodynamical framework [3,27] that the tilted initial source parametrized by means of Eq. (12) allow to correctly reproduce the η dependence of the directed flow of charged particle in peripheral relativistic collisions, once a suitable value of the parameter η m is chosen. As every initialization that generates a vorticity in the system, this kind of initial conditions breaks the boost invariance on the longitudinal direction since the distribution of particles according to Eqs. (7)- (12) is not uniform in space-time rapidity η s even in the central range. reaction plane. The middle and bottom panels show the evolution of the proper density profile at two later times, t = 2 fm/c and t = 4 fm/c respectively. The information on how this spatial anisotropy is transferred to the momentum space is encoded on the rapidity-odd directed flow of charged particle v ch 1 , which we show as a function of pseudorapidity in Fig. 4; our result (magenta line) fairly matches the STAR experimental data for 5-40% central collisions [89] (maroon circles) in the range |η| < 1. It is important to note that the bulk QGP medium does not rotate like a rigid body but the angular momentum of the system constituted by the two initial nuclei is transferred to the plasma as a longitudinal shear flow, differently with respect to what happens in heavy-ion collisions at lower energy where a significant space rotation of the participant region or the strong bouncing-off in the transverse direction determine a quite sizeable transverse flow [90,91]. As mentioned in the previous section, the charm and anticharm quark distributions are initialized in coordinate space according to hard binary scatterings between nucleons of the initial colliding nuclei. Hence, their production is symmetric with respect to the reflection η s → −η s and does not follow the profile introduced in Eqs. (11)-(12) for characterizing the distribution of gluons and light quarks. In Ref. [28] this shift between the spatial profiles of bulk matter and binary collisions was introduced. It implies that at η s = 0 the area of heavy quark production points in the transverse plane is shifted with respect to the transverse section of the fireball, and the shift increases at higher absolute values of the spacetime rapidity. The asymmetrically-distributed bulk matter drags the heavy quarks, whose directed flow results to be amplified and about one order of magnitude larger than that of charged particles. For a more sophisticated calculation we should include, from one hand, a dependence on N coll of the wounded nucleon weight function in Eq. (8) for the energy deposition of the bulk matter in order to take into account the contribution from binary collisions and, from the other hand, the heavy quark production profile may present a mild asymmetry due to fluctuations in nucleon positions in the colliding nuclei; for the sake of simplicity, we neglect both effects. These are however expected to be subdominant with respect to the main mechanism generating the large v 1 of heavy quarks.

IV. VORTICITY OF THE QGP
Due to the initial condition explained in the previous section the fireball owns a nonzero angular momentum and a strong vortical flow. In this section we focus on Au+Au collisions at top RHIC energy √ σ N N = 200 GeV with impact parameter b = 7 fm; in all plots we show quantities averaged over about 1000 events. The total orbital angular momentum J of the QGP can be computed summing the contributions J i of all particles: where r i is the position vector of the particle respect to the centre of the overlap region and p i is its momentum. We find that J is mainly directed perpendicularly to the reaction plane; due to our convention of the reference frame its dominant component J y lies along the negative y axis. Regarding the event-averaged values, for RHIC collisions at √ σ N N = 200 GeV with b = 7 fm we find J y ≈ −10400 ( units) while the other components are four order of magnitude smaller, in agreement with the picture expected and discussed in previous works [1][2][3][4][5]. However, in each event the x and z components can reach values of few hundred . For what concern the time dependence, we have checked that the total angular momentum stay almost constant during the whole evolution, being conserved within 0.2 %. Our result for J y is in agreement with calculations based on the HIJING model [4]. However, a strict comparison of the magnitude of the QGP angular momentum computed through different approaches is not very meaningful; indeed, the precise value of J depends on several factors, first of all on the size of the fireball along the longitudinal direction, and a small difference in that length could produce a very different amount of angular momentum. A more suitable quantity for comparing different models is the local vorticity distribution. The classical vorticity is defined as the curl of the velocity field v: being the velocity of a fluid cell identified with the velocity of the particle flow inside the cell and hence computed as v = ( i p i )/( i E i ), where the sum runs over all particles in a given cell and p i and E i are momentum and energy of the i particle.
Since we are dealing with a relativistic system, it would be more appropriate to go beyond the non-relativistic computation, where several vorticities can be defined [3]. However, for the goals of this paper we do not need to use the values of the vorticity and the main aim of this section is to show that within our model we obtain a realistic magnitude and time behaviour of the vortical flow in the fireball in agreements with the ongoing research on the Λ polarization.
In Fig. 5 we plot the distribution of the component of the classical vorticity perpendicular to the reaction plane ω y on the η s − x plane for |y| < 0.25 fm at the times t = 0.6, 1, 4 fm/c. In agreement with other works [3,5], ω y shows a quadrupole pattern with nearly vanishing values along the axis x = 0 and η s = 0 and increasing absolute values for increasing x and η s . Within our convention for the reference frame (nucleus with centre at positive x moving towards positive z) ω y is negative/positive in the regions where x and η s have the same/opposite sign. However, as pointed out in Ref. [5], a locally nonzero vorticity is not directly related to a global angular momentum and these vortical patterns could be mainly due to the radial flow of the fireball; the global rotational motion can be pinpointed by performing an average over the fireball in order to cancel out the radial flow contributions to the local vorticity and obtain a quantity that better reflects the strength of the vorticity acting on the entire interaction area [4,5]. The space-averaged vorticity is defined by where g(x ⊥ , η s ) is a weighting function. In Fig. 6 we show the time evolution of the dominant component of the classical vorticity ω y averaged according to Eq. (15) over the whole transverse plane and space-time rapidities |η s | < 1; we consider three weighting functions: the particle density n (dashed green line), the energy density ε (dot-dashed red line) and a sort of moment-of-inertia density I = ρ 2 ε (solid blue line), where ρ is the distance of the cell from the y axis. At early times | ω y | increases and reaches a peak at about 0.5 fm/c; then it begins to decrease with time approaching zero at the late stage of the evolution. We can compare our solid blue line with the same quantity computed in Ref. [5] with the AMPT model: in both cases | ω y | weighted with I evolves after 1 fm/c roughly as ∼ exp(−t/2) but the peak value in our model is three times higher than that obtained in Ref. [5]. However, the values of the local non-relativistic vorticity on the reaction plane at t = 1 fm/c depicted in the middle panel of Fig. 5 agree to a great extent in the two models, with maximum absolute values of about 2 fm −1 in the area shown in the plot. Despite the different initial conditions and models, our results for the angular momentum, the longitudinal velocity profile (Fig. 2) and the space-averaged vorticity (Fig. 6) are qualitatively in agreement with previous works [4,5].
It is worth to comment on the relation between the spatial profiles of vorticity and density and the consequences on final flow observables. Looking at the quadrupolar pattern of the local distribution of vorticity on the reaction plane shown in Fig. 5 we observe that the system is rotating clockwise in the regions where ω y < 0 and counterclockwise where ω y > 0. As we can see from the upper panels of Fig. 3 the regions with ω y < 0 correspond to the areas with higher proper density. This is the reason why, as evident from Fig. 6, the vorticity weighted over the fireball is negative and the fireball is globally rotating clockwise. As a consequence the directed flow generated in the central rapidity region has a negative slope, see Fig. 4.

V. DIRECTED FLOW OF D MESONS AT TOP RHIC ENERGY
One of the main goals of this paper is to show how the rapidity dependence of the directed flow v 1 of neutral D mesons is originated by the vorticity of the bulk matter coming from a a tilted longitudinal distribution and the non-perturbative interaction of HQ that induce a transverse pressure push that drags the HQs in the transverse direction; on top of this there is a charge motion induced by the electromagnetic fields produced in the collision mainly by spectator protons. We follow the ideas of Refs. [28,38,40]. In this section we refer to Au+Au collisions at top RHIC energy √ σ N N = 200 GeV, whose experimental data for the v 1 (y) of D 0 and D 0 mesons have been measured by the STAR Collaboration in the centrality bin 10-80% [64]. We performed simulation at b = 9 fm in order to compute the v 1 of charmed mesons and disentangle the effects of the electromagnetic fields and the vorticity. As explained in the previous sections, angular momentum and vorticity is created in our model by means of an asymmetric initial distribution of the QGP density in the x − η s plane and the parameters in Eqs. (11)-(12) are fixed to reproduce the experimental data of the directed flow of charged particles. We use the STAR data in the centrality classes 5-40% and 40-80% in order to have a reference band for our simulations at b = 9 fm. In Fig. 7 we compare the experimental data of the  Fig. 7 we can also get information on how the initial density asymmetry introduced with Eqs. (7)-(12) is translated to the final azimuthal asymmetry measured by the directed flow. Indeed, a larger value of η m corresponds to a milder asymmetry between the positive and negative η s sides of the reaction plane and this leads in turn to a lower value of the particle directed flow. It can be noticed that the dependence of the v 1 of the bulk is significantly dependent of the value of η m and allow to tune it with good accuracy. However for the aim of the paper it will be even more interesting to see that the HQ v 1 appears to be nearly independent on it.
In Fig. 8 we plot the rapidity dependence of the averaged directed flow of D 0 and D 0 mesons. In the top panel we show the result for two different values of the parameter η m of Eqs. (12) as in Fig. 7: the solid blue line and the dashed red line are obtained respectively with η m = 1.1 and η m = 1.2. We notice that, even though the charged particle v 1 is very sensitive to the parameter η m , the directed flow of D mesons is influenced only slightly.
A far big influence on v D 1 is instead exerted by the initial condition for the space distribution of charmed quarks, namely if c and c are also distributed along η s according to Eqs. (7)-(12) as done for the light quarks, hence giving them since the beginning of the collision evolution an Our result for the directed flow of D mesons is few times larger than that found in previous works [28,92] and is in the right ballpark with respect to the STAR data [64] labelled with black diamonds; however, the experimental uncertainties are still quite large for a very clear comparison. Moreover, the values of v D 1 in the central rapidity range |y| < 1 are about 20-30 times larger than those of charged particle shown in Fig. 7. As anticipated in Sec. III, this is a result of the shift in the transverse profiles of the tilted QGP fireball and the symmetricallydistributed emission points of heavy quarks from initial hard scatterings [28]; this shift leads to an amplification of directed flow of D mesons due to the collisions undergone with the QGP and the drag exerted from the bulk matter. However, if the c − c pairs are initialized according with a tilted source, as done for gluons and light quarks, the transverse spatial distributions of heavy quarks and quark-gluon plasma are superimposed. We show this in Fig. 9 where the density profile of both the bulk and the charm quarks are shown as a function of the transverse coordinate x at the space rapidity η s = 1. Once the charm are distributed as the bulk matter the large v 1 of the HQs disappears. In fact in this case the small v 1 comes only by the weak rotation as for the bulk matter and hence this leads to a v 1 even smaller then the one of the bulk, as shown by the dashed line in the bottom plot of Fig. 8. This allows us to understand a fundamental aspect of the directed flow. Its origin appears to be different from the one of light particles. The large v 1 is due to a pressure push of the bulk that having a different space distribution pushes the heavy quarks toward the negative x−direction at positive η s .
The effect of such a pressure gradient in the transverse direction is however effective only because the interaction of the HQ with the bulk is largely non-perturbative. To show this we plot in Fig. 10 the v 1 (y) if one keeps the non-tilted space distribution but assume a drag and diffusion of the charm quarks according to pQCD. It is clear that in this case the bulk is not able to transmit the transverse push because the drag is much smaller. Therefore in turn the large v 1 of D meson is a signature of the non-perturbative interaction of the HQ with the bulk hot QCD matter that is however visible only because the bulk matter is tilted in the longitudinal direction. We get a reasonably good prediction of v D 1 because we assume for heavy quarks the drag and diffusion extrapolated from the phenomenological analysis of R AA (p T ) and v 2 (p T ) as summarized in Ref. [58]. GeV with impact parameter b = 9 fm. The experimental data from STAR Collaboration [64] are also shown.
In Ref. [38] some of us pointed out for the first time that the electromagnetic field generated in the early stage of the heavy-ion collisions can induce a splitting in the v 1 between charm and anti-charm. The calculation were performed in a Langevin framework discarding the longitudinal tilted distribution of the bulk. We present here the first results incorporating the electromagnetic field dynamics in our Boltzmann transport approach, as described in Sec. II, including the vorticity and the tilted bulk distribution. In Fig. 11 it is shown the influence of the electromagnetic fields on the rapidity dependence of the directed flow of D 0 and D 0 mesons for Au+Au collisions at top RHIC energy with impact parameter b = 9 fm. The simulations include the Lorentz force in the Boltzmann equation and therefore the effect of the electromagnetic fields on particle dynamics. For the computation of the time evolution of the electric and magnetic fields the value of the electric conductivity used is σ el = 0.023 fm −1 , a central value from lattice QCD calculation at T = 2 T c . We have performed also simulations without the electromagnetic fields and found that the two lines of the v 1 of D 0 and D 0 lie one on top of each other. We see clearly from Fig. 11 that a splitting of the two curves is generated when the electromagnetic fields are considered [38,40]: if we focus on the region y > 0 the v 1 of D 0 is pushed downwards and that of D 0 is pushed upwards. Since we are assuming that D 0 and D 0 are produced by fragmentation of c and c quarks respectively, the v 1 of neutral D mesons is driven by the v 1 of charm quarks and antiquarks; from Fig. 11 we deduce that at forward rapidity positively-charged particles (c) get a negative contribution from the electromagnetic fields to the v 1 and negatively-charged particles (c) get a positive contribution -and vice versa at backward rapidity. This means that the total effect coming from B y and E x which push particles in opposite directions is slightly dominated by the electric field. This is also consistent with the results obtained in previous studies for light mesons, both in hydrodynamic calculations of RHIC Au+Au and LHC Pb+Pb collisions [35,39] and in transport simulations of asymmetric collisions [36,37] and small systems [41] at top RHIC energy. In addition to the v 1 splitting induced by the electromagnetic field, there may be a further mismatch of the directed flow of D mesons conveyed by the initial orbital angular momentum of the two colliding nuclei, which contributes to the directed flow of u and d quarks [41,93,94]: as a consequence, the D 0 (cu system) gets a larger v 1 at forward rapidity and a smaller v 1 at backward rapidity, following the directed flow pattern of spectators. This effect, that is not included in our model, would increase the splitting between D 0 and D 0 especially at higher y, while in the central rapidity region the force exerted by the electromagnetic field dominates [41]. However it can be expected to give correction of the order of 10 −3 on v 1 given the small transverse flow of the bulk. We have also investigated how the electric conductivity affect the evolution of the electromagnetic fields and, consequently, the v 1 splitting of D mesons. The last one can be easily quantified by the difference in the directed flow between one particle and the corresponding antiparticle, i.e., for D mesons ∆v D This quantity is plotted in Fig. 12 for different values of the electric conductivity: the result for σ el = 0.023 fm −1 (solid blue line), which is the same of the solid lines in Fig. 11, is compared with those obtained with σ el = 0.0115 fm −1 (dotted magenta line) and σ el = 0.46 fm −1 (dashed brown line), which correspond respectively to halving and doubling the value of electric conductivity firstly considered. Only in the latter case, larger σ el , we see a mild change in the ∆v D 1 for higher rapidities |y| > 1. This can be understood by looking at the middle and right panels of Fig. 1 in which we show the temporal behaviour at forward space-time rapidity of the fields for the considered values of electrical conductivity: for σ el = 0.46 fm −1 the magnetic field has a milder decrease with respect to smaller conductivities and this in turn means that the produced electric field is lower and then its effect on the propagation of particle is smaller; this explains why the splitting ∆v D 1 is smaller with respect to that obtained with lower σ el .
In Fig. 13 we plot the time evolution of the v 1 of D 0 and D 0 versus rapidity for simulations with σ el = 0.023 fm −1 . We notice that at midrapidity most of the v 1 of D 0 and D 0 is generated within the first 3 fm/c; this holds for both the average and the difference of the D meson v 1 , which give information on the temporal behaviour of the imprint of the vorticity and the electromagnetic field respectively. In particular, in Fig. 14 we show that for the y−dependence of both v D 1 and ∆v D 1 the slope in the central rapidity region |y| < 0.5 reaches half of its maximum (absolute) value in the first 1 fm/c; then, the averaged D meson directed flow completely saturates at 3 fm/c while for the splitting there is still a mild effect of about a 20% at later times. plotted in Fig. 16 in comparison with the experimental data from ALICE Collaboration [65]. Thin lines correspond to the simulation with the standard initialization (which generate a zero directed flow of charged particles) whereas thick curves are the results with initial conditions given by Eqs. (7)- (12). Both simulations includes the electromagnetic fields, which generate a splitting in the v 1 (η) of D mesons: at forward rapidity D 0 (solid green line) has a lower directed flow than its antiparticle (dashed orange line) and the opposite trend is seen at backward rapidity. This pattern is the same of that found for RHIC collisions, indicating again that the contribution of the electric field E x in the Lorentz force dominate over that produced by the magnetic field B y . Though with big uncertainties, the ALICE data indicate a positive slope for the combined directed flow of D 0 and D 0 that is smaller then the one observed at RHIC and in agreement with the predictions of our approach. However it is certainly required to have experimental data with much higher precision to draw any conclusions.
Besides the slope of the v 1 (η) of the two particles that can be connected also to the vorticity of the fireball, the measurements of ALICE Collaboration for the split-  Fig. 17 for simulations with impact parameter b = 7.5 fm: we found that ∆v D 1 is negative at forward rapidity and positive at backward rapidity while the ALICE experimental data have the opposite trend. The same disagreement in sign has been found by other theoretical calculations for the charge-odd directed flow of light [39] and heavy mesons [40]. Not only the sign of the splitting, but also its magnitude in absolute value measured by ALICE cannot be explained by our model: we found a |∆v D 1 | that is two order of magnitude smaller than the experimental data. Moreover, from Fig. 17 we see a comparable splitting in both simulations with and without vorticity, indicating that the ∆v D 1 is not affected by the initialization used for describing the effect of the angular momentum of the system.
The understanding of the ALICE measurements of the charge-dependent directed flow of D mesons represents a challenge from the theoretical point of view, being them is contrast with model calculations but also with theoretical and experimental results at lower center-of-mass energy. Possible explanations may deal with the very early time dynamics of the collision, that is the less understood stage especially at very high energy. Even the evolution of the electromagnetic themselves is not clear in the first fractions of fm/c when the system could be constituted mainly of color fields or strings and the QGP matter that sustain the electromagnetic field evolution with its electrical conductivity may be not present yet. But charm quark and antiquarks, produced in initial hard scatterings, would be witnesses of the early dynamics and of the vortical and electromagnetic fields pertinent to it. In order to see that even a mild modification of the latter could give an insight into how the splitting of D 0 and D 0 directed flow could change sign, we modified by hand the B y component of the electromagnetic field, rising by 25% its value in each space-time cell, and keeping unchanged the values of the electric field. The outcome of this trial is represented by the dot-dashed brown line in Fig. 17 for the splitting of the D meson directed flow at LHC energy. The effect of the increasing the magnetic field is not only to flip the sign of ∆v D 1 , going in the right direction to reproduce the trend indicated by the ALICE data, but also to generated a |∆v D 1 | even larger. This means that a much detailed investigation of the spacetime profile of the electromagnetic field has to be performed. Recently in Ref. [42] in a schematic calculation within the Langevin approach that neglected the vortical dynamics and the tilted longitudinal distribution it has been shown that a magnetic field evolving slowly with respect to to modellings of a medium in equilibrium at constant electric conductivity can lead to a sign change of ∆v D 1 and strength comparable to the one observed at ALICE.

VII. CONCLUSIONS
We have presented a study of the build-up of the directed flow of charm quarks at both RHIC and LHC energy. The predictions have been obtained by means of realistic simulation of the uRHICs based on relativistic Boltzmann transport equation evolving on initial conditions that take into account both the vorticity of the created matter and its tilted longitudinal distribution as well as the action of the electromagnetic fields.
The origin of the directed flow for the D meson appears to be different with respect to the one of the bulk matter. It is, in fact, not coming from the quite weak effect of rotation induced by the initial angular momentum of the spectators transferred to the created bulk matter. The surprising large v 1 (y) observed for D 0 and D 0 mesons, about 30 times larger than the one of the light charged mesons, originates from the pressure push of the bulk matter in the transverse direction. This is possible only if the peak of the bulk matter is tilted with respect to the x−axis and to the charm distribution.We have pointed out that such a mechanism is effective only if the interaction of HQ in the hot QCD matter is non-perturbative, otherwise the bulk matter would not be able to transfer such a push in the transverse direction to the charm quarks. In particular, we have show that assuming perturbative QCD drag and diffusion for the charm would lead to a v 1 (y) at least one order of magnitude smaller than the one observed by STAR. We find the important result to have a quite good prediction for the v 1 (y) of D mesons once the in medium HQ quark interaction able to correctly predict both R AA and v 2,3 (p T ) at RHIC and LHC is considered and the bulk longitudinal distribution is tuned to correctly predict the v 1 (y) for charged particles. The absence of only one of these two features, tilted bulk and non-perturbative HQ interaction, would lead to a v 1 (y) much smaller than the one observed and comparable with the one of the charged particles. Therefore, the very large v 1 (y) is indeed a further confirmation of the non-perturbative dynamics and in the future more precise measurement may contribute to the determination of the space-diffusion coefficient, in particular by Bayesian analysis [57]. We have also pointed out that the formation time of v 1 (y) for the charm quarks is quite short being of about only 2 fm/c. This should imply that it is more sensitive to the HQ transport coefficient at the initial high temperature phase providing complementary information with respect to v 2 that is known to be more sensitive to the interaction at lower temperature close to T c , as discussed in [95]. The other aspect studied is the splitting of v 1 (y) for D 0 and D 0 due to the electromagnetic fields. We have found that this is consistent with the prediction at RHIC energy, but due to the current experimental data error bars would also be consistent with zero. On the other hand the splitting at LHC energy appears to be quite large experimentally even if with still large error bars, but at present the theoretical models based on electromagnetic field profile assuming a medium in equilibrium at constant electric conductivity is not able to account for it. However, it is likely that at LHC energies such an assumption is significantly inappropriate. In this respect it has to be considered that modellings like the one employed here and in most of the current literature predicts a quite similar magnetic fields at both RHIC and LHC, see Fig. 1. On the other hand at LHC energy the initial maximum value expected in the vacuum is about a factor 50 larger than the one in a medium at the assumed σ el , while such a gap is instead only about a factor of 3 − 4 at RHIC energy. It is likely that at least a study of the electromagnetic field under a medium with a quickly varying electric conductivity is necessary, as well as, a study of the possible effects of anomalous chiral conductivity.