Ballistic diffusion of heavy quarks in the early stage of relativistic heavy ion collisions at RHIC and LHC

We study the diffusion of charm quarks in the early stage of high energy nuclear collisions at the RHIC and the LHC. The main novelty of the present study is the introduction of the color current carried by the heavy quarks that propagate in the evolving Glasma (Ev-Glasma), that is responsible of the energy loss via polarization of the medium. We compute the transverse momentum broadening, $\sigma_p$, of charm in the pre-thermalization stage, and the impact of the diffusion on the nuclear modification factor in nucleus-nucleus collisions. The net effect of energy loss is marginal in the pre-thermalization stage. The study is completed by the calculation of coordinate spreading, $\sigma_x$, and by a comparison with Langevin dynamics. $\sigma_p$ in Ev-Glasma overshoots the result of standard Langevin dynamics at the end of the pre-hydro regime. We interpret this as a result of memory of the color force acting on the charm quarks that implies $\sigma_p\propto t^2$. Moreover, $\sigma_x\propto t^2 $ in the pre-hydro stage shows that the charm quark in the Ev-Glasma is in the regime of ballistic diffusion.

Heavy quarks, charm and beauty, are good probes of the system created in high energy nuclear collisions, both for the pre-equilibrium part and for the thermalized quark-gluon plasma (QGP), see  and references therein. In fact, their formation time is τ form ≈ 1/(2m) with m the quark mass which gives τ form ≤ 0.07 fm/c, which is shorter than the formation time of light quarks; thus, charm and beauty propagate in the Ev-Glasma and probe its evolution. The large mass, the early production and the low concentration of charm and beauty makes heavy quarks the perfect probes of the medium produced * Electronic address: ruggieri@lzu.edu.cn in the collisions.
The evolution of heavy quarks in the Ev-Glasma has attracted a lot of interest recently, [45][46][47][48][49][50][51]. The study of diffusion and energy loss was started in [45], while in [46,49,51] it was shown for the first time how the evolution of charm and beauty in the Ev-Glasma can affect the nuclear modification factor, R AA , and the elliptic flow of these quarks both in proton-nucleus and in nucleusnucleus collisions. More studies have been devoted to the diffusion and momentum broadening [47,48,50].
In previous studies, the energy loss of the heavy quarks has been neglected: this is due to the polarization of the gluon medium induced by the color current of the heavy quarks [56][57][58][59], see also [60] for a treatment of the problem within a classical model. Adding this current results in a drag force acting on the heavy quarks: it is thus a back-reaction. Neglecting this sounds as a reasonable approximation: in fact, the energy density of the Ev-Glasma is much larger than that of the QGP phase, therefore the momentum broadening due to diffusion is expected to be more important than the energy loss due to polarization. Nevertheless, adding the color current is a well defined procedure, therefore it is possible to add it and compute its effect on the observables. This is one of the main purpose of the present study.
We compute the effect of the color current on the nuclear modification factor of charm in the Ev-Glasma, on the momentum broadening and on the transverse coordinates diffusion. This is achieved by adding this current to the CYM equations, that we solve consistently with the kinetic equations of motion of the heavy quarks. We find that the effect on the R AA is present but not large: this is due both to the small magnitude of the current of charm, arXiv:2011.05818v1 [hep-ph] 11 Nov 2020 as well as to the short lifetime of the pre-thermalization stage. The effect of the drag force is more visible on the momentum and transverse plane coordinates diffusion: however, even for these quantities the net effect of the drag is to slow down the diffusion of at most 20%.
We estimate both the diffusion coefficient in momentum space, D, and the drag coefficient, γ: in particular, we find that γ is quite small, certainly smaller than the value it is usually used in the QGP phase. This means that the equilibration time of the charm in the Ev-Glasma stage is much larger than the lifetime of the Ev-Glasma itself. Thus, the motion of charm in the Ev-Glasma is dominated by diffusion because the equilibration time, τ therm = 1/γ, is much larger than the lifetime of the Ev-Glasma, as it happens for the standard Brownian motion.
We remark that this addition solves the classical problem completely and consistently: this procedure adds the classical radiation produced by the moving heavy quarks in a consistent way, and is qualitatively similar to what one should do in classical electrodynamics for the problem of the propagation of a classical electric charge in a classical electromagnetic field. In fact, it is well known that approaching this classical electrodynamics problem leads to the production of a near and a far fields, the latter being responsible of energy loss by radiation. In solving the classical problem for the heavy quarks, we clearly ignore the quantum processes and in particular the hard gluon emission by bremsstrahlung: these processes might be introduced by adding a random force in the equations of motion of the heavy quarks, but it is known that these would contribute only in a range of p T way larger than the one that we consider here.
We also study in detail the momentum broadening, σ p = (p(t) − p 0 ) 2 , where p 0 denotes the initial value of momentum of charm. For a standard Brownian motion with uncorrelated noise σ p = 2Dt for t 1/γ, while the drag affects later evolution. For charm we find that σ p ∝ t 2 in the very early part of the evolution: we interpret this as the effect of the memory in the correlator of the force exerted by the gluon fields on the charm; in fact, for a Brownian motion with a nontrivial memory kernel σ p ∝ t 2 . In the case of charm in the Ev-Glasma this can be interpreted as a time correlation of the force, F , acting on the charm, , where x(τ ) denotes the position of the charm at time τ 1 . Although we do not compute the correlator of the force since it would require a different approach to the solution of the CYM equations, we have computed the correlator of the electric field at different times and found that this is characterized by a finite time decay, suggesting finite time correlation of the force. Comparison with the σ p of a standard, uncorrelated Brownian motion we find that the effect of memory is to slow down the diffusion in the very early stage, but after a short transient σ p in the Ev-Glasma overshoots the one of the uncorrelated motion.
We complete the study by computing diffusion in the transverse coordinate space. We find that σ x = (x(t) − x 0 ) 2 follows the qualitative path of a Brownian motion in its early stage, σ x = at 2 + bt 3 . In this stage memory plays a less relevant role for σ x because it would affect only terms of order O(t 4 ) which are smaller than the O(t 2 ). The σ x ∝ t 2 shows that the diffusion of charm in Ev-Glasma is in a ballistic regime [61,62].
The plan of the article is as follows: in Section II we review the theoretical setup of the work; in Section III we review briefly the solution of the Langevin equations for the Brownian motion, emphasizing the effect of a memory kernel on the early evolution of σ p ; in Section IV we present our results on σ p , σ x and R AA of charm; in Section V we compare σ p of charm in the Ev-Glasma and in a thermal medium and present a qualitative comparison of the R AA in the two cases; finally in Section VI we summarize our work and discuss possible future improvements.

