Topological charge fluctuations in the Glasma

TL;DR

This paper analytically calculates the early axial charge fluctuations in the Glasma, revealing significant differences from Gaussian approximations at larger distances due to non-linear QCD dynamics.

Contribution

It provides the first-principles analytical computation of the Chern-Simons current correlator in the Glasma, highlighting non-Gaussian effects at larger transverse separations.

Findings

Exact correlator decays slower than Gaussian approximation at large distances.

Non-linear dynamics cause color screening at relatively large transverse scales.

Results support previous indications of large-distance correlation suppression.

Abstract

The early-time evolution of the system generated in ultra-relativistic heavy ion collisions is dominated by the presence of strong color fields known as Glasma fields. These can be described following the classical approach embodied in the Color Glass Condensate effective theory, which approximates QCD in the high gluon density regime. In this framework we perform an analytical first-principles calculation of the two-point correlator of the divergence of the Chern-Simons current at proper time $\tau\!=\!0^+$, which characterizes the early fluctuations of axial charge density in the plane transverse to the collision axis. This object plays a crucial role in the description of anomalous transport phenomena such as the Chiral Magnetic Effect. We compare our results to those obtained under the Glasma Graph approximation, which assumes gluon field correlators to obey Gaussian statistics. While this approach proves to be equivalent to the exact calculation in the limit of short transverse separations, important differences arise at larger distances, where our expression displays a remarkably slower fall-off than the Glasma Graph result ($1/r^4$ vs.\ $1/r^8$ power-law decay). This discrepancy emerges from the non-linear dynamics mapping the Gaussianly-distributed color source densities onto the Glasma fields, encoded in the classical Yang-Mills equations. Our results support the conclusions reached in a previous work, where we found indications that the color screening of correlations in the transverse plane occurs at relatively large distances.

Equations

$$Q_{\text{w}}=\frac{g^2}{16{\pi }^2}\int d^{4}x\text{Tr}\left\{F_{\mu \nu }(x){\tilde{F}}^{\mu \nu }(x)\right\},$$

$$\frac{d{N}_5}{dt}=\int d^{3}x{\partial }_{\mu }{j}_5^{\mu }(x)=-\frac{g^2{N}_f}{8{\pi }^2}\int d^{3}x\text{Tr}\left\{F_{\mu \nu }(x){\tilde{F}}^{\mu \nu }(x)\right\}=\frac{g^2{N}_f}{2{\pi }^2}\int d^{3}x\dot{\nu }(x),$$

$$\left[D_{\mu },F^{\mu \nu }\right]={\delta }^{\nu +}{\rho }^a(x^{-},x_{\text{vbox}\perp })t^a,$$

$$⟨\mathcal{O}[\rho ]⟩=\frac{1}{\mathcal{N}}\int [d\rho ]W[\rho ]\mathcal{O}[\rho ],$$

$$⟨{\rho }^a(x^{-},x_{\text{vbox}\perp }){\rho }^b(y^{-},y_{\text{vbox}\perp }){⟩}_{\text{scaleto}\text{MV}0.12cm}={\mu }^2(x^{-}){\delta }^{ab}\delta (x^{-}-y^{-}){\delta }^2(x_{\text{vbox}\perp }-y_{\text{vbox}\perp }).$$

$$⟨{\rho }^a(x^{-},x_{\text{vbox}\perp }){\rho }^b(y^{-},y_{\text{vbox}\perp })⟩$$

$$\equiv \lambda (x^{-},b_{\text{vbox}\perp }){\delta }^{ab}\delta (x^{-}-y^{-})f(x_{\text{vbox}\perp }-y_{\text{vbox}\perp }),$$

$$A_{\text{scaleto}14.2pt,\text{scaleto}24.2pt}^{\pm }(x^{\mp },x_{\text{vbox}\perp })$$

$$A_{\text{scaleto}14.2pt,\text{scaleto}24.2pt}^i(x^{\mp },x_{\text{vbox}\perp })$$

$$\equiv \theta (x^{\mp }){\alpha }_{\text{scaleto}14.2pt,\text{scaleto}24.2pt}^i(x_{\text{vbox}\perp }),$$

$$U_{\text{scaleto}14.2pt,\text{scaleto}24.2pt}(x^{\mp },x_{\text{vbox}\perp })={\text{P}}^{\mp }\exp \left\{-ig{\int }_{\text{scaleto}-\infty 5pt}^{x^{\mp }}d{z}^{\mp }\int d{z}_{\text{vbox}\perp }^{2}G(z_{\text{vbox}\perp }-x_{\text{vbox}\perp }){\tilde{\rho }}_{\text{scaleto}14pt,\text{scaleto}24pt}(z^{\mp },z_{\text{vbox}\perp })\right\},$$

$$A^{\pm }(x^{\mp },x_{\text{vbox}\perp })$$

$$A^i(x^{\mp },x_{\text{vbox}\perp })$$

$${\alpha }^i(\tau =0^{+},x_{\text{vbox}\perp })$$

$$\alpha (\tau =0^{+},x_{\text{vbox}\perp })$$

$$E^z(\tau =0^{+},x_{\text{vbox}\perp })=$$

$$B^z(\tau =0^{+},x_{\text{vbox}\perp })=$$

$${\dot{\nu }}_{\text{scaleto}04.2pt}(x_{\text{vbox}\perp })$$

$$=-g^2{\delta }^{ij}{ϵ}^{kl}{\alpha }_{\text{scaleto}14.2ptx}^{i,a}{\alpha }_{\text{scaleto}24.2ptx}^{j,b}{\alpha }_{\text{scaleto}14.2ptx}^{k,c}{\alpha }_{\text{scaleto}24.2ptx}^{l,d}\text{Tr}\{[t^a,t^b][t^c,t^d]\}$$

$$=\frac{g^2}{2}{\delta }^{ij}{ϵ}^{kl}{f}^{abn}{f}^{cdn}{\alpha }_{\text{scaleto}14.2ptx}^{i,a}{\alpha }_{\text{scaleto}24.2ptx}^{j,b}{\alpha }_{\text{scaleto}14.2ptx}^{k,c}{\alpha }_{\text{scaleto}24.2ptx}^{l,d},$$

$${\alpha }^i(x_{\text{vbox}\perp })$$

$$={\int }_{-\infty }^{\infty }d{z}^{-}\int d{z}_{\text{vbox}\perp }^{2}G(z_{\text{vbox}\perp }-x_{\text{vbox}\perp }){\partial }^i{\tilde{\rho }}^a(z^{-},z_{\text{vbox}\perp })U^{ab}(z^{-},x_{\text{vbox}\perp })t^b\equiv {\alpha }^{i,b}(x_{\text{vbox}\perp })t^b,$$

$$⟨{\dot{\nu }}_{\text{scaleto}04.2pt}(x_{\text{vbox}\perp })⟩=\frac{g^2}{2}{\delta }^{ij}{ϵ}^{kl}{f}^{abn}{f}^{cdn}⟨{\alpha }_{\text{scaleto}14.2pt}^{i,a}(x_{\text{vbox}\perp }){\alpha }_{\text{scaleto}14.2pt}^{k,c}(x_{\text{vbox}\perp })⟩⟨{\alpha }_{\text{scaleto}24.2pt}^{j,b}(x_{\text{vbox}\perp }){\alpha }_{\text{scaleto}24.2pt}^{l,d}(x_{\text{vbox}\perp })⟩,$$

$$⟨{\alpha }^{i,a}(x_{\text{vbox}\perp }){\alpha }^{j,b}(x_{\text{vbox}\perp })⟩$$

$$L(r_{\text{vbox}\perp })\equiv \int d{z}_{\text{vbox}\perp }^{2}d{u}_{\text{vbox}\perp }^{2}G(z_{\text{vbox}\perp }-x_{\text{vbox}\perp })G(u_{\text{vbox}\perp }-y_{\text{vbox}\perp })f(z_{\text{vbox}\perp }-u_{\text{vbox}\perp }).$$

$${\partial }_x^i{\partial }_y^{j}L(r_{\text{vbox}\perp })=A(r_{\text{vbox}\perp }){\delta }^{ij}+B(r_{\text{vbox}\perp })\left(\frac{{\delta }^{ij}}{2}-\frac{r^i{r}^j}{r^2}\right).$$

$$A(r_{\text{vbox}\perp })$$

$$B(r_{\text{vbox}\perp })$$

$$⟨U^{a^{\prime }a}(x^{-},x_{\text{vbox}\perp })U^{a^{\prime }b}(x^{-},y_{\text{vbox}\perp })⟩$$

$$\equiv {\delta }^{ab}{C}_{\text{adj}}^{(2)}(x^{-};x_{\text{vbox}\perp },y_{\text{vbox}\perp }).$$

$$\Gamma (r_{\text{vbox}\perp })=2(L(0_{\text{vbox}\perp })-L(r_{\text{vbox}\perp })).$$

$$⟨{\alpha }^{i,a}(x_{\text{vbox}\perp }){\alpha }^{j,b}(x_{\text{vbox}\perp })⟩=-\frac{1}{2}{\delta }^{ab}{\delta }^{ij}{\bar{\mu }}^{2}h(x_{\text{vbox}\perp }){\partial }^{2}L(0_{\text{vbox}\perp })=-\frac{1}{2}{\delta }^{ab}{\delta }^{ij}\bar{\lambda }(x_{\text{vbox}\perp }){\partial }^{2}L(0_{\text{vbox}\perp }).$$

$$\underset{r\to 0}{\lim }{\partial }_x^i{\partial }_y^{j}L(r_{\text{vbox}\perp })=\frac{{\delta }^{ij}}{2}\int \frac{d^2{q}_{\text{vbox}\perp }}{(2\pi {)}^2}\hat{f}(q_{\text{vbox}\perp })\frac{1}{q^2}\equiv -\frac{1}{2}{\delta }^{ij}{\partial }^{2}L(0_{\text{vbox}\perp }).$$

$$⟨{\dot{\nu }}_{\text{scaleto}04.2pt}(x_{\text{vbox}\perp })⟩=\frac{g^2}{8}({\partial }^{2}L(0_{\text{vbox}\perp }){)}^2{\bar{\lambda }}_{\text{scaleto}14.2pt}(x_{\text{vbox}\perp }){\bar{\lambda }}_{\text{scaleto}24.2pt}(x_{\text{vbox}\perp }){\delta }^{ij}{ϵ}^{kl}{f}^{abn}{f}^{cdn}{\delta }^{ac}{\delta }^{ik}{\delta }^{bd}{\delta }^{jl}=0,$$

$${\dot{\nu }}_{\text{scaleto}04.2pt}(x_{\text{vbox}\perp }){\dot{\nu }}_{\text{scaleto}04.2pt}(y_{\text{vbox}\perp })=\frac{g^4}{4}{\delta }^{ij}{ϵ}^{kl}{\delta }^{i^{\prime }{j}^{\prime }}{ϵ}^{k^{\prime }{l}^{\prime }}{f}^{abn}{f}^{cdn}{f}^{a^{\prime }{b}^{\prime }m}{f}^{c^{\prime }{d}^{\prime }m}{\alpha }_{1x}^{i,a}{\alpha }_{2x}^{j,b}{\alpha }_{1x}^{k,c}{\alpha }_{2x}^{l,d}{\alpha }_{1y}^{i^{\prime },a^{\prime }}{\alpha }_{2y}^{j^{\prime },b^{\prime }}{\alpha }_{1y}^{k^{\prime },c^{\prime }}{\alpha }_{2y}^{l^{\prime },d^{\prime }},$$

$$⟨{\dot{\nu }}_{\text{scaleto}04.2pt}(x_{\text{vbox}\perp }){\dot{\nu }}_{\text{scaleto}04.2pt}(y_{\text{vbox}\perp })⟩=\frac{g^4}{4}{ϵ}^{kl}{ϵ}^{k^{\prime }{l}^{\prime }}{f}^{abn}{f}^{cdn}{f}^{a^{\prime }{b}^{\prime }m}{f}^{c^{\prime }{d}^{\prime }m}⟨{\alpha }_x^{i,a}{\alpha }_x^{k,c}{\alpha }_y^{i^{\prime },a^{\prime }}{\alpha }_y^{k^{\prime },c^{\prime }}{⟩}_{\text{scaleto}14.2pt}⟨{\alpha }_x^{i,b}{\alpha }_x^{l,d}{\alpha }_y^{i^{\prime },b^{\prime }}{\alpha }_y^{l^{\prime },d^{\prime }}{⟩}_{\text{scaleto}24.2pt}.$$

Full text

1 Introduction

High-energy heavy ion collisions (HICs from now on) generate a hot, chiral-symmetric medium –the Quark Gluon Plasma (QGP)– where parity (P) and charge-parity (CP)-violating fluctuations are expected to happen with relatively high probability due to the chiral anomaly of QCD Fukushima:2008xe . The anomalous term weights gauge field configurations according to their associated value of the topological invariant known as winding number (or topological charge):

[TABLE]

where $F^{\mu\nu}$ and $\tilde{F}^{\mu\nu}\!=\!\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$ correspond respectively to the field strength tensor and its dual. $Q_{\text{w}}$ labels degenerate (but topologically inequivalent) vacuum states separated by potential barriers that greatly suppress the transition probability except at high temperatures, such as the ones reached in the QGP phase (for a review of the topological aspects of gauge field theories, see Lenz:2001me ). The tunneling transitions, mediated by localized field configurations called sphalerons, induce a transformation of left- into right-handed quarks (and viceversa) whose rate can be computed in the massless fermion limit as the spatial integral of the chiral anomaly:

[TABLE]

where $N_{5}$ and $j_{5}^{\mu}$ correspond respectively to the axial charge density and current, $N_{f}$ is the number of flavors and $\dot{\nu}$ is the divergence of the Chern-Simons current.

HICs thus provide an appropriate environment –the QGP– for the manifestation of the non-trivial topological structure of QCD. Specifically, off-central HICs yield favourable conditions for the search of observable signs of such features. These collisions give rise to large background electromagnetic fields, which in the presence of deconfined chirally-imbalanced matter may induce a separation of positive and negative charges along the direction of angular momentum Kharzeev:2007jp . This effect –known as the Chiral Magnetic Effect (CME)– thus creates a preferential direction for the emission of charged particles that would in turn translate into non-trivial azimuthal correlations in the hadronic spectrum. The search for such signatures of this and other transport phenomena connected to the chiral anomaly (generically called anomalous transport phenomena) has been carried at the Relativistic Heavy-Ion Collider (RHIC) and at the Large Hadron Collider (LHC) Kharzeev:2015znc . Although these experiments have provided numerous measurements that are indeed compatible with said phenomena Abelev:2009ac ; Abelev:2012pa ; Adamczyk:2014mzf ; Adamczyk:2015eqo ; Tribedy:2017hwn , the presence of large background effects (e.g. transverse momentum Bzdak:2010fd and local charge Schlichting:2010qia conservation, which also give rise to intrinsic back-to-back correlations; and final state interactions Ma:2011uma ) prevent from drawing definite conclusions. Hence, there is a strong interest from the high-energy QCD community in reducing this uncertainty. Significant advances have been achieved on the experimental side, including the development of different detection techniques Voloshin:2004vk ; Bzdak:2012ia ; Wen:2016zic ; Xu:2017qfs and, most recently, the implementation of an isobaric collision program at RHIC aimed at the isolation of CME backgrounds Huang:2018duz . Still, a thorough approach to this task demands for better theoretical constraints on the dynamical origin of correlations between detected particles.

Non-trivial final state correlations in heavy ion experiments have been shown to reflect not only the collective dynamics of QGP, but also the event-by-event fluctuations that characterize the initial phase of the collisions Luzum:2013yya . These early fluctuations provide a natural source of topological charge that competes with sphaleron transitions. Although the latter are known to dominate axial charge production in the QGP phase, throughout the pre-equilibrium stage both mechanisms are likely to yield a significant contribution111It has been argued that in the early stage of the collision sphaleron transitions are suppressed due to the boost invariance of the generated fields Kharzeev:2001ev . However, as the system evolves towards thermalization and boost invariance wears down, they would be significantly enhanced Kharzeev:2001ev ; Mace:2016svc . Whether or not event-by-event fluctuations dominate over this or other mechanisms of axial charge production –like thermal fluctuations of the field strength– is out of the scope of this paper. For more exhaustive discussions on this topic, the reader is referred to the aforementioned studies.. It is thus essential to understand and quantitatively constrain the influence of each source in the experimentally observed correlations.

However, there is currently no consensus about either the origin or the practical description of early event-by-event fluctuations. Different prescriptions arise from a wide variety of phenomenological models whose goal is to provide initial conditions for the quasi-ideal relativistic hydrodynamical evolution of QGP. The numerical values of the physical quantities required as input by these models are constrained by agreement with data, sometimes varying greatly from one model to another Luzum:2013yya ; Retinskaya:2013gca . Such discrepancy introduces a significant amount of uncertainty in both the precision and physical interpretation of most phenomenological studies of the expansion and cooling of QGP. This general issue is of particular importance for those hydrodynamical descriptions that mimic the effects induced by the chiral anomaly Newman:2005hd ; Son:2009tf ; Sadofyev:2010pr since event-by-event fluctuations are expected to contribute a significant fraction of the initial axial charge densities. Therefore, a better theoretical control of the early stages of the collision is essential to properly characterize the origin and description of anomalous transport phenomena.

The Color Glass Condensate (CGC) effective theory (see e.g. Weigert:2005us ; Gelis:2010nm for a review) is arguably the most promising framework for the description of the early phase of HICs. CGC describes the high density of small-$x$ gluons carried by nuclei as strong color fields whose dynamics obey the classical Yang-Mills equations. The classical approximation is founded in the fact that for very large occupation numbers the quantum fluctuations represent a negligible correction to the strong background field. This condition is at the base of the McLerran-Venugopalan model (MV model from now on) McLerran:1993ni ; McLerran:1993ka ; McLerran:1994vd , where nuclei are represented by an ensemble of SU($N_{c}$) color charges that act as sources of the classical fields. The MV model thus assumes a separation of degrees of freedom that is performed at an arbitrary light-cone momentum $\Lambda^{+}$: particles with $p^{+}\!>\!\Lambda^{+}$ are taken as hard color charges (which represent the valence quarks), and these generate the small-$x$ dynamical modes (the strong color fields), which satisfy $p^{+}\!\!<\!\Lambda^{+}$. CGC incorporates the means to compute the quantum corrections to the MV model via the B-JIMWLK equations, which describe the evolution of the theory with $\Lambda^{+}$. This framework has been extensively applied in the description of the early, non-equilibrium stage of HICs known as Glasma phase Lappi:2006fp ; Lappi:2006hq ; Dumitru:2008wn ; Fukushima:2011nq .

In a previous work Albacete:2018bbv we provided a first-principles calculation quantifying the size and extent of the transverse correlations of the energy-momentum tensor of Glasma at early times (and in a subsequent paper Giacalone:2019kgg the results of said work were applied in the description of anisotropic flow harmonic coefficients to excellent agreement with data measured at both RHIC and LHC). In this follow-up paper we extend the classical treatment to the divergence of the Chern-Simons current: where $E^{i}\!=\!-F^{0i}$ and $B^{k}\!=\!\frac{1}{2}\epsilon^{ijk}F^{ij}$ are, respectively, the Glasma chromo-electric and -magnetic fields, evaluated at an infinitesimal positive proper time $\tau\!=\!0^{+}$ at a point x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}} of the plane transverse to the collision axis. This object is proportional to the strong CP-violating term of the QCD Lagrangian, which is the source for local production of axial charge (Eq.​ (2)). In the present work we evaluate the correlation function \langle\dot{\nu}_{\scaleto{0}{4.2pt}}(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\dot{\nu}_{\scaleto{0}{4.2pt}}(y_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\rangle, which characterizes the early event-by-event fluctuations of the $\dot{\nu}$ distribution. In evaluating this object we find relatively long-range correlations that contrast with both naive expectations one could have from the MV model –where we assume local correlations at the level of color source distributions– and previously obtained results from the Glasma Graph approximation Lappi:2017skr –which assumes a linear mapping of the statistics followed by the color sources onto the Glasma fields. Although such a discrepancy mirrors the results found in Albacete:2018bbv , it is worth remarking that the calculations presented in this work yield an even larger difference with respect to those obtained under the Glasma Graph approximation.

This paper is organized as follows. In Sec.​ 2 we introduce a generalization of the MV model with explicit impact parameter dependence and relaxed transversal locality. In this framework we outline the solution to the Yang-Mills equations with two sources at an infinitesimal proper time after the collision $\tau\!=\!0^{+}$, which acts as boundary condition for the ensuing evolution in the future light-cone. In Sec.​ 3 we calculate the expectation value of the divergence of the Chern-Simons current in the previously presented framework. In Sec.​ 4 we compute its two-point correlator, comparing our results with the aforementioned calculation performed under the Glasma Graph approximation Lappi:2017skr . We also provide the first orders of the $N_{c}$-expansion, as well as the strict MV model limit of our expression. Note that a big part of the technical challenges faced during the calculations presented in this paper –such as the computation of non-trivial projections of the correlator of four Wilson lines in the adjoint representation and the decomposition of correlators of $m$ valence color sources and $n$ Wilson lines– were analyzed in depth in a previous work Albacete:2018bbv , being the reader referred to said paper for detailed derivations. Finally, in Sec.​ 5 we discuss the physical implications and potential phenomenological applications of our result, as well as its role in future studies.

2 Setup: The classical approach to gluon production in high-energy heavy ion collisions

In the MV model we represent the high density of small-$x$ gluons carried by each nucleus with gauge fields $A^{\mu}_{\scaleto{1,2}{5.8pt}}$ whose dynamics follow from the classical Yang-Mills equations:

[TABLE]

where $\rho^{a}$ represents the density of valence (large-$x$) partons and $t^{a}$ is the generator of $SU(N_{c})$ in the fundamental representation. The $\delta^{\nu+}$ factor indicates that the source generates a color current only in the + direction. This suggests a physical picture of the interaction where the valence partons do not recoil from their light-cone trajectory as the gluons they continuously exchange with the medium are too soft to affect their motion (eikonal approximation).

The MV model accounts for the event-by-event fluctuations of color charge configurations by taking $\rho^{a}$ as an stochastic quantity with a certain probability distribution $W[\rho]$ associated as weight function. Thus, the observables are obtained as ensemble averages over the background classical fields:

[TABLE]

where ${\cal N}$ is a normalization constant equal to $\int[d\rho]W[\rho]$. The main assumption adopted in the MV model is that in nuclei with large mass numbers the valence parton configurations emerge from a large number of separate nucleons and therefore are uncorrelated. Thus, invoking the central limit theorem, this model approximates $W[\rho]$ with a Gaussian distribution, yielding the following fundamental result:

[TABLE]

Here, $\mu^{2}$ is a parameter proportional to the color source number density that enters as the variance of the Gaussian weight. However, as we intend to apply a more general approach, we choose to relax some of the approximations implied in Eq.​ (5) by considering the following, more general, two-point function of the sources:

[TABLE]

where we allow the possibility of finite, non-homogeneous nuclei by explicitly introducing an impact parameter (b_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}\!\equiv\!(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}\!+y_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})/2) dependence as previously done in Chen:2015wia . Also, we drop the assumption that interactions are local in the transversal plane by introducing an undetermined function f(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}\!-y_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) instead of a Dirac delta. These extensions of the original MV model might prove especially useful in subsequent phenomenological applications of our results.

Although there is no general solution for the Yang-Mills equations with two sources, the MV model provides an analytical expression of the gauge fields at proper time $\tau\!=\!0^{+}$ (i.e. an infinitesimal positive proper time after the collision). These are obtained in terms of the gauge fields that characterize each nucleus before the collision, which can be computed independently as:

[TABLE]

by solving the Yang-Mills equations with one source in the Fock-Schwinger gauge. Here $\tilde{\rho}$ is the color charge density in the covariant gauge and $U$ is the Wilson line, an SU($N_{c}$) element that represents the effect of the interaction with the classical gluon field over the fast valence partons in the eikonal approximation, i.e. a rotation in color space. It is defined as a path-ordered exponential:

[TABLE]

where G(z_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}\!\!-\!x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) is the Green’s function for the two-dimensional Laplace operator. The solution at $\tau\!=\!0^{+}$ is found by proposing the following ansatz:

[TABLE]

where we adopted the comoving coordinate system, defined by proper time $\tau\!=\!\sqrt{2x^{+}x^{-}}$ and rapidity $\eta\!=\!\frac{1}{2}\ln(x^{+}/x^{-})$. Then, we invoke a physical ‘matching condition’ that requires Yang-Mills equations to be regular in the limit $\tau\!\rightarrow\!0$ . In doing so, the following relations are obtained:

[TABLE]

which act as boundary conditions of the subsequent $\tau$-evolution222Several approaches of both analytical and numerical nature have been applied for this computation in the literature. For instance, in Fries:2006pv an analytical approximation based on an expansion of the Yang-Mills equations in powers of $\tau$ is proposed. However, this is out of the scope of the work presented in this paper.. From these expressions we can compute the chromo-electric and -magnetic fields at $\tau\!=\!0^{+}$ as:

[TABLE]

being these the only non-vanishing components. This peculiar configuration of boost-invariant longitudinal fields motivates the Glasma flux tube picture, which predicts short-range transverse correlations Iancu:2002aq .

3 One-point correlator of the divergence of the Chern-Simons current

Before evaluating the two-point function, in this section we will show that the expectation value of the divergence of the Chern-Simons current over the classical background fields is 0, indicating that there is no overall CP violation in the process. Although this result has been obtained previously in the literature Lappi:2006fp ; Lappi:2017skr , we deem it convenient to include this preface as it allows us to establish the notation used in the rest of the paper. Let us first write $\dot{\nu}_{\scaleto{0}{4.2pt}}$ in terms of the gluon fields:

[TABLE]

where we adopted the shorthand notation \alpha^{i,a}(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\!\equiv\!\alpha^{i,a}_{x}. As an intermediate step we expand the color structure of the gluon fields as:

[TABLE]

where we used the relation between Wilson lines in the fundamental and adjoint representations $U^{\dagger}t^{a}U\!=\!U^{ab}t^{b}$. The correlator of Eq.​ (16) factorizes like:

[TABLE]

as in the MV model the color source fluctuations of each nucleus are assumed to be independent. Thus, the building block of $\langle\dot{\nu}_{\scaleto{0}{4.2pt}}\rangle$ is the correlator of two gauge fields evaluated in the same transverse position, \langle\alpha^{i,a}(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\alpha^{k,c}(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\rangle. We calculate this object as a limit:

[TABLE]

where r\!=\!|r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}|\!=\!|x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}-y_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}|. Here we introduced the following function:

[TABLE]

From its symmetries and dimension, the double derivative of L(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) featured in Eq.​ (19) can be parameterized as:

[TABLE]

This formula accounts for an explicit separation of the contributions of the unpolarized A(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) and linearly polarized B(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) parts of the gluon distribution. We can express these coefficients in terms of f(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) explicitly as:

[TABLE]

where \hat{f}(q_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) is the Fourier transform of f(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}). As for the last factor in Eq.​ (19), it corresponds to the dipole function in the adjoint representation, which stems from the following correlator:

[TABLE]

Here we introduced the notation \bar{\lambda}(x^{-},x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\!=\!\int^{x^{-}}_{-\infty}dz^{-}\lambda(z^{-},x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}), as well as the following function:

[TABLE]

Computing the indicated limit explicitly, Eq.​ (19) yields:

[TABLE]

(For a more detailed calculation of this correlator we refer the reader to Albacete:2018bbv ). The factor \bar{\lambda}(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) corresponds to $\lambda$ integrated from $-\infty$ to $\infty$ in the longitudinal direction (in general, we will identify functions integrated in the longitudinal direction from $-\infty$ to $\infty$ by simply omitting their longitudinal dependence). The factor \partial^{2}L(0_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) is a model-dependent constant that emerges from the following limit:

[TABLE]

Some assumptions had to be made about the functions h(b_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}), f(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) introduced in Eq.​ (6) in order to arrive at the expressions presented in this section. Specifically, we take h(b_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) as a slowly varying function over lengths of the order of an infrared length scale $1/m$ (or smaller), which can be understood as a cut-off that imposes color neutrality at the nucleon size. In addition, we assume that f(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) behaves in such a way that its Fourier transform \hat{f}(k_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) tends to unity in the infrared limit (see Albacete:2018bbv for details).

When finally substituting our building block Eq.​ (26) into Eq.​ (18), we get:

[TABLE]

which vanishes due to the antisymmetric property of the Levi-Civita tensor. As mentioned earlier, this null average accounts for the Glasma state being globally CP-symmetric. However, as will be shown below, local axial charge fluctuations are expected to happen on an event-by-event basis. Our object of interest is therefore the two-point correlator of $\dot{\nu}_{\scaleto{0}{4.2pt}}$, whose computation we outline in the following section.

4 Two-point correlator of the divergence of the Chern-Simons current

In this section we describe the calculation of \langle\dot{\nu}_{\scaleto{0}{4.2pt}}(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\dot{\nu}_{\scaleto{0}{4.2pt}}(y_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\rangle, which characterizes the early fluctuations of the divergence of the Chern-Simons current in the transverse plane. As we did in the previous section, we start by expanding $\dot{\nu}_{\scaleto{0}{4.2pt}}$ in terms of the gluon fields:

[TABLE]

then, the correlator reads:

[TABLE]

Color algebra-wise, this expression presents the same level of complexity than the two-point correlator of the energy density, previously computed in Albacete:2018bbv . Happily, it features a much simpler transverse index structure, which we can rewrite as:

[TABLE]

yielding:

[TABLE]

The building block for this computation is the correlator of four gluon fields in two different transverse positions:

[TABLE]

where $1/\nabla^{2}$ is the shorthand notation we adopt to denote a convolution with the Laplacian Green’s function. The above correlator is a highly non-trivial object whose calculation poses a number of outstanding technical challenges such as the computation of non-trivial projections of the correlator of four Wilson lines in the adjoint representation and the decomposition of correlators of four valence color sources and four Wilson lines. For a detailed calculation the reader is referred to Albacete:2018bbv . Here we will simply indicate the result, for which we need to make some definitions. We start with the ‘connected’ function:

[TABLE]

which accounts for the contribution of correlations between the ‘external’ color source densities and those arranged inside the Wilson lines. The remaining terms stem from a complete factorization of the correlations of color source densities and Wilson lines:

[TABLE]

These contributions are described in terms of the two ‘disconnected’ functions:

[TABLE]

and:

[TABLE]

where:

[TABLE]

The remarkable length of this function is a consequence of an specific step of its derivation process, namely the propagation of a non-trivial color state via a $6\!\times\!6$ matrix that represents the correlator of four Wilson lines in the adjoint representation in color space (see Albacete:2018bbv for details). Having defined these functions, we can write our building block compactly as:

[TABLE]

Substituting this result in Eq.​ (32) and performing the ensuing index contractions (for which we use the Mathematica package FeynCalc Mertig:1990an ; Shtabovenko:2016sxi ), we obtain the main result of this work:

[TABLE]

where the dependencies have been omitted for readability. Here we defined the saturation scale characterizing each nucleus as:

[TABLE]

A(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) and B(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) were introduced in Eq.​ (21) as coefficients of the unpolarized and linearly polarized parts of the gluon distribution, respectively. Although we provided their general expressions in terms of f(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) (Eq.​ (22) and Eq.​ (23)), in order to analyze our result quantitatively we need to adopt an specific model. In the particular case of the strict MV model, where f(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) is taken as a Dirac delta, we can compute these functions as:

[TABLE]

where $K_{\scaleto{0}{4.2pt}}$ is a modified Bessel function and $m$ is the infrared cut-off mentioned in the previous section, now employed to regularize the divergent Green’s function G(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}):

[TABLE]

Note that although we could have used an unrelated regulator mass, for the sake of simplicity we choose it to be the same one introduced before. In our calculation we will consider only the leading behavior in the $m\!\rightarrow\!0$ limit, which is:

[TABLE]

$B_{\scaleto{\text{MV}}{0.12cm}}$, being a constant, yields a negligible correction to this logarithm. In the same limit, the leading behavior of \Gamma(r_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}) corresponds to the following expression:

[TABLE]

Except for $B_{\scaleto{\text{MV}}{0.12cm}}$, all these factors exhibit logarithmic divergences. However, as all logarithms stemming from $A$ and $\Gamma$ are cancelled in the prefactor of Eq.​ (44), the only divergences that we need to deal with are the ones included in the saturation scale $Q_{s}^{2}$ (Eq.​ (45)), which diverges in both infrared $m\!\rightarrow\!0$ and ultraviolet $r\!\rightarrow\!0$ limits. Different prescriptions with a varying level of sophistication are available in the literature to treat this issue. In order to give a general idea of the magnitude and analytical features of our solution, we will adopt the GBW model, which in practice consists simply in neglecting all logarithmic dependencies. In this framework, on Fig.​ 1 we draw the ratio of Eq.​ (44) to the square of the energy density average:

[TABLE]

(whose computation can be found in Albacete:2018bbv ) as a function of the dimensionless product $rQ_{s}$ for $Q_{s\scaleto{1}{4pt}}\!=\!Q_{s\scaleto{2}{4pt}}$. Note that we are also taking h(b_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\!=\!1 (strict MV model).

Although Eq.​ (44) is somewhat lengthy, the following simplifying limits provide us with remarkably more compact formulas. In the limit of small transverse separations $r\!\to\!0$ the two-point function tends to:

[TABLE]

The ratio with the product of the energy density averages at each transverse position reads:

[TABLE]

which displays the characteristic $1/(N_{c}^{2}-1)$ suppression factor of non-trivial color correlations. In the opposite limit, $rQ_{s}\gg 1$, we obtain:

[TABLE]

The previous expressions, being more ‘user-friendly’ than our complete result, greatly simplify the potential application of this work to further analytical or phenomenological calculations. Also, they allow for a straightforward comparison of our approach to other analytical frameworks available in the literature, which is the main subject of the next subsection.

4.1 The Glasma Graph approximation

The correlators presented in this paper, along with the energy density two-point function, were previously calculated in Lappi:2017skr under the so-called Glasma Graph approximation. In this approach it is assumed that the four-point correlation functions of the gluon fields factorize into products of two-point correlation functions such that:

[TABLE]

This Wick theorem-like decomposition is equivalent to assuming that the gluon fields conserve the Gaussian character of the color source distributions. This is not generally correct, as the dynamical generation of the former by the latter (encoded in the Yang-Mills equations) is non-linear. However, as observed in a previous work Albacete:2018bbv , this assumption yields a good approximation of the exact result in the limit of small transverse separations $r\!\to\!0$. In this limit an effective linearization of the fields’ dynamics takes place, as the connected function Eq.​ (34) vanishes and the disconnected contributions become equivalent to the two-point function of gluon fields. This results in a mapping of the Gaussian statistics followed by the color source distributions onto the gluon fields.

A comment is in order with respect to the designation of this approach. The original Glasma Graph approximation combines Gaussian statistics with the assumption that the valence quarks interact with the classical field by exchanging only two gluons, being applicable in the dilute limit Lappi:2015vta . This results in a factorization of double parton distributions into all possible products of single parton distributions, which yields great simplification in the context of the calculation of di-hadron correlators Lappi:2017skr . In the same spirit, the decomposition defined in Eq.​ (55) proposes a similar approach to the calculation of \langle\dot{\nu}(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\dot{\nu}(y_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\rangle, which is thus expressed in terms of the building block defined for \langle\dot{\nu}(x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})\rangle (namely the two-point correlator of gluon fields). Ignoring conceptual differences, in this paper we will give the name ‘Glasma Graph approximation’ to the approach based on said decomposition.

In Fig.​ 1 and Fig.​ 2 we compare our result (in the strict MV model and with $Q_{s\scaleto{1}{4pt}}\!=\!Q_{s\scaleto{2}{4pt}}$) with the one computed according to the Glasma Graph approximation. As can be seen in Fig.​ 1, although both results agree exactly in the small transverse separation limit $r\!\rightarrow\!0$, in the rest of the spectrum (approximately for $r\!>\!1/Q_{s}$) our computation yields a significantly harder curve. As shown in Fig.​ 2, this is also the case for the two-point function of the energy density, computed in the exact analytical approach in Albacete:2018bbv . Another major difference observed in said paper –and one of its main results– consists in a comparatively slow vanishing behavior in the infrared limit, where the covariance of the energy density decreases following a $\sim\!1/r^{2}$ asymptotic curve, in stark contrast with the much steeper $\sim\!1/r^{4}$ decreasing behavior displayed by the Glasma Graph result. Remarkably, in the calculation presented here this difference becomes larger (as can be seen in the right panel of Fig.​ 2). In the present case, while our result shows a $\sim\!1/r^{4}$ decreasing behavior (see Eq.​ (54)), the Glasma Graph approximation yields a much steeper $\sim\!1/r^{8}$ tail:

[TABLE]

This discrepancy potentially implies much different numerical results and physical interpretations for any observable built from this quantity, in particular those based in fluctuations of the global amount of axial charge per unit rapidity (as they are proportional to the two-dimensional transverse integral of Eq.​ (54)).

4.2 $N_{c}$-expansion

In order to complete the analysis of our final expression Eq.​ (44), in this subsection we will show its expansion in orders of $N_{c}$. The leading order term, of order $N_{c}^{0}$, reads:

[TABLE]

and the next term, of order $N_{c}^{-2}$, reads:

[TABLE]

As it is also the case for the covariance of the energy density $\epsilon_{0}$, the first two orders of the $N_{c}$-expansion of Eq.​ (44) yield a neat approximation of the complete result (see Fig.​ 3), but not a significant improvement regarding the practicality of the formulas.

5 Discussion and outlook

In this paper we performed a first-principles analytical calculation of the two-point correlator of the divergence of the Chern-Simons current. This object characterizes a source of fluctuations of axial charge density in the Glasma state produced in the initial stage of an ultra-relativistic HIC. With this calculation we expand on the results presented in a previous work Albacete:2018bbv , where we computed the covariance of the Glasma energy-momentum tensor. We performed both calculations following a classical approach based on the CGC effective theory, which we introduced by summarizing the computation of the gluon fields produced at $\tau\!=\!0^{+}$. Our framework features an explicit impact parameter dependence in the two-point correlator of color source densities (first introduced in Chen:2015wia ), as well as a generalization of the transverse profile of the interaction. These modifications were also incorporated in Albacete:2018bbv with the aim of expanding the potential phenomenological applications of our results. In the present work, however, we limited our analysis to the GBW prescription within the original MV model for simplicity.

With this setup we compare our result for the two-point correlator with the one obtained under the Glasma Graph approximation Lappi:2017skr . As was also the case for the energy-momentum tensor Albacete:2018bbv , the exact computation shows complete agreement with the Glasma Graph result in the ultraviolet $r\!\rightarrow\!0$ limit. However, a strong discrepancy emerges in the infrared $rQ_{s}\!\gg\!1$ limit: the exact two-point correlator of $\dot{\nu}$ (normalized to the product of the average energy densities at each transverse position) decays following a $\propto\!1/r^{4}$ power-law tail, whereas the Glasma Graph expression exhibits a much more rapidly decaying $\propto\!1/r^{8}$ behavior. Remarkably, the gap between both results is even larger than the one showed by the two-point correlator of the energy density, which in this limit disagrees with the Glasma Graph result by ‘only’ a $1/r^{2}$ factor. This suggests that the non-linear dynamics followed by the gluon fields have an even greater effect over the long-range transverse fluctuations of axial charge density than they do over those of the deposited energy. On the other hand, the results show that for both calculations the fact that the Glasma field correlators do not obey Gaussian dynamics can be overlooked in the ultraviolet limit, or to a good approximation for correlation distances shorter than $1/Q_{s}$. This outcome agrees with the expected validity range of the Glasma Graph approximation Lappi:2015vta .

One feature of our previous work Albacete:2018bbv that is not reproduced by the results of the present paper is the logarithmic enhancement exhibited by the correlation length. The computation of the two-dimensional transverse integral of Eq.​ (54) is dominated by the lower bound $r\!\sim\!1/Q_{s}$, as opposed to the case of the corresponding energy density correlator, which is sensitive to the infrared cut-off $r\!\sim\!1/m$ via a logarithmic factor $\ln(Q_{s}/m)$. This result thus seems somewhat more consistent with the conjectured Glasma flux tube picture, which predicts the range of the transverse color screening of correlations to be of size $1/Q_{s}$ Iancu:2002aq . Nevertheless, Eq.​ (54) still displays a remarkably slow fall-off that contrasts with the behavior one could naively expect from correlations between Gaussianly-distributed color charges.

The results of this work add further evidence on the importance of the non-linear dynamics relating color source densities and gauge field correlators beyond the validity range of the Glasma Graph approximation (thus supporting the conclusions reached in Albacete:2018bbv in this regard). In addition, the expressions presented in this paper provide analytical insight on the early-time local fluctuations of axial charge density in the transverse plane. By following the practical steps first presented in Lappi:2017skr , our formulas can be directly applied in phenomenological studies of anomalous transport phenomena. From Eq.​ (2) (rewritten using the covariant coordinate system, (\tau,\eta,x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}})), we obtain:

[TABLE]

Taking the first order of the $\tau$-expansion of $\dot{\nu}$ and integrating, we get to the differential distribution of axial charge density at early times:

[TABLE]

From this expression one can straightforwardly relate the two-point function of the divergence of the Chern-Simons current computed here to the correlation function \left\langle\frac{dN_{5}}{d^{2}x_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}d\eta}\frac{dN_{5}}{d^{2}y_{\mathbin{\mathchoice{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}{\vbox{\hbox{\scriptscriptstyle\perp}}}}}d\eta}\right\rangle. This object serves as the fundamental input for the Monte-Carlo modelization of initial conditions of axial charge density Lappi:2017skr , required by those hydrodynamical simulations that aim at describing anomalous transport phenomena. The results based on the expressions presented here would of course be subject to higher order corrections in $\tau$. The computation of said terms, as well as the calculation of observables relevant to QGP phenomenology, are left for future work.

Acknowledgements.

The author thanks Cyrille Marquet, Javier L. Albacete, Giuliano Giacalone and Tuomas Lappi for their useful feedback during the development of this work. Partial funding by a FP7-PEOPLE-2013-CIG Grant of the European Commission, reference QCDense/631558, the MINECO project FPA2016-78220 of the Spanish Government, and the ‘La Caixa’ Banking Foundation is gratefully acknowledged by the author.