Linearly polarized gluons and axial charge fluctuations in the Glasma

We calculate of the one- and two-point correlation functions of the energy density and the divergence of the Chern-Simons current in the nonequilibrium Glasma state formed in a high-energy nuclear collision. We show that the latter depends on the difference of the total and linearly polarized gluon transverse momentum distributions. Since the divergence of the Chern-Simons current provides the source of axial charge, we infer information about the statistical properties of axial charge production at early times. We further develop a simple phenomenological model to characterize axial charge distributions in terms of distributions of the energy density.


Introduction
Novel transport phenomena associated with the Chiral magnetic [1][2][3] and related effects have recently caused an excitement across different fields of physics. In the high-energy QCD context, experimental measurements at RHIC and LHC have provided intriguing hints at possible signatures of such anomalous transport phenomena [4][5][6][7]. However, the interpretation of these experimental results remains inconclusive [8,9] due to the presence of large background effects [10][11][12][13][14]. Despite significant progress on the theory side in developing different microscopic [15][16][17][18][19][20] and macroscopic [21,22] descriptions of the coupled dynamics of vector and axial charges, a first principles description of the effects in highenergy heavy-ion collisions remains an outstanding challenge. Present phenomenological predictions [23][24][25][26] have to rely to a varying extent on modeling assumptions. Most importantly, all phenomenological descriptions based e.g. on anomalous hydrodynamics [21,22] require information about the early time dynamics as an initial condition for the subsequent space-time evolution. Even though significant progress has been achieved in understanding the early time dynamics of the conserved energy-momentum tensor, both from a theoretical perspective [27,28] as well as through sophisticated model/data comparisons [29,30], achieving a similar level of understanding of the space-time dynamics of axial charge production and anomalous transport processes during the very early pre-equilibrium stages ( 1fm/c) remains a key challenge.
One important difference between the dynamics of vector and axial charges, is the fact that the density of axial charge is not conserved. This is due to the axial anomaly, which for N f flavors of (approximately) massless fermions takes the form where F µν denotes the field strength andF µν = 1 2 ε µνρσ F ρσ its dual. Hence understanding the dynamics of axial charges and currents in a QCD plasma inevitably requires some knowledge about the structure of non-abelian gauge fields entering on the right-hand side of Eq. (1.1). Even though it is understood that in the long time and long wave-length limit, topological (sphaleron) transitions dominate the production/dissociation of axial charge (see e.g. [31] and references therein), it is not clear to what extent these considerations apply to the typical time and length scales relevant during the early stages of high-energy heavy-ion collisions [32]. Despite the fact that the rate of topological transitions can be significantly enhanced during the early time pre-equilibrium stage [32], various kinds of short distance field strength fluctuations can also contribute significantly to axial charge production at early times. Consequently, it is of crucial importance to understand different mechanisms of axial charge production in order to estimate the magnitude and features and isolate the most relevant effects.
One more direct way to study the strong gauge fields that dominate the initial stages of heavy ion collisions is to probe them with a dilute probe, such as in high energy deep inelastic scattering. The experimental program at a future Electron-Ion Collider [33] will be able to characterize the spacetime structure of partons inside nucleons and nuclei in a variety of ways. Out of these the linearly polarized gluon transverse momentum distribution [34,35] has recently been of particular interest to the small-x community. Based on the Color-Glass-Condensate (CGC) picture it has been shown [36][37][38] that the linearly polarized gluon distribution can be related to, and ultimately calculated from, the same Wilson line correlators that characterize unpolarized gluon distributions. We will show in this paper, that the correlation structure of the gauge fields at the earliest times after a heavy ion collision is sensitive to both the linearly polarized and unpolarized gluon distributions. It turns out that the correlations and fluctuations of axial charge are particularly sensitive to the polarized distributions. Whereas for energy density fluctuations the polarized and unpolarized contributions add up, for the axial charge they appear with a different sign. This observation opens up a fascinating new connection between correlation studies in deep inelastic scattering and local CP-violating fluctuations in hadronic collisions.
The aim of this paper is to calculate the statistical properties of axial charge production at the earliest stages of a high energy heavy ion collision. The calculation is based on the description of the early time dynamics in the Color Glass Condensate framework [39][40][41][42], which leads to the presence of longitudinal chromo-electric and chromo-magnetic fields at very early times after the collision. We start with brief discussion of the space-time struc-ture of chromo-electric and chromo-magnetic fields at very early times in Sec. 2. We then review in Sec. 3 the structure of the linearly polarized and unpolarized Weiszäcker-Williams (WW) gluon distributions in the CGC framework. With a Gaussian approximation for the field correlators (the "glasma graph" approximation) we then perform an analytic calculation of the one and two-point correlation functions of the energy density ε(x) and the divergence of the Chern-Simons currentν(x) in terms of the WW correlators. Here our calculation generalizes the closely related earlier work of [43]. We then in Sec. 4 relate our calculation to works studying two-gluon correlations using the glasma graph approximation. We finally discuss the implications of our results for the basic phenomenological properties of axial charge production in the Glasma in Sec. 5, developing a simple algorithm for using our results in anomalous hydrodynamical calculations. We conclude in Sec. 6 with a summary of our results and perspectives for future studies.

