Fluctuations of topological charge and chiral density in the early stage of high energy nuclear collisions

We study systematically the topological charge density and the chiral density correlations in the early stage of high energy nuclear collisions: the intial condition is given by the McLerran-Venugopalan model and the evolution of the gluon fields is studied via the Classical Yang-Mills equations up to proper time $\tau\approx 1$ fm/c for an $SU(2)$ evolving Glasma. Topological charge is related to the gauge invariant $\bm E \cdot \bm B$ where $\bm E$ and $\bm B$ denote the color-electric and color-magnetic fields, while the chiral density is produced via the chiral anomaly of Quantum Chromodynamics. We study how the correlation lengths are related to the collision energy, and how the correlated domains grow up with proper time in the transverse plane for a boost invariant longitudinal expansion. We estimate the correlation lengths of both quantities, that after a short transient results of the order of the typical energy scale of the model, namely the inverse of the saturation scale. We estimate the proper time for the formation of a steady state in which the production of the chiral density in the transverse plane per unit rapidity slows down, as well as the amount of chiral density that would be present at the switch time between the Classical Yang-Mills evolution and the relativistic transport or hydro for the quark-gluon plasma phase.


I. INTRODUCTION
The study of the initial condition in high energy collisions is a difficult but interesting problem related to the physics of relativistic heavy ion collisions (RHICs), as well as to that of high energy proton-proton (pp) and proton-nucleus (pA) collisions. The dynamics of the central rapidity region is determined by the small Bjorken x gluons before the collision where saturation takes place [1][2][3][4]. In the saturation region the gluon occupation number is large enough that classical effective field theories (EFTs) can be used [5][6][7][8][9]: the Lorentz-contracted colliding nuclei are idealized to fly along the light cone, with the large-x partons behaving as static sources of the small-x modes that make the color-glass condensate (CGC) fields inside the two nuclei, see Refs. [2][3][4][10][11][12][13] for reviews. Right after the collision, color sources form on the two collision sheets as a result of the interaction of the two CGC colliding sheets, in such a way longitudinal color electric, E, and color magnetic, B, fields are formed: this particular field configuration is named the Glasma [14] and it serves as the initial condition for the evolution of the classical gluon field after the collision that can be studied by the Classical Yang-Mills (CYM) equations [15][16][17][18][19][20][21] up to the formation of the quark-gluon plasma.
At the initial time E · B = 0 in Glasma, that results in a nonvanishing topological charge density, ρ T : the initial condition then consists of approximately N = πR 2 A Q 2 s uncorrelated domains of topological charge density [22], where R A is the nucleus radius and Q s the saturation scale, in which the amplitude of the fluctuations of the topological charge density is given by √ N ≈ R A Q s . It has been discussed that the high energy collisions might be described by the decay of instanton like configurations, the sphalerons [23][24][25]. The decay of these longitudinal electric and magnetic fields is in fact the decay of topological charges. Since ρ T is charge conjugation and parity (CP) odd, it might be the source of large CP violating fluctuations in heavy ion collisions, see Refs. [26][27][28][29], and [30][31][32][33][34][35] for recent reviews.
A nonzero ρ T naturally induces the production of chirality imbalance, N 5 , by virtue of the chiral anomaly of Quantum Chromodynamics (QCD). In this article, we focus on the topological charge density ρ T and the associated chiral charge density N 5 produced in the evolving Glasma. Differently from [15,36] where the evolution has been followed up to late times, we will focus on the very early stage, that is the proper time range in which the description based on CYM has some phenomenological interest for high energy nuclear collisions. In particular, we are interested to follow the evolution of the correlators of ρ T and of the chiral density in the early stage, computing the size of the correlation domains as a function of the center of mass energy of the collision; we also compute the expectation value of the chiral density in the steady state, namely in the range or proper time in which the fields are diluted enough that further chiral density is not produced, and study how this depends on the collision energy. The latter estimate can have some phenomenological interest as it turns out that the proper time needed to form the steady state is in the same ballpark of the thermalization time, namely the time at which the CYM evolution is switched to that of the quark-gluon plasma via relativistic hydro or transport [19,20,[37][38][39][40][41]: the amount of chiral density that we compute will therefore be present at this stage and it is likely to persist for the full evolution of the system up to hadronization time since the relaxation time for N 5 in the hot quark-gluon plasma due to the helicity flipping processes turns out to be larger than the typical lifetime of the fireball produced in heavy ion collisions [42]. This study is a natural continuation of previous works on the correlators in the evolving Glasma [43][44][45], and paves the way to future studies on the anomalous particle production by the chiral magnetic effect (CME) [46][47][48][49][50][51].
The article is organized as follows. In Sec.II, we briefly review the McLerran-Venugopalan (MV) model and the classical Yang-Mills (CYM) equation. In Sec.III, we introduce the topological charge density and chiral charge density via Adler-Bell-Jackiw anomaly equation, then we define the associated correlators. In Sec.IV, we show our results for topological charge density and chiral charge density correlations at different energy scale. Then, we estimate the size of the correlated domains for both quantities. At last, we analyze their time dependence as well as energy scale dependence. In Sec.V, we draw our conclusion and make some outlooks.

II. GLASMA AND CLASSICAL YANG-MILLS EQUATION
In this section, we briefly review the McLerran-Venugopalan (MV) model, by which we describe the initial condition of the classical gluon field produced after the collision [2][3][4]6], which is then evolved by the classical Yang-Mills (CYM) equations . We remark that in our notation the gauge fields have been rescaled by the QCD coupling A µ → A µ /g therefore g does not appear explicitly in the equations.