A. Glasma and classical Yang-Mills equations
In this section, we review the Glasma and the McLerran-Venugopalan (MV) model [1][2][3]63]. In this work we scale the gauge fields as A µ → A µ /g where g is the QCD coupling. In the MV model for the collision of two nuclei labeled as A and B the static color charge densities ρ a on A and B are assumed to be random variables that are normally distributed with zero mean and variance given by here, a and b denote the adjoint color index; in this work we consider the case of the SU (2) color group therefore a, b = 1, 2, 3. In Eq. (1) g 2 µ A denotes the color charge density, g 2 µ = O(Q s ) [64]. For protons, the dependence of Q s of the average x = p T / √ s can be estimated via the GBW fit [65][66][67][68], Q 2 s = Q 2 0 (x 0 /x) λ , with λ = 0.277, Q 0 = 1 GeV and x 0 = 4.1 × 10 −5 . For nuclei we borrow the modification of the GBW fit obtained within the IP-Sat model [69], namely Other forms of the generalized GBW fit are possible [70,71], but these do not lead to significant changes of Q s . Using the numerical result Q s /g 2 µ = 0.57 of [64] we find g 2 µ Pb = 3.4 GeV for the Pb nucleus for collisions at the LHC energy, or Q s = 1.9 GeV. The static color sources {ρ} generate pure gauge fields outside and on the light cone, which in the forward light cone combine and give the initial Glasma fields. In order to determine these fields we solve the Poisson equations for the gauge potentials generated by ρ A and ρ B , namely Wilson lines are computed as V † (x T ) = e iΛ (A) (x T ) , W † (x T ) = e iΛ (B) (x T ) , and the pure gauge fields of the two colliding nuclei are given by α In terms of these fields the solution of the CYM in the forward light cone at initial time, namely the Glasma gauge potential, can be written as for i = x, y and A z = 0, and the Glasma fields are [8,9] where z is the direction of the collision and the transverse fields vanish. The evolution of the initial condition is achieved via the classical Yang-Mills (CYM) equations, namely we have put where f abc = ε abc with ε 123 = +1, and the standard summation convention has been used. We name the evolving field as the Ev-Glasma, leaving the name Glasma to the initial condition. At variance with previous calculations [45][46][47][48][49][50][51] we include the color current, j a i , carried by the heavy quarks. This is essential to describe the energy loss of the colored particles interacting with the evolving Glasma [56][57][58][59], due to the polarization of the medium.
The lack j a i in previous calculations gives a purely diffusive motion of heavy quarks, and interaction with the gluon fields leads merely to momentum broadening. Instead, adding the color current and solving consistently the field equations and the kinetic equations of the heavy quarks, see below, we take into account both momentum broadening and energy loss. Our solution of the problem is purely numerical, however we do not rely on any assumption on equilibration and thermalization of both the gluon medium and the heavy quarks, as well as we do not assume linear response theory and do not make any assumption on the trajectories and velocities of the heavy quarks. This approach solves the problem of propagation of heavy quarks in the Ev-Glasma completely, as far as classical field theory can do; this solution is similar to Electrodynamics, in which solving consistently the Maxwell equations with the kinetic equations for the charges gives both near and far fields produced by the charges themselves; in particular, the far fields are responsible of electromagnetic radiation.