Glasma flux tubes and axial charge production
The CGC effective theory description of a high energy nucleon or nucleus is based on a a separation of scales. Degrees of freedom carrying a large fraction of the energy of the projectile/target are described as a color charge, which acts as a source for the small-x gluons. The color field of a single nucleus can be expressed analytically in terms of the color charges. When transformed to light cone gauge, these fields (which we denote here by α i and β i for the two nuclei) are "transverse pure gauge" fields [44,45] Here U x and V x are light-like Wilson lines, which are scattering amplitudes for the eikonal interaction of a color charge passing through the color field. Based on this picture, a high-energy heavy ion collision is realized when two such systems pass though each other. In this case the color fields of the projectile and target interact with each other, leading to the formation of a non-equilibirum "Glasma" [46] state. By requiring that the fields be continuous over the future light cone one obtains [47,48] the gauge fields immediately after the collision at τ = 0 + as In terms of the field strength tensor these correspond to longitudinal color-electric and color-magnetic fields [27,[46][47][48][49] We use the sign convention D µ = ∂ µ +igA µ for the covariant derivative, and take the electric and magnetic fields in usual Minkowski coordinates to be E i ≡ F 0i and B i = − 1 2 ε ijk F jk . For the components in proper-time rapidity coordinates we define The subsequent dynamics at very early times has been studied in great detail analytically e.g. based on small proper time expansions [50,51], as well as numerically through real-time lattice simulations [52][53][54][55][56][57][58]. On a time scale τ ∼ 1/Q s the classical Yang-Mills dynamics leads to the decoherence of the longitudinal fields building up transverse field strengths E i and B i . Eventually the longitudinal expansion leads to a significant reduction of the field strength, where the semi-classical description becomes inapplicable [59][60][61] and the system undergoes a kinetic regime before approaching local thermal equilibrium [62,63].
Even though the structure of the boost-invariant fields in Eq. (2.3) is topologically trivial [64], the strong longitudinal chromo electric and chromo electric fields at early times can still contribute significantly to axial charge production. Despite the fact that the axial charge is of course carried by the fermionic degrees of freedom, an estimate of this effect can be immediately deduced from the axial anomaly relation. In this spirit, a first estimate of the fluctuations of the net axial charge density per unit rapidity was provided in Ref. [64] based on explicit numerical simulations of the early time dynamics (see also [23] for a parametric estimate used in phenomenological studies). We will follow a different approach and estimate the fluctuations of the axial charge directly from the analytic expressions for the initial fields in Eq. (2.3), including also the structure of fluctuations of the axial charge density in the transverse plane. Starting from the explicit form of the axial anomaly relation (1.1) in Bjorken coordinates 1 and defining a shorter notationν(x) ≡ trE · B for the divergence of the Chern-Simons current we first note that the term ∂ η j η (5)x vanishes by virtue of the boost invariance assumption. Next we note that -at sufficiently early times -we can neglect the effect of the axial currents ∂ i j i (5) (x), such that the source term on the right-hand sidė leads to local production of axial charge imbalance before axial charge starts to diffuse in the transverse plane. Based on this approximation one can then estimate the local density of axial charge per unit rapidity at each point in the transverse plane according to 1 Note that the transformation to co-moving coordinates can performed by expressing the left-hand side This allows us to compute axial charge production directly from the correlation functions of light-like Wilson lines. As we will discuss shortly the expectation value of the quantitẏ ν(x) is vanishes in accordance with the fact that there is no CP violation in the process. Nevertheless, there can be sizeable fluctuations on an event-by-event basis, which are characterized by the correlation function ν(x)ν(y) at two different points x, y in the transverse plane. Sinceν(x) is a dimension four operator, it is most naturally compared to the energy density ε(x) of the system, and we will also compute the correlation functions of the energy-density ε(x)ε(y) for comparison.

Energy density and Chern-Simons currents in the Glasma
Before we turn to the evaluation of correlation functions of the energy density and the divergence of the Chern-Simons current, we will briefly review the calculation of the corresponding one-point functions. Even though the results are well established in the literature [50,51,65] this exercise is nevertheless useful to illustrate the procedure and fix our notations.

Expectation values of one-point functions
Based on the analytic expressions for the color-electric and color-magnetic fields at τ = 0 + we can immediately compute the expectation value of the local energy density ε(x) and the divergence of the Chern-Simons currentν(x) as Evaluating the color structures by decomposing α, β over the Lie Algebra, noting that tr[t c t c ] = δ cc 2 and separating the averages over projectile and target fields, we obtain Since only color-singlet expectation values are non-vanishing, such that we can evaluate the color structure as f abc f abc = N c (N 2 c − 1). Upon factorization of the averages of the projectile and target fields we obtain where W ik (U/V ) (xy) are the Weizsäcker-Williams gluon distributions of the two nuclei (3.7) Generally speaking, the Weizsäcker-Williams distribution can be further decomposed into various different tensor structures. We start from the usual momentum space decomposition into unpolarized G (1) and linearly polarized h (1) ⊥ gluon distributions in an unpolarized hadron 2 3 The corresponding tensor decomposition in coordinate space takes the form where the coordinate space functions G (1) (x, y) and h (1) ⊥ (x, y) are given by Note that due to the angular structure of the integration, h  The Weizsäcker-Williams distribution W ij is a 2 × 2 matrix with eigenvalues G (1) ± h (1) ⊥ . As expectation values of positive definite operators, both W ij (U ) (x, x) (at the same coordinate x = y) and the impact parameter averaged W ij (k) (for general k) should be positive definite. This leads to positivity constraints [34] in both coordinate and momentum space, which in our notation read G (1) (x, x) ≥ |h   ⊥ (|k|) in our normalization) as is expected at high transverse momentum, this is not true for the coordinate space functions due to the behavior of the Bessel functions near the origin, which will be important for our discussion in the following.
Collecting everything and expressing the result in terms of the G (1) and h (1) we obtain the following expression for the local operator expectation values 2 Note that taking into account impact parameter b dependence the most general decomposition requires additional tensor structures involving b as well as combinations of b and k. However, since we are only interested in the application to the collisions of large nuclei, we will ignore these subtleties and proceed as usual. We refer the interested reader to Ref. [50,51] for a detailed discussion on the implications for energy density correlators in the glasma. 3 Normalization conventions for the Weizsäcker-Williams distributions vary, for example the one introduced in [66] is related to the one here by where the index contractions lead to a vanishing result for the expectation value of the CP odd operatorν. We see that the linearly polarized distribution does not contribute to the expectation value 4 .