A. Glasma
In the MV model, the color charge densities ρ a act as static sources of the transverse CGC fields in two colliding nuclei: they are assumed to be random variables that, for each of the two colliding nuclei, are normally distributed with zero average and variance specified by the equation (1) where a and b denote the adjoint color indices, x ⊥ and y ⊥ denote transverse plane coordinates and η 1 , η 2 the spacetime rapidities. In this work, we limit ourselves to a SU(2) glasma, therefore, a, b = 1, 2, 3. In the MV model, g 2 µ is the only energy scale which is related to the saturation momentum Q s : rather, we simply refer to the estimate of [52] namely that Q s /g 2 µ ≈ 0.57. Due to the δ−functions in Eq. (1) the color charge densities are uncorrelated in transverse plane as well as in rapidity. To specify the initial condition, namely the Glasma fields, it is convenient to work in Bjorken coordinates (τ, η) in the radial gauge, where In order to compute Glasma fields we firstly solve the Poisson equations, namely where A and B denote the two colliding nuclei. The solutions of these equations are while the Wilson line is defined as U (x ⊥ ) ≡ Pexp dz µ α µ (z(x ⊥ )) , with P is the path order operator and z(x ⊥ ) is a trajectory. In terms of these fields, the glasma gauge potential at τ = 0 + can be written as [53,54]: Solving the Yang-Mills equation near the light cone, one finds that the transverse color electric and color magnetic fields vanish as τ → 0, but the longitudinal electric and magnetic fields are non-vanishing [55]: In all the discussion above we have neglected the possibility of fluctuations that, among other things, would break the longitudinal boost invariance: while these are relavant for the onset of the hydrodynamical flow [15-17, 21, 56-59], they are not crucial for the production of chiral density and for the evolution of the topological charge density; in any case, we plan to add these fluctuations in near future works. We also assume that g 2 µ is the same in the two colliding nuclei and it has no dependence on the transverse plane coordinates: we will remove this assumption in the future, to mimic the energy density profile that would be produced in realistic collisions [19,20].

