Topological charge fluctuations in the Glasma
Pablo Guerrero-Rodríguez
 
A review of swarmalators and their potential in bio-inspired computing
    Christian Bettstetter
    December 30, 2023
======================================================================







§ 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 <cit.>. The anomalous term weights gauge field configurations according to their associated value of the topological invariant known as winding number (or topological charge):
Q_w=g^2/16π^2∫d^4x Tr{F_μν(x)F̃^μν(x)},
where F^μν and F̃^μν=1/2ϵ^μνρσF_ρσ correspond respectively to the field strength tensor and its dual. Q_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 <cit.>). 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:
dN_5/dt=∫d^3x ∂_μj^μ_5(x)=-g^2N_f/8π^2∫d^3x Tr{F_μν(x)F̃^μν(x)}=g^2N_f/2π^2∫d^3x ν̇(x),
where N_5 and j_5^μ correspond respectively to the axial charge density and current, N_f is the number of flavors and ν̇ 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 <cit.>. 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) <cit.>. Although these experiments have provided numerous measurements that are indeed compatible with said phenomena <cit.>, the presence of large background effects (e.g. transverse momentum <cit.> and local charge <cit.> conservation, which also give rise to intrinsic back-to-back correlations; and final state interactions <cit.>) 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 <cit.> and, most recently, the implementation of an isobaric collision program at RHIC aimed at the isolation of CME backgrounds <cit.>. 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 <cit.>. 

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 contribution[It has been argued that in the early stage of the collision sphaleron transitions are suppressed due to the boost invariance of the generated fields <cit.>. However, as the system evolves towards thermalization and boost invariance wears down, they would be significantly enhanced <cit.>. 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 <cit.>. 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 <cit.> 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. <cit.> 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) <cit.>, 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 Λ^+: particles with p^+>Λ^+ 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^+<Λ^+. 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 Λ^+. This framework has been extensively applied in the description of the early, non-equilibrium stage of HICs known as Glasma phase <cit.>.




In a previous work <cit.> 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 <cit.> 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:
ν̇(τ=0^+,x_)≡ν̇_04.2pt(x_)=Tr{E(τ=0^+,x_)B(τ=0^+,x_)},
where E^i=-F^0i and B^k=1/2ϵ^ijkF^ij are, respectively, the Glasma chromo-electric and -magnetic fields, evaluated at an infinitesimal positive proper time τ=0^+ at a point x_ 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 ((<ref>)). In the present work we evaluate the correlation function ⟨ν̇_04.2pt(x_)ν̇_04.2pt(y_)⟩, which characterizes the early event-by-event fluctuations of the ν̇ 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 <cit.> –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 <cit.>, 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:setup 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 τ=0^+, which acts as boundary condition for the ensuing evolution in the future light-cone. In sec:cs1 we calculate the expectation value of the divergence of the Chern-Simons current in the previously presented framework. In sec:cs2 we compute its two-point correlator, comparing our results with the aforementioned calculation performed under the Glasma Graph approximation <cit.>. 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 <cit.>, being the reader referred to said paper for detailed derivations. Finally, in sec:end we discuss the physical implications and potential phenomenological applications of our result, as well as its role in future studies.



§ 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^μ_1,25.8pt whose dynamics follow from the classical Yang-Mills equations:

[ D_μ,F^μν]=δ^ν+ρ^a(x^-,x_) t^a,
where ρ^a represents the density of valence (large-x) partons and t^a is the generator of SU(N_c) in the fundamental representation. The δ^ν+ 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 ρ^a as an stochastic quantity with a certain probability distribution W[ρ] associated as weight function. Thus, the observables are obtained as ensemble averages over the background classical fields:

⟨O[ρ] ⟩=1/N∫[dρ] W[ρ]O[ρ],
where N is a normalization constant equal to ∫[dρ] W[ρ]. 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[ρ] with a Gaussian distribution, yielding the following fundamental result:

⟨ρ^a(x^-,x_)ρ^b(y^-,y_)⟩_MV0.12cm = μ^2(x^-)δ^abδ(x^–y^-)δ^2(x_-y_).
Here, μ^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 (<ref>) by considering the following, more general, two-point function of the sources:

    ⟨ρ^a(x^-,x_)ρ^b(y^-,y_)⟩   = μ^2(x^-)h(b_)δ^abδ(x^–y^-)f(x_-y_) 
       ≡λ(x^-,b_)δ^abδ(x^–y^-)f(x_-y_),

where we allow the possibility of finite, non-homogeneous nuclei by explicitly introducing an impact parameter (b_≡(x_+y_)/2) dependence as previously done in <cit.>. Also, we drop the assumption that interactions are local in the transversal plane by introducing an undetermined function f(x_-y_) 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 τ=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:

    A^±_14.2pt,24.2pt(x^∓,x_)   =0 
    
    A^i_14.2pt,24.2pt(x^∓,x_)   =θ(x^∓)∫^∞_-∞dz^∓∫ dz^2_ G(z_-x_)U_14.2pt,24.2pt^†(z^∓,x_)∂^iρ̃_14pt,24.2pt(z^∓,z_)U_14.2pt,24.2pt(z^∓,x_)
       ≡θ(x^∓)α_14.2pt,24.2pt^i(x_),


by solving the Yang-Mills equations with one source in the Fock-Schwinger gauge. Here ρ̃ 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:

    U_14.2pt,24.2pt(x^∓,x_)=P^∓exp{ -ig∫^x^∓_-∞5ptdz^∓∫ dz^2_G(z_-x_)ρ̃_14pt,24pt(z^∓,z_)},

