Impact of Glasma on heavy quark observables in nucleus-nucleus collisions at LHC

In the pre-thermal equilibrium stage of relativistic heavy-ion collisions, a strong quasi-classical transverse gluon field emerges at about $\tau_0 \simeq 0.1 \, \rm fm/c$ and evolves together with their longitudinal counterparts according to the classical Yang-Mills (CYM) equations. Recently it has been shown that these fields induce a diffusion of charm quarks in momentum space resulting in a tilt of their spectrum without a significant drag. We find that in nucleus-nucleus collisions at LHC such a novel dynamics of charm quarks leads to an initial enhancement of the nuclear modification factor ($R_{AA}$) at $p_T$ larger than 2 GeV$/c$ contrary to the standard lore. Moreover, the same dynamics leads to a larger final elliptic flow ($v_2$) inducing a relation between $R_{AA}$ and $v_2$ that is quite close to the experimental measurements. Our study also shows that such an initial pre-thermal stage is unlikely to be described in terms of a standard drag and diffusion dynamics, because even if one tune such coefficients to reproduce the same $R_{AA}(p_T)$ this would imply a significantly smaller $v_2$.

In the pre-thermal equilibrium stage of relativistic heavy-ion collisions, a strong quasi-classical transverse gluon field emerges at about τ0 ≃ 0.1 fm/c and evolves together with their longitudinal counterparts according to the classical Yang-Mills (CYM) equations. Recently it has been shown that these fields induce a diffusion of charm quarks in momentum space resulting in a tilt of their spectrum without a significant drag. We find that in nucleus-nucleus collisions at LHC such a novel dynamics of charm quarks leads to an initial enhancement of the nuclear modification factor (RAA) at pT larger than 2 GeV/c contrary to the standard lore. Moreover, the same dynamics leads to a larger final elliptic flow (v2) inducing a relation between RAA and v2 that is quite close to the experimental measurements. Our study also shows that such an initial pre-thermal stage is unlikely to be described in terms of a standard drag and diffusion dynamics, because even if one tune such coefficients to reproduce the same RAA(pT ) this would imply a significantly smaller v2.