B. The classical Yang-Mills equation
After preparing the initial condition of CYM equations, we solve these numerically. In the gauge A τ = 0, the Lagrangian density reads and the canonical momenta are given by As a consequence, the Hamiltonian density is The CYM equations in Bjorken coordinates are where ij F ij as the i and η components of the color magnetic field. In the next section, these equations will be solved in a 4 fm × 4 fm box with periodic boundary conditions on the transverse plane. The lattice spacing that we use is a = 0.04 fm. We solve the CYM equations via the fourth order Runge-Kutta method as in [21,36]: results with the given approach are in good agreement with the Yang-Mills solver based on the gauge links [16,17], and we leave the more rigorous implementation of the numerical problem based on gauge links to a future work.

III. THE TOPOLOGICAL CHARGE AND THE CHIRAL DENSITIES
As mentioned above, immediately after the collision nonzero longitudinal components of the color electric, E = E a T a , and color magnetic, B = B a T a fields are formed (here and in the following we use boldface to denote vectors in color space): because of the Adler-Bell-Jackiw anomaly equation [60,61], these lead to the nonconservation of the U (1) axial symmetry of QCD, that in Bjorken coordinates is written as where, in the same coordinate system, In writing Eq. (18) we have used the boost-invariance assumption as well as that the divergence of the axial current in the transverse plane vanishes, the latter requirement implying that there is no net flux of the axial current across the transverse plane. The gauge invariant quantity on the right hand side of Eq. (18) is called the topological charge density since its integral over the full volume and on the time history of the system gives the change in the Chern-Simons number of the gluon configuration, Therefore, here we call the invariant on the right hand side of Eq. (18) as the topological charge density, where we have made explicit the dependence of this quantity on τ and x ⊥ that comes from that of E and B. It is possible to use j τ 5 to define a chiral charge per unit volum, namely where N 5 corresponds to the chiral charge and τ dηd 2 x ⊥ is the 3-dimensional volume of a cell with extension d 2 x ⊥ in the transverse plane and τ dη in the longitudinal direction; from this we can define the chiral charge distribution in the transverse plane per unit rapidity as In this study, we solve the CYM equations to get ρ T by means of Eq. (21) and we use it in the right hand side of Eq. (18) that can be formulated in terms of n 5 by writing j τ 5 = n 5 /τ then solving for n 5 , namely Here we study how topological charge density and chiral density correlation domains form with time, and we estimate their size on the transverse plane. According to Eq. (23) quarks should form as soon as the chiral anomaly starts to play its game: in principle, these quarks bring additional sources in the Yang-Mills equations (16) and (17), but we neglect these because we limit ourselves to study the early stages after the collision, and in this stage the number of quarks per unit volume is small leaving the dynamics dominated by the gluons.
Next we turn to the definition of the gauge invariant correlators that we are interested to; in this study, we consider only equal time correlators, leaving the different times ones to a future study. In continuum limit, the correlator of ρ T we analyze is defined as where we have made explicit the dependence on the transverse plane coordinate, x ⊥ and denotes the ensemble average. On the lattice, this average is done by running a finite amount of events, N events , then summing over all these and dividing by N events , where each event is initialized with a different random condition. In addition to this average, we introduce an average over the volume of the box in order to improve the statistics. Therefore, we compute the average of an observable O as where x j denotes the jth lattice point and N is the number of lattice points involved in the summation: concretely, for a given  Similarly to Eq. (25) we define the correlator of the chiral density, From this correlator we can also define n 2 5 ≡ n 5 (|x ⊥ |)n 5 (|x ⊥ |) that estimates the amount of fluctuations of chiral density in the transverse plane and that has phenomenological interest like in Refs. [31,33]. The study of the correlators (27) and (25) is useful because these allow to estimate the size of the correlation domains of chiral and topological densities in the early stage of high energy nuclear collisions, as we will discuss in the next section.

IV. RESULTS
In this section we present our results: we firstly discuss the correlation of the topological charge, then we turn to that of the chiral density. Finally, we estimate the size of the correlated domains for both quantities. In particular, we analyze the time dependence as well as the dependence on the energy scale g 2 µ of the correlators.