B. Wong equations for heavy quarks
The dynamics of charm quarks in the Ev-Glasma is studied by the Wong equations [46,47,49,51,72,73], that for a single quark can be written as where i = x, y, z; these correspond to the Hamilton equations of motion for the coordinate and its conjugate momentum, while the third equation corresponds to a classical description of the conservation of the color current.
In the third Wong equations, Q a with a = 1, . . . , N 2 c −1 corresponds to an effective color charge of quarks; this has not be confused with the standard QCD color charge of quarks, because quarks sit in the fundamental representation of the color group SU (N c ) so they carry N c colors, while Q a has the adjoint color index. This effective charge can be understood as a function that allows to describe classically the color current carried by the heavy quarks, namely [26,27,74] in the continuum limit, or in the discretized version, where the sum is understood over all particles present in a given lattice cell. Notice that the current is multiplied by the squared of the QCD coupling, g 2 : this is simply due to the scaling of the gluon fields used for the CYM equations mentioned in Section II and can be verified easily starting from the QCD lagrangian. We initialize Q a by a uniform distribution with support in the range (−1, +1). For each heavy quark we produce an antiquark as well: for this, we assume the same initial position of the companion quark, opposite momentum and a random color charge. Solving the Wong equations is equivalent to solve the Boltzmann-Vlasov equations for a collisionless plasma made of heavy quarks, which propagate in the Ev-Glasma We enlarge the number of heavy quarks by N p test particles to improve statistics: this amounts to replace g 2 → g 2 /N p in Eq. (13). Heavy quarks are initialized at time τ form = 1/(2m). In calculations based on relativistic transport the heavy quarks are assumed to do a free streaming between their formation time and the initialization of the quark-gluon plasma phase, see [75] and references therein; we do not have a free streaming period and the heavy quarks are formed exactly at their formation time and interact immediately with the gluon background.

