Evolution of initial stage fluctuations in the Glasma
EEvolution of initial stage fluctuations in the Glasma
Pablo Guerrero-Rodr´ıguez
1, 2 and Tuomas Lappi
1, 2 Department of Physics, University of Jyv¨askyl¨a P.O. Box 35, 40014 University of Jyv¨askyl¨a, Finland Helsinki Institute of Physics, P.O. Box 64, 00014 University of Helsinki, Finland
We perform a calculation of the one- and two-point correlation functions of energy density andaxial charge deposited in the Glasma in the initial stage of a heavy ion collision at finite propertime. We do this by describing the initial stage of heavy ion collisions in terms of freely-evolvingclassical fields whose dynamics obey the linearized Yang-Mills equations. Our approach allows usto systematically resum the contributions of high momentum modes that would make a powerseries expansion in proper time divergent. We evaluate the field correlators in the MV model usingthe Glasma Graph approximation, but our approach for the time dependence can be applied to ageneral four-point function of the initial color fields. Our results provide analytical insight into thepre-equilibrium phase of heavy ion collisions without requiring a numerical solution to the Yang-Millsequations.
I. INTRODUCTION
Heavy ion collisions (HICs) open an experimental win-dow to the most extreme regimes of Quantum Chromody-namics (QCD). More specifically, they give access to theQuark Gluon Plasma (QGP), an extremely hot and densephase of matter made up of deconfined partons. Theheavy ion physics programs at the Relativistic Heavy IonCollider (RHIC) and the Large Hadron Collider (LHC)are devoted to making precision measurements of QGPproperties, such as its collective fluid-like behavior. Thisfeature manifests itself experimentally through nontrivialcorrelations in the final state of the collision. However,such correlations are reflective not only of the propertiesof the QGP, but also of the fluctuations of energy and mo-mentum deposited during the earliest stages of HICs—before the QGP is formed (see e.g. [1]). The descriptionof these fluctuations requires a degree of phenomenolog-ical modeling to incorporate contributions from differentdynamical mechanisms (randomly changing nucleon po-sitions, fluctuations of partonic degrees of freedom, colorcharge density fluctuations) whose relative importance isnot yet very clear. A further source of uncertainty isunderstanding the time evolution of these fluctuations inthe very early pre-equilibrium stage, before hydrodynam-ics or even kinetic theory is applicable. Addressing thistime development is the primary purpose of this paper.Despite the inherent non-perturbative character of thedynamics in the early stages of a heavy ion collision, theColor Glass Condensate (CGC) effective theory (see e.g.[2–4] for reviews) provides a formalism to analyze it inthe weak coupling limit. In this framework, the high den-sity of soft (small- x ) partons carried by the colliding nu-clei is described in terms of gluon fields whose dynamicsobey the classical Yang-Mills equations. This approach isvalid at high energies, where the occupation numbers ofsmall- x gluons (more specifically, those with transversemomentum of the order of the saturation scale Q s ) arelarge. This justifies the use of a classical approximation,which makes it possible to both include nonlinear satu- ration effects, and address the time development of thesystem in an explicit calculation. The source of the gluonfields are the large- x partons, represented as a collectionof SU( N c ) charges. The fluctuations are implemented inthis context by describing the density of color chargesas a stochastic quantity that varies on an event-by-eventbasis. In the McLerran-Venugopalan (MV) model [5–7]it is assumed that these random variations obey Gaus-sian statistics. In this work we will adopt the generalizedGaussian approach, where correlators of the color chargesare expressed in terms of one two-point function, whichneed not, however, have precisely the form of the two-point function of the original MV model.In the CGC picture, the multiple interactions ensuingafter a HIC give rise to a coherent, highly dense substanceknown as Glasma [8]. Within a short time, this transientstate is believed to undergo an evolution process thatleads the system to local thermal equilibrium. This tran-sition is known as thermalization, and its precise theoreti-cal characterization remains one of the most fundamentalopen problems of the field [9–12]. The CGC formalismis well suited to describe the earliest phase of this evolu-tion, during which the gluon density is large enough towarrant a classical treatment. However, as the system ex-pands and becomes dilute, the classical field descriptionstarts to break down. A series of recent works [13–16]argues that an effective kinetic theory might provide theintermediate step that matches the classical descriptionwith hydrodynamics, which govern the subsequent evo-lution of the QGP phase. Although this discussion is offundamental interest in the field, it is out of the scope ofthis paper.Rather, in this work we focus on the evolution of theGlasma state at even earlier times and its impact on theprimordial fluctuations. Quite generally, the magnitudeof such fluctuations is encoded in the following differenceof correlators: S f ( x ⊥ , y ⊥ ) = (cid:104) f ( x ⊥ ) f ( y ⊥ ) (cid:105) − (cid:104) f ( x ⊥ ) (cid:105)(cid:104) f ( y ⊥ ) (cid:105) , (1)where f ( x ⊥ ) denotes the value of a property of the a r X i v : . [ h e p - ph ] F e b Glasma at a point x ⊥ of the transverse plane and thenotation (cid:104) ... (cid:105) represents an average over the backgroundfields. Such correlations have been computed analyti-cally in previous works for both the energy density [17–19] and the axial charge [8, 18, 20] deposited by the nucleiright after the collision (i.e. for an infinitesimal positiveproper time τ = 0 + ). In the case of the energy den-sity, this calculation has found a practical application inthe construction of a model of initial conditions for HICsthat was used in the description of eccentricity fluctu-ations in Refs. [20–23]. Similarly, it was proposed thatthe two-point function of the divergence of the Chern-Simons current might be applied in the modeling of ini-tial conditions of anomalous hydrodynamical simulations[18]. These and other potential applications are foundedon the assumption that, at least up until proper times τ ∼ /Q s , the thermalization process does not induce asignificant modification of the fluctuations of the rele-vant properties. While this assumption is supported bystudies in the case of eccentricities [24], it is uncertainwhether the same can be said about their fluctuations,whose computation implies the integration of the energydensity two-point function over the transverse plane. Itis thus essential to have a theoretically robust notion ofthe evolution of this object at later times.The evolution of the Glasma is encoded in the clas-sical Yang-Mills equations, for which no analytical so-lution has been found so far. A full nonperturbativesolution of the time dependence, such as in the IP-Glasma model [24–26], relies on a numerical lattice cal-culation [27, 28] to evolve the system from τ = 0 + (whereanalytical solutions can be obtained and used as bound-ary conditions to the evolution) up to a time when thesolution is matched to kinetic theory or hydrodynamics.As will be detailed in Sec.III, some analytical approachesare based on approximations such as the weak field limit,where one performs a systematic expansion in orders ofthe small color source densities [18, 29–32]. The downsideof the weak field limit is that it makes the results appli-cable only in the dilute-dense regime. Another strategy,first proposed in [33] and further developed in [34–36],is based on a systematic expansion of the gluon fields inpowers of τ . This expansion turns the Yang-Mills equa-tions into an infinite system of differential equations thatcan be solved recursively. However, this is a series aroundthe point τ = 0 where the solution in the MV model isnot analytical [17], and the resulting terms are found tobe plagued by UV-divergences.Here we will follow the approach also proposed inRef. [34], where it is argued that using the Abelian(linearized) version of the Yang-Mills effectively resumsthe singularities, making the series finite. Remarkably,this resummation ansatz achieves the same functional τ -dependence found in previous works without resorting tothe weak field limit. This would in principle come atthe cost of the validity range in τ , which according toRef. [34] would be reduced to very short proper times after the collision ( τ (cid:28) /Q s ). However, we argue thatthis limitation affects only the low momentum modes,which are less important for the discussion presented inthis work (we focus our study on largely UV-dominatedobjects—the energy density and the divergence of theChern-Simons current). Treating the dense-dense caseby means of the linearized Yang-Mills equations allowsus to propagate the full initial conditions of the colli-sion (i.e. at τ = 0 + ) instead of taking the first ordersof an expansion of the Wilson lines in powers of colorcharge densities. As we will discuss (see also [37]) thisapproach introduces some gauge dependence that is notthere for a systematical expansion in a small source den-sity. However, for a relatively UV dominated quantitysuch as the energy density both this gauge dependenceand the screening effects caused by nonlinearities in thefinal state [38, 39] should be small. In combining thenonlinear initial condition with a linearized time evolu-tion, we preserve the full gluon saturation features of thesystem at the initial condition.The paper is organized as follows. In Sec. II, we brieflydiscuss the basic elements of the Glasma fields at τ = 0 + ,as well as the specific approximations applied in thiswork. We also devote a part of this section to reviewthe calculation of the one- and two-point functions of theenergy density and the divergence of the Chern-Simonscurrent of the Glasma at τ = 0 + . In Sec. III we detailour calculation of the τ -evolution of these correlators. Inthe final section we present our conclusions, as well aspotential applications and continuations of this work. Ashort description of the divergence of the Chern-Simonscurrent and its role in the search for Chiral MagneticEffect signals is included in Appendix A, and a morecomplete result for the energy density correlator in Ap-pendix B. Lastly, in Appendix C we briefly discuss theresults obtained under the MV model with a fixed cou-pling constant. II. INITIAL CONDITIONS
In the following section we briefly revisit the calcula-tion of one- and two-point correlators at τ = 0 + . Theresults described here are well established in the liter-ature [8, 17–20]. However, reviewing them serves thetwofold purpose of showing a fuller picture of Glasmaevolution and fixing notations that will be used through-out the paper. Firstly, we outline the basic elements ofour framework. A. The MV model
This model, presented in [5–7], provides a descriptionof the parton content of a nucleus in the Infinite Mo- In Ref. [34] it is also proposed that this situation can be improvedby considering nonlinear corrections to the linearized Yang-Millsequations. mentum Frame (IMF), where one sees it moving in thepositive x direction with a very large light-cone momen-tum P + (cid:29) Λ QCD . The IMF naturally gives rise to twogroups of partons with vastly different dynamical fea-tures: the ‘hard’ modes, which carry a large momen-tum fraction p + = xP + and which could be thought ofas valence-like degrees of freedom; and the ‘soft’ modes,which would correspond to the small- x gluons . Due toLorentz time dilation, the latter perceive the former tobe ‘frozen’ in light-cone time x + . Also, the uncertaintyprinciple tells us that in the IMF the hard modes appearto be sharply localized around the light cone (within adistance ∆ x − ∼ /p + ). These kinematic considerationsmotivate the MV model to represent the valence quarksas static color currents: J µ = δ µ + ρ ( x ⊥ , x − ) , (2)where the color source densities ρ are assumed to be veryclose to a delta function in x − . As for the soft modes,at high energies they are characterized by large occupa-tion numbers (of the order of 1 /α s for modes with trans-verse momenta under the saturation scale Q s ), whichmakes them amenable to a description in terms of clas-sical fields. The dynamics of this system are encoded inthe classical Yang-Mills equations:[ D µ , F µν ] = J ν , (3)where the valence quarks enter as the source of the gluonfields. Solving these equations in the light cone gaugewe obtain the fields carried by the nucleus as pure gaugefields: α i ( x ⊥ ) = − ig U ( x ⊥ ) ∂ i U † ( x ⊥ ) . (4)Here U ( x ⊥ ) are the Wilson lines, SU( N c ) elements de-fined as path-ordered exponentials: U ( x ⊥ ) = P − exp (cid:26) − ig (cid:90) ∞−∞ dz − Φ( z − , x ⊥ ) (cid:27) , (5)which characterize the interaction of an external probewith the gluon field. The fields Φ satisfy: −∇ ⊥ Φ( x − , x ⊥ ) = ˜ ρ ( x − , x ⊥ ) , (6)where ˜ ρ is the color charge density in the covariant gauge.The field of Eq. (4) represents the small- x gluons carriedby one nucleus moving close to the speed of light in thepositive x direction. Similar expressions apply in thecase of a nucleus moving in the opposite direction (up toa x − → x + change). The separation between both groups is performed at an arbitrarymomentum Λ + . The evolution of the theory with Λ + is givenby the renormalization equation known as B-JIMWLK, whichcompletes the CGC framework [40–46]. In the case of a collision between two nuclei, first an-alyzed in Ref. [29], the current that enters in the Yang-Mills equations has the following form: J µ = δ µ + ρ ( x ⊥ , x − ) + δ µ − ρ ( x ⊥ , x + ) , (7)where the subscripts 1, 2 label the nuclei moving in thepositive and negative directions, respectively. In orderto solve these equations one has to separately considerdifferent regions in spacetime. In the forward light-cone( x ± >
0) one can parametrize the solutions in the form: A i ( τ, x ⊥ ) = α i ( τ, x ⊥ ) (8) A ± ( τ, x ⊥ ) = ± x ± α ( τ, x ⊥ ) . (9)This form enforces the Fock-Schwinger gauge condition x + A − + x − A + = 0, which acts as a sort of interpolationof the light-cone gauges of each nuclei. Combining Eqs.(8) and (9) with the fields carried by the individual nucleibefore the collision (i.e. x ± < A ± = ± θ ( x + ) θ ( x − ) x ± α ( τ, x ⊥ ) (10) A i = θ ( x − ) θ ( − x + ) α i ( x ⊥ )+ θ ( x + ) θ ( − x − ) α i ( x ⊥ ) + θ ( x + ) θ ( x − ) α i ( τ, x ⊥ ) . (11)No analytical solution has been found thus far for the sys-tem defined above at finite τ . It is possible, however, tofind exact expressions of the fields generated an infinites-imal proper time after the collision ( τ = 0 + ) in termsof α , . This is done by requiring the Yang-Mills equa-tions to be regular at τ = 0, which leads to the followingboundary conditions: α i ( τ = 0 + , x ⊥ ) = α i ( x ⊥ ) + α i ( x ⊥ ) (12) α ( τ = 0 + , x ⊥ ) = ig (cid:2) α i ( x ⊥ ) , α i ( x ⊥ ) (cid:3) . (13)These expressions, along with the Yang-Mills equations,constitute a system that yields a unique solution for τ -evolution once the color source densities ρ , are given.In the MV model the color charges are assumed to bestochastic variables that obey Gaussian statistics. Allthe information on the correlation functions of the colorcharges is contained in the two-point function, which inthe MV model is taken to be (cid:104) ρ a ( x ⊥ ) ρ b ( y ⊥ ) (cid:105) = g µ δ ab δ (2) ( x ⊥ − y ⊥ ) , (14)where the variance of the distribution, µ , is a parameterproportional to the color source number density and re-lated to the saturation scale Q s . This is the basic build-ing block for the calculation of observables in the MVmodel, which are computed as ensemble averages overthe color charge distributions.In this paper we will need higher point correlators ofgluon fields in order to calculate the energy density two-point function. We will calculate these in a Gaussianapproximation, by which we retain the feature that ev-ery higher point correlator can be expressed in terms of atwo-point function. However, we will use different func-tional forms of the relevant two-point function, gener-alizing the local correlator (14). The building block ofour calculations is the correlator of two gluon fields incoordinate space: (cid:104) α a,i ( x ⊥ ) α b,j ( y ⊥ ) (cid:105) = δ ab (cid:18) δ ij G ( x ⊥ , y ⊥ )+ (cid:18) δ ij − r i r j r (cid:19) h ( x ⊥ , y ⊥ ) (cid:19) , (15)where r ⊥ = x ⊥ − y ⊥ and r = | r ⊥ | . The functions G ( x ⊥ , y ⊥ ), h ( x ⊥ , y ⊥ ) are related to the unpolarizedˆ G ( b ⊥ , k ⊥ ) and linearly polarized ˆ h ( b ⊥ , k ⊥ ) gluon distri-butions by the following Hankel transforms [18]: G ( x ⊥ , y ⊥ ) = 12 π (cid:90) dkkJ ( kr ) ˆ G (cid:18) x ⊥ + y ⊥ , k (cid:19) (16) h ( x ⊥ , y ⊥ ) = 12 π (cid:90) dkkJ ( kr ) ˆ h (cid:18) x ⊥ + y ⊥ , k (cid:19) . (17)As in this paper we work mainly in coordinate space, wewill henceforth refer to G ( x ⊥ , y ⊥ ) and h ( x ⊥ , y ⊥ ) simplyas the unpolarized and linearly polarized gluon distri-butions, respectively. In Gaussian models, the polarizedand unpolarized gluon distributions are not independent,but can both be explicitly related to the gluon dipole am-plitude [18, 47], which contains all the information aboutthe strong interactions involved in Deep Inelastic Scat-tering (DIS) processes. Within the CGC framework onecan use different parameterizations of this object. Wewill in this work compare the GBW dipole distributionto a variant of the MV model that is regularized in aspecific way, as detailed below.The Golec-Biernat Wusthoff (GBW) model [48] wasoriginally introduced as a purely phenomenologicalparametrization to fit DIS data. The GBW dipole am-plitude reads: D GBW ( r ⊥ ) = exp (cid:26) − N c − N c Q s r (cid:27) (18)and the corresponding gluon distributions become: G GBW ( r ⊥ ) = Q s g N c − exp (cid:110) − Q s r (cid:111) Q s r / h GBW ( r ⊥ ) = 0 . (20)The MV model dipole amplitude, on the other hand, isobtained by starting from the color charge correlator (14)that is local in coordinate space. In order to obtain thedipole amplitude and the gluon field correlator, one needsto spread out the delta function-like color charges in the -6 -5 -4 -3 -2 -1
0 5 10 15 20
GBWMVMV (no running coupling)
0) the unpolarized gluon distribution G MV diverges logarithmically, while the polarized oneapproaches a constant. In momentum space this corre-sponds to the unintegrated gluon distribution in Eq. (16)behaving as ∼ /k at large momenta. In itself, this isone of the physically attractive features of the MV model,since it makes the integrated gluon distribution grow log-arithmically with the hard scale Q , corresponding to arudimentary approximation of DGLAP evolution. How-ever, it is problematic for the calculation of the τ = 0 + Glasma energy density, which is proportional to G ( r = 0).As discussed in detail in [17], the Glasma energy densityat finite τ is finite; in the numerical CYM calculationsthis UV divergence at precisely τ = 0 + is regularized bythe lattice spacing and would in fact be divergent in thecontinuum limit. Here we will regularize this divergenceby using a running coupling prescription. The choice wemake is to assume that the coupling g appearing as anoverall coefficient (but not the coupling in the exponent)in Eqs. (22) and (23) runs with the distance r accordingto: g ( r ) = g (¯ µ )ln (cid:16) e − γe − m r + e (cid:17) . (24)In the UV or small r limit, the factor on the second lineof Eqs. (22) and (23) approaches unity. The leading be-havior of the unpolarized gluon distribution G MV ( r ⊥ ), isthen dictated by the K ( mr ) function. Taking the cou-pling in the prefactor to be the running coupling (24) anddeveloping the Bessel functions the unpolarized distribu-tion becomes G MV ( r ⊥ ) −→ r → g ( r )¯ µ πN c ln (cid:18) e − γ e − m r (cid:19) = g (¯ µ )4 π ¯ µ N c ln (cid:16) e − γe − m r (cid:17) ln (cid:16) e − γe − m r + e (cid:17) −→ r → g (¯ µ )4 π ¯ µ N c . (25)This has a nice physical interpretation: the total gluondistribution as probed by a small probe r → µ times a coupling evaluated at the scale corresponding tothis number density g (¯ µ ).In the dilute limit the number density of color charges¯ µ is the variable that has a clear physical meaning. Inthe saturation regime, on the other hand, we should notbe discussing in terms of the number density, but in termsof the saturation scale Q s , which is defined as the inverseof the probe size r for which the nonlinear behavior of theexponential on the second line of Eq. (22) and (23) startsto matter. Thus we eliminate the explicit dependence on¯ µ in favor of the saturation scale, defined as: Q s = α s (¯ µ )¯ µ N c . (26)After this procedure our explicit expression for the gluon distributions is: G MV ( r ⊥ ) = Q s g (¯ µ ) N c ( mrK ( mr ) − K ( mr ))ln (cid:16) e − γe − m r + e (cid:17) (27) × − exp (cid:110) Q s m ( mrK ( mr ) − (cid:111) Q s m ( mrK ( mr ) − h MV ( r ⊥ ) = − Q s g (¯ µ ) N c mrK ( mr )ln (cid:16) e − γe − m r + e (cid:17) (28) × − exp (cid:110) Q s m ( mrK ( mr ) − (cid:111) Q s m ( mrK ( mr ) − . The mass regulator m is parametrically a confinementscale quantity, and in this work we will take its value as m = 0 . Q s . With such a prescription, the gluon distri-butions G MV and h MV tend to the corresponding GBWmodel expressions in the UV limit:lim r → G MV ( r ⊥ ) = Q s g N c (29)lim r → h MV ( r ⊥ ) = 0 . (30)A related procedure was previously applied in [18] inthe calculation of correlators at τ = 0 + by writing theMV model results for the gluon distributions in such away that there is an explicit factor of g in the exponent(using the fact that Q s ∼ g µ ). By taking this couplingto run one eliminates all the logarithms of mr at small r from the Bessel functions, and arrives at a form that inpractice is similar to GBW. In contrast, here we applythe running coupling prescription defined by Eqs. (24),(26) only to the prefactors on the first lines of Eqs. (22),(23) and not to the exponent.The difference between the running coupling ap-proaches in different sources in the literature amountsto some extent to the question of whether the couplingshould run such that the saturation scale Q s , the colorcharge density g ¯ µ , or the number density of colorcharged particles ¯ µ stays fixed. Ultimately different ap-proaches used in the literature, including ours, are some-what ad hoc. They should be judged on their effect onphysics, where our current approach is constructed topreserve a bit more the features of the fixed couplingMV model. The effect of this modification is better ob-served in the momentum space distributions, shown inFig. 1. Here we can see that our prescription indeed pre-serves a power law tail in the unpolarized gluon distri-bution. However, the power law is visibly steeper thanwith the fixed coupling MV model, which distinctly yieldslim k →∞ ˆ G ( k ⊥ ) ∼ /k . Such a steeper power is required forour calculation here, in order to obtain finite results incoordinate space.Another difference between our approach and the oneused in [18] lies in the exact form of Eq. (24), specificallyin the e − γ E − factor contained in the argument of thelogarithm in the denominator. With such a prescription,the subleading terms in the mr → (cid:104) α a α b α c α d (cid:105) .This is a highly complex object with a rich color struc-ture whose full calculation is discussed in depth in [19]for the simpler case of two and three transverse posi-tions. It could, in principle, be fully evaluated in theGaussian approximation. Doing so would, however, berather complicated. For simplicity’s sake, we shall thusin the present work assume that the averaging of gluonfields obeys Gaussian dynamics: (cid:104) α i,ax α j,by α k,cu α l,dv (cid:105) = (cid:104) α i,ax α j,by (cid:105)(cid:104) α k,cu α l,dv (cid:105) + (cid:104) α i,ax α k,cu (cid:105)(cid:104) α j,by α l,dv (cid:105) + (cid:104) α i,ax α l,dv (cid:105)(cid:104) α j,by α k,cu (cid:105) , (31)just like the color source densities do. Here we haveadopted the shorthand notation α i,ax ≡ α i,a ( x ⊥ ). TheWick theorem-like decomposition featured in Eq. (31) isknown in the literature as the Glasma Graph approxi-mation, and it has been used in many phenomenologi-cal studies of semihard two-particle correlations in thedense-dense collision regime [51–55]. It will also greatlysimplify the calculations presented in this paper. Thisassumption has been shown to yield exact results in theUV limit of Eq. (1), where the connected, nonlinear con-tributions computed in Ref. [19] vanish. In this regime,the dynamics linearize in such a way that a Gaussiandistribution for the color source densities is effectivelymapped onto another one for the gauge fields. Althoughthis approximation breaks down for correlation distanceslarger than 1 /Q s , it provides valuable analytical insightfor the quantities computed in the present work. Ourlinearized approach for the time evolution would, in it-self, be straightforwardly generalizable to a full nonlinearGaussian calculation of the gauge field four-point func-tion at τ = 0 + . B. Glasma correlators at τ = 0 + The structure that emerges right after a HIC in theCGC picture is characterized by the presence of purelylongitudinal chromo-electric and chromo-magnetic fields.In the Fock-Schwinger gauge used for the matching at τ = 0, these fields are given by: E η ( τ = 0 + , x ⊥ ) ≡ E η ( x ⊥ ) = − ig δ ij [ α i ( x ⊥ ) , α j ( x ⊥ )](32) B η ( τ = 0 + , x ⊥ ) ≡ B η ( x ⊥ ) = − ig (cid:15) ij [ α i ( x ⊥ ) , α j ( x ⊥ )] , (33)in terms of the gluon fields carried by each nucleus priorto the collision, α i , . As Eqs. (32), (33) are the onlynon-vanishing components of the field strength tensor at τ = 0 + , the gauge invariant energy density ε and thedivergence of the Chern-Simons current ˙ ν have the fol-lowing simple forms: ε ( τ = 0 + , x ⊥ ) ≡ ε ( x ⊥ ) = Tr { E η E η + B η B η } (34)˙ ν ( τ = 0 + , x ⊥ ) ≡ ˙ ν ( x ⊥ ) = Tr { E η B η } . (35)The expectation values of these quantities are well knownresults in the CGC literature [8, 17]. Writing the pre-vious expressions in terms of the gluon fields and thensubstituting the corresponding correlators (Eq. (15)) weobtain [17]: (cid:104) ε (cid:105) = g N c ( N c − G ( x ⊥ , x ⊥ ) G ( x ⊥ , x ⊥ ) (36) (cid:104) ˙ ν (cid:105) = 0 , (37)where G , correspond to the unpolarized gluon distri-butions of nuclei 1 and 2, respectively. In this resultwe can see that linearly polarized gluon distributions donot contribute to these expectation values. In the GBWmodel, Eq.(36) simply reads (cid:104) ε (cid:105) = C F Q s Q s /g , where C F = ( N c − / (2 N c ) is the Casimir of the fundamentalrepresentation. Note that in the MV model this objectcontains a logarithmic UV divergence, which we regular-ized by means of our running coupling prescription, but isvisible in the numerical CYM calculation as a divergencein the continuum limit of the energy density at τ = 0 + .Let us now focus on the two-point functions of ε and˙ ν . Writing them in terms of the gluon fields carried byeach nucleus, we obtain the following expressions: (cid:104) ε ( x ⊥ ) ε ( y ⊥ ) (cid:105) = g δ ij δ kl + (cid:15) ij (cid:15) kl )( δ i (cid:48) j (cid:48) δ k (cid:48) l (cid:48) + (cid:15) i (cid:48) j (cid:48) (cid:15) k (cid:48) l (cid:48) ) f abn f cdn f a (cid:48) b (cid:48) m f c (cid:48) d (cid:48) m (cid:104) α i,ax α k,cx α i (cid:48) ,a (cid:48) y α k (cid:48) ,c (cid:48) y (cid:105) (cid:104) α j,bx α l,dx α j (cid:48) ,b (cid:48) y α l (cid:48) ,d (cid:48) y (cid:105) (38) (cid:104) ˙ ν ( x ⊥ ) ˙ ν ( y ⊥ ) (cid:105) = g δ ij (cid:15) kl δ i (cid:48) j (cid:48) (cid:15) k (cid:48) l (cid:48) f abn f cdn f a (cid:48) b (cid:48) m f c (cid:48) d (cid:48) m (cid:104) α i,ax α k,cx α i (cid:48) ,a (cid:48) y α k (cid:48) ,c (cid:48) y (cid:105) (cid:104) α j,bx α l,dx α j (cid:48) ,b (cid:48) y α l (cid:48) ,d (cid:48) y (cid:105) . (39)These objects have an intricate color algebra composi- tion, featuring four SU( N c ) structure constants whose in-dices are contracted with two four-point correlators (oneper nucleus). Furthermore, the indices labeling the trans-verse plane vectors also contribute to the complexity ofthe formulas (and even more in the finite τ case, as willbe discussed in Sec. III). Evaluating such expressions ischallenging even when employing the simplest approachavailable (this is, adopting both the Glasma Graph ap-proximation and the GBW model), which is why we usedthe Mathematica package FeynCalc [56, 57] to performthe algebraic manipulations required throughout the cal-culation process. By applying the Glasma Graph ap-proximation and adopting the GBW model one is ableto obtain the following formulas: (cid:104) ε ( x ⊥ ) ε ( y ⊥ ) (cid:105)(cid:104) ε ( x ⊥ ) (cid:105)(cid:104) ε ( y ⊥ ) (cid:105) − N c − (cid:34) (cid:32) − e − Q s r / Q s r / (cid:33) + 23 (cid:32) − e − Q s r / Q s r / (cid:33) (cid:35) (40) (cid:104) ˙ ν ( x ⊥ ) ˙ ν ( y ⊥ ) (cid:105)(cid:104) ε ( x ⊥ ) (cid:105)(cid:104) ε ( y ⊥ ) (cid:105) = 38( N c − (cid:32) − e − Q s r / Q s r / (cid:33) (41)(first presented in this compact form in [18]).By performing the full calculation of the four-pointfunction of the fields (i.e. without assuming a Gaussian-like decomposition) one arrives at much lengthier expres-sions that, remarkably, feature different asymptotic be-haviors in the Q s r (cid:29) ∝ /r in the case of the en-ergy density and ∝ /r for the divergence of the Chern-Simons current) [19]. We expect that, since the consid-ered quantities are largely UV-dominated, the breakdownof the Glasma Graph approximation should not have abig impact on the qualitative conclusions of the presentwork. Conversely, the discrepancy at large correlationdistances will likely be much more relevant when consid-ering quantities that are sensitive to the infrared. An im-portant example of this is given by the integral of Eq.(40)over the transverse plane, which is a key quantity in thecalculation of eccentricity fluctuations. In fact, the in-frared sensitivity enters due to the 1 /r fall-off featuredin the full calculation. This asymptotic behavior givesrise to a correlation length that depends on the infraredscale 1 /m through a logarithmic factor ln ( Q s /m ). Thestudy of the τ -dependence of this property is left for fu-ture work. III. PROPER TIME EVOLUTION
Our goal in this section is to generalize the previouscorrelators to finite proper times, τ >
0. As the respec-tive supports of the color charge sources are Lorentz-contracted almost to delta functions ( ρ , ∼ δ ( x ± )), inthe space-time region τ > D µ , F µν ] = 0. The separate com-ponents of these equations in the comoving coordinate system ( τ, η, i ) read: igτ [ α, ∂ τ α ] − τ (cid:2) D i , ∂ τ α i (cid:3) = 0 (42)1 τ ∂ τ τ ∂ τ ( τ α ) − (cid:2) D i , (cid:2) D i , α (cid:3)(cid:3) = 0 (43)1 τ ∂ τ ( τ ∂ τ α i ) − igτ (cid:2) α, (cid:2) D j , α (cid:3)(cid:3) − (cid:2) D j , F ji (cid:3) = 0 . (44)We now want to solve these equations by linearizingthem, i.e. neglecting all the terms that are higher or-der in the gauge potentials α, α i . This linearization un-avoidably introduces a gauge dependence in the calcula-tion, since any gauge transformation changes the valueof the gauge potentials, and thus the value of the ne-glected higher order terms. However, there are severalphysical constraints that allow one to choose the propergauge. Firstly, we are solving the equations of motionin a spacetime region without any external charges thatwould introduce a static Coulomb field. It is thereforenatural to maintain the Fock-Schwinger gauge condition A τ = 0, which is achieved by refraining from τ -dependentgauge transformations. Secondly, our physical situationis boost-invariant, and it is natural to maintain this boostinvariance at the level of the gauge potentials. Thismeans that we do not allow our gauge transformationto depend on the spacetime rapidity η . These conditionsrestrict us to performing gauge transformations that onlydepend on the transverse coordinate. The general suchtransformation can be written as α ( τ, x ⊥ ) = U ( x ⊥ ) β ( τ, x ⊥ ) U † ( x ⊥ ) (45) α i ( τ, x ⊥ ) = U ( x ⊥ ) (cid:18) β i ( τ, x ⊥ ) − ig ∂ i (cid:19) U † ( x ⊥ ) . (46)The linearized equations of motion for the transformedgauge potential β read: ∂ τ ∂ i β i = 0 (47)1 τ ∂ τ τ ∂ τ ( τ β ) − ∂ i ∂ i β = 0 (48)1 τ ∂ τ ( τ ∂ τ β i ) − ( ∂ k ∂ k δ ij − ∂ i ∂ j ) β j = 0 . (49)A natural choice for fixing the last degree of freedomin the gauge transformation is to choose the transverseCoulomb gauge defined by the condition ∂ i β i = 0. Thischoice minimizes the amplitudes of the transverse com-ponents of the fields, and one could thus argue that it isthe choice where a linearized approximation works best.Since we refrain from time dependent gauge transforma-tions, we can in general only enforce the Coulomb condi-tion at one specific value of τ . With the general nonlinearequations of motion the condition ∂ i β i = 0 is not con-served in τ [58]. However, for the linearized equations theCoulomb condition is explicitly conserved in time, (47).Thus with linearized evolution we can work completelyin Coulomb gauge at all values of τ , and this is what weshall do in the following. In this case the equations ofmotion reduce to ∂ τ τ ∂ τ ( τ β ) = τ ∂ i ∂ i β (50) ∂ τ ( τ ∂ τ β i ) = τ ∂ k ∂ k β i , (51)which can be straightforwardly solved by Fourier-transforming. The general plane wave solutions for thetime dependence of the momentum modes are β ( τ, k ⊥ ) = β ( k ⊥ ) 2 J ( ωτ ) ωτβ i ( τ, k ⊥ ) = β i ( k ⊥ ) J ( ωτ ) , (52)with the free dispersion relation ω ( k ⊥ ) = | k ⊥ | .We must now relate the initial conditions of the gaugepotentials, β ( k ⊥ ) and β i ( k ⊥ ) to the gauge fields of theincoming nuclei. In principle, we would do this by takingthe initial conditions for the gauge fields from Eqs. (12)and (13), and performing the gauge transformation toCoulomb gauge. This procedure was done numericallyin Ref. [37]. In this study it was found that a lineartime evolution starting from such a Coulomb gauge initialcondition is a very good approximation of the full CYMevolution, apart from the very soft modes whose evolu-tion is genuinely nonlinear due to screening. Finding therequired gauge transformation (46) to Coulomb gauge isalso possible analytically in the dilute-dense case [18, 31],leading to nice k ⊥ -factorized expressions for the gluonspectrum in proton-nucleus collisions that are often usedin phenomenology [59, 60] (a calculation recently ex-tended one order higher in the color charge density inthe proton [61]). However, in the full dense-dense casethat we want to address here, finding an explicit form forthe gauge transformation U ( x ⊥ ) from Eq.(46) is not pos-sible analytically. A way out of this conundrum can befound by fully exploiting the fact that we have resolved touse a linearized approximation for the time evolution inthe region τ >
0, and that the physics we are interestedin really depend on the electric and magnetic fields, notthe gauge potentials themselves. Thus, in fact, we do notneed an explicit form of the gauge transformation U ( x ⊥ ).Instead, it is enough to find a gauge potential that:1. Satisfies the Coulomb gauge linearized equation ofmotion, i.e. is of the form (52), and2. Reproduces the initial field strength (32), (33) at τ = 0 + , when calculating the electric and magneticstrength from the gauge potential using the lin-earized approximation, consistently with the lineartime evolution.Thus, instead of finding an explicit form for the gaugetransformation, we obtain the gauge potential at τ > τ = 0 + . Alternatively, one couldformulate our approach in terms of linear equations of motion for the electric and magnetic fields themselves,as the Abelian Maxwell equations are usually written.This equation can then be matched to an initial condi-tion that includes the full nonlinear dependence in thecolor charges of the colliding projectiles. This procedureenables us to use include the full nonlinear structure ofgluon saturation in both colliding projectiles at the initialtime, but nevertheless obtain an analytical expression forthe time evolution.Taking this approach, it is now straightforward to fixthe initial conditions in the general solution (52) to ob-tain: β ( τ, k ⊥ ) = τk E η ( k ⊥ ) J ( kτ ) β i ( τ, k ⊥ ) = − i (cid:15) ij k j k B η ( k ⊥ ) J ( kτ ) ., (53)with E η ( k ⊥ ) and B η ( k ⊥ ) given by the Fourier-transformsof Eqs. (32) and (33). The evolution is thus describedsimply as the propagation of the initial conditions as afree plane wave. This is a reasonable approximation forhigh momentum modes only, and therefore appropriatefor largely UV-dominated quantities such as the energydensity or the divergence of the Chern-Simons current.Correspondingly, the full electric and magnetic fields asfunctions of τ are given by: E η ( τ, k ⊥ ) = E η ( k ⊥ ) J ( kτ ) (54) E i ( τ, k ⊥ ) = − i(cid:15) ij k j k B η ( k ⊥ ) J ( kτ ) B η ( τ, k ⊥ ) = B η ( k ⊥ ) J ( kτ ) B i ( τ, k ⊥ ) = − i(cid:15) ij k j k E η ( k ⊥ ) J ( kτ ) . In order to infer the τ -dependence in coordinate spacewe rewrite the initial condition Eq. (32) as the followingFourier transform: E η ( x ⊥ ) = − igδ ij (cid:90) d k ⊥ (2 π ) (cid:90) d u ⊥ [ α i ( u ⊥ ) , α j ( u ⊥ )] e ik ⊥ ( x − u ) ⊥ ≡ (cid:90) d k ⊥ (2 π ) E η ( k ⊥ ) e ik ⊥ x ⊥ , (55)and thus, from Eq. (54) we get to: E η ( τ, x ⊥ ) = − igδ ij (cid:90) d k ⊥ (2 π ) (cid:90) d u ⊥ [ α i ( u ⊥ ) , α j ( u ⊥ )] × J ( kτ ) e ik ⊥ ( x − u ) ⊥ . (56)Similar transformations of the initial conditions are ob-tained for E i and the chromo-magnetic fields.Let us briefly discuss the relation of our approach tothat of Refs. [33–36], where one performs a systematicpower expansion of the fields in orders of τ . In this ap-proach, the Yang-Mills equations are re-formulated as aninfinite system of equations for the coefficients of the se-ries. Taking the initial conditions Eqs. (32), (33) as thelowest order τ , the remaining coefficients are obtainedrecursively. The main drawback of this approach comesin the form of UV divergences, which are carried by eachterm of the series. These singularities originate in thespatial derivatives and give rise to the dominant termsof the expansion. The interpretation of this UV diver-gence in terms of the Bessel functions J , ( kτ ) that givethe (free) time dependence of the momentum modes isevident. Dimensionally, every power of τ is associatedwith an additional power of momentum that is integratedover, making the higher coefficients ever more UV diver-gent. Resumming the whole series into a Bessel function,on the other hand, actually improves the UV behaviorat any non-zero value of τ , because of the suppression J (0 , ( kτ ) ∼ / √ k for k (cid:29) /τ . This property is whatmakes the energy density in the MV model (even with-out any running coupling) UV finite at τ >
0, even if itis UV divergent at τ = 0 . The same asymptotic be- havior J (0 , ( kτ ) ∼ / √ τ for τ (cid:29) /Q s ∼ /k directlyresults in the large time ∼ /τ behavior of the energydensity. This is naturally interpreted as the effect ofthe boost-invariant expansion of the whole system, andwould be difficult to reach in a power series expansionin τ . This was also pointed out in Ref. [34], where thehighest-derivative terms of the Yang-Mills equation areresummed into a free field evolution. This ansatz is tech-nically identical to considering the linearized evolutionequations we use. Therefore, Eq. (53) can be understoodas a resummation of the high momentum modes of a τ -expansion of the forward light cone fields. A. One-point functions
Using the formulas presented above, the τ -dependentexpectation values of the energy density and the diver-gence of the Chern-Simons current read: (cid:104) ε ( τ, x ⊥ ) (cid:105) = (cid:104) Tr { E η E η + B η B η + E i E i + B i B i }(cid:105) = g δ ij δ kl + (cid:15) ij (cid:15) kl ) f abn f cdn (cid:90) p,k (cid:90) u,v (cid:104) α i,au α k,cv (cid:105) (cid:104) α j,bu α l,dv (cid:105) × (cid:18) J ( pτ ) J ( kτ ) − p ⊥ · k ⊥ p k J ( pτ ) J ( kτ ) (cid:19) e ip ⊥ ( x − u ) ⊥ e ik ⊥ ( x − v ) ⊥ (57) (cid:104) ˙ ν ( τ, x ⊥ ) (cid:105) = (cid:104) Tr { E η B η + E i B i }(cid:105) = g f abn f cdn (cid:90) p,k (cid:90) u,v (cid:104) α i,au α k,cv (cid:105) (cid:104) α j,bu α l,dv (cid:105) × (cid:18) δ ij (cid:15) kl J ( pτ ) J ( kτ ) − (cid:15) ij δ kl p ⊥ · k ⊥ p k J ( pτ ) J ( kτ ) (cid:19) e ip ⊥ ( x − u ) ⊥ e ik ⊥ ( x − v ) ⊥ , (58)where (cid:82) p stands for the integration over momentum (cid:82) d p ⊥ (2 π ) , while (cid:82) u corresponds to (cid:82) d u ⊥ . There aresome operations that can be performed on these gen-eral expressions before adopting a specific dipole model.For instance, by substituting the general two-point func-tion (15) in Eq. (58) one directly obtains that (cid:104) ˙ ν (cid:105) = 0 forany value of τ . This is a trivial result that reflects thatour field averages are performed over a CP-even ensem-ble. Let us now focus on the more interesting case ofthe average energy density. This object contains Fouriertransforms of Bessel functions of zeroth and first order,which can be computed analytically. Let us show this cal- For a finite UV cutoff Λ the MV model energy density has finitelimit Λ → ∞ at τ >
0, but diverges as ln Λ at τ = 0. In thelattice calculation the UV cutoff Λ corresponds to 1 /a , where a is the lattice spacing. Thus one refers to the limit Λ → ∞ asthe “continuum limit”. From this behavior one can deduce thatif one takes the limit Λ → ∞ first and then approaches τ = 0 + from above, the continuum limit energy density will diverge asln /τ . Thus, in the MV model at fixed coupling a power seriesexpansion in τ is a series expansion around a point where thecontinuum limit energy density is nonanalytic [17]. culation explicitly. Integrating over the angular variablesof the momenta θ p and θ k , the second line of Eq. (57) be-comes: (cid:90) dp p (2 π ) dk k (2 π ) (cid:18) J ( | x − u | p ) J ( | x − v | k ) J ( pτ ) J ( kτ )+cos ( θ x − u − θ x − v ) J ( | x − u | p ) J ( | x − v | k ) J ( pτ ) J ( kτ ) (cid:19) . (59)Then, applying the orthogonality condition of the Besselfunctions: (cid:90) ∞ J ν ( kr ) J ν ( sr ) rdr = δ ( k − s ) s , (60)we obtain: (cid:104) ε ( τ, x ⊥ ) (cid:105) = g δ ij δ kl + (cid:15) ij (cid:15) kl ) f abn f cdn (cid:90) u,v (cid:104) α i,au α k,cv (cid:105) (cid:104) α j,bu α l,dv (cid:105) × δ ( | x ⊥ − u ⊥ |− τ )2 πτ δ ( | x ⊥ − v ⊥ |− τ )2 πτ (1 + cos( θ x − u − θ x − v )) . (61)0 u ?
Let us now focus on the calculation of two-point func-tions. By expanding the τ -dependent chromo-electricand chromo-magnetic fields (54) in terms of the initialconditions, we obtain the following expressions:1 Qs tau
GBWMVno logs MV Q s = Q s m = Q s /
0, thus plotting the ratio φ = (cid:104) ε ( τ ) (cid:105) / (cid:104) ε (cid:105) would not make sense. (cid:104) ε ( τ, x ⊥ ) ε ( τ, y ⊥ ) (cid:105) = g δ ij δ kl + (cid:15) ij (cid:15) kl )( δ i (cid:48) j (cid:48) δ k (cid:48) l (cid:48) + (cid:15) i (cid:48) j (cid:48) (cid:15) k (cid:48) l (cid:48) ) f abn f cdn f a (cid:48) b (cid:48) m f c (cid:48) d (cid:48) m (cid:90) p,k (cid:90) ¯ p, ¯ k (cid:90) u,v (cid:90) ¯ u, ¯ v (cid:104) α i,au α k,c ¯ u α i (cid:48) ,a (cid:48) v α k (cid:48) ,c (cid:48) ¯ v (cid:105)(cid:104) α j,bu α l,d ¯ u α j (cid:48) ,b (cid:48) v α l (cid:48) ,d (cid:48) ¯ v (cid:105)× (cid:18) J ( pτ ) J (¯ pτ ) − p ⊥ · ¯ p ⊥ p ¯ p J ( pτ ) J (¯ pτ ) (cid:19)(cid:18) J ( kτ ) J (¯ kτ ) − k ⊥ · ¯ k ⊥ k ¯ k J ( kτ ) J (¯ kτ ) (cid:19) × e ip ⊥ ( x − u ) ⊥ e ik ⊥ ( y − v ) ⊥ e i ¯ p ⊥ ( x − ¯ u ) ⊥ e i ¯ k ⊥ ( y − ¯ v ) ⊥ = g δ ij δ kl + (cid:15) ij (cid:15) kl )( δ i (cid:48) j (cid:48) δ k (cid:48) l (cid:48) + (cid:15) i (cid:48) j (cid:48) (cid:15) k (cid:48) l (cid:48) ) f abn f cdn f a (cid:48) b (cid:48) m f c (cid:48) d (cid:48) m (cid:90) u,v (cid:90) ¯ u, ¯ v (cid:104) α i,au α k,c ¯ u α i (cid:48) ,a (cid:48) v α k (cid:48) ,c (cid:48) ¯ v (cid:105)(cid:104) α j,bu α l,d ¯ u α j (cid:48) ,b (cid:48) v α l (cid:48) ,d (cid:48) ¯ v (cid:105)× δ ( | x ⊥ − u ⊥ |− τ )2 πτ δ ( | x ⊥ − ¯ u ⊥ |− τ )2 πτ δ ( | y ⊥ − v ⊥ |− τ )2 πτ δ ( | y ⊥ − ¯ v ⊥ |− τ )2 πτ × (1 + cos( θ x − u − θ x − ¯ u ))(1 + cos( θ y − v − θ y − ¯ v )) (64) (cid:104) ˙ ν ( τ, x ⊥ ) ˙ ν ( τ, y ⊥ ) (cid:105) = g f abn f cdn f a (cid:48) b (cid:48) m f c (cid:48) d (cid:48) m (cid:90) p,k (cid:90) ¯ p, ¯ k (cid:90) u,v (cid:90) ¯ u, ¯ v (cid:104) α i,au α k,c ¯ u α i (cid:48) ,a (cid:48) v α k (cid:48) ,c (cid:48) ¯ v (cid:105)(cid:104) α j,bu α l,d ¯ u α j (cid:48) ,b (cid:48) v α l (cid:48) ,d (cid:48) ¯ v (cid:105)× (cid:18) δ ij (cid:15) kl J ( pτ ) J (¯ pτ ) − (cid:15) ij δ kl p ⊥ · ¯ p ⊥ p ¯ p J ( pτ ) J (¯ pτ ) (cid:19)(cid:18) δ i (cid:48) j (cid:48) (cid:15) k (cid:48) l (cid:48) J ( kτ ) J (¯ kτ ) − (cid:15) i (cid:48) j (cid:48) δ k (cid:48) l (cid:48) k ⊥ · ¯ k ⊥ k ¯ k J ( kτ ) J (¯ kτ ) (cid:19) × e ip ⊥ ( x − u ) ⊥ e ik ⊥ ( y − v ) ⊥ e i ¯ p ⊥ ( x − ¯ u ) ⊥ e i ¯ k ⊥ ( y − ¯ v ) ⊥ = g f abn f cdn f a (cid:48) b (cid:48) m f c (cid:48) d (cid:48) m (cid:90) u,v (cid:90) ¯ u, ¯ v (cid:104) α i,au α k,c ¯ u α i (cid:48) ,a (cid:48) v α k (cid:48) ,c (cid:48) ¯ v (cid:105)(cid:104) α j,bu α l,d ¯ u α j (cid:48) ,b (cid:48) v α l (cid:48) ,d (cid:48) ¯ v (cid:105) δ ( | x ⊥ − u ⊥ |− τ )2 πτ δ ( | x ⊥ − ¯ u ⊥ |− τ )2 πτ × δ ( | y ⊥ − v ⊥ |− τ )2 πτ δ ( | y ⊥ − ¯ v ⊥ |− τ )2 πτ ( δ ij (cid:15) kl + (cid:15) ij δ kl cos( θ x − u − θ x − ¯ u ))( δ i (cid:48) j (cid:48) (cid:15) k (cid:48) l (cid:48) + (cid:15) i (cid:48) j (cid:48) δ k (cid:48) l (cid:48) cos( θ y − v − θ y − ¯ v )) . (65)The physical interpretation of the results (64) and (65)is quite clear in light of our discussion of the one pointfunction in Sec. III A. It consists of energy densities atpoints x ⊥ , y ⊥ resulting from spherically expanding wavesstarting from points u ⊥ , ¯ u ⊥ , v ⊥ , ¯ v ⊥ . As discussed in Sec. II, the building block of these cor-relators is the four-point function of the gluon fields, (cid:104) α i,au α k,c ¯ u α i (cid:48) ,a (cid:48) v α k (cid:48) ,c (cid:48) ¯ v (cid:105) . In the present work we apply theGlasma Graph approximation in order to circumvent thecomplexity that this object presents. Moreover, in the2formulas presented below we also adopt the GBW model,which significantly simplifies the resulting expressions.However, even in the simplest setup available we obtainrelatively lengthy formulas. After performing the vari- able changes s ⊥ = x ⊥ − u ⊥ , ¯ s ⊥ = x ⊥ − ¯ u ⊥ , t ⊥ = y ⊥ − v ⊥ ,¯ t ⊥ = y ⊥ − ¯ v ⊥ and computing the analytically doable in-tegrals, we obtain: (cid:104) ε ( τ, x ⊥ ) ε ( τ, y ⊥ ) (cid:105) = g N c ( N c − (cid:90) π dθ s π dθ ¯ s π dθ t π dθ ¯ t π (1 + cos( θ s − θ ¯ s ))(1 + cos( θ t − θ ¯ t )) × (cid:20)(cid:18)
12 ( N c − G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ ) G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ )+ 12 G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ ) G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ )+ 12 G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ ) G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ )+ 14 G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ ) G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ )+ 14 G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ ) G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ ) (cid:19) + (cid:18) ↔ (cid:19)(cid:21) (66) (cid:104) ˙ ν ( τ, x ⊥ ) ˙ ν ( τ, y ⊥ ) (cid:105) = g N c ( N c − (cid:90) π dθ s π dθ ¯ s π dθ t π dθ ¯ t π (cid:20)(cid:18) G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ ) G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ )(1 − cos( θ s − θ ¯ s ))(1 − cos( θ t − θ ¯ t )) − G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ ) G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ )(1 − cos( θ s − θ ¯ s ))(1 − cos( θ t − θ ¯ t ))+ 12 G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ ) G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ )(1 + cos( θ s − θ ¯ s ) cos( θ t − θ ¯ t ))+ 12 G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ ) G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ )(cos( θ s − θ ¯ s ) + cos( θ t − θ ¯ t ))+ 14 G (( s − t ) τ ) G ((¯ s − ¯ t ) τ ) G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ )(1 + cos( θ s − θ ¯ s ))(1 + cos( θ t − θ ¯ t )) (cid:19) + (cid:18) ↔ (cid:19)(cid:21) , (67)where G , are the unpolarized gluon distributions de-fined by Eq. (15). The general Glasma Graph expressionincluding the linearly polarized distribution is shown inAppendix B. Here we have defined ( s − t ) τ ≡ τ (ˆ s ⊥ − ˆ t ⊥ ).This is a subtraction of two vectors s and t whose mod-uli have been replaced by τ as a result of the integrationof the delta functions. Note that this does not repre-sent a transverse shift of length τ ; instead, the value of | ( s − t ) τ | depends on the integration variables θ s , θ t . Fig-ure 4 provides a graphical interpretation of the transversecoordinate structure corresponding to this expression.The Glasma Graph approximation allows us to distin-guish between contributions stemming from “connected” or “disconnected” correlators, defined as: (cid:104) α i,ax α j,bx α k,cy α l,dy (cid:105) = disconnected (cid:122) (cid:125)(cid:124) (cid:123) (cid:104) α i,ax α j,bx (cid:105)(cid:104) α k,cy α l,dy (cid:105) + (cid:104) α i,ax α k,cy (cid:105)(cid:104) α j,bx α l,dy (cid:105) + (cid:104) α i,ax α l,dy (cid:105)(cid:104) α j,bx α k,cy (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) connected . (68)This is the usual terminology employed in previous cal-culations of correlation functions at τ = 0 + . Althoughwe have not explicitly made this distinction here, in Eqs.(66) and (67) one can identify the connected contribu-tions as those that depend on the correlation distance r ⊥ . Another remarkable aspect about the previous for-mulas lies in the interference factors carried by each term.In the case of the energy density Eq. (66) we have anoverall factor that results simply from squaring that ofthe one-point function (cid:104) ε (cid:105) . However, the Chern-Simonscharge Eq. (67) features a much more involved structure.3 y ?
In this work we have presented an analytical calcu-lation of one- and two-point correlators of energy den-sity and axial charge at finite proper times. These ob-jects characterize the average and fluctuations of energydensity and axial charge deposited throughout the ini-tial stage of HICs, during which a classical descriptionof the system is appropriate. In our calculation we as-sume a free-field evolution of the gluon fields. This setuphas been previously applied in the literature for the de-scription of dilute-dense collisions [18, 31, 32, 62, 63]. Inthese works the initial conditions are expanded in pow-ers of the weak field describing the gluon content of thedilute nucleus. To the lowest order in this expansion, thedynamics of the system are described by the linearizedYang-Mills equations. In the present paper we proposethat the same equations can be applied to the evolutionof the full initial conditions, which encode the saturationfeatures of the system. Our claim is based on the dom-inance of the high-momentum modes of a power seriesin proper time, observed in [33–35, 64]. The UV diver-gences carried by all terms of this series are effectivelyresummed by considering the higher order derivatives ofthe Yang-Mills equations [34]. We argue that nonlinearcorrections to this approach are negligible when focusing4 -4 -3 -2 -1
0 1 2 3 4 5
Qs r
Qs tau = 0Qs tau = 1Qs tau = 2Qs tau = 4 Q s r
0 1 2 3 4 5 C o v a r i an c e / a v e r age Qs r
Qs tau = 0Qs tau = 1Qs tau = 2Qs tau = 4 Q s ⌧ = 0 Q s ⌧ = 1 Q s ⌧ = 2 Q s ⌧ = 4
0 1 2 3 4 5
Qs r
Qs tau = 0Qs tau = 1Qs tau = 2Qs tau = 4 Q s r
FIG. 5. Correlation functions of the energy density and the divergence of the Chern-Simons current at different values of τ for the GBW model (left plot) and the MV model (right plot). The thin black lines in the left plot correspond to Eqs. (40)and (41). on quantities dominated by large momenta, such as theenergy density and the divergence of the Chern-Simonscurrent.In this paper we evaluate the initial conditions inthe Glasma Graph approximation, which assumes aGaussian-like decomposition of four-point correlatorsinto two-point correlators for the gluon fields. In or-der to obtain quantitative results, we adopt the GBWand MV models, including running coupling correctionsin the last case. Our running coupling regularization ofthe MV model yields a finite result in the UV limit, whilepreserving a perturbative power-law tail at large trans-verse momenta of the unintegrated gluon distribution (asignature feature of the MV model). This running cou-pling approach differs from that of previous works [18],hence yielding slightly different curves at τ = 0 + . In theGBW model, in contrast, the unintegrated gluon distri-bution falls off as a Gaussian at large transverse momen-tum. Nevertheless, for the energy density correlator thetwo parametrizations give qualitatively similar results.Yang-Mills evolution was analytically computed for allcorrelators up until a proper time τ = 4 /Q s . The evolu-tion of the average energy density deposited on a givenpoint can be described through a dimensionless dilutionfactor, which accounts for both the attenuation and theinterference of the expanding plane waves that describethe gluon propagation in the plane transverse to the col-lision. Noticeably, we observe that in the MV model thisdilution effect is more rapid than for the GBW dipole cor-relator. This can be naturally understood as the effect ofa larger relative contribution of fast-evolving large trans-verse momentum modes in the MV model, where theyare only suppressed as a power law and not a Gaussian.The effect of time evolution over the considered two-pointcorrelators can be described as an elongation of the cor-relation length, a trend that is suggestive of a transitiontowards the hydrodynamical regime.The Glasma Graph approximation has been adopted here merely for calculational convenience, and we expectit to be a relatively well satisfied apart from the longestcoordinate separations r > ∼ /Q s . However, our approachfor computing the time dependence is much more general,and can be applied to any initial color fields, as long asthe four-point function of the Weizs¨acker-Williams fieldsof the colliding nuclei can be computed. A reliable ap-plication of our results into the study of eccentricity har-monics will require doing this and re-evaluating our τ -dependent correlators in the full MV model, i.e. beyondthe Glasma Graph approximation. This is due to thesensitivity of the mean-squared eccentricity fluctuationsto large-scale correlations ( Q s r > V. ACKNOWLEDGEMENTS
The authors thank S¨oren Schlichting for providinguseful insight on the running coupling prescription pre-sented in [18]. PGR thanks Cyrille Marquet and HeikkiM¨antysaari for their useful feedback during the develop-ment of this work. TL thanks Matt Sievert for pertinentquestions on the time dependence of the Glasma energydensity. This work was supported by the Academy ofFinland, project 321840 and under the European Union’sHorizon 2020 research and innovation programme bythe European Research Council (ERC, grant agreementNo. ERC-2015-CoG-681707) and by the STRONG-2020project (grant agreement No 824093). The content ofthis article does not reflect the official opinion of the Eu-ropean Union and responsibility for the information andviews expressed therein lies entirely with the authors.5
A. THE DIVERGENCE OF THE CHERN-SIMONS CURRENT
The chiral anomaly of QCD induces a transformationof left- into right-handed quarks, or equivalently, a gen-eration of axial charge N : dN dt = d ( N R − N L ) dt = − g N f π (cid:90) d x Tr (cid:110) F µν ( x ) ˜ F µν ( x ) (cid:111) . (A1)Here, N f represents the number of flavors. In this for-mula we can see how the production rate of axial charge isdirectly related to the properties of the gauge fields enter-ing the right-hand side of the equation. This term implic-itly includes various contributions to axial charge produc-tion, one of them stemming directly from the topologicalstructure of QCD.We can define distinct topological classes of gauge fieldconfigurations labeled by the following quantity: Q w = g π (cid:90) d x Tr (cid:110) F µν ( x ) ˜ F µν ( x ) (cid:111) , (A2)a topological invariant known as winding number. Theseclasses are separated by potential barriers with heightsof order Λ QCD , which suppress the Q w fluctuations atlow temperature. However, that is not the case in a hightemperature medium such as the QGP. At T ∼ Λ QCD , Q w fluctuations happen with a rate that can be directlyrelated to dN /dt . This relation is expressed in terms of the Chern-Simons current, defined as: K µ = (cid:15) µνρσ A aν (cid:16) F aρσ + g A bρ A cσ (cid:17) , (A3)and whose divergence reads: ∂ µ K µ = −
14 Tr (cid:110) F µν ( x ) ˜ F µν ( x ) (cid:111) ≡ ˙ ν ( x ) . (A4)It is convenient to rewrite Eq. (A1) in terms of the diver-gence of the Chern-Simons current: dN dt = g N f π (cid:90) d x ˙ ν ( x ) . (A5)Based on this relation, in the present work we take ˙ ν as the fundamental object controlling axial charge gener-ation. This effect may manifest experimentally in off-central HICs, where large background electromagneticfields are generated [65]. These fields, when in the pres-ence of deconfined chirally-unbalanced matter, induce aseparation of charges that gives rise to an electric dipolealong the direction of angular momentum. This is knownas the Chiral Magnetic Effect (CME), and it translatesinto nontrivial azimuthal correlations in the hadron spec-trum [66, 67]. This signal, however, is obscured by largebackground effects that produce similar back-to-back cor-relations. Reducing this uncertainty is currently a majorgoal of the high energy QCD community, both on theexperimental [68–72] and the theoretical side [18, 73–75].In this context, it is of paramount importance to con-strain the dynamical origin of CME signatures. The cal-culation of ˙ ν correlators presented in this paper accountsfor the τ -evolution of potential contributions emergingfrom event-by-event color charge fluctuations during theGlasma phase. This source, although unrelated to thetopological structure of QCD, can have a sizeable effecton final state observables related to CP violation. The ex-ploration of this possibility is left for future phenomeno-logical studies.6 -4 -3 -2 -1
0 1 2 3 4 5
Qs r
MV tau=2MV (no logs) tau=2MV (no running coupling) tau=2 Q s r
0 1 2 3 4 5
Qs r
MV tau=2MV (no logs) tau=2MV (no running coupling) tau=2 Q s r
B. TWO-POINT FUNCTION OF ENERGY DENSITY IN THE GLASMA GRAPH APPROXI-MATION
In the following appendix we display the general expression of the correlator of the energy density deposited in twopoints x ⊥ and y ⊥ of the transverse plane at a proper time τ within the Glasma Graph approximation: (cid:104) ε ( τ, x ⊥ ) ε ( τ, y ⊥ ) (cid:105) = g N c ( N c − (cid:90) π dθ s π dθ ¯ s π dθ t π dθ ¯ t π (1 + cos( θ s − θ ¯ s ))(1 + cos( θ t − θ ¯ t )) × (cid:20)(cid:18)
12 ( N c − G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ ) G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ )+ 12 G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ ) (cid:16) G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ ) + h (( s − t ) τ + r ⊥ ) h ((¯ s − ¯ t ) τ + r ⊥ ) (cid:17) + 12 G (( s − ¯ s ) τ ) G (( t − ¯ t ) τ ) (cid:16) G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ ) + h (( s − ¯ t ) τ + r ⊥ ) h ((¯ s − t ) τ + r ⊥ ) (cid:17) + 14 G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ ) G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ )+ 14 h (( s − t ) τ + r ⊥ ) h ((¯ s − ¯ t ) τ + r ⊥ ) h (( s − t ) τ + r ⊥ ) h ((¯ s − ¯ t ) τ + r ⊥ )+ 14 G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ ) G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ )+ 14 h (( s − ¯ t ) τ + r ⊥ ) h ((¯ s − t ) τ + r ⊥ ) h (( s − ¯ t ) τ + r ⊥ ) h ((¯ s − t ) τ + r ⊥ )+ 14 G (( s − t ) τ + r ⊥ ) G ((¯ s − ¯ t ) τ + r ⊥ ) h (( s − ¯ t ) τ + r ⊥ ) h ((¯ s − t ) τ + r ⊥ )+ 14 G (( s − ¯ t ) τ + r ⊥ ) G ((¯ s − t ) τ + r ⊥ ) h (( s − t ) τ + r ⊥ ) h ((¯ s − ¯ t ) τ + r ⊥ ) (cid:19) + (cid:18) ↔ (cid:19)(cid:21) . (B1)The indicated integrals were performed numerically with Mathematica to obtain the curves shown in Figure 5. C. TWO-POINT FUNCTION OF ENERGYDENSITY IN THE MV MODEL WITHFIXED COUPLING
The energy density of the Glasma fields in the MVmodel is finite for all τ >
0, which allows to computethe τ -evolution of its correlators without resorting to anyUV regularization scheme. Under this approach we couldcompare results in the GBW and fixed coupling MV mod- els exclusively at finite τ values, or alternatively, regular-ize the UV divergence only at τ = 0 + . Instead of that,in this work we adopt the running coupling prescriptiondescribed in Sec. II A, which defines a variant of the MVmodel that we consistently apply throughout the whole τ range considered. This ansatz also gives rise to lessnumerically demanding integrals.However, the fixed coupling MV model yields some7nontrivial results that are worth remarking. In Fig. 6one can see that the correlator corresponding to the MVmodel with no running coupling (thick curve) featuresa significantly steeper fall-off at large transverse separa-tions. The correlator enters this fast decay region some-what abruptly (at r ≈ /Q s for τ = 1 /Q s and at r ≈ /Q s for τ = 2 /Q s ), giving rise to a nontrivial profile that hintsat the formation of a relatively flat plateau at short corre-lation distances. Remarkably, this characteristic shape isalso visible in the correlators obtained under our runningcoupling prescription, although it appears less promi-nently and at later proper times (around τ = 2 /Q s , ascan be seen in Fig. 6(b)). In Fig. 6 we also display the correlator obtained byapplying the running coupling prescription Eq. (24) toall the coupling constants of Eqs. (22) and (23) (greendashed curve). This curve, labeled as “MV (no logs)”,shows none of the features mentioned above, being qual-itatively identical to the GBW result. The comparisonbetween these correlators supports our particular choicefor the running coupling prescription, as it was madewith the aim of preserving (to the extent possible) thenontrivial aspects of the fixed coupling MV model whilesimultaneously providing finite results at τ = 0 + . [1] M. Luzum and H. Petersen, Initial State Fluctuationsand Final State Correlations in Relativistic Heavy-IonCollisions , J. Phys. G (2014) 063102[ arXiv:1312.5503 [nucl-th] ].[2] F. Gelis, E. Iancu, J. Jalilian-Marian andR. Venugopalan, The Color Glass Condensate , Ann.Rev. Nucl. Part. Sci. (2010) 463 [ arXiv:1002.0333[hep-ph] ].[3] T. LAPPI, Small x physics and RHIC data , Int. J. Mod.Phys. E (2011) 1 [ arXiv:1003.1852 [hep-ph] ].[4] J. L. Albacete and C. Marquet, Gluon saturation andinitial conditions for relativistic heavy ion collisions , Prog. Part. Nucl. Phys. (2014) 1 [ arXiv:1401.4866[hep-ph] ].[5] L. D. McLerran and R. Venugopalan, Computing quarkand gluon distribution functions for very large nuclei , Phys. Rev. D (1994) 2233 [ arXiv:hep-ph/9309289 ].[6] L. D. McLerran and R. Venugopalan, Gluon distributionfunctions for very large nuclei at small transversemomentum , Phys. Rev. D (1994) 3352[ arXiv:hep-ph/9311205 ].[7] L. D. McLerran and R. Venugopalan, Green’s functionsin the color field of a large nucleus , Phys. Rev. D (1994) 2225 [ arXiv:hep-ph/9402335 ].[8] T. Lappi and L. McLerran, Some features of the glasma , Nucl. Phys. A (2006) 200 [ arXiv:hep-ph/0602189 ].[9] K. Dusling, T. Epelbaum, F. Gelis and R. Venugopalan,
Role of quantum fluctuations in a system with strongfields: Onset of hydrodynamical flow , Nucl. Phys. A (2011) 69 [ arXiv:1009.4363 [hep-ph] ].[10] F. Gelis,
Initial state and thermalization in the ColorGlass Condensate framework , Int. J. Mod. Phys. E (2015) 1530008 [ arXiv:1508.07974 [hep-ph] ].[11] J. Berges, M. P. Heller, A. Mazeliauskas andR. Venugopalan, Thermalization in QCD: theoreticalapproaches, phenomenological applications, andinterdisciplinary connections , arXiv:2005.12299[hep-th] .[12] A. Kurkela, Initial state of Heavy-Ion Collisions:Isotropization and thermalization , Nucl. Phys. A (2016) 136 [ arXiv:1601.03283 [hep-ph] ].[13] A. Kurkela and E. Lu,
Approach to Equilibrium inWeakly Coupled Non-Abelian Plasmas , Phys. Rev. Lett. (2014) 182301 [ arXiv:1405.6318 [hep-ph] ].[14] A. Kurkela and Y. Zhu,
Isotropization and hydrodynamization in weakly coupled heavy-ioncollisions , Phys. Rev. Lett. (2015) 182301[ arXiv:1506.06647 [hep-ph] ].[15] L. Keegan, A. Kurkela, A. Mazeliauskas and D. Teaney,
Initial conditions for hydrodynamics from weaklycoupled pre-equilibrium evolution , JHEP (2016) 171[ arXiv:1605.04287 [hep-ph] ].[16] A. Kurkela, A. Mazeliauskas, J.-F. Paquet,S. Schlichting and D. Teaney, Effective kineticdescription of event-by-event pre-equilibrium dynamicsin high-energy heavy-ion collisions , Phys. Rev. C (2019) 034910 [ arXiv:1805.00961 [hep-ph] ].[17] T. Lappi, Energy density of the glasma , Phys. Lett. B (2006) 11 [ arXiv:hep-ph/0606207 ].[18] T. Lappi and S. Schlichting,
Linearly polarized gluonsand axial charge fluctuations in the Glasma , Phys. Rev.D (2018) 034034 [ arXiv:1708.08625 [hep-ph] ].[19] J. L. Albacete, P. Guerrero-Rodr´ıguez and C. Marquet, Initial correlations of the Glasma energy-momentumtensor , JHEP (2019) 073 [ arXiv:1808.00795[hep-ph] ].[20] R. S. Bhalerao, G. Giacalone, P. Guerrero-Rodr´ıguez,M. Luzum, C. Marquet and J.-Y. Ollitrault, Relatingeccentricity fluctuations to density fluctuations inheavy-ion collisions , Acta Phys. Polon. B (2019)1165 [ arXiv:1903.06366 [nucl-th] ].[21] G. Giacalone, P. Guerrero-Rodr´ıguez, M. Luzum,C. Marquet and J.-Y. Ollitrault, Fluctuations inheavy-ion collisions generated by QCD interactions inthe color glass condensate effective theory , Phys. Rev. C (2019) 024905 [ arXiv:1902.07168 [nucl-th] ].[22] F. Gelis, G. Giacalone, P. Guerrero-Rodr´ıguez,C. Marquet and J.-Y. Ollitrault,
Primordial fluctuationsin heavy-ion collisions , arXiv:1907.10948 [nucl-th] .[23] G. Giacalone, F. Gelis, P. Guerrero-Rodr´ıguez,M. Luzum, C. Marquet and J.-Y. Ollitrault, Evolutionof fluctuations in the initial state of heavy-ion collisionsfrom RHIC to LHC , Springer Proc. Phys. (2020)453 [ arXiv:1911.04720 [nucl-th] ].[24] B. Schenke, P. Tribedy and R. Venugopalan,
Event-by-event gluon multiplicity, energy density, andeccentricities in ultrarelativistic heavy-ion collisions , Phys. Rev. C (2012) 034908 [ arXiv:1206.6805[hep-ph] ].[25] B. Schenke, P. Tribedy and R. Venugopalan, Fluctuating Glasma initial conditions and flow in heavyion collisions , Phys. Rev. Lett. (2012) 252301[ arXiv:1202.6646 [nucl-th] ].[26] H. M¨antysaari, B. Schenke, C. Shen and P. Tribedy,
Imprints of fluctuating proton shapes on flow inproton-lead collisions at the LHC , Phys. Lett. B (2017) 681 [ arXiv:1705.03177 [nucl-th] ].[27] A. Krasnitz and R. Venugopalan,
Nonperturbativecomputation of gluon minijet production in nuclearcollisions at very high-energies , Nucl. Phys. B (1999) 237 [ arXiv:hep-ph/9809433 ].[28] A. Krasnitz, Y. Nara and R. Venugopalan,
Coherentgluon production in very high-energy heavy ioncollisions , Phys. Rev. Lett. (2001) 192302[ arXiv:hep-ph/0108092 ].[29] A. Kovner, L. D. McLerran and H. Weigert, Gluonproduction at high transverse momentum in theMcLerran-Venugopalan model of nuclear structurefunctions , Phys. Rev. D (1995) 3809[ arXiv:hep-ph/9505320 ].[30] Y. V. Kovchegov and D. H. Rischke, Classical gluonradiation in ultrarelativistic nucleus-nucleus collisions , Phys. Rev. C (1997) 1084 [ arXiv:hep-ph/9704201 ].[31] A. Dumitru and L. D. McLerran, How protons shattercolored glass , Nucl. Phys. A (2002) 492[ arXiv:hep-ph/0105268 ].[32] L. McLerran and V. Skokov,
Odd Azimuthal Anisotropyof the Glasma for pA Scattering , Nucl. Phys. A (2017) 83 [ arXiv:1611.09870 [hep-ph] ].[33] R. J. Fries, J. I. Kapusta and Y. Li,
Near-fields andinitial energy density in the color glass condensatemodel , arXiv:nucl-th/0604054 .[34] H. Fujii, K. Fukushima and Y. Hidaka, Initial energydensity and gluon distribution from the Glasma inheavy-ion collisions , Phys. Rev. C (2009) 024909[ arXiv:0811.0437 [hep-ph] ].[35] G. Chen, R. J. Fries, J. I. Kapusta and Y. Li, EarlyTime Dynamics of Gluon Fields in High Energy NuclearCollisions , Phys. Rev. C (2015) 064912[ arXiv:1507.03524 [nucl-th] ].[36] M. E. Carrington, A. Czajka and S. Mrowczynski, Theenergy-momentum tensor at the earliest stage ofrelativistic heavy ion collisions – formalism , arXiv:2012.03042 [hep-ph] .[37] J. P. Blaizot, T. Lappi and Y. Mehtar-Tani, On thegluon spectrum in the glasma , Nucl. Phys. A (2010)63 [ arXiv:1005.0955 [hep-ph] ].[38] K. Boguslavski, A. Kurkela, T. Lappi and J. Peuron,
Broad excitations in a 2+1D overoccupied gluon plasma , arXiv:2101.02715 [hep-ph] .[39] K. Boguslavski, A. Kurkela, T. Lappi and J. Peuron, Highly occupied gauge theories in 2+1 dimensions: Aself-similar attractor , Phys. Rev. D (2019) 094022[ arXiv:1907.05892 [hep-ph] ].[40] J. Jalilian-Marian, A. Kovner, A. Leonidov andH. Weigert,
The Wilson renormalization group for low xphysics: Towards the high density regime , Phys. Rev. D (1998) 014014 [ arXiv:hep-ph/9706377 ].[41] J. Jalilian-Marian, A. Kovner and H. Weigert, TheWilson renormalization group for low x physics: Gluonevolution at finite parton density , Phys. Rev. D (1998) 014015 [ arXiv:hep-ph/9709432 ].[42] H. Weigert, Unitarity at small Bjorken x , Nucl. Phys. A (2002) 823 [ arXiv:hep-ph/0004044 ]. [43] E. Iancu, A. Leonidov and L. D. McLerran,
Nonlineargluon evolution in the color glass condensate. 1. , Nucl.Phys. A (2001) 583 [ arXiv:hep-ph/0011241 ].[44] E. Ferreiro, E. Iancu, A. Leonidov and L. McLerran,
Nonlinear gluon evolution in the color glass condensate.2. , Nucl. Phys. A (2002) 489[ arXiv:hep-ph/0109115 ].[45] I. Balitsky,
Operator expansion for high-energyscattering , Nucl. Phys. B (1996) 99[ arXiv:hep-ph/9509348 ].[46] Y. V. Kovchegov,
Small x F(2) structure function of anucleus including multiple pomeron exchanges , Phys.Rev. D (1999) 034008 [ arXiv:hep-ph/9901281 ].[47] A. Metz and J. Zhou, Distribution of linearly polarizedgluons inside a large nucleus , Phys. Rev. D (2011)051503 [ arXiv:1105.1991 [hep-ph] ].[48] K. J. Golec-Biernat and M. Wusthoff, Saturation effectsin deep inelastic scattering at low Q**2 and itsimplications on diffraction , Phys. Rev. D (1998)014017 [ arXiv:hep-ph/9807513 ].[49] T. Lappi, Wilson line correlator in the MV model:Relating the glasma to deep inelastic scattering , Eur.Phys. J. C (2008) 285 [ arXiv:0711.3039 [hep-ph] ].[50] J. Jalilian-Marian, A. Kovner, L. D. McLerran andH. Weigert, The Intrinsic glue distribution at very smallx , Phys. Rev. D (1997) 5414[ arXiv:hep-ph/9606337 ].[51] K. Dusling, F. Gelis, T. Lappi and R. Venugopalan, Long range two-particle rapidity correlations in A+Acollisions from high energy QCD evolution , Nucl. Phys.A (2010) 159 [ arXiv:0911.2720 [hep-ph] ].[52] A. Dumitru, K. Dusling, F. Gelis, J. Jalilian-Marian,T. Lappi and R. Venugopalan,
The Ridge inproton-proton collisions at the LHC , Phys. Lett. B (2011) 21 [ arXiv:1009.5295 [hep-ph] ].[53] K. Dusling and R. Venugopalan,
Azimuthal collimationof long range rapidity correlations by strong color fieldsin high multiplicity hadron-hadron collisions , Phys. Rev.Lett. (2012) 262001 [ arXiv:1201.2658 [hep-ph] ].[54] K. Dusling and R. Venugopalan,
Evidence for BFKLand saturation dynamics from dihadron spectra at theLHC , Phys. Rev. D (2013) 051502 [ arXiv:1210.3890[hep-ph] ].[55] K. Dusling and R. Venugopalan, Explanation ofsystematics of CMS p+Pb high multiplicity di-hadrondata at √ s NN = 5 . TeV , Phys. Rev. D (2013)054014 [ arXiv:1211.3701 [hep-ph] ].[56] R. Mertig, M. Bohm and A. Denner, FEYN CALC:Computer algebraic calculation of Feynman amplitudes , Comput. Phys. Commun. (1991) 345.[57] V. Shtabovenko, R. Mertig and F. Orellana, NewDevelopments in FeynCalc 9.0 , Comput. Phys.Commun. (2016) 432 [ arXiv:1601.01167[hep-ph] ].[58] K. Boguslavski, A. Kurkela, T. Lappi and J. Peuron,
Spectral function for overoccupied gluodynamics fromreal-time lattice simulations , Phys. Rev. D (2018)014006 [ arXiv:1804.01966 [hep-ph] ].[59] Y. V. Kovchegov and K. Tuchin, Inclusive gluonproduction in DIS at high parton density , Phys. Rev. D (2002) 074026 [ arXiv:hep-ph/0111362 ].[60] J. P. Blaizot, F. Gelis and R. Venugopalan, High-energypA collisions in the color glass condensate approach. 1.Gluon production and the Cronin effect , Nucl. Phys. A (2004) 13 [ arXiv:hep-ph/0402256 ].[61] M. Li and V. V. Skokov, First Saturation Correction inHigh Energy Proton-Nucleus Collisions: I. Timeevolution of classical Yang-Mills fields beyond leadingorder , arXiv:2102.01594 [hep-ph] .[62] Y. V. Kovchegov, Can thermalization in heavy ioncollisions be described by QCD diagrams? , Nucl. Phys.A (2005) 298 [ arXiv:hep-ph/0503038 ].[63] Y. V. Kovchegov,
Thoughts on non-perturbativethermalization and jet quenching in heavy ion collisions , Nucl. Phys. A (2006) 476 [ arXiv:hep-ph/0507134 ].[64] R. J. Fries,
Early Time Evolution of High Energy HeavyIon Collisions , J. Phys. G (2007) S851[ arXiv:nucl-th/0702026 ].[65] V. Skokov, A. Y. Illarionov and V. Toneev, Estimate ofthe magnetic field strength in heavy-ion collisions , Int.J. Mod. Phys. A (2009) 5925 [ arXiv:0907.1396[nucl-th] ].[66] D. Kharzeev, A. Krasnitz and R. Venugopalan, Anomalous chirality fluctuations in the initial stage ofheavy ion collisions and parity odd bubbles , Phys. Lett.B (2002) 298 [ arXiv:hep-ph/0109253 ].[67] D. E. Kharzeev, L. D. McLerran and H. J. Warringa,
The Effects of topological charge change in heavy ioncollisions: ’Event by event P and CP violation’ , Nucl.Phys. A (2008) 227 [ arXiv:0711.0950 [hep-ph] ].[68] S. A. Voloshin,
Parity violation in hot QCD: How todetect it , Phys. Rev. C (2004) 057901 [ arXiv:hep-ph/0406311 ].[69] A. Bzdak, V. Koch and J. Liao, Charge-DependentCorrelations in Relativistic Heavy Ion Collisions andthe Chiral Magnetic Effect , Lect. Notes Phys. (2013) 503 [ arXiv:1207.7327 [nucl-th] ].[70] F. Wen, J. Bryon, L. Wen and G. Wang,
Event-shape-engineering study of charge separation inheavy-ion collisions , Chin. Phys. C (2018) 014001[ arXiv:1608.03205 [nucl-th] ].[71] H.-j. Xu, J. Zhao, X. Wang, H. Li, Z.-W. Lin, C. Shenand F. Wang, Varying the chiral magnetic effect relativeto flow in a single nucleus-nucleus collision , Chin. Phys.C (2018) 084103 [ arXiv:1710.07265 [nucl-th] ].[72] H. Huang et. al. in , p. TUPAF006, 2018.[73] S. Schlichting and S. Pratt, Charge conservation atenergies available at the BNL Relativistic Heavy IonCollider and contributions to local parity violationobservables , Phys. Rev. C (2011) 014913[ arXiv:1009.4283 [nucl-th] ].[74] P. Guerrero-Rodr´ıguez, Topological charge fluctuationsin the Glasma , JHEP (2019) 026 [ arXiv:1903.11602[hep-ph] ].[75] J. Hammelmann, A. Soto-Ontoso, M. Alvioli, H. Elfnerand M. Strikman, Influence of the neutron-skin effect onnuclear isobar collisions at energies available at theBNL Relativistic Heavy Ion Collider , Phys. Rev. C (2020) 061901 [ arXiv:1908.10231 [nucl-th]arXiv:1908.10231 [nucl-th]