A. Correlation of the topological charge density
In Fig. 1, we plot the correlator of the topological charge density, Eq. (25), versus the transverse plane coordinate at different proper times; the correlators have been normalized to their values at x ⊥ = 0 (an overall constant in each correlator is not relevant at all if we focus on the correlation length), and have been computed for g 2 µ = 1 GeV. The calculations have been performed on a lattice with transverse size 4 fm × 4 fm with a lattice spacing a = 0.04 fm giving g 2 µa = 0.2. The results show that the correlation of ρ T decays quickly in transverse plane: in fact, for all the cases shown in the figure already for x = 0.16 fm, corresponding do a dimensionless distance of g 2 µx = 0.8, correlation reduces to approximately the 40% of the initial value. Although the qualitative behavior of the correlators is the same at each of the finite time we analyze, we notice some fluctuation of the shape as time evolution goes on: compare for example the cases at g 2 µτ = 2, 3 and 4 in Fig. 1: these fluctuations are due to the continuous exchange of energy between the longitudinal and the transverse fields that happens in the evolution, but they lead to almost no quantitative effect on the correlation length, see below. We also notice that some modest anticorrelation shows up in time: an anticorrelation was also obtained for the color electric and color magnetic fields [21,44] and in this case shows the tendency of the topological charge density to flip its sign on length scales much larger than the correlation domains.
In Fig. 2, we plot the correlator of the topological charge density versus the transverse plane coordinate for three values of g 2 µ: this piece of information never appeared in the literature before, but it is important because changing g 2 µ amounts to change the energy of the collision and a larger g 2 µ corresponds to a larger collision energy. In this paper, we have worked up to g 2 µ = 3.4 GeV that corresponds roughly to the estimate for this quantity for a Au-Au collision at the maximum RHIC energy. We find that increasing the g 2 µ results in a broader correlator when this is studied as a function of the dimensionless length g 2 µ|x ⊥ |. In the lower panel of Fig.  2 we plot the same quantity versus the physical length: increasing the g 2 µ results in a quicker decay of the correlator, which suggests that the correlation length of the topological charge density becomes smaller.
The results in Fig. 2 allow to compute a correlation length, λ T , of ρ T in the transverse plane: we define this length by the value of |x ⊥ | such that the correlator decays to 1/e of its initial value. Since the correlators depend on time, this definition in principle depends on time as well: therefore, we perform an average over the full time history of the system in addition to the ensemble average. Concretely, what we do is that firstly we compute the correlators at any time by taking the ensemble average of the proper quantity; then, for each time we find the value of |x ⊥ | such that C(x ⊥ )/C(0) = 1/e, and we average over these values, keeping also the maximum and the minimum of the range to define the uncertainty. We show the result of this definition in Fig. 3 in which we plot the correlation length versus g 2 µ. In the upper panel we plot the dimensionless correlation length, while in the lower panel we show the correlation length in units of fm. In the figures, the blue dots correspond to the average value while the error bar denotes the maximum and the minimum value of λ T achieved during time evolution. We find that on average g 2 µλ T lies within about 0.5 − 2.5 in the range of g 2 µ analyzed, which is consistent with the expectation that the correlation domains of ρ T are of the same transverse size of the correlation domains of color electric and color magnetic fields, see for example [21]. Changing to physical units, we find that λ T lies in the range (0.1 fm, 0.2 fm), namely the domains remain microscopic with respect to the transverse size of nucleons. Altogether, these results show that increasing the energy of the collision, the number of correlation domains of ρ T in the transverse plane increases as ≈ S/λ 2 T ≈ S(g 2 µ) 2 where S denotes the transverse area of the nucleon, so the density of these domains in the transverse plane behaves as ≈ (g 2 µ) 2 : increasing the energy results in more and finer domains of topological charge density in the transverse plane. Physical correlation length of ρT versus g 2 µ. The parameters are the same as in Fig. 1.