III. A QUICK REMINDER ON THE DIFFUSION IN THE BROWNIAN MOTION
In this section, we review briefly the classical Brownian motion in one spatial dimension: this reminder is useful to fix a few key results, that allow to understand better those that we obtain for the motion of the heavy quarks in the Ev-Glasma. For simplicity we study only the case of a nonrelativistic particle: results about momentum broadening are valid also in the relativistic case. In order to highlight the most important characteristics of the Brownian motion we make several assumptions along the way: these assumptions have illustrative purposes and are not done in the full solution of the problem for charm presented in the next section.
Brownian motion is the motion of a heavy particle, with mass M , in a bath that interacts with the particle via a time-dependent random force, ξ(t), plus a viscous force, The viscous force is necessary for the equilibration of the heavy particle with the medium. The random force ξ is assumed to satisfy ξ(t) = 0 and in the simplest case, the time correlations of the random force are neglected so it is assumed that namely the motion of the heavy particle is a Markov process. On the other hand, for the very early stage of the propagation of heavy quarks in the Ev-Glasma it is useful to introduce memory effects, namely assume that the time correlator of the random force driven by the gluon fields vanishes only if |t 1 − t 2 | τ mem . The evolution of momentum of the heavy particle is governed by the equation For the discussion we assume τ mem τ therm , which is satisfied by heavy quarks in the Ev-Glasma fields. Under this assumption, it is legitimate to consider Eq. (17) in three limits, namely t τ mem that we call the very early stage, τ mem t τ therm that we call the preequilibrium stage and τ therm t that we call the equilibrium stage.
Firstly, we focus on the very early stage. The general solution of Eq. (17) can be written in terms of Laplace transforms [60] and will be presented elsewhere; for the purpose of the present study it suffices to say that for t τ mem τ therm the drag force can be neglected and, putting σ p ≡ (p(t) − P 0 ) 2 with P 0 = p(t = 0) we have where we have used ξ(t) = 0, P 0 ξ(t) = 0 as well as Eq. (15). In the very early stage we can expand the right hand side of the above equation in powers of t/τ mem , namely F (t) ≈ F (0)t + F (0)t 2 /2 with each prime denoting a derivative with respect to t, and We notice that if f (t − t 1 ) has not singularities then F (0) = 0. The only case in which F (0) = 0 is when f (t − t 1 ) is singular, for example for t 1 = t: this happens in particular for the Markov processes. For the case of heavy quarks in Ev-Glasma the random force is related to the correlators of the electric and magnetic color fields that are not singular [45] so it is meaningful to study this case.
For the sake of illustration we assume normalized as We notice that σ p ∝ t 2 in the very early stage. Clearly, changing Eq. (21) will not change Eq. (23) modulo an overall constant factor as long as the correlator of the force is a regular function. Assuming such regular correlator we can write, in general, We will use the above result in Section IV A. For the pre-thermalization and equilibrium stages the time range is much larger than the decay time of the memory and we can effectively assume Eq. (16) instead of Eq. (21) to evaluate σ p ; moreover, in this case, assuming for simplicity γ(t) = 2γδ(t), Eq. (17) is replaced by The solution of Eq. (17) is given by where P 0 = p(t = 0). After a straightforward calculation we get which represents the standard momentum broadening of a particle subject to a Brownian motion without memory. From the above equation we get, in the prethermalization and equilibrium stages, and where means that the quantity on the left tends asymptotically to the one on the right in the large time limit. Notice that Eqs. (27) and (28) are valid also in the relativistic limit since no assumption has been done on the relation between energy and momentum. In particular, Eq. (28) implies that p 2 (t) D/γ. It can easily be verified that the evolution of σ p at small and large times in Eqs. (27) and (28) agrees with that obtained by the solution of the one-dimensional Fokker-Planck equation solved with a δ−function initial condition [52].
The physical meaning of Eqs. (23), (27) and (28) is that for times much smaller than the memory time, momentum spreads quadratically with time until the heavy particle enters the pre-thermalization stage, with linear momentum spreading. In both these regimes the motion is dominated by the diffusion, while the drag force appears only to higher orders in time and affects the later motion of the heavy particle. On the other hand, for times much larger than the thermalization time, the heavy particle equilibrates with the medium, as a resut of the balance of the momentum spreading given by ξ(t) and the energy loss driven by f drag . From we notice the last addendum on the right hand side of the above equation, that gives the characteristic linear growing of σ x at large times. Once again, it is convenient to consider the limits of small and large times; for the former we get and for the large times which is the classic result of the Brownian motion. Considering the very early stage in which memory is effective adds a term O(t 4 ) to Eq. (30) which is less important than the ballistic term O(t 2 ). In plain words, Eqs. (30) and (31) state that for times much smaller than the thermalization time the heavy particle experiences an accelerated motion due to the random force; this motion is gradually slowed down by the friction and waiting enough time, the balance between friction and random force will lead to the linear spreading of the position of the particle.