I. INTRODUCTION
The relativistic heavy-ion program provides the possibility to scrutiny Quantum Chromodynamics (QCD) phase diagram at finite temperature and density in the region where a transition of nuclear matter into a plasma of quarks and gluons (QGP) has been predicted. In particular, in the last decade the main properties of such a matter have been studied both at the Relativistic Heavy Ion Collider (RHIC) and at the Large Hadron Collider (LHC), clarifying that the high temperature medium created in high energy nuclear collisions is characterized by a large scattering rate and thus a low shear viscosity over entropy density ratio, in disagreement with the naive expectation of a weakly coupled plasma.
Heavy quarks (Charm and Beauty) have been considered to have a unique role in such a study since they are generated in the early stage τ 0 ≃ m −1 HQ < O(10 −1 ) fm/c according to next-to leading order perturbative QCD and hence are witness of the entire evolution of the QGP; furthermore because M HQ >> T c they still preserve their identity at hadronization by picking-up a light quark or undergoing an independent fragmentation. These properties, together with a thermalization time τ therm that is comparable to the lifetime of the QGP phase [1,2], makes them a probe able to preserve key information about the time evolution of their interaction in the hot QCD medium. In addition, it is in perspective possible to perform a direct comparison of the transport properties of the heavy quarks with lattice QCD (lQCD) calculations. Recently, it has been shown that the determination of the space-diffusion transport coefficient D s from the phenomenology [1,3] is in agreement with first lQCD calculations in quenched approximation [4,5], even if within still significant uncertainty both in the phenomenology as well as in the lQCD approach.
The study of the HQ physics in AA collisions has been successful and also has clearly shown that in the low momentum regime p T < 10 GeV the interaction is strongly non-perturbative and implies a space diffusion coefficient 2πT D s ∼ 2 − 5 around T c [1]. The determination of such a coefficient is mainly driven by the phemonemological prediction of two main observables: the nuclear modification factor R AA and the elliptic flow v 2 . However, there has been always a tension between these two observables that are hard to be correctly predicted simultaneously [1,[6][7][8][9][10][11]. A significant part of such a tension is reduced when an increasing temperature dependence of D s and an hadronization by coalescence plus fragmentation are included, as discussed in [12]. However, especially at the LHC energy such a tension persists and is currently partially moderated by the still significant large error bars in the experimental data, especially for the elliptic flow [3,13].
We notice that while the dynamics of the HQs in the QGP phase has been thoroughly studied and also the possible impact of the later stage of hadronic re-scattering has been discussed in several approaches [14][15][16][17][18], the dynamics of HQs in the early stage has been addressed only recently [19]. The early stage dynamics can be quite rele-vant especially for HQ's considering their short formation time and the fact that HQ thermalization time is comparable to the lifetime of the QGP. This, in fact, implies a larger sensitivity to the early time evolution of the HQ during the fireball expansion, hence potentially keeping memory of the initial dynamics. Indeed, recently, it has been shown that HQs dynamics should be particularly sensitive to both the initial electromagnetic field and the tilted initial condition of the medium [20,21], a prediction that appears to be confirmed by early experimental results at both RHIC and LHC energies [22].
According to the the color glass condensate (CGC) effective theory [23][24][25][26], the dense gluon system produced by the interaction of two colored glasses immediately after the collision can be described in terms of classical longitudinal fields named the Glasma, whose evolution in the early stage can be described by means of the classical Yang-Mills (CYM) equations. Recent works have investigated the impact of the propagation of charm and beauty in the Glasma [19,[27][28][29] with particular reference to the diffusion in momentum space. Here we study for the first time how the interaction of the heavy quarks with the Glasma affects the dynamics of HQs in AA collisions; we achieve this by implementing a simulation of the HQs dynamics including both the glasma and the QGP stages. In this preliminary study we have not implemented the longitudinal expansion in the YM evolution: in order to overcome this problem we have used a small value of the saturation scale, in agreement with the lowest bound for the estimate of this scale for AA collisions at the LHC energy [19]. We show that as a result of the diffusion of HQs in the Glasma, there is a significant modification of the relation between R AA and v 2 ; we expect a significant impact also on several other observables, both in pA and AA collisions. Although a systematic study is feasible in perspective, here we report specifically on Pb-Pb collisions at √ s N N = 5.02 TeV in order to emphasize the novel idea.
In this Letter, we firstly discuss how we set up the initial condition and the evolution of HQs in the Glasma by means of the Wong equations [30]; then we discuss how this initial stage evolution is embedded in a modified Fokker-Planck equation that is able to reproduce a very similar dynamics, but in addition allows for the implementation of the subsequent standard HQs evolution in the QGP. We then focus on the impact of the initial stage dynamics on R AA and v 2 estimating the uncertainty on these observables due to the lifetime of the Glasma stage; in addition to this, we also show the difference with standard approaches that neglect the diffusion of HQs in the glasma phase. Finally we draw our conclusions with an outlook for upcoming studies.

II. MODEL SETUP
Within the CGC effective theory the pre-thermal equilibrium stage of the Pb-Pb collisions can be described in terms of strong gluon fields, namely the Glasma, that evolves according to the CYM equations. In the temporal gauge A a 0 = 0, the Hamiltonian density takes the following form [31,32] where As in [19,28] we deal with the SU(2) gauge theory for simplicity, we thus have f abc = ε abc with ε 123 = +1. The CYM equations are which we solve in a static box in three spatial dimensions [19,[31][32][33]. We have adopted periodic boundary conditions at finite time. We neglect the longitudinal expansion in the YM stage, so the natural coordinate system to use is the (t, x, y, z) with t the lab time, (x, y) the transverse plane and z the longitudinal coordinate; invariance along the longitudinal direction is assumed, analogously to the rapidity invariance of the expanding system. The box size is 4 fm×4 fm in the transverse plane with a lattice spacing a=0.04 fm. The initialization of the classical gluon field is given by the standard Mclerran-Venugopalan (MV) model [23][24][25], in which the transverse color charge densities ρ a on the nucleus A (same for B in an AB collision at high energy) are assumed to be distributed with zero mean and variance specified by with a, b = 1, 2, 3 for SU(2). We set g 2 µ A = 3 GeV in this study, in agreement with the estimate of the saturation scale obtained within the IP-Sat model for the relevant LHC energy [19,34,35]. In order to determine these fields we firstly solve the Poisson equations for the gauge potentials generated by the color charge distributions of the nuclei A and B, namely (a similar equation holds for the distribution belonging to B). Wilson lines are computed as , 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 initial longitudinal Glasma fields are while the transverse fields are vanishing. It is useful to remind that the initial color charges give no contribution to the gluon fields at finite time. This is the standard assumption of any calculation based on the effective theory of the color-glass-condensate. These details have been discussed many times in the literature, see for example [36]. Charm and anti-charm quarks are produced in the prethermal equilibrium stage of relativistic heavy ion collisions; the formation time is The equations of motion of these heavy color probes in the pre-equilibrium stage are the Wong equations [30] with i = x, y, z and E = p 2 + m 2 . The conservation of the color charge is guaranteed by the additional equation with the color charge Q a (a =1,2,3) of charm quarks that we initialize randomly in the range (−1, +1) with uniform distribution for each color charge. We study the Pb+Pb collisions at LHC √ s N N = 5.02 TeV at fixed impact parameter b = 9.25 fm in order to simulate the collisions at centrality range 30% − 50%. At the initial time τ f orm = 0.1 fm/c, the position coordinate of charm quarks are distributed according to the binary nucleon-nucleon collisions profile derived from the standard Glauber model, while the initial momentum of the charm quarks are distributed from the one obtained within Fixed Order+Next-to-Leading Log(FONLL) QCD, which we parametrize it as: where the parameters are x 0 = 20.28, x 1 = 1.950, x 2 = 3.136 and x 3 = 0.0751 respectively. HQs are distributed with p z =0 at the initial time, because in our work we are considering the rapidity range around y=0. Starting with this initial conditions for the charm distribution and with a background of evolving Glasma as described above, we follow the evolution of the charm quarks by means of the Wong's equation, Eqs. (10) and (11) until the formation of QGP at τ glasma = 0.3 − 0.5fm/c (which mean a glasma lifetime in the range 0.2 − 0.4 fm/c), where the label "glasma" stands for the time at which this phase ends and the standard hydrodynamical QGP evolution starts.
We have considered such a range of values that represents a typical range for the beginning of the hydrodynamical expansion of the QGP phase at LHC energies. In Fig. 1, we show the results for the R AA (p T ) when we assume that the glasma phase ends at τ glasma = 0.3 fm/c and at τ glasma = 0.5 fm/c. The results shown have been obtained employing 75 configurations of the glasma fields. The glasma is seen to push charm quarks from low to high p T , which is shown by the solid green and solid red lines of Fig. 1 and the R AA increases with p T up to 1.2 for τ glasma = 0.3 fm/c and 1.5 for τ glasma = 0.5 fm/c. Therefore, we find a significant enhancement of the R AA (p T ) at intermediate p T associated to a depletion at low momenta. Such a dynamics is the result of the diffusion of the heavy quarks in the evolving gluon fields as already studied in the context of pA collisions in [28]. Certainly in AA collisions we have never observed such a behavior of R AA (p T ) because the initial stage is followed by a long phase of in-medium charm quark scattering. The aim of the present work is however to point out the impact on final observable of this initial glasma dynamics. To this end it is necessary to merge this initial glasma phase with the subsequent standard evolution of charm quarks in the QGP. The standard method applied for heavy quark (HQ) dynamics in QGP phase is to follow their evolution by means of Fokker-Planck equation, which is usually solved stochastically by the Langevin equations: where dx i and dp i are the change of the coordinate and momentum of each HQ in each time step dt. Γ and C ij are the drag and the covariance matrix related to the diffusion tensor by where P ⊥ ij = δ ij − pipj p 2 and P ij = pipj p 2 are the transverse and longitudinal projector operators respectively. We employ the standard assumption B 0 = B 1 = D for all momenta in this study as Refs. [8,9,[37][38][39], though this is strictly valid only for p → 0. Under this assumption, Eq. (15) becomes simply C ij = √ 2Dδ ij . This is the approach that most of the groups have used to study the heavy quark dynamics in AA collisions [8,[39][40][41][42][43][44][45][46].
At present an approach that starting from the colored glasma field evolves into the locally thermalized quarkgluon plasma is missing. However we notice that the effect of the glasma fields resembles that of an anomalous Fokker-Planck diffusion without a sizeable drag that degrades the charm momentum, as can be seen in Fig.3 of Ref. [19]. Because of this, we can simulate the dynamics of HQs in the pre-thermal equilibrium phase of heavy ion collisions by using the Fokker-Planck equation. To this end we modified the initial Langevin dynamics to mimic the glasma effect on charm quark transport by gauging a momentum dependent diffusion coefficient to generate the same R AA (p T ) obtained within the Wong equations dynamics. The results are shown in Fig. 1 by dashed lines for the two cases considered: a glasma dynamics up to τ glasma = 0.3 fm/c (green dashed line) and τ glasma = 0.5 fm/c (red dashed line) and can be compared to the corresponding results from the Wong equations shown by solid lines. The results with the Langevin simulation are obtained with Γ = 0 and a diffusion coefficient parametrized as D(p) = D 0 (1 + cp) 2 , where c affects only R AA at large p T regulating its increase by increasing the parameter c. To reproduce the R AA (p T ) due to the effect of the glasma, solid lines in Fig.1, we find that D 0 = 3.6 GeV 2 /fm, and c = 0.02 GeV −1 for the case τ glasma = 0.3 fm/c and a slightly increased c = 0.025 GeV −1 to have a nearly perfect matching with the results with Wong's equation also for τ glasma = 0.5 fm/c .
Given that we can gauge the Fokker-Planck equation to mimic the initial stage dynamics up to τ Glasma by mean of D(p) as determined above, we continue the charm dynamical evolution at τ > τ Glasma accordin to the standard charm evolution given by the Fokker-Planck evolution with a drag Γ and a diffusion D related by the fluctuation-dissipation theorem. We note that the diffusion coefficient D(p) in the pre-thermal equilibrium glasma stage is larger than the one in thermally equilibrated QGP stage, which means that there exists a transition between these two stages. As there is a lack of model that describes this smooth transition, we thus simply switch the diffusion coefficient to the one in QGP phase at τ > τ Glasma . The background medium is given by the expanding bulk according to viscous hydrodynamics by a relativistic transport Boltzmann solved at fixed shear viscosity to entropy density η/s [47][48][49]. More specifically, we determine the heavy flavor in medium scattering employing the drag Γ and the diffusion D derived from a quasi-particle model (QPM) [50][51][52].
The QPM approach accounts for the non-perturbative dynamics by T −dependent quasi-particle masses, with m 2 q = 1/3g 2 (T )T 2 and m 2 g = 3/4g 2 (T )T 2 , plus a T −dependent background field known as a bag constant. with g(T ) tuned to fit the thermodynamics of the lattice QCD [10,53].
This approach has been shown to lead to a good description of the experimental data for the R AA (p T ), both at RHIC and LHC employing an enhancement factor K ∼ 2 of the drag and diffusion coefficient. This and similar approaches have allowed first estimates of the spacediffusion coefficient, as discussed in [1,13,54]. However all these approaches do not consider an evolving Glasma stage and have an initialization time τ 0 ≃ 0.3 − 0.6 fm/c.
We choose three cases in order to see the effect of the glasma field: a) free streaming up to τ = 0.3 fm/c followed by a Langevin dynamics in a hydro-like bulk for τ > 0.3 fm/c; b) only diffusion mimicking the glasma evolution up to τ glasma = 0.3 fm/c followed by standard Langevin dynamics with drag and diffusion in a hydrolike but at τ > 0.3 fm/c; c) drag and diffusion from QPM in a Langevin dynamics starting at τ = 0.1 fm/c. This third case is included to see what happens if one simulate the initial stage just starting with an initial time τ 0 = 0.1 fm/c letting charm quarks to undergo scattering in the QGP medium. In all these case we have tuned by a K factor the Γ and D in the QGP phase to reproduce a very similar final R AA (p T ) that reasonably describes the experimental data, as can be seen from Fig. 3 (right  panel).
The aim of this Letter is to point out that heavy-quark interaction with glasma fields induces a dynamics that is opposite to the standard heavy-quark in medium scattering with the bulk medium. We can see in Fig. 2 the transverse momentum dependence of R AA of charm quarks at time τ = 0.3 (left panel) and 1 fm/c (right panel). In the left panel of Fig. 2, we can see that at time 0.3 fm/c, the glasma leads to the enhancement of R AA at high p T (red solid line), while charm propagation in a hydro bulk produces a significant suppression at p T > 2 GeV shown by the blue solid line. The right panel of Fig. 2 shows that drag and diffusion leads to the suppression of R AA of high p T in an efficient way, which can be seen by the significant dropping of the differences compared to the left panel between the three cases already at τ ≃ 1 fm/c. From this we can already understand that to consider the evolution under the gluonic fields at early time is quite different from moving the initial time evolution of the QGP phase down to τ 0 = 0.1 fm/c. After the evolution of charm quarks for the three different cases mentioned in the above section, we hadronize charm quarks to D mesons by a standard hadronization by fragmentation as done in [13,39].
In the left panel of Fig. 3, we show R AA for the three cases, and compare them to the experimental results for the same collision system and centrality [55,56]. It is seen that we can nicely reproduce almost the same R AA in all three cases, but the associated v 2 is quite different (right panel). It should be noted by magenta dotted line that, if the system evolves after the glasma phase at τ = 0.3 fm/c with K factor same as the case without the glasma phase, v 2 does not change so much while R AA will increase. In Fig. 3, the dashed green line shows v 2 of standard evolution in a hydrodynamic bulk starting at τ 0 = 0.3 fm/c, which determines values of v 2 below the experimental results. However, with the inclusion of the Glasma field before the formation of QGP, v 2 is significantly larger, about a 15 − 20% at p T larger than 2 GeV/c and is in the lower limit of experimental results, which is shown by the solid red line of right panel of Fig. 3. If we start the hydrodynamic evolution earlier at τ 0 = 0.1 fm/c, which is shown by the solid blue line, v 2 decreases further by about a 20% at p T > 2 GeV/c. This shows that the initial glasma dynamics cannot be simply simulated by decreasing the intial time τ 0 as done for example in [57]. In fact this last case would lead to estimate about a 30% smaller ellitptic flow of the D mesons which also means a v 2 (p T ) quite smaller with respect to the experimental data.
The generation of a larger v 2 , when a glasma phase is taken into account, is related with the initial enhancement of R AA (p T ). In fact as discussed in Ref. [12] R AA and v 2 are correlated in such a way that a lower nuclear modification factor is associated to a larger v 2 . However the size of this anti-correlation depends on the time evolution of R AA . The initial enhancement due to the diffusion of charm quarks in the evolving Glasma implies that during the QGP phase more interaction is needed in order to get a R AA (p T ) that agrees with the experimental data. This means that charm quarks interact more with the medium when the bulk has a larger ellipticity, so that the bulk itself can transfer it to the charm quark more efficiently. As the lifetime of the glasma field is not yet clearly defined, we explore also the impact of its uncertainty varying it in the range τ glasma = 0.3 − 0.5 fm/c that is the usual time interval within which a hydrodynam-ical expansion of the bulk QGP matter is considered. In Fig.4 (right panel) the band indicates the uncertainty in v 2 (p T ) associated to the uncertainty in the glasma lifetime that appears comparable with the current error bars of the data with the upper bound associated to the largest lifetime considered. It appears a quite good agreement with both the R AA and v 2 of the experimental data. We also emphasize that usually to have a reasonable agreement with both the R AA and v 2 in the phenomenological model one prefers to chooce an initial time τ 0 for the evolution in the range 0.3 − 0.6fm/c with a preference on larger times because otherwise the associated v 2 is smaller, in fact this is shown by the blue solid line in Fig. 3 (left panel). Indeed assuming that nothing relevant happens in the HQ dynamical evolution, is not in general really justified considered also the short HQ formation time τ f orm ≤ O(10 −1 ) fm/c. The mechanism of the initial R AA enhancement allows to not discard the initial phase while acquiring an even better description of the experimental observables.

III. CONCLUSIONS AND DISCUSSIONS
In this study we have pointed out the effect of the initial gluon field on the evolution of charm quarks in high energy Pb+Pb collisions at 5.02 TeV. We have found that the gluon field can push the charm quarks from low to high p T , and this effect can be seen as an initial diffusion of charm quarks in momentum space without a degrading drag force. The effect discussed here is qualitatively different w.r.t. the standard dynamics of HQ in the hot QCD medium, because it implies a rise and fall of the R AA (p T ) during its time evolution. Using the Fokker-Planck equation extended to include both the effect of the glasma and QGP, we have found that the correlation between the R AA and v 2 of D mesons is modified leading to a final large v 2 at the same R AA when the evolution in the glasma is considered. This points towards a better description of both observables and we have discussed that this can be understood as an effect due to the delay in the formation of R AA . We have also found that a longer lifetime of the glasma can lead to a larger effect on the enhancement of v 2 , but already for τ glasma ∼ 0.3 fm/c the effect can be quite large, about a 20%. We remark that the longitudinal expansion of the colliding systems in the glasma phase is not included in this study. This may affect quantitatively the effect on R AA and v 2 , and should be considered in future studies. However, the qualitative result presented here will not be affected by the presence of the longitudinal expansion, as this will have the effect to dilute the energy density but it will not affect the physical mechanisms described. We note that there is a discontinuity when we switch the diffusion coefficient from glasma phase to QGP phase. A complete modeling should be developed to describe the evolution of the glasma phase in a QGP one through a mix stage. This is certainly a challenge for future developments. We did not include the cold nuclear matter effect which affect R AA especially at p T < 2 GeV, even if it would enhance the nuclear modification factor at intermediate p T further contributing to the main effect discussed here. We have not included for simplicity the back-reaction on the HQs in the glasma phase. This has to be certainly studied in future developments. However, considering that the only energy scale is g 2 µ which is in the range 3-5 GeV, we may expect that the effect of the gluon radiation will be not substantial at least up to this p T scale.
In this Letter we have presented the first results on the two main observables in HIC, R AA and v 2 . However a thorough understanding of the initial stage dynamics is a timely fundamental task and should affect at least two other observables like the triggered D −D angular correlation and the splitting in the directed flow v 1 that recently has been shown to probe the initial strong electromagnetic field that reaches its maximum value in the same time range of the glasma stage [20]. In our opinion, studies in this direction have to be pursued because they will also allow a unified description between the dynamics in pA and in AA collisions. Furthermore such a direction joins the current lively activity in the study of the initial stage of ultra-relativistic collisions to the HQ physics.