B. Production of chiral density
In Fig. 4 we plot the averaged n 2 5 , shown in a logscale, versus proper time for three values of g 2 µ. In the upper panel we show the time measured in units of g 2 µ while in the lower panel we plot the quantity versus time measured in fm/c. We notice that in all the cases considered here, for g 2 µτ ≈ 1.0 the bulk of chiral density is already formed, since a steady state is reached within this time range: this fixed point is due to the dilution of the fields with the expansion that lowers the value of the topological charge density and stops the production of n 5 via the chiral anomaly. The actual value of the average n 2 5 in the steady state depends on g 2 µ, which is obvious because increasing the latter results in a higher average of the gluon fields in the initial condition.
In both panels of Fig. 4 we have highlighted by filled triangles the values of the chiral charge density at g 2 µτ = 1 that corresponds roughly to the proper time range necessary to reach the steady state; we denote this time by τ s . In calculations based on CYM the g 2 µ is a param- eter: it is interesting to give numerical estimates of the chiral charge density in the steady state and of its dependence on g 2 µ. From Fig. 4 we calculate the squared chiral charge per unit volume and unit rapidity at τ = τ s , where the physical volume per unit rapidity is V = A T × τ s with A T denoting the transverse area and τ s = 1/g 2 µ, and denotes ensemble and transverse plane average [we remind that our definition of n 5 corresponds to chiral charge per unit transverse area and unit rapidity, see Eqs. (23) and (24)]. Our findings in the steady state are N 2 5 /V 2 ≈ 1.79 × 10 −1 fm −3 at g 2 µ = 1.0 GeV, N 2 5 /V 2 ≈ 8.85 × 10 −1 fm −3 at g 2 µ = 1.7 GeV, N 2 5 /V 2 ≈ 6.48 fm −3 at g 2 µ = 3.4 GeV.

C. Correlation of the chiral density
In the upper panel of Fig. 5 we plot the correlator of the chiral density (normalized to x ⊥ = 0 in each case) versus g 2 µ|x ⊥ |, for several values of the proper time, for g 2 µ = 1 GeV. Increasing time results in a progressive broadening of the correlator. In the lower panel of Fig. 5 we plot D(x ⊥ )/D(0) at fixed time g 2 µτ = 1 for three different values of g 2 µ. In agreement with the results for ρ T presented in the previous subsection, increasing g 2 µ affects mildly the correlator of n 5 .
From the results in Fig. 5 we can estimate the correlation length, λ 5 , of n 5 by looking for the value of |x ⊥ | such that the correlator decays to 1/e of its value at x ⊥ = 0; we follow the same procedure depicted in the previous subsection for the definition and computation of λ T . We summarize the results of this estimate in Fig. 6 where we plot λ 5 versus g 2 µ, both as a dimensionless quantity (upper panel) and in physical units (lower panel). The correlation length sits in the range (0.5, 2.5) in units of g 2 µ, in agreement with the results we have found for the topological charge: the correlations of ρ T are transmitted to those of n 5 almost unaffected.