Saturation models for Weiszsäcker-Williams distribution
In order to provide explicit results for the one-and two-point correlation functions, we need to specify a model for the Weiszsäcker-Williams gluon distribution. We follow previous works and exploit the fact [67,68] that in Gaussian models the Weiszsäcker-Williams gluon distribution can be related to the Dipole gluon distribution for which a number of phenomenologically useful parametrizations exist. Based on this standard procedure, described for completeness in Appendix A, we obtain is the expectation value of the dipole operator. From this it is relatively easy, assuming that the dipole distribution only depends on the distance r ≡ |x − y|, to extract the individual distributions as We can explicitly evaluate the correlation function in a number of simple models. One finds for instance that in the Golec-Biernat Wusthoff (GBW) model [69] for the dipole amplitude the linearly polarized gluon distribution vanishes identically h ⊥,GBW (x, y) = 0 and the un-polarized gluon distribution is simply given by (3.17) such that in the limit y → x, relevant for the expectation value of the local energy density Conversely, in the (screened) McLerran-Venugopalan (MV) model the dipole amplitude is given by such that the unpolarized and linearly polarized distributions become At fixed coupling g the function G (1) .

(3.23)
We will take the suggestion of some previous works (e.g. [70]) and regulate the logarithmic divergence by introducing a running of the coupling via the following replacement in Eqs. (3.20), (3.21): with the coordinate space running coupling We then absorb the superfluous parameters into physical ones by expressing the correlators in terms of the physical momentum scale of the problem, the saturation scale Q s . This can be done by taking the limit Λ ∼ m µ ∼ Q s , in which the unpolarized distribution becomes We then require that this has the same normalization in term of the saturation scale Q s as in the GBW parametrization, see Eq. (3.18). This can be achieved by setting Expressed in terms of Q s we now have the same short distance behavior as in the GBW model and we will employ this prescription in the following when presenting numerical results.
It is interesting to note that a for a dipole parametrization that has an UV anomalous dimension, i.e. ln D xy ∼ −(x − y) 2γ , corresponding to h Fits to HERA data using the BK equation favor values γ 1 for the initial condition, which evolves to γ 1 during the evolution. All of these are in the region γ ≥ 1/2 required by the positivity bound . At the limiting value γ = 1 of the MV model the analytical structure changes: h (1) ⊥ changes sign and for γ ≥ 1 the convergence of the Fourier-integral for the coordinate space distribution G (1) in terms of the momentum space one starts to require regularization; see the related discussion in [71].

Expectation values of two-point correlation functions
We now turn to the evaluation of correlation functions of ε(x)ε(y) andν(x)ν(y) characterizing local fluctuations of the energy density and divergence of the Chern-Simons current in the transverse plane. By performing the same steps as outlined above, we obtain for the correlation functions the following expressions We now have to evaluate correlation functions of the gluon field Even though it is in principle possible to evaluate such objects numerically in Gaussian models as discussed e.g. in [72], we will follow a different approach in order to obtain some analytic insight. Namely, we will assume that the four-point correlation functions of the gluon fields can be factorized into products of two-point correlation functions such that and similarly for the second nucleus Stated differently, this procedure corresponds to a factorization of the relevant double parton distribution into all possible products of single parton distributions. We note that this approximation scheme has been frequently employed in the literature e.g. in the context of di-hadron correlations (Glasma graphs) [73][74][75][76][77][78][79] and the quality of approximation has been investigated e.g. in [80]. Based on the above expression for the four-point correlation functions of the gluon fields, we can then proceed to evaluate the color structures in the expressions. Distinguishing the terms by connected and disconnected contractions as indicated in Eqs.
where we used the identities tr T a adj T b adj = N c δ ab and tr T a adj T b adj T a adj T c adj = 1 2 N 2 c δ bc to evaluate the final expressions. Collecting all the different terms we then obtain for the correlation function where we note that -as usual -all terms involving a connected contraction are suppressed by a factor 1/(N 2 c − 1) relative to the fully disconnected contribution. By performing also all of the contractions of the transverse tensors, we obtain after some algebra our result for the two point correlation function of the energy density which receives three distinct contributions, related to the disconnected-connected, connecteddisconnected and connected-connected contributions 5 . In contrast for the Chern-Simons correlator, all disconnected contractions vanish identically and only the connected-connected contractions give rise to a non-vanishing contribution. Our final result and the central result of this paper reads . 5 We note that the above expression corrects the earlier result of [43], where the connected-connected term in the last two lines was given incorrectly. We have checked the calculation of Ref. [43] step-by-step.
In the notation of the reference, we find that the pre-factor of the fully connected contribution to M1 should be 3/16 instead of 3/8 and that M5 + M6 ] featuring a relative minus sign between the unpolarized and linearly polarized contributions.
We note that the unpolarized and linearly polarized distributions contribute with different relative signs. However, as discussed in Sec. 3.2, at sufficiently small distances |x − y| the unpolarized contribution dominates and the correlation function is manifestly positive. Of course, the calculation outlined above can also be performed more or less entirely using modern computer algebra tools such as FeynCalc [81,82] or Form [83] and we have cross-checked our results in this way.

Correlators in momentum space
Even though we have performed the entire calculation above in coordinate space, our calculation is in fact closely related to the Glasma graph analysis of double inclusive particle production. In order to illustrate the similarities and differences it is useful to generalize our previous expressions to finite time by approximating the dynamics in the forward lightcone (τ > 0) in terms of the free field evolution. Based on the linearized evolution equations for abelian gauge fields the dynamics of the two independent polarizations corresponding to non-zero E η and respectively B η at τ = 0 + decouples from each other. By matching the general solution of Eq.

2)
Staying consistently at lowest order in the abelian approximation to the dynamics in the forward light-cone, the non-abelian field strength can be determined as One subtle issue is that the quality of the abelian approximation for the dynamics in the forward light-cone depends on the gauge choice. Even though the above expressions show that dynamics in the abelian approximation can be entirely formulated in terms of correlation functions of chromo-electric and chromo-magnetic fields, objects such as E a (τ = 0 + , x)E b (τ = 0 + , y) are in fact not gauge invariant. One natural gauge choice is the transverse Coulomb gauge ∂ i A i (τ = 0 + , x) = 0 which minimizes the transverse gauge field amplitudes, and it has been established from numerical simulations in [85] that the effects of final state interactions at τ > 0 become small in this gauge.
It is not generally known how to find the gauge transformation to Coulomb gauge analytically. However, the problem becomes considerably simpler in the case where either the projectile or target can be considered as dilute [86,87]. Specifically, if this is the case for the second nucleus (V x = 1 + igA (V ) (x)), a gauge transformation with V † U † yields the desired result to leading order in the dilute expansion. One finds that in this case, the non-vanishing components of the field strength tensor are given by and performing the transformation to Fourier space, the field-strength bi-linears ε(τ, x) andν(τ, x) can be compactly expressed as where k stands for Similarly, the two-point correlation functions of interest take the following form ε(τ, x)ε(τ, y) = (ig) 4 p,p k,k qq l,l × e i(p+k+p+k)x e i(q+l+q+l)y .

(4.13)
At early times τ 1/Q s the products of Bessel functions are dominated by J 2 0 ≈ 1, corresponding to the limit discussed in Sec. 3. Beyond early times only the disconnected contribution has a delta function setting p = −p, k = −k etc in such a way that the Bessel functions are arranged into combinations J 2 0 (x) + J 2 1 (x) with the same argument x. Based on the approximate relation J 2 0 (x) + J 2 1 (x) ≈ 2/(πx) one then obtains the usual behavior of the energy density as ε(τ ) ∼ 1/τ . On the other hand, simplifications of this nature do not occur for the disconnected contributions, and the Bessel functions oscillate out of phase. Hence, we expect the correlation signal for the energy density and the divergence of the Chern-Simons current in coordinate space to vanish for τ 1/Q s .

Glasma graph two gluon correlation
While the coordinate space correlation can be argued to vanish at τ 1/Q s , the situation is different for particle production, which is measured in momentum space. Here one integrates over the coordinates x, y and the corresponding x, y for the conjugate amplitude. This gives an additional momentum conservation delta function, which always sets p + k = −(p + k) and q + l = −(q + l) also for the connected contributions in the expressions analogous to Eqs. (4.12), (4.13) (for an illustration of the momentum flow see Fig. 1). This leads to the "glasma graph" [77] momentum space correlation structure. Even though this has not been the main focus of our paper, it is illustrative to derive this momentum space correlation signal here. This will clarify the relation of the calculation of Sec. 3 to the earlier literature on these "Glasma graph" correlations [73][74][75][76][77][78][79].
In order to obtain single and double inclusive particle spectra at leading order accuracy in the LSZ formalism one usually considers the limit τ → ∞ and projects gauge fixed equaltime correlation functions onto plane wave modes ξ µ,(λ) k (τ ) according to (4.14) By use of the explicit form of the plane wave solutions in transverse Coulomb gauge [88] ξ |k|τ H 1 (|k|τ ) .
the above expression evaluates to where in the dilute-dense regime the correlation functions in Coulomb gauge are given by The single inclusive gluon spectrum is obtained by evaluating the expectation value of (4.17) directly in momentum space as where we defined the Dipole gluon distribution 7   Figure 1. Examples of (a) completely disconnected diagram and examples of (b) a disconnectedconnected correlation (i.e. a "rainbow diagram" in the terminology of [74]) and (c) a connectedconnected correlation in the Glasma graph approximation. For particle production, the coordinates in the amplitude x, y are different from those in the conjugate amplitude, x, y, and are related by the momenta of the produced gluon. For the energy density and axial charge correlators, on the other hand, we integrate over momenta of the final state gluons, setting x = x, y = y.
Note that the dipole distribution is explicitly proportional to the momentum. Thus in a decomposition into polarization states similarly as for the Weiszäcker-Williams distribution, the unpolarized and polarized distributions are equal and there is only one scalar Using these expressions we obtain the following result for the single inclusive distribution We are now in a situation to repeat the calculation of double inclusive gluon production in [75] in our notations. Since the Glasma graph contribution to the double inclusive spectrum is simply given by the square of the single inclusive spectrum, we obtain where we inserted the explicit expressions for the chromo-electric and chromo-magnetic fields in order to make the similarities and differences with Eq. (4.13) most apparent. One immediately observes that both Eq. (4.13) and Eq. (4.23) involve the same correlation function of the gluon fields in momentum-space, allowing for the same interpretation in terms of a diagrammatic analysis. Specifically, the various different contractions in the projectile and target fields can be associated with the usual Glasma graphs as illustrated in Fig. (1). Even though the diagrammatics is essentially the same for double inclusive production and two-point correlation functions of local operators, there are of course some crucial differences in the calculation. Besides the appearance of a different operator structure in the middle of Eq. (4.23), another key difference is that for the local operator correlation function ν(x)ν(y) (and similarly for ε(x)ε(y) ) all expressions are to be evaluated at the same coordinates x = x and y = y in the amplitude and complex conjugate amplitude. Moreover, for the double inclusive spectrum the relevant correlation function of Wilson lines is given by can be decomposed into a complete set of color singlet structures [89]. Evaluating the full color structure is however quite challenging, and following [73][74][75] one usually resorts to an approximation of the full color structure in terms of the leading components in a dilute expansion. Specifically one expands the adjoint Wilson lines in Eq. (4.25) to lowest order in the target fields and performs a Gaussian averaging in terms of the fields A (U ) according to (4.32) where we dropped the contributions of all terms proportional to delta functions of a single momentum, i. e. δ(p), δ(p), δ(q), δ(q), as these do not contribute to particle production. We note that the approximation of the adjoint four point function in Eqns. (4.29), (4.30), (4.31), (4.32) is equivalent to the approximation used for the four-point function of the Weiszäcker-Williams field in Eq. (3.43) to leading order in the dilute limit. However, they correspond to different selective resummations of higher order terms away from the dilute limit. Ultimately this difference originates in the approximations used for the higher point functions of Wilson lines in Eqs. (3.32) and (4.26). We also stress that one cannot perform a naive decomposition of the four-point function of adjoint Wilson lines (4.25) into pairwise contractions of nonsinglet 2-point functions. Such a procedure would, for example, not reproduce the correct N c counting in the high momentum dilute limit, which we can check using the dilute approximation. Thus the "glasma graph" appproximation must be used with care, since it only really works in the dilute limit. In particular, we have not been able to find a k T -factorized expression for the 2-particle correlation function in the dilute-dense pA-case although one, involving the dipole distribution, does exist for the single gluon cross section.
Evaluating the individual terms, we obtain the following result for the contributions to the double-inclusive spectrum where the "disconnected-connected" (DC), "connected-disconnected" (CD) and (symme-try/asymmetric) "connected-connected" (CC S /CC A ) contributions are given by (4.37) One interesting feature of Eqs. (4.33),(4.34), (4.35), (4.36) and (4.37) -which is also visible in Eq. (3.44) -is the polarization structure on the proton side, which we maintained in full generality. One sees that for the "disconnected" contribution on the proton side (the CD-term) only the unpolarized gluon distribution appears, whereas on the "connected" side (DC and CC-terms) one is sensitive to the sum of the squares of the unpolarized and linearly polarized distributions. This is a subtle effect of a full treatment of a nontrivial linear polarization structure on the two-gluon correlations in momentum space.

Discussion
We now return to the central objective of this paper -to characterize axial charge production in the Glasma. Based on our calculation in Sec. 3.3, we find that the expectation value of the divergence of the Chern-Simons current ν(x) = 0 vanishes identically, such that on average no imbalance axial charge imbalance is created. However, the variance ν(x)ν(y) is finite, such that local fluctuations of the axial charge density should be expected on an event-by-event basis. Specifically, in the GBW saturation model, we obtain the following result for the correlation functions which is depicted in the left panel of Fig. 2. We note that except for the 1/(N 2 c − 1) suppression factor characteristic for fluctuations, there is no parametric suppression of ν(x)ν(y) compared to the energy density ε(x) , indicating that locally Glasma flux tubes can induce a significant imbalance of the axial charge density. However, it is also evident from Fig. 2 that the correlation length of these Glasma flux tubes in the transverse plane is microscopically small ∼ 1/Q s -such that a large number of uncorrelated domains should be expected in a realistic event. Besides the analytic results obtained in the GBW saturation model, the right panel of Fig. 2 shows the same quantities calculated in the MV model (see Sec. 3.2 for details).
Based on the above results for the source for axial charge production, we can also estimate local fluctuations of the axial charge density. Using our approximate treatment in Eq. (2.7) we find that for times τ 1/Q s the local fluctuations can be estimated as whereas fluctuations of the global amount of axial charge are suppressed by the overall number of Glasma flux tubes 1/Q 2 s S ⊥ and approximately given by where κ = π(44 ln(2) − 27 ln (3)) ≈ 2.6262 and S ⊥ denotes the transverse size of the overlap area.
We also note that our result in Eq. (5.3) can directly be used to model the initial conditions of the axial charge density dN 5 d 2 xdη in anomalous hydrodynamics or other calculations that attempt to relate anomalous transport phenomena to experimental measurements. We start from the assumption that on an event-by-event basis one knows the average energy density profile (x) as a function of the transverse coordinates, e.g. from a Monte Carlo Glauber model. This energy density profile should be thought of as the "average" energy density in the the sense that color charge fluctuations at the scale Q s are not included. The fluctuations at longer length scales, such as those resulting from fluctuations of the positions of nucleons inside the nucleus, should be averaged over separately as an external input to our calculation. Assuming that the energy density profile is sampled at a discrete set of points x in the transverse plane, one straightforward way to generate configurations  By sampling individual configurations of the axial charge distribution according to where ξ(z) are uncorrelated random numbers with zero mean ξ(z) = 0 and unit variance ξ(z)ξ(z ) = δ z,z , it is then straightforward to verify that the correlation function is correctly reproduced on average. Similarly, our result in Eq. (5.1) can also be used to include additional sub-nucleonic Q s -scale fluctuations of the energy density (x) (y) − (x) (y) on top of the average energy density profile (x) by following the same procedure outlined above. This provides a simplistic way to include the kind of Q s -scale energy density fluctuations that are present in the IPglasma model [27,49], although the analytic expressions used here are just approximations of the full numerical result. We emphasize that the procedure can be applied to any model or parametrization for the initial energy density at very early times τ 1/Q s . Even if the initial average energy density does not come from an explicit saturation model calculation, one can estimate the corresponding saturation scale by solving for Q s from the initial energy density ε ≈ 1 s . We illustrate this procedure in Fig. 3, with the example of a peripheral Pb+Pb event. Based on the average energy density profile obtained from the T R ENTO event generator [90] shown in the first panel of Fig. 3, we include fluctuations of the energy density and axial charge distribution following the above procedure. Despite the fact that average energy density profile is rather smooth, with typical variations on size scales ∼ fm, sub-nucleonic fluctuations give rise to fluctuations of the energy density at length scales ∼ 1/Q s as can be seen from the central panel of Fig. 3. Similarly, variations of the axial charge distribution due to Glasma flux tubes occur on microscopic length scales with a characteristic size ∼ 1/Q s . However, due to the approximate boost invariant nature of the Glasma fields, these structures are elongated in rapidity. It will be interesting to see from phenomenological calculations whether such small structures can have a sizeable effect on hadronic observables. In order to facilitate the use of our result in this context, we provide the source code for generating axial charge distributions as in Fig. 3 as supplementary material.

Conclusions & Perspectives
Based on known analytic solutions for the Glasma fields we calculated energy and axial charge fluctuations at early times τ 1/Q s after the collision of heavy nuclei at high energies. Our calculation generalizes the earlier work of [43] to a more general structure for the gluon distribution and, more importantly, to derive an expression for the Chern-Simons correlator. Generally, we find that the expressions for energy and axial charge fluctuations in Eqs. (3.30) involve the correlation function of two Weiszäcker-Williams (WW) gluon distributions, represented as a correlator of eight light-like Wilson lines for each nucleus. We evaluated this correlation function in the "Glasma graph" approximation, where the relevant double parton distribution is factorized into a product of single parton distributions. Based on previous calculations [80], we expect the Glasma graph approximation to be quite close to the full result. Extending this calculation to the full nonlinear Gaussian treatment would require working out an eight-point function of Wilson lines in the similar way as the four-point function in Appendix A.2. Based on the complexity of the expressions it appears unlikely that this could be done analytically, but a numerical evaluation similar to the recent one in [72] should certainly be feasible. We also note that, based on our primary interest of applications to the collision of large nuclei, we neglected some more subtle effects related to position-momentum correlations in the gluon distribution (see e.g. [91,92]), which may be interesting to investigate in further applications to small systems.
Our result in Eqs. In view of possible phenomenological applications of our result, we provided a practical algorithm to use our result to implement quantitatively the axial-charge density fluctuations in the Glasma. Of course, this relation relies on a rough treatment of the time dependence of the Chern-Simons charge in the Glasma (c.f. Sec. 2) and ultimately a full classical Yang-Mills calculation including dynamical fermions along the lines of [19,93,94] as well as a realistic geometry will be needed. Even with this approximation our result should, however, enable a better control of the initial conditions for anomalous hydrodynamics simulations or other calculations that are needed to relate these ideas to experimental measurements. We caution, however, that axial charge changing processes e.g. due to sphaleron transitions or thermal fluctuations of the field strength tensor continue to take place throughout the entire space-time evolution of the Quark-Gluon Plasma. Clearly such effects should also be included in realistic model calculations and further theoretical progress will be required.
We finally note that several calculations similar to ours have been performed for momentum space gluon correlations based on the Glasma graph approximation [73][74][75][76][77][78][79]. However, these calculations are performed in the dilute limit and do not give access to the linear polarization structure of the gluon distribution. Focusing only on the coordinate space correlation structures of the fields at τ 1/Q s enables us to do a calculation in a manifestly gauge invariant way and more cleanly elucidate the role of the gluon polarization. The relation between our present work and the Glasma graph calculations of ridge correlations is explained in more detail in Sec. 4.

A Evaluation of Weiszäcker-Williams distribution in Gaussian models
We start by decomposing the gluon fields α i x over the Lie Algebra such that the Weizsäcker-Williams distribution is given by By re-expressing the derivatives in terms of new coordinates x, y and making use of the SU (N c ) Fierz identity t a ij t a kl = the relevant correlation function of Wilson lines then take the form Clearly the second term vanishes upon taking the derivative and setting coordinates x = x and y = y equal to each, as U x ∂ i U † x are elements of the Lie algebra and thus traceless. We are then left with evaluating the first term involving the quadrupole correlator.

A.1 Evaluation of the Wilson line correlators in Gaussian model
We perform the Gaussian averaging of the correlators of Wilson line, by expressing the usual Gaussian integral over color charges in terms of a stochastic process in the evolution variable z ∈ [0, 1] such that the Wilson lines at z = 0 are given by V x (z = 0) = 1 and in each step where ξ a x are stochastic variables with ξ a x (z)ξ b y (z ) = 1 g 2 C F δ ab λ xy δ(z − z ) , (A. 6) where C F = (N 2 c − 1)/(2N c ) is the fundamental Casimir. Starting for simplicity with the dipole operator, we can then evaluate where we introduced the correlation functions such that the dipole correlator is simply given by

A.2 Quadrupole & Dipole-Dipole correlators
Similarly for the quadrupole, we obtain upon use of the SU (N c ) Fierz identity the evolution equation where the transition function T xx,yy given by T xx,yy = λ xy + λ yx − λ xy − λ xy = G xy + G yx − G xy − G xy (A.11) This has to be supplemented with the evolution equation for the dipole-dipole correlator T xy,yx tr U x U † y tr U † x U y (A.12) with the transition function T xy,yx given by We obtain the coupled set of evolution equations where the evolution operator M (x, x, y, y) takes the form M (x, x, y, y) = G xy + G yx − 1 Of course, for this simple example we could easily calculate the full expression as done in [95]. However, for our purpose it is more useful to first take the derivatives and set the coordinates x = x and y = y equal to each other, such that the relevant expression Similarly, using the relations we obtain the derivative of the evolution operator as