where G(z_-x_) is the Green's function for the two-dimensional Laplace operator. The solution at τ=0^+ is found by proposing the following ansatz:

    A^±(x^∓,x_)   =± x^±α(τ=0^+,x_)
    
    A^i (x^∓,x_)   =α^i(τ=0^+,x_),

where we adopted the comoving coordinate system, defined by proper time τ=√(2x^+x^-) and rapidity η=1/2ln (x^+/x^-). Then, we invoke a physical `matching condition' that requires Yang-Mills equations to be regular in the limit τ→0

. In doing so, the following relations are obtained:

    α^i (τ=0^+,x_)   = α_14.2pt^i(x_)+α_24.2pt^i(x_)
    α (τ=0^+,x_)   = ig/2[ α_14.2pt^i(x_) , α_24.2pt^i(x_) ],

which act as boundary conditions of the subsequent τ-evolution[Several approaches of both analytical and numerical nature have been applied for this computation in the literature. For instance, in <cit.> an analytical approximation based on an expansion of the Yang-Mills equations in powers of τ 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 τ=0^+ as:
E^z(τ=0^+,x_)=   -ig δ^ij[α^i_14.2pt(x_),α^j_24.2pt(x_)]

B^z(τ=0^+,x_)=   -ig ϵ^ij[α^i_14.2pt(x_),α^j_24.2pt(x_)],
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 <cit.>.



§ 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 <cit.>, 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 ν̇_04.2pt in terms of the gluon fields:
ν̇_04.2pt(x_)   =-g^2δ^ijϵ^klTr{[α^i_14.2pt x,α^j_24.2pt x][α^k_14.2pt x,α^l_24.2pt x]}

   =-g^2δ^ijϵ^klα^i,a_14.2pt xα^j,b_24.2pt xα^k,c_14.2pt xα^l,d_24.2pt xTr{[t^a,t^b][t^c,t^d]}

   =g^2/2δ^ijϵ^klf^abnf^cdnα^i,a_14.2pt xα^j,b_24.2pt xα^k,c_14.2pt xα^l,d_24.2pt x,
where we adopted the shorthand notation α^i,a(x_)≡α^i,a_x. As an intermediate step we expand the color structure of the gluon fields as:
α^i(x_)   =∫^∞_-∞dz^-∫dz^2_ G(z_-x_)∂^iρ̃^a(z^-,z_)U^†(z^-,x_)t^aU(z^-,x_)

   =∫^∞_-∞dz^-∫dz^2_ G(z_-x_)∂^iρ̃^a(z^-,z_)U^ab(z^-,x_)t^b≡α^i,b(x_) t^b,
where we used the relation between Wilson lines in the fundamental and adjoint representations U^†t^aU=U^abt^b.
The correlator of (<ref>) factorizes like:
⟨ν̇_04.2pt(x_)⟩=g^2/2δ^ijϵ^klf^abnf^cdn⟨α^i,a_14.2pt(x_)α^k,c_14.2pt(x_)⟩⟨α^j,b_24.2pt(x_)α^l,d_24.2pt(x_)⟩,
as in the MV model the color source fluctuations of each nucleus are assumed to be independent. 
Thus, the building block of ⟨ν̇_04.2pt⟩ is the correlator of two gauge fields evaluated in the same transverse position, ⟨α^i,a(x_)α^k,c(x_)⟩. We calculate this object as a limit:
⟨α^i,a(x_)α^j,b(x_) ⟩   =δ^ablim_r→0∫^∞_-∞dz^-λ(z^-,b_)∂^i_x∂^j_yL(r_)C^(2)_adj(z^-;x_,y_),
where r=|r_|=|x_-y_|. Here we introduced the following function:
L(r_)≡∫ dz^2_du^2_G(z_-x_)G(u_-y_)f(z_-u_).
From its symmetries and dimension, the double derivative of L(r_) featured in (<ref>) can be parameterized as:
∂^i_x∂^j_yL(r_)=A(r_)δ^ij+B(r_)( δ^ij/2-r^ir^j/r^2).
This formula accounts for an explicit separation of the contributions of the unpolarized A(r_) and linearly polarized B(r_) parts of the gluon distribution. We can express these coefficients in terms of f(r_) explicitly as:
A(r_)   =1/2∫d^2q_/(2π)^2f̂(q_)e^iq_·r_/q^2

B(r_)   =-∫d^2q_/(2π)^2f̂(q_)e^iq rcosθ/q^2cos(2θ),
where f̂(q_) is the Fourier transform of f(r_).
As for the last factor in (<ref>), it corresponds to the dipole function in the adjoint representation, which stems from the following correlator:
⟨U^a^'a(x^-,x_)U^a^'b(x^-,y_) ⟩   =δ^abexp{ -g^2N_c/2Γ(r_)λ̅(x^-,b_)}

   ≡δ^abC^(2)_adj(x^-;x_,y_).
Here we introduced the notation λ̅(x^-,x_)=∫^x^-_-∞dz^-λ(z^-,x_), as well as the following function:
Γ(r_)=2(L(0_)-L(r_)).
Computing the indicated limit explicitly, (<ref>) yields:
⟨α^i,a(x_)α^j,b(x_) ⟩=-1/2δ^abδ^ijμ̅^2h(x_)∂^2L(0_)=-1/2δ^abδ^ijλ̅(x_)∂^2L(0_).
(For a more detailed calculation of this correlator we refer the reader to <cit.>). The factor λ̅(x_) corresponds to λ integrated from -∞ to ∞ in the longitudinal direction (in general, we will identify functions integrated in the longitudinal direction from -∞ to ∞ by simply omitting their longitudinal dependence). The factor ∂^2L(0_) is a model-dependent constant that emerges from the following limit:
lim_r→0∂^i_x∂^j_yL(r_)=δ^ij/2∫d^2q_/(2π)^2f̂(q_)1/q^2≡-1/2δ^ij∂^2L(0_).
Some assumptions had to be made about the functions h(b_), f(r_) introduced in (<ref>) in order to arrive at the expressions presented in this section. Specifically, we take h(b_) 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_) behaves in such a way that its Fourier transform f̂(k_) tends to unity in the infrared limit (see <cit.> for details).

When finally substituting our building block (<ref>) into (<ref>), we get:
⟨ν̇_04.2pt(x_)⟩=g^2/8(∂^2L(0_))^2λ̅_14.2pt(x_)λ̅_24.2pt(x_)δ^ijϵ^klf^abnf^cdnδ^acδ^ikδ^bdδ^jl=0,
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 ν̇_04.2pt, whose computation we outline in the following section.



§ TWO-POINT CORRELATOR OF THE DIVERGENCE OF THE CHERN-SIMONS CURRENT
 


In this section we describe the calculation of ⟨ν̇_04.2pt(x_)ν̇_04.2pt(y_)⟩, 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 ν̇_04.2pt in terms of the gluon fields:
ν̇_04.2pt(x_)ν̇_04.2pt(y_)=g^4/4δ^ijϵ^klδ^i'j'ϵ^k'l'f^abnf^cdnf^a'b'mf^c'd'mα^i,a_1xα^j,b_2xα^k,c_1xα^l,d_2xα^i',a'_1yα^j',b'_2yα^k',c'_1yα^l',d'_2y,
then, the correlator reads:
⟨ν̇_04.2pt(x_)ν̇_04.2pt(y_)⟩=g^4/4ϵ^klϵ^k'l'f^abnf^cdnf^a'b'mf^c'd'm⟨α^i,a_xα^k,c_xα^i',a'_yα^k',c'_y⟩_14.2pt⟨α^i,b_xα^l,d_xα^i',b'_yα^l',d'_y⟩_24.2pt.
Color algebra-wise, this expression presents the same level of complexity than the two-point correlator of the energy density, previously computed in <cit.>. Happily, it features a much simpler transverse index structure, which we can rewrite as:

ϵ^klϵ^k'l'=δ^kk'δ^ll'-δ^kl'δ^lk',
yielding:
⟨ν̇_04.2pt(x_)ν̇_04.2pt(y_)⟩=g^4/4f^abnf^cdnf^a'b'mf^c'd'm(   ⟨α^i,a_xα^k,c_xα^i',a'_yα^k,c'_y⟩_14.2pt⟨α^i,b_xα^l,d_xα^i',b'_yα^l,d'_y⟩_24.2pt

   -⟨α^i,a_xα^k,c_xα^i',a'_yα^l,c'_y⟩_14.2pt⟨α^i,b_xα^l,d_xα^i',b'_yα^k,d'_y⟩_24.2pt).
The building block for this computation is the correlator of four gluon fields in two different transverse positions:
⟨α^i,a(x_) α^k,c(x_) α^i^',a^'(y_)α^k^',c^'(y_) ⟩=∫^∞_-∞dz^-dw^-dz^-'dw^-'
⟨∂^iρ̃^e(z^-,x_)/∇^2U^ea(z^-,x_)

∂^kρ̃^f(w^-,x_)/∇^2U^fc(w^-,x_)∂^i^'ρ̃^e^'(z^-',y_)/∇^2U^e'a'(z^-',y_)∂^k^'ρ̃^f^'(w^-',y_)/∇^2U^f'c'(w^-',y_) ⟩,
where 1/∇^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 <cit.>. Here we will simply indicate the result, for which we need to make some definitions. We start with the `connected' function:
C^ij;kl_ab;cd(x_,    y_,x_,y_)=f^acef^bde∂^i_x∂^j_yL(x_-y_)∂^k_xΓ(x_-y_)∂^l_yΓ(y_-x_)

   ×( 4/Γ^3 g^4 N_c^3-(λ̅^2(b_) /2ΓN_c+4/Γ^3 g^4 N_c^3+2 λ̅(b_)/Γ^2 g^2 N_c^2)C_adj^(2)(x_,y_)),
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:
 ∫^∞_-∞dz^-dw^-dz^-'dw^-'   ⟨∂^iρ̃^e(z^-,x_)/∇^2∂^kρ̃^f(w^-,x_)/∇^2∂^i^'ρ̃^e^'(z^-',y_)/∇^2∂^k^'ρ̃^f^'(w^-',y_)/∇^2⟩

   ×⟨U^ea(z^-,x_)U^fc(w^-,x_)U^e'a'(z^-',y_)U^f'c'(w^-',y_) ⟩.
 These contributions are described in terms of the two `disconnected' functions:
D^ik;i'k'_ac;a'c'(x_,x_,y_,y_)=1/4δ^ikδ^i'k'(∂^2L(0_))^2δ^acδ^a^'c^'λ̅^2(b_),
and:
D^ij;kl_ab;cd(x_,y_,x_,y_)=2(δ^abδ^cd[N_c^2-4/2N_c^2f_14pt+2/N_c^2f_24pt+N_c+2/4N_cf_34pt+N_c-2/4N_cf_44pt].

+δ^acδ^bd[1/N_c^2-1f_54pt-N_c+2/2N_c(N_c+1)f_34pt+N_c-2/2N_c(N_c-1)f_44pt]

+δ^adδ^bc[-N_c^2-4/2N_c^2f_14pt-2/N_c^2f_24pt+N_c+2/4N_cf_34pt+N_c-2/4N_cf_44pt]

+d^abmd^cdm[-1/N_cf_14pt+1/N_cf_24pt+1/4f_34pt-1/4f_44pt]+d^admd^cbm[ 1/N_cf_14pt-1/N_cf_24pt+1/4f_34pt-1/4f_44pt]

.+d^acmd^bdm[N_c/N_c^2-4f_24pt-N_c+4/4(N_c+2)f_34pt+N_c-4/4(N_c-2)f_44pt] )T^ ij;kl(x_,y_;x_,y_),
where:
f_14pt=   2/(N_cg^2Γ)^2(1-C^(2)_adj(x_,y_))^2

    f_24pt =   2/N_cg^2Γ( 2/N_cg^2Γ(1-C^(2)_adj(x_,y_))-λ̅(b_)C^(2)_adj(x_,y_))

    f_34pt =    ( 4/N_c(N_c+2)g^4Γ^2(1-C^(2)_adj(x_,y_)).

      .-2/(N_c+2)(N_c+1)g^4Γ^2(1-(C^(2)_adj(x_,y_))^2exp{-g^2Γλ̅(b_) }))

    f_44pt =    ( 4/N_c(N_c-2)g^4Γ^2(1-C^(2)_adj(x_,y_)).

       .-2/(N_c-2)(N_c-1)g^4Γ^2(1-(C^(2)_adj(x_,y_))^2exp{g^2Γλ̅(b_) }))

    f_54pt =    2/N_cg^2Γ( λ̅(b_)-2/N_cg^2Γ(1-C^(2)_adj(x_,y_))).
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×6 matrix that represents the correlator of four Wilson lines in the adjoint representation in color space (see <cit.> for details).
Having defined these functions, we can write our building block compactly as:
⟨α^i a(x_) α^k c(x_) α^i^'a^'(y_)α^k^'c^'(y_) ⟩= D^ik;i'k'_ac;a'c'(x_,x_,y_,y_)+D^ii';kk'_aa';cc'(x_,y_,x_,y_)

+D^ik';ki'_ac';ca'(x_,y_,x_,y_)+C^ii';kk'_aa';cc'(x_,y_,x_,y_)+C^ik';ki'_ac';ca'(x_,y_,x_,y_)

+C^kk';ii'_cc';aa'(x_,y_,x_,y_)+C^ki';ik'_ca';ac'(x_,y_,x_,y_).
Substituting this result in (<ref>) and performing the ensuing index contractions (for which we use the Mathematica package FeynCalc <cit.>), we obtain the main result of this work:


⟨   ν̇_04.2pt(x_)ν̇_04.2pt(y_)⟩=16A^4-B^4/g^4Γ^4N^2_c([ N_c^6+2 N_c^4-19 N_c^2+8/2(N_c^2-1)^2-2N_c^6-3 N_c^4-26 N_c^2+16/N_c^4-5 N_c^2+4e^-Q_s14pt^2r^2/4 

   +(N_c^2-1)(1-e^-Q_s14pt^2r^2/4(1+Q_s14pt^2r^2/4))(1-e^-Q_s24pt^2r^2/4(1+Q_s24pt^2r^2/4))

    +r^4/4Q_s14pt^2Q_s24pt^2-2r^2Q_s14pt^2(1-e^-Q_s24pt^2r^2/4)+ 2 (N_c^2-8) (N_c^2-1) (N_c^2+4)/(N_c^2-4)^2e^-(Q_s14pt^2+Q_s24pt^2)r^2/4

    +(N_c-1)(N_c+3)N^3/2(N_c+1)^2(N_c+2)^2(N_c/2e^-(N_c+1) r^2 Q_s24pt^2/2 N_c+(N_c+2)-2(N_c+1)e^-Q_s24pt^2r^2/4)e^-(N_c+1) r^2 Q_s14pt^2/2 N_c 

    +(N_c+1)(N_c-3)N^3/2(N_c-1)^2(N_c-2)^2(N_c/2e^-(N_c-1) r^2 Q_s24pt^2/2 N_c+(N_c-2)-2(N_c-1)e^-Q_s24pt^2r^2/4)e^-(N_c-1) r^2 Q_s14pt^2/2 N_c] 

    + [ 1↔2]),
where the dependencies have been omitted for readability. Here we defined the saturation scale characterizing each nucleus as:

    r^2Q_s^2/4=g^2N_c/2Γ(r_)λ̅(b_).

A(r_) and B(r_) were introduced in (<ref>) 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_) ((<ref>) and (<ref>)), 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_) is taken as a Dirac delta, we can compute these functions as:
A(r_)_MV0.12cm   =-1/2G(r_)=1/4πK_04.2pt(mr)

B(r_)_MV0.12cm   =1/4π,
where K_04.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_):
G(r_)=-∫d^2q_/(2π)^2e^iq_·r_/q^2+m^2=-1/2πK_04.2pt(m r).
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→0 limit, which is:
A(r_)_MV0.12cm≈1/8πln(4/m^2r^2).
B_MV0.12cm, being a constant, yields a negligible correction to this logarithm. In the same limit, the leading behavior of Γ(r_) corresponds to the following expression:
Γ(r_)_MV0.12cm=1/2πm^2-r/2πmK_14.2pt(mr)≈r^2/8πln(4/m^2r^2).

Except for B_MV0.12cm, all these factors exhibit logarithmic divergences. However, as all logarithms stemming from A and Γ are cancelled in the prefactor of (<ref>), the only divergences that we need to deal with are the ones included in the saturation scale Q_s^2 ((<ref>)), which diverges in both infrared m→0 and ultraviolet r→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 Result1 we draw the ratio of (<ref>) to the square of the energy density average:
⟨ϵ_03.8pt(x_)⟩_MV0.12cm=C_F0.4em/g^2Q^2_s14ptQ^2_s24pt ,
(whose computation can be found in <cit.>) as a function of the dimensionless product rQ_s for Q_s14pt=Q_s24pt. Note that we are also taking h(b_)=1 (strict MV model).



Although (<ref>) is somewhat lengthy, the following simplifying limits provide us with remarkably more compact formulas. In the limit of small transverse separations r→0 the two-point function tends to:

    lim_r→0 ⟨ν̇_04.2pt(x_)ν̇_04.2pt(y_)⟩_ MV0.12cm=3(N^2_c-1)/32g^4N_c^2Q_s14pt^4Q_s24pt^4.

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

    lim_r→ 0(⟨ν̇_04.2pt(x_)ν̇_04.2pt(y_)⟩/⟨ϵ_03.8pt(x_)⟩⟨ϵ_03.8pt(y_)⟩)_ MV0.12cm   =3/8(N^2_c-1),

which displays the characteristic 1/(N_c^2-1) suppression factor of non-trivial color correlations.
In the opposite limit, rQ_s≫1, we obtain:





    lim_r Q_s≫1(⟨ν̇_04.2pt(x_)ν̇_04.2pt(y_)⟩/⟨ϵ_03.8pt(x_)⟩⟨ϵ_03.8pt(y_)⟩)_ MV0.12cm   =32/(N^2_c-1)^2r^4Q_s14pt^2Q_s24pt^2.

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.



 §.§ The Glasma Graph approximation
 

The correlators presented in this paper, along with the energy density two-point function, were previously calculated in <cit.> 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:

    ⟨α^i,a(x_) α^k,c(x_) α^i^',a^'(y_)α^k^',c^'(y_) ⟩_GG0.12cm=    ⟨α^i,a(x_) α^k,c(x_)⟩⟨α^i^',a^'(y_)α^k^',c^'(y_) ⟩
    
    +    ⟨α^i,a(x_)α^i^',a^'(y_)⟩⟨α^k,c(x_)α^k^',c^'(y_)⟩
    
    +    ⟨α^i,a(x_)α^k^',c^'(y_)⟩⟨α^k,c(x_)α^i^',a^'(y_)⟩.

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 <cit.>, this assumption yields a good approximation of the exact result in the limit of small transverse separations r→0. In this limit an effective linearization of the fields' dynamics takes place, as the connected function (<ref>) 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 <cit.>. 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 <cit.>. In the same spirit, the decomposition defined in (<ref>) proposes a similar approach to the calculation of ⟨ν̇(x_)ν̇(y_)⟩, which is thus expressed in terms of the building block defined for ⟨ν̇(x_)⟩ (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 Result1 and Result2 we compare our result (in the strict MV model and with Q_s14pt=Q_s24pt) with the one computed according to the Glasma Graph approximation. As can be seen in Result1, although both results agree exactly in the small transverse separation limit r→0, in the rest of the spectrum (approximately for r>1/Q_s) our computation yields a significantly harder curve. As shown in Result2, this is also the case for the two-point function of the energy density, computed in the exact analytical approach in <cit.>. 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 ∼1/r^2 asymptotic curve, in stark contrast with the much steeper ∼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 Result2). In the present case, while our result shows a ∼1/r^4 decreasing behavior (see (<ref>)), the Glasma Graph approximation yields a much steeper ∼1/r^8 tail:

    lim_rQ_s≫1(⟨ν̇_03.8pt(x_)ν̇_03.8pt(y_)⟩/⟨ϵ_03.8pt(x_)⟩⟨ϵ_03.8pt(y_)⟩)_GG0.3em=96/(N_c^2-1)r^8Q_s14pt^4Q_s24pt^4.

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 (<ref>)).



 §.§ N_c-expansion
 

In order to complete the analysis of our final expression (<ref>), in this subsection we will show its expansion in orders of N_c. The leading order term, of order N_c^0, reads:

    [⟨ν̇_03.8pt(x_)ν̇_03.8pt(y_)⟩]_N_c^0=[1/g^4 r^8e^-r^2/2(Q_s14pt^2+Q_s24pt^2)(8+16 e^Q_s14pt^2 r^2/2-32e^Q_s14pt^2 r^2/4+ 24 e^r^2/2(Q_s14pt^2+Q_s24pt^2)..
    ..-8 e^r^2/4(2 Q_s14pt^2+Q_s24pt^2)( 8+Q_s24pt^2r^2)+e^r^2/4(Q_s14pt^2+Q_s24pt^2)(Q_s14pt^2 Q_s24pt^2 r^4+4r^2(Q_s14pt^2+Q_s24pt^2)+48))]+[1↔2],

and the next term, of order N_c^-2, reads:

    [⟨ν̇_03.8pt(x_)ν̇_03.8pt(y_)⟩]_N_c^-2=[1/N_c^2g^4 r^8e^-r^2/2(Q_s14pt^2+Q_s24pt^2)(2 Q_s14pt^2r^2(8+Q_s14pt^2r^2)e^Q_s24pt^2 r^2/2.
    
    +8  e^r^2/4(2Q_s14pt^2+Q_s24pt^2)( 4Q_s14pt^2 r^2+Q_s24pt^2 r^2-4)+4e^r^2/2(Q_s14pt^2+Q_s24pt^2)(Q_s14pt^2 Q_s24pt^2 r^4-4r^2(Q_s14pt^2+Q_s24pt^2)+4)
    
    -4(Q_s14pt^4 r^4+12Q_s14pt^2 r^2+32)e^Q_s24pt^2 r^2/4-e^r^2/4(Q_s14pt^2+Q_s24pt^2)(Q_s14pt^2 Q_s24pt^2 r^4+4r^2(Q_s14pt^2+Q_s24pt^2)-80)
    
    +(Q_s14pt^2r^2+Q_s24pt^2r^2+8)^2)]+[1↔ 2].

As it is also the case for the covariance of the energy density ϵ_0, the first two orders of the N_c-expansion of (<ref>) yield a neat approximation of the complete result (see Result3), but not a significant improvement regarding the practicality of the formulas.



§ 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 <cit.>, 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 τ=0^+.  Our framework features an explicit impact parameter dependence in the two-point correlator of color source densities (first introduced in <cit.>), as well as a generalization of the transverse profile of the interaction. These modifications were also incorporated in <cit.> 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 <cit.>. As was also the case for the energy-momentum tensor <cit.>, the exact computation shows complete agreement with the Glasma Graph result in the ultraviolet r→0 limit. However, a strong discrepancy emerges in the infrared rQ_s≫1 limit: the exact two-point correlator of ν̇ (normalized to the product of the average energy densities at each transverse position) decays following a ∝1/r^4 power-law tail, whereas the Glasma Graph expression exhibits a much more rapidly decaying ∝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 <cit.>.

One feature of our previous work <cit.> 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 (<ref>) is dominated by the lower bound r∼1/Q_s, as opposed to the case of the corresponding energy density correlator, which is sensitive to the infrared cut-off r∼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 <cit.>. Nevertheless, (<ref>) 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 <cit.> 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 <cit.>, our formulas can be directly applied in phenomenological studies of anomalous transport phenomena. From (<ref>) (rewritten using the covariant coordinate system, (τ,η,x_)), we obtain:
dN_5/d^2x_dη=∫dτ τ∂_μj^μ_5=g^2N_f/2π^2∫dτ τν̇(x).





Taking the first order of the τ-expansion of ν̇ and integrating, we get to the differential distribution of axial charge density at early times:
dN_5/d^2x_dη|_τ=0^+=τ^2/2g^2N_f/2π^2ν̇_03.8pt(x_),
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 ⟨dN_5/d^2x_dηdN_5/d^2y_dη⟩. This object serves as the fundamental input for the Monte-Carlo modelization of initial conditions of axial charge density <cit.>, 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 τ. The computation of said terms, as well as the calculation of observables relevant to QGP phenomenology, are left for future work.








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.


JHEP-2modM

10

Fukushima:2008xe
K. Fukushima, D. E. Kharzeev and H. J. Warringa,  
  http://dx.doi.org/10.1103/PhysRevD.78.074033Phys. Rev. D78 (2008) 074033 [http://arXiv.org/abs/0808.3382arXiv:0808.3382 [hep-ph]]



Lenz:2001me
F. Lenz,    http://dx.doi.org/10.1007/978-3-540-31532-2_2Lect. Notes Phys. 659 (2005) 7
  [http://arXiv.org/abs/hep-th/0403286arXiv:hep-th/0403286
  [hep-th]]. [,7(2001)]



Kharzeev:2007jp
D. E. Kharzeev, L. D. McLerran and H. J. Warringa,  
  http://dx.doi.org/10.1016/j.nuclphysa.2008.02.298Nucl. Phys.
  A803 (2008) 227 [http://arXiv.org/abs/0711.0950arXiv:0711.0950 [hep-ph]]



Kharzeev:2015znc
D. E. Kharzeev, J. Liao, S. A. Voloshin and G. Wang,  
  http://dx.doi.org/10.1016/j.ppnp.2016.01.001Prog. Part. Nucl.
  Phys. 88 (2016) 1 [http://arXiv.org/abs/1511.04050arXiv:1511.04050 [hep-ph]]





Abelev:2009ac
STAR collaboration, B. I. Abelev et. al.,  
  http://dx.doi.org/10.1103/PhysRevLett.103.251601Phys. Rev.
  Lett. 103 (2009) 251601 [http://arXiv.org/abs/0909.1739arXiv:0909.1739 [nucl-ex]]



Abelev:2012pa
ALICE collaboration, B. Abelev et. al.,  
  http://dx.doi.org/10.1103/PhysRevLett.110.012301Phys. Rev.
  Lett. 110 (2013)no. 1 012301
  [http://arXiv.org/abs/1207.0900arXiv:1207.0900
  [nucl-ex]]



Adamczyk:2014mzf
STAR collaboration, L. Adamczyk et. al.,  
  http://dx.doi.org/10.1103/PhysRevLett.113.052302Phys. Rev.
  Lett. 113 (2014) 052302 [http://arXiv.org/abs/1404.1433arXiv:1404.1433 [nucl-ex]]



Adamczyk:2015eqo
STAR collaboration, L. Adamczyk et. al.,  
  http://dx.doi.org/10.1103/PhysRevLett.114.252302Phys. Rev.
  Lett. 114 (2015)no. 25 252302
  [http://arXiv.org/abs/1504.02175arXiv:1504.02175
  [nucl-ex]]



Tribedy:2017hwn
STAR collaboration, P. Tribedy,  
  http://dx.doi.org/10.1016/j.nuclphysa.2017.05.078Nucl. Phys.
  A967 (2017) 740 [http://arXiv.org/abs/1704.03845arXiv:1704.03845 [nucl-ex]]



Bzdak:2010fd
A. Bzdak, V. Koch and J. Liao,  
  http://dx.doi.org/10.1103/PhysRevC.83.014905Phys. Rev. C83 (2011) 014905 [http://arXiv.org/abs/1008.4919arXiv:1008.4919 [nucl-th]]



Schlichting:2010qia
S. Schlichting and S. Pratt,  
  http://dx.doi.org/10.1103/PhysRevC.83.014913Phys. Rev. C83 (2011) 014913 [http://arXiv.org/abs/1009.4283arXiv:1009.4283 [nucl-th]]



Ma:2011uma
G.-L. Ma and B. Zhang,  
  http://dx.doi.org/10.1016/j.physletb.2011.04.057Phys. Lett.
  B700 (2011) 39 [http://arXiv.org/abs/1101.1701arXiv:1101.1701 [nucl-th]]



Voloshin:2004vk
S. A. Voloshin,  
  http://dx.doi.org/10.1103/PhysRevC.70.057901Phys. Rev. C70 (2004) 057901 [http://arXiv.org/abs/hep-ph/0406311arXiv:hep-ph/0406311 [hep-ph]]



Bzdak:2012ia
A. Bzdak, V. Koch and J. Liao,  
  http://dx.doi.org/10.1007/978-3-642-37305-3_19Lect. Notes Phys.
  871 (2013) 503 [http://arXiv.org/abs/1207.7327arXiv:1207.7327 [nucl-th]]



Wen:2016zic
F. Wen, J. Bryon, L. Wen and G. Wang,  
  http://dx.doi.org/10.1088/1674-1137/42/1/014001Chin. Phys. C42 (2018)no. 1 014001 [http://arXiv.org/abs/1608.03205arXiv:1608.03205 [nucl-th]]



Xu:2017qfs
H.-j. Xu, J. Zhao, X. Wang, H. Li, Z.-W. Lin, C. Shen and F. Wang,  
  http://dx.doi.org/10.1088/1674-1137/42/8/084103Chin. Phys. C42 (2018)no. 8 084103 [http://arXiv.org/abs/1710.07265arXiv:1710.07265 [nucl-th]]



Huang:2018duz
H. Huang et. al. in Proceedings, 9th International Particle
  Accelerator Conference (IPAC 2018): Vancouver, BC Canada, p. TUPAF006,
  2018



Luzum:2013yya
M. Luzum and H. Petersen,  
  http://dx.doi.org/10.1088/0954-3899/41/6/063102J. Phys. G41 (2014) 063102 [http://arXiv.org/abs/1312.5503arXiv:1312.5503 [nucl-th]]



Kharzeev:2001ev
D. Kharzeev, A. Krasnitz and R. Venugopalan,  
  http://dx.doi.org/10.1016/S0370-2693(02)02630-8Phys. Lett. B545 (2002) 298 [http://arXiv.org/abs/hep-ph/0109253arXiv:hep-ph/0109253 [hep-ph]]



Mace:2016svc
M. Mace, S. Schlichting and R. Venugopalan,  
  http://dx.doi.org/10.1103/PhysRevD.93.074036Phys. Rev. D93 (2016)no. 7 074036 [http://arXiv.org/abs/1601.07342arXiv:1601.07342 [hep-ph]]



Retinskaya:2013gca
E. Retinskaya, M. Luzum and J.-Y. Ollitrault,  
  http://dx.doi.org/10.1103/PhysRevC.89.014902Phys. Rev. C89 (2014)no. 1 014902 [http://arXiv.org/abs/1311.5339arXiv:1311.5339 [nucl-th]]



Newman:2005hd
G. M. Newman,  
  http://dx.doi.org/10.1088/1126-6708/2006/01/158JHEP 01
  (2006) 158 [http://arXiv.org/abs/hep-ph/0511236arXiv:hep-ph/0511236 [hep-ph]]



Son:2009tf
D. T. Son and P. Surowka,  
  http://dx.doi.org/10.1103/PhysRevLett.103.191601Phys. Rev.
  Lett. 103 (2009) 191601 [http://arXiv.org/abs/0906.5044arXiv:0906.5044 [hep-th]]



Sadofyev:2010pr
A. V. Sadofyev and M. V. Isachenkov,  
  http://dx.doi.org/10.1016/j.physletb.2011.02.041Phys. Lett.
  B697 (2011) 404 [http://arXiv.org/abs/1010.1550arXiv:1010.1550 [hep-th]]



Weigert:2005us
H. Weigert,    http://dx.doi.org/10.1016/j.ppnp.2005.01.029Prog. Part. Nucl. Phys. 55 (2005) 461
  [http://arXiv.org/abs/hep-ph/0501087arXiv:hep-ph/0501087
  [hep-ph]]



Gelis:2010nm
F. Gelis, E. Iancu, J. Jalilian-Marian and R. Venugopalan,  
  http://dx.doi.org/10.1146/annurev.nucl.010909.083629Ann. Rev.
  Nucl. Part. Sci. 60 (2010) 463
  [http://arXiv.org/abs/1002.0333arXiv:1002.0333 [hep-ph]]



McLerran:1993ni
L. D. McLerran and R. Venugopalan,  
  http://dx.doi.org/10.1103/PhysRevD.49.2233Phys. Rev. D49
  (1994) 2233 [http://arXiv.org/abs/hep-ph/9309289arXiv:hep-ph/9309289 [hep-ph]]



McLerran:1993ka
L. D. McLerran and R. Venugopalan,  
  http://dx.doi.org/10.1103/PhysRevD.49.3352Phys. Rev. D49
  (1994) 3352 [http://arXiv.org/abs/hep-ph/9311205arXiv:hep-ph/9311205 [hep-ph]]



McLerran:1994vd
L. D. McLerran and R. Venugopalan,  
  http://dx.doi.org/10.1103/PhysRevD.50.2225Phys. Rev. D50
  (1994) 2225 [http://arXiv.org/abs/hep-ph/9402335arXiv:hep-ph/9402335 [hep-ph]]



Lappi:2006fp
T. Lappi and L. McLerran,  
  http://dx.doi.org/10.1016/j.nuclphysa.2006.04.001Nucl. Phys.
  A772 (2006) 200 [http://arXiv.org/abs/hep-ph/0602189arXiv:hep-ph/0602189 [hep-ph]]



Lappi:2006hq
T. Lappi,  
  http://dx.doi.org/10.1016/j.physletb.2006.10.017Phys. Lett.
  B643 (2006) 11 [http://arXiv.org/abs/hep-ph/0606207arXiv:hep-ph/0606207 [hep-ph]]



Dumitru:2008wn
A. Dumitru, F. Gelis, L. McLerran and R. Venugopalan,  
  http://dx.doi.org/10.1016/j.nuclphysa.2008.06.012Nucl. Phys.
  A810 (2008) 91 [http://arXiv.org/abs/0804.3858arXiv:0804.3858 [hep-ph]]



Fukushima:2011nq
K. Fukushima and F. Gelis,  
  http://dx.doi.org/10.1016/j.nuclphysa.2011.11.003Nucl. Phys.
  A874 (2012) 108 [http://arXiv.org/abs/1106.1396arXiv:1106.1396 [hep-ph]]



Albacete:2018bbv
J. L. Albacete, P. Guerrero-Rodríguez and C. Marquet,  
  http://dx.doi.org/10.1007/JHEP01(2019)073JHEP 01 (2019)
  073 [http://arXiv.org/abs/1808.00795arXiv:1808.00795
  [hep-ph]]



Giacalone:2019kgg
G. Giacalone, P. Guerrero-Rodríguez, M. Luzum, C. Marquet and J.-Y.
  Ollitrault,    http://arXiv.org/abs/1902.07168arXiv:1902.07168 [nucl-th]



Lappi:2017skr
T. Lappi and S. Schlichting,  
  http://dx.doi.org/10.1103/PhysRevD.97.034034Phys. Rev. D97 (2018)no. 3 034034 [http://arXiv.org/abs/1708.08625arXiv:1708.08625 [hep-ph]]



Chen:2015wia
G. Chen, R. J. Fries, J. I. Kapusta and Y. Li,  
  http://dx.doi.org/10.1103/PhysRevC.92.064912Phys. Rev. C92 (2015)no. 6 064912 [http://arXiv.org/abs/1507.03524arXiv:1507.03524 [nucl-th]]



Fries:2006pv
R. J. Fries, J. I. Kapusta and Y. Li,  
  http://arXiv.org/abs/nucl-th/0604054arXiv:nucl-th/0604054
  [nucl-th]



Iancu:2002aq
E. Iancu, K. Itakura and L. McLerran,  
  http://dx.doi.org/10.1016/S0375-9474(03)01477-5Nucl. Phys. A724 (2003) 181 [http://arXiv.org/abs/hep-ph/0212123arXiv:hep-ph/0212123 [hep-ph]]



Mertig:1990an
R. Mertig, M. Bohm and A. Denner,  
  http://dx.doi.org/10.1016/0010-4655(91)90130-DComput. Phys.
  Commun. 64 (1991) 345



Shtabovenko:2016sxi
V. Shtabovenko, R. Mertig and F. Orellana,  
  http://dx.doi.org/10.1016/j.cpc.2016.06.008Comput. Phys.
  Commun. 207 (2016) 432 [http://arXiv.org/abs/1601.01167arXiv:1601.01167 [hep-ph]]



Lappi:2015vta
T. Lappi, B. Schenke, S. Schlichting and R. Venugopalan,  
  http://dx.doi.org/10.1007/JHEP01(2016)061JHEP 01 (2016)
  061 [http://arXiv.org/abs/1509.03499arXiv:1509.03499
  [hep-ph]]