V. CONCLUSION AND OUTLOOK
We have studied the correlations of the topological charge density, ρ T , carried by the strong gluon fields in the early stage of high energy nuclear collisions; besides, we have analyzed the production and the correlations of the chiral density per unit rapidity in the transverse plane, n 5 , produced by the chiral anomaly of QCD. In fact in the early stage of the collisions, assuming the McLerran-Venugopalan initialization and the Glasma picture, ρ T ∝ E ·B = 0 and dn 5 /dτ = τ ρ T , where E and B denote the color electric and color magnetic fields and τ is the proper time. We have described the evolution of the gluon fields by means of the Classical Yang-Mills (CYM) equations, that we have solved numerically for the case of SU (2): we have followed the evolution up to τ = 1 fm/c, in agreement with the switch time from the CYM evolution to relativistic hydro or transport in nucleus-nucleus collisions at the RHIC energy [19,20,[37][38][39][40]. We have analyzed the case of a longitudinal boost invariant expansion that is more relevant for the modeling of the early stage of the collisions. Despite some amount of information available in the literature on the formation of topological charge domains in the early stage of high energy nuclear collisions, see for example [14,30,32], a concrete calculation that takes into account the full evolution of the gluon field with the initial condition given by the MV model is still missing: we aim to cover this lack of information here. Besides its theoretical own interest, this study is relevant for the phenomenology of high energy nuclear collisions: for example, just recently it has been pointed out that the production of chiral density in the early stage can affect the polarization of the Λ particles [31,33] and puts constraints on the strength of the magnetic field produced in the collisions: it is therefore important to have both a qualitative and a quantitative understanding of the n 5 that is produced in the early stage where classical gluon fields dominate the dynamics.
We have computed the correlator of ρ T in the transverse plane, studying how the correlation develop and how the correlation length depends on g 2 µ, namely on the density of color charges in the transverse plane: this number is the only energy scale in the initial condition, and increasing it amounts to model a collision with a higher center of mass energy. Examining the range (1, 3.4) GeV for g 2 µ, we have found, for the correlation length λ T , that λ T ≈ O(1) × 1/g 2 µ: increasing the collision energy amounts to fit more correlated domains within the transverse plane, the density of these being given approximately by O(1) × (g 2 µ) 2 . These domains remain microscopic with respect to the nucleon, in the sense that in physical units the correlation length sits in between the 10% and the 20% of the proton radius. The correlators and correlation lengths that we have found agree with previous studies [21,30,32,43,44].
We have also studied the production of the chiral density in the transverse plane per unit of rapidity, n 5 ≡ dN 5 /dηd 2 x ⊥ , that is formed by the chiral anomaly of QCD. For this quantity, we have analyzed the production rate and how it depends on g 2 µ, giving some estimate of this value when the system reaches a steady state. The n 2 5 depends on g 2 µ, which is obvious since increasing the latter amounts to increase the magnitude of the gluon fields in the initial condition. We have estimated the proper time to reach a steady state, τ s , in the production of n 5 to be τ s ≈ 1/g 2 µ. To give concrete numbers, we have found N 2 5 /V 2 ≈ 1.79 × 10 −1 fm −3 at g 2 µ = 1.0 GeV, N 2 5 /V 2 ≈ 8.85 × 10 −1 fm −3 at g 2 µ = 1.7 GeV, N 2 5 /V 2 ≈ 6.48 fm −3 at g 2 µ = 3.4 GeV. The average chiral density per unit rapidity can be quite large for the largest value of g 2 µ that we have used, at least until the steady state is formed: the longitudinal ex-pansion will dilute N 2 5 /V 2 as 1/τ 2 in the steady state, when the topological charge density is low enough that no substantial new n 5 is produced.
We have analyzed the correlation domains in the transverse plane and how their size is affected by g 2 µ: as for ρ T we have found that the correlation length of n 5 , namely λ 5 ≈ O(1) × 1/g 2 µ, and within numerical uncertainty λ 5 ≈ λ T .
The qualitative picture that arises from our study is that increasing the collision energy, the transverse plane gets populated by a larger amount of correlated domains of ρ T and n 5 , the size of these becoming smaller as 1/g 2 µ and their density increasing as (g 2 µ) 2 ; chiral density per unit rapidity forms immediately after the collision, and a steady state is achieved after a short proper time range τ s ≈ (0.1 fm/c, 0.2 fm/c), the lowest value corresponding to g 2 µ = 3.4 GeV and the highest to g 2 µ = 1 GeV. We mention that we have also also found some anticorrelation at the initial, time, in agreement with previous studies of the gauge invariant correlators of the gluon fields [21,43], but the evolution cancels this anticorrelation in a time range g 2 µτ = O(1).
A natural extension of the work reported here is the study of the correlations in rapidity, introducing fluctu-ations along the longitudinal direction in the initial condition; besides, the systematic study we have presented here can be enriched by modeling the realistic transverse plane geometries of nucleus-nucleus and proton-nucleus collisions. In addition to this, the diffusion of the chiral density in the transverse plane in the early stage might be an interesting topic to investigate. Even more, the production of photons in the early stage due to coexistence of n 5 and an electromagnetic field in the early stage is worth of being studied. We plan to address these topics in the near future.