IV. RESULTS
A. Diffusion in the transverse momentum space: toy model initial condition To begin with, we prepare a δ−like initialization in transverse momentum, p T , and study the evolution of the distribution function, dN/dp T , and of momentum and energy of the charm quark. For this toy model initalization we use n c = 15 heavy quarks, that roughly corresponds to the number of charm quarks produced in Pb-Pb collisions at midrapidity at the LHC energies [76]. We show the results for charm quarks only, since they look very similar for the case of beauty quarks.
In Fig. 1 we plot dN/dp T at t = 0.2 fm/c (blue lines), t = 0.6 fm/c (orange lines) and at t = 1 fm/c (green lines); the solid lines correspond to calculations with the color current while the dashed line to those without the current. The upper panel corresponds for an initial p T = 0.5 GeV while in the lower panel the initial p T = 5 GeV; in all calculations we initialize the quarks with momentum p z = 0. The effect of the backreaction, due to the color current, on the charm quark is evident in the small p T case: in fact, the evolution of dN/dp T when the current is taken into account is slower with respect to the case in which the current is not considered, as expected by a drag force; this is seen by the naked eye by examining both the evolution of the peak value and the broadening of dN/dp T . When we consider larger values of p T , we find that the effect of the drag force is not strong.
In Fig. 2 we plot the momentum variance of charm quarks, σ p , versus time, for three values of the initial p T , where we have put with p 0x , p 0y denoting the x, y components of the initial transverse momentum and p 2 T = p 2 x + p 2 y . The solid lines in Fig. 2 denote the results with current taken in to account, dashed lines correspond to calculations without current. Results correspond to g 2 µ = 3.4 GeV. We notice that at small p T the effect of the drag force is quite large, lowering the momentum broadening of ≈ 40% after t = 1 fm/c of evolution in the gluon field; the effect of the current becomes smaller for larger values of p T . However, for the typical lifetime of the pre-hydro stage in nuclear collisions at the LHC energy, τ ≈ 0.3 fm/c [76], we find that even for small p T including the drag force does not affect the σ p substantially: for example, for p 0 = 0.5 GeV we find that the effect of the current is to lower σ p of ≈ 13%.
We notice that the evolution of σ p (t) is not linear in the whole time range considered, see in particular the early time behavior of p 0 = 0.5 GeV in the lower panel of Fig. 2 which is the one more relevant for the role of Ev-Glasma in the early stage of relativistic heavy ion collisions. This nonlinearity be understood as the effect of memory in the very early stage of the evolution of charm in the gluon field, similarly to the early stage of the non-Markovian Brownian motion discussed in Section III for which σ p ∝ t 2 . The calculation of the correlators of the force is beyond the purpose of this article, however to check the plausibility of memory effects in the Ev-Glasma we have computed the correlator of the electric field at different times and found a decay time ≈ 0.06 fm/c ≈ 1/g 2 µ. Using Eq. (18) to fit the data in the lower panel of Fig. 2 we estimate τ mem ≈ 0.07 fm/c, in agreement with the result of the correlator.
We estimate the diffusion and drag coefficients of charm by fitting the data in Fig. 2 with Eq. (26) up to t = 2 fm/c, starting from t = 0.2 fm/c to remove the early stage that is dominated by the memory. We get D = 3.37 GeV 2 /fm and γ = 0.026 fm −1 for p 0 = 0.5 GeV; we use γ to estimate the thermalization time of the charm in the Glasma, namely τ therm = 1/γ ≈ 38 fm/c.
In Fig. 2 we compare the results with those obtained by solving the standard Langevin equations without memory Eq. (26), with the values of D and γ that we get for charm in the Ev-Glasma; we represent the data with blue dot-dashed lines. We notice that initially the ∝ t 2 of charm in the gluon fields makes momentum broadening slower than the corresponding Markovian dynamics. On the other hand, for t τ mem the broadening in the gluon fields overshoots the Langevin results and gives a faster diffusion in momentum space. We present more comparisons with the Langevin dynamics in Section V.
We can summarize our findings by writing that if we had to follow the diffusion of charm in the Ev-Glasma for times up to ≈ 1fm/c, then this would appear largely as a standard Brownian motion with drag and diffusion, however mostly dominated by diffusion since equilibration time is quite larger than 1 fm/c. On the other hand, limiting to consider the timescales ≈ 0.3 fm/c which are those relevant for the relativistic heavy ion collisions, the memory is qualitatively important as it slows down the momentum broadening of the charm quarks in the very early stage, then gives a boost and puts σ p above the result we would measure if the diffusion was a Markov process. Having added the drag force by the color current, we have found that the net effect of the drag in this short time range is quite modest.

B. Diffusion in the transverse momentum space: realistic initial condition
Next we turn to a realistic initialization of heavy quarks in transverse momentum space. To this end, at the formation time we assume the prompt spectrum obtained within Fixed Order + Next-to-Leading Log (FONLL) QCD that reproduces the D-mesons spectra in pp collisions after fragmentation [77][78][79] the parameters that we use in the calculations are x 0 = 20.2837, x 1 = 1.95061, x 2 = 3.13695 and x 3 = 0.0751663 for charm quarks; the slope of the spectrum has been calibrated to a collision at √ s = 5.02 TeV. We use n c = 15 charm quarks as in the δ−function initializations, corresponding to the estimated number of charm quarks produced in Pb-Pb collisions at midrapidity at the LHC energies [76]. Moreover, we assume that the initial longitudinal momentum vanishes. In coordinate space, for the setup of AA collisions, we simulate the most central interaction region in which we assume that width of the random color density fluctuations, given by g 2 µ, is constant: as a consequence, we assume that the probability of formation of heavy quarks in this region is also uniform. We can quantify the effect of the propagation of charm quarks in the Ev-Glasma by introducing the modification factor, R AA , defined as where the prompt spectrum is given by Eq. (33) and dN/d 2 p T evolved corresponds to the spectrum after the evolution in the Glasma fields: this is a time dependent quantity so in general R AA depends on time as well. If R AA = 1 for all values of p T it means that the spectrum after the evolution is the one computed from hard scatterings in pQCD; on the other hand, R AA = 1 signals the interaction of the charm with the medium.
In Fig. 3 we plot the nuclear modification factor for c−quarks at three different times, computed with and without the color current in the YM equations. The tilting of the spectrum discussed above naturally results in R AA smaller than one at low p T , as a result of the diffusion of these charms to higher p T ; the larger the time of the evolution in the gluon field, the larger the effect on R AA . The drag force induced by the polarization of the medium slows down the evolution of the spectrum, and it naturally results in a slower evolution of R AA as well.
The results in Fig. 3 agree qualitatively with those presented in [46,49]: the main novelty of the present work is to upgrade those results taking into account the drag force that results from the polarization of the medium induced by the color current carried by the quarks. Quantitatively, we find that at up to t = 0.3 fm/c the effect of the drag force is negligible, while it is substantial at t = 1 fm/c. The typical lifetime of the pre-hydro stage is τ ≈ 0.3 fm/c for collisions at the LHC energies [76], therefore the results of the present study suggest that the inclusion of the color current of the charm quarks will not affect drastically the evolution of the spectrum at the LHC energies in comparison to the results published in [46,49].
The qualitative shape of R AA that we get at the end of the pre-thermalization stage is different from the one that is usually found after the evolution in the quark-gluon plasma, see [38,80] and references therein: there, it is evident the diffusion of the large p T charm quarks to lower p T states due to energy loss. It should be noted that in the present calculation the energy density of the bulk is quite larger than the one in the quark-gluon plasma phase, therefore in the latter the effect of the drag force will be larger than the one we have found here. Moreover, the drag coefficient of charm in the quark-gluon plasma phase is larger than the one we have found in the Ev-Glasma. Both these factor make the motion of charm in the pre-thermalization stage gluon fields dominated by diffusion and low-p T flow to higher p T .

C. Diffusion of charm in the transverse coordinate space
In Fig. 4 we plot the transverse coordinate variance of charm quarks, σ x , versus time, where for the cases with and without the color current in the YM equations. The calculation setup corresponds to that of Fig. 4. Energy loss slows down the diffusion since σ x is bent downwards when the color current is introduced. Once again, the effect of energy loss is quite modest for the very early times up to t ≈ 0.3 fm/c, while it becomes more substantial for larger times.
Combining the results of this and previous subsections, we conclude that the motion of charm in the prethermalization stage is that of a ballistic diffusion. In fact, equilibration time is much larger than the lifetime of the pre-hydro stage, which gives σ x ∝ t 2 in the time range of interest, and σ p ∝ t 2 in the same time range overshooting the standard diffusion for which σ p ∝ t.

V. COMPARISON WITH LANGEVIN DYNAMICS
In this section, we compare the evolution of charm quarks in the Ev-Glasma, with that obtained by solving standard Langevin equations without memory. To facilitate the comparison we use the same diffusion coefficient in all calculations, namely D = 3.37 GeV 2 /fm that matches what we estimated in Sections IV A and IV B.
In Fig. 5 we plot σ p versus time for the initialization p 0 = 0.5 GeV: Langevin 1 uses the same γ of Ev-Glasma namely γ = 0.026 fm −1 , while in Langevin 2 and 3 we have used the γ that would be required by the Fluctuation-Dissipation theorem, γ = D/ET , for T = 1 GeV and T = 1.5 GeV respectively: these are γ = 1.72 fm −1 in Langevin 2 and γ = 2.59 fm −1 in Langevin 3. In the cases Langevin 2 and 3 it is obvious that charm quarks equilibrate with the thermal medium within ≈ 1 fm/c. In the lower panel of Fig. 5 we zoom on the very early stage of the evolution, up to t = 0.2 fm/c: we notice that σ p in all cases is different both qualitatively and quantitatively from the one in the Ev-Glasma. In particular, the effect of memory is clearly visible in the Ev-Glasma for σ p ∝ t 2 rather than ∝ t up to t ≈ τ mem , and σ p in the Ev-Glasma overshoots that in the Langevin dynamics for t τ mem .

VI. CONCLUSIONS AND OUTLOOK
We have studied the diffusion of charm quarks in the Ev-Glasma produced in high energy nucleus-nucleus collisions. We have solved consistently the Yang-Mills equations for the evolution of the gluon field and the Wong equations for the heavy quarks. The main novelty of this study is the inclusion of the color current carried by heavy quarks in the classical Yang-Mills equations of the gluon field, that is necessary to describe the energy loss of heavy quarks due to the polarization of the medium [56][57][58][59]. This study concludes the one started in [46,49,51], in which the phenomenological impact of the diffusion of heavy quarks in Glasma have been studied for the first time. There, the idea that diffusion in the early stage affects the nuclear modification factor as well as the elliptic flow of heavy quarks, despite the short lifetime of the pre-hydro stage, was investigated, but the calculations neglected the color current carried by the heavy quarks and the subsequent backreaction on the motion of heavy quarks themselves. We fill this gap here.
Qualitatively, our results agree with the expectation that the energy loss slows down the momentum broadening of charm and beauty in the Ev-Glasma. This affects the evolution of the nuclear modification factor in the pre-hydro stage: however, taking into account energy loss we get R AA that is consistent, within the 10%, with results previously computed without the current. This modest effect can be understood easily because the current carried by charm and beauty is very tiny due to the low density of these quarks in the initial stage. Therefore, we confirm the previous findings [46,49,51] that the diffusion of the heavy quarks in Glasma is responsible of a tilt of their spectrum, effectively moving low p T quarks to higher p T . We have achieved this conclusion by studying δ−function initializations as well as realistic p T −initializations, looking at both the momentum broadening and the R AA of charm.
We have investigated more closely the motion of the heavy quarks. We have found that overall the diffusionwith-drag can be interpreted in terms of the Brownian motion at late times, plus a motion with memory effects in the very early stage of the evolution. We achieve this by studying the p T −broadening versus time, σ p (t), and identify an initial range in which σ p ∝ t 2 , interpreting this as an effect of memory related to the finite time width of the correlators of the electric and magnetic color fields. This non-Markovian regime lasts in the very early stage of the evolution and is dominated by diffusion, and is followed by a standard Brownian motion regime with drag and diffusion; the net effect of the drag is however small because the lifetime of the prehydro stage is smaller than the thermalization time of heavy quarks and in this case the leading contribution to momentum broadening comes from diffusion. Significant effects of drag have been found at later times, t ≈ 1 fm/c; however, these times are well beyond the lifetime of the pre-equilibrium stage of high energy nucleus-nucleus collisions, therefore it cannot affect any observable. For the time range in which the Ev-Glasma can play a role in collisions the diffusion in the early stage is the relevant one, in which the effect of memory is important. In the very early stage we have found some quantitative difference with the diffusion of a standard Brownian motion studied via Langevin equations without a memory kernel. In particular, momentum broadening in the Ev-Glasma proceeds slower than the linear increase of the standard Brownian motion, then overshoots the latter: at the end of the Ev-Glasma evolution, τ ≈ 0.4 fm/c, the σ p that we get from Ev-Glasma is larger than the one we would obtain by Langevin dynamics with same drag and diffusion coefficients.
We have also studied the diffusion of charm in coordinate space. We have found that coordinate broadening, σ x , evolves as σ x ∝ at 2 + bt 3 hence faster than the steady state Brownian motion result σ x ∝ t. This faster diffusion in coordinate space means that the diffusion of charm in the Ev-Glasma is in the ballistic regime.
In conclusions, our findings support the diffusion-withno-drag advertised in [46,49,51], and allow to understand it easily: the time range relevant for the propagation of the heavy quarks in the Ev-Glasma is much smaller than the equilibration time, and in this regime the motion is diffusion-dominated; a short transient where memory is important is replaced by a standard Brownian motion at later times. The drag, that we have computed self-consistently in this study, does not affect in a considerable way the observables that we have studied, in particular the nuclear modification factor, because substantial energy loss is effective only on time scales comparable with the thermalization time.
While this study answers the question whether energy loss is an important ingredient to study the diffusion of heavy quarks in the Ev-Glasma, it opens up other questions. The role of fluctuations in the initial stage should be considered: it is well known that fluctuations enhance isotropization already in the initial condition of Glasma [24,33], therefore it is interesting to study how the a larger amount of isotropization affects the evolution of the heavy quarks. In additon to this, it is important to focus on phenomenological calculations aimed to compute the impact of the early stage diffusion on observables, mostly hadron spectra, two-bodies correlations and collective flows, both in proton-nucleus and nucleus-nucleus collisions. Protons can be initialized similarly to nuclei according to the constituent quark model, see [81][82][83][84][85] and references therein. We have not included the cold nuclear matter effects [54,55,[86][87][88][89][90][91] in the initialization of the charm quarks, therefore adding them is a further improvement of the present work. We will report on these subjects in the near future.