Initial correlations of the Glasma energy-momentum tensor

We present an analytical calculation of the covariance of the energy-momentum tensor associated to the gluon field produced in ultra-relativistic heavy ion collisions at early times, the Glasma. This object involves the two-point and single-point correlators of the energy-momentum tensor (〈Tμν (x⊥)Tσρ(y⊥)〉 and 〈Tμν (x⊥)〉, respectively) at proper time τ = 0+. Our approach is based on the Color Glass Condensate effective theory, which allows us to map the fluctuations of the valence color sources in the colliding nuclei to those of the energy-momentum tensor of the produced gluon fields via the solution of the classical equations of motion in the presence of external currents. The color sources in the two colliding nuclei are characterized by Gaussian correlations, albeit in more generality than in the McLerran-Venugopalan model, allowing for non-trivial impact parameter and transverse dependence of the two-point correlator. We compare our results to those obtained under the Glasma Graph approximation, finding agreement in the limit of short transverse separations. However, important differences arise at larger transverse separations, where our result displays a slower fall-off than the Glasma Graph result (1/r2 vs. 1/r4 power-law decay), indicating that the color screening of the correlations in the transverse plane occurs at distances larger than 1/Qs by a logarithmic factor sensitive to the infrared. In the Glasma flux tube picture, this implies that the color domains are larger than originally estimated.


Introduction
Understanding the dynamical features of the matter produced in the early stages of heavy ion collisions and its eventual thermalization into a Quark Gluon Plasma (QGP from now on) is one of the most pressing questions in the field of heavy ion collisions, both at the experimental and theoretical levels.
The study of the correlations between the detected particles plays a main role for the understanding of the problem of QGP formation -and its characterization -in heavy ion collisions, since non-trivial correlations provide a clear indication of the collective behavior of the produced medium. However, it has also become clear over the last years that the observed correlations reflect as much the collective dynamics of the produced medium as they do the initial state correlations, namely those dynamically generated during the early stages of the collision (before an eventual thermalization of the system) or already built in the wave function of the colliding nuclei (see e.g. [1]). This observation relates to the very small ratio of viscosity over entropy density extracted from hydrodynamical simulations or, equivalently, to the low dissipation of the dynamics mapping early and late times of the collision [2]. Therefore, a detailed and theoretically robust characterization of initial JHEP01(2019)073 state correlations is mandatory for a proper understanding of the medium transport and dynamical properties.
Indeed, the existence in the literature of a broad variety of phenomenological models for the description of the initial stages of heavy ion collisions reflects the importance of this kind of studies (for a review see e.g. [3]). The main practical use of such models is to generate initial conditions of the energy density and velocity profiles for further evolution of the system, typically described by quasi-ideal relativistic hydrodynamics during the QGP phase followed by kinetic transport during the hadronic afterburner. All such models allow for fluctuations of the energy and momentum deposited in the collision area. The dynamical origin and practical description of such fluctuations vary largely from model to model -from the positions of the nucleons in the transverse plane at the collision time to fluctuations of the sub-nucleonic degrees of freedom and their inelasticity density -and are, in general, subject to a large degree of phenomenological modeling. Clearly, a higher degree of theoretical control on the description of the initial collision profile is desirable, and adding to it is precisely the main goal of this work.
The classical approach that we shall follow is embodied in the Color Glass Condensate (CGC) effective theory (see e.g. [4,5] for a review), arguably the most complete theoretical framework for the description of the early time dynamics in heavy ion collisions. The 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 based on the fact that for very large occupation numbers the quantum fluctuations represent a negligible correction to the strong background field. Quantum corrections are incorporated in the CGC framework via the JIMWLK renormalization group equations [6][7][8][9][10][11][12][13]. They ensure that the physical observables are independent of the arbitrary longitudinal momentum scale at which the separation between slow (dynamical) and fast (static sources) degrees of freedom, on which the CGC effective theory is build up, is performed.
The properties of the medium produced in heavy ion collisions at early times, dubbed as Glasma, have been extensively studied in a series of works in the CGC framework [14][15][16][17]. This kind of studies start by solving the classical equations of motion for the produced gluon field in the presence of two external color sources -the valence degrees of freedom of the two colliding nuclei. The picture that emerges is that of the Glasma as a strongly correlated, maximally anisotropic system dominated by strong classical fields. The fact that the chromo-electric and magnetic fields are parallel to the collision axis immediately after the collision leads to a very peculiar form for the energy-momentum tensor [14], where 0 is the initial average energy density. The most striking feature is that the longitudinal pressure is negative, reminiscent of strings, or flux tubes, stretching in the longitudinal direction. This picture is reinforced by the observation that the correlations of the classical fields extend over long rapidities in the longitudinal direction. In turn, correlations on the plane transverse to the collision axis are expected to be short-range -parametrically of the order of the inverse of the saturation scale ∼ 1/Q s , much smaller than the nucleon size -since color charges in the projectile (or target) are correlated only over this typical distance, which effectively plays the role of the scale for color neutrality in the nuclear wave function. In this work we shall provide explicit results

JHEP01(2019)073
to quantify the size and extent of the transverse correlations, not entirely supporting the qualitative expectations on its short-range character. How the above-described coherent ensemble of classical flux tubes decays and whether it eventually thermalizes into a QGP is a subject of open debate and intense investigation over the last years and is beyond the scope of this work. For a review we refer the reader to [18] or to the more recent works on the matching of the CGC description with effective kinetic theory as an intermediate dynamical step before the hydrodynamization of the system (e.g. [19]).
Rather, our goal in this work is to further explore the properties of the strictly classical Glasma dynamics by presenting a first analytical calculation of the two-point correlator of the energy-momentum tensor of the Glasma fields right after the collision time. We start from the assumption that the relevant correlations among the fast color charges in the wave function of the colliding nuclei are known. Then we calculate how the collision dynamics, described under the classical approximation, maps such correlations onto correlations of the energy-momentum tensor of the produced gluon field right after the collision. Hence, the only source of fluctuations in our approach is that of the incoming color sources, since the collision dynamics are fully deterministic in the leading-order classical approach. Detailed knowledge about them or, more generally, about the wave function of the colliding nuclei, can be obtained either at dedicated experiments like the proposed Electron Ion Collider [20] or, in the absence of direct empiric data, via theoretical modeling sustained by the abundant empiric information on the proton partonic structure provided by the HERA experiment (see [3]).
Specifically, we perform the analytical calculation of the following covariance: where T µν 0 (x ⊥ ) is the energy-momentum tensor (EMT) associated to the gluon field produced over an infinitesimal positive proper time after two heavy ion nuclei with mass numbers A 1 , A 2 collide at relativistic speed. We rely on an extended version of the McLerran-Venugopalan (MV) model for the description of the valence color sources of the colliding nuclei [21], whereby we assume that they obey Gaussian statistics, as in the original MV model, but we allow for a more general form of the two-point correlator, ρ a (x − , x ⊥ )ρ b (y − , y ⊥ ) , in order to expand the possibilities for phenomenological applications. Our specific modifications consist of relaxing the assumption of local transverse interactions, as well as including an explicit impact parameter dependence that allows the possibility of describing finite, non-homogeneous nuclei. However, for the sake of simplicity, in some sections we shall discuss our results in terms of the original MV model.
Following the approach outlined above, and despite the complexity of the calculation and of the full result, we obtain a remarkably compact expression for the covariance of the EMT in the limit of large transverse separations, rQ s 1 with r ≡ |x ⊥ − y ⊥ |: (1. 2) The factors Q s1,2 (r ⊥ , b ⊥ ) andQ s1,2 (b ⊥ ) -two definitions of the saturation scales characterizing each nuclei -will be introduced later along with the factor L(0 ⊥ ). Eq. (1.2) is one

JHEP01(2019)073
of the most important results of the paper, as it could challenge the conjectured physical picture of Glasma flux tubes or color field domains [16] -which basically states that when two sheets of CGC pass through each other, color flux tubes of transverse size 1/Q s are created. In our result, although the transverse correlation length is parametrically of the order of 1/Q s , the correlations decrease only according to a 1/r 2 power-law tail at large distances, extending further in the transverse plane than what was indicated by previous calculations. Such slowly vanishing covariance could potentially have an impact in both physical interpretations and numerical results for any observable built from this quantity. For instance, the 2-dimensional transverse integral of eq. (1.2) will be dominated by the upper bound (the infrared cut-off r ∼ 1/m) rather than the lower bound r ∼ 1/Q s , which is what happens under the Glasma Graph approximation [22] (that features a 1/r 4 fall-off as we will discuss later), or even in the case of a more naive exponential fall-off. This indicates that the range of the transverse color screening of the correlations, which determines the size of the color domains in the interaction region, is actually bigger than 1/Q s , as it features a logarithmic enhancement ln(Q s /m) sensitive to the infrared. Similar observations were made in [23] in the context of two-particle correlations: a sensitivity of the color domain size to the infrared was observed numerically, with it getting larger as the infrared cut-off was decreased. In the case of EMT correlations, our qualitative discussion also remains to be quantified with numerical calculations.
As an input to hydrodynamical simulations, eq. (1.2) also has important implications. Indeed, neglecting logarithmic dependencies, we can write In Monte Carlo Glauber models, where eccentricity fluctuations are created by uncorrelated, small-scale fluctuations in the transverse plane, this quantity is taken as a constant proportional to d 2 x ⊥ T 00 0 (x ⊥ ) [24]. In our calculation at τ = 0 + , which takes into account sub-nucleonic degrees of freedom (but nevertheless give rise to long-range correlations), that quantity is not a constant. The dimensionless ratio of eq. (1.3) to the integrated energy density characterizes the strength of the eccentricity fluctuations [25], and in our calculation is given by (8π/(N 2 Therefore, we find this ratio bigger in the middle of the overlap region than near the edge, which brings new insight for the characterization of the initial stage of heavy-ion collisions. This ratio also displays the usual 1/(N 2 c − 1) suppression characteristic of non-trivial color correlations.
In more general terms, the calculation presented in this work provides further analytical insight to the dynamics of the classical fields produced in relativistic heavy ion collisions, otherwise also accounted for in the well-known IP-Glasma model [26,27] and related numerical methods, where the classical equations of motion that we discuss here are solved numerically to higher proper times τ > 0 + . However, counting with exact analytical expressions for the description of the initial state could simplify to a large extent the phenomenological analyses of data by reducing the amount of numerical work. Our result could be directly applied, for instance, in the multi-parametric fits based on Bayesian JHEP01(2019)073 statistics aimed to determine the medium properties [28]. Also, upon the proper spectral decomposition, they may allow to perform mode-by-mode studies of the hydrodynamical propagation of the initial fluctuations as was proposed in [29,30], or be used to determine the initial eccentricities fluctuations as proposed in [25].
Our paper is organized as follows. In section 2 we introduce a generalization of the MV model with relaxed transversal locality and explicit impact parameter dependence. In this framework we outline the solution to the Yang-Mills equations with one and 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 section 3 we calculate the EMT correlator in the previously presented framework. In section 4 we compute the correlator of two EMTs. Using the results of these two sections, we calculate the covariance of the EMT and show the first orders of its N c -expansion, as well as the strict MV model limit. Our final expression for the EMT covariance is presented in eq. (4.36), the main result of this work. We also compare our results with the previously mentioned computation, performed under the Glasma Graph approximation [22]. Remarkably, throughout this calculation we face a number of outstanding technical challenges such as the calculation of non-trivial projections of the correlator of four Wilson lines in the adjoint representation and the decomposition of correlators of m color sources and n Wilson lines. We analyze these problems in depth on appendices B and C. Finally, in section 5 we discuss the physical implications and phenomenological applications of our result, as well as its role in future works.

The classical approach to gluon production in heavy ion collisions
In the following section we compute the gluon field generated in ultra-relativistic heavy ion collisions. Although this calculation has been done previously in the literature, we deem it convenient to include this preface as it allows us to introduce our modifications to the MV model and establish the notation used in the rest of the paper. We will follow the derivation steps first presented in [31].
In the MV model we represent the high density of small-x gluons carried by each nuclei with gauge fields A µ 1,2 (x) whose dynamics follow from the classical Yang-Mills equations: The source of the fields is a color current J ν,a that represents the flow of large-x valence partons. If we assume a nucleus moving in the positive x 3 direction with a large light-cone momentum p + , we can fix the initial form of J ν,a based on kinematic considerations: where ρ a is the color charge density. The δ ν+ factor indicates that the source generates a color current only in the + direction. This suggests a physical picture of the interaction where the fast 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). As for ρ a , we might factorize its x − dependence by assuming that the JHEP01(2019)073 currents are shaped as a Dirac delta on x − = 0 (last approximate equality). This approximation is motivated by the Lorentz contraction experienced by the relativistic nuclei. However, we choose not to make any assumptions about the longitudinal structure of the nuclei, thus leaving it undetermined for now. In the MV model the calculation of gauge-invariant quantities requires performing an average over the background classical field sources, operation that we denote by . . . . The physical picture for this average emerges from the process-to-process fluctuations of color charges observed inside the nuclei. We take the spatial configuration of color sources ρ a as an stochastic quantity with a certain probability distribution W [ρ] associated as weight function. Thus, observables are obtained as expectation values: . (2.4) Here µ 2 (x − ) is a parameter proportional to the color source number density that acts as the variance of the Gaussian weight. The main implication of the MV model is the following two-point correlator: However, as we intend to apply a more general approach, we choose to relax some of the approximations implied in eq. (2.5) by considering the following, more general, two-point correlator: 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 [32]. 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. This allows to implement the JIMWLK evolution of W [ρ] within the so-called Gaussian truncation [33][34][35][36]. These extensions of the original MV model might prove especially useful in subsequent phenomenological applications. In section 3 we will go into detail about the specific behavior assumed for both h(b ⊥ ) and f (x ⊥ − y ⊥ ). When attempting to describe the medium generated in the collision of two nuclei in the framework outlined above, we encounter a crucial problem: there is no general analytical solution for the Yang-Mills equations with two sources. Thus, we need to turn to either analytical or numerical approximations. A good starting point for these methods is the inner surface of the light-cone, τ = 0 + (i.e. an infinitesimal positive proper time after the collision), as in this region it is possible to find an analytical expression of the gauge fields. In order to do so, it is convenient to divide the space-time into four quadrants as indicated in figure 1. The MV model provides the appropriate framework to calculate the gauge fields that characterize each nuclei before the collision (quadrants 1 and 2). These fields define the boundary conditions for the solution in the future light-cone (quadrant 3). As for τ < 0 the nuclei are located in causally disconnected regions of space-time, we can compute each gauge field independently. Let us take, for instance, a nucleus moving in the positive x 3 direction (which we indicate with the label 1). By solving eq. (2.1) in the light-cone gauge (see e.g. [37] for a detailed resolution), we obtain: which is a non-abelian Weizsäcker-Williams (WW) field. Hereρ is the color charge density in the covariant gauge 1 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. U (x − , x ⊥ ) is defined as a path-ordered

JHEP01(2019)073
exponential: where G(z ⊥ −x ⊥ ) is the Green's function for the 2-dimensional Laplace operator. In the previous expression we show explicitly the definition of the differential operator 1/∇ 2 , which is the notation we adopt to denote a convolution with G(z ⊥ −x ⊥ ). The choice of the integration lower limit x − 0 is arbitrary, with different choices giving us solutions A i connected by residual, 2-dimensional gauge transformations. We shall adopt x − 0 = −∞, which implies that the fields vanish in the remote past (retarded boundary conditions). Likewise, for the nucleus moving in the opposite direction 2 (indicated with the label 2), we have: where: Thus, the total gauge field outside the light-cone reads: (2.14) The gluon field sources vanish everywhere except at the very light-cone (τ = 0), and thus at τ = 0 + the Yang-Mills equations become homogeneous. In order to solve them the following ansatz is proposed in [31]: where we adopted the comoving coordinate system, defined by proper time τ = √ 2x + x − and rapidity η = 1 2 ln(x + /x − ). Substituting the above expressions, the separate components of the homogeneous Yang-Mills equations [D µ , F µν ] = 0 take the following form [38]:

JHEP01(2019)073
This system provides the initial conditions for the τ -evolution of the gluon fields in the future light-cone, be it computed via analytical or numerical methods. In order to relate them to the fields prior to the collision (eq. (2.14)) 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: 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 [39] an analytical approximation based on an expansion of the previous solution in powers of τ is proposed. However, this is out of the scope of the work presented in this paper.

The EMT one-point correlator in the classical approximation
Using eq. (2.15), eq. (2.16) along with the boundary conditions of eq. (2.20), eq. (2.21) we obtain the following expression for the EMT at τ = 0 + [14]: ] µν (in terms of the Cartesian coordinate system), F αβ is the field strength tensor, and 0 (x ⊥ ) is the energy density at proper time τ = 0 + in a point x ⊥ of the transverse plane. Note that the characteristic diagonal structure of this tensor is a feature of the specific proper time at which we are setting our calculation. The ensuing time evolution brings non-trivial off-diagonal corrections that largely modify this initial form, as indicated by the higher order terms of the τ -expansion proposed in [39]. At τ = 0 + , however, the classical approximation yields a remarkably simple, diagonal EMT even at the event-by-event level, prior to the computation of its average over the background fields.
As mentioned earlier, another remarkable aspect of this tensor is the maximum pressure anisotropy denoted by the negative value in the longitudinal direction. The negative pressure slows down the longitudinal expansion of the system, while the remaining components force it to expand in the transverse directions. However, prior to the interpretation of this object we must compute its average over the background fields. We have [15]:

JHEP01(2019)073
As the trace in this expression is performed over color space, in order to compute it we need to expand the color structure of our fields: Here we used the relation between Wilson lines in the fundamental and adjoint represen- Substituting in eq. (3.2) we get: In the last step we factorize the average over color source densitiesρ 1 andρ 2 , since in the MV model we assume the source fluctuations in each nuclei to be independent of each other: Thus, the building block of 0 is the average of two gauge fields evaluated in the same transverse coordinate: Nevertheless, it will prove useful to perform this calculation for different transverse coordinates x ⊥ , y ⊥ and eventually take the limit y ⊥ → x ⊥ : (3.6) The average in the right hand side of this expression contains, for each transverse coordinate, an infinite product ofρ factors: one external and the rest arranged inside the Wilson lines. Since we are assuming that the color sources obey Gaussian statistics, we can apply Wick's theorem, which states that any correlator can be expressed in terms of products of two-point functions. In our particular case, the only nonvanishing terms of the infinite possibilities available are the ones that correspond to a factorization of the external sources from those inside the Wilson lines (see appendix B for a general analysis of the decomposition of the correlator of n Wilson lines and m external sources): (3.7) As the differential operators 1/∇ 2 , ∂ i commute with the average operation, the factor involving the external sources can be calculated via an almost direct application of the JHEP01(2019)073 two-point correlator. In the MV model (eq. (2.5)) this yields a quite simple expression: However, using our generalized version (eq. (2.6)), we have: and then: yielding: In the same spirit than [32], in eq. (3.10) we implicitly make the assumption that the impact parameter profile h(b ⊥ ) introduced earlier is a slowly varying function over lengths of the order of an infrared length scale 1/m, or smaller. Therefore we take 1/m to be an intermediate scale between the inverse saturation scale and the nuclear radius R A : One can think of 1/m as a cut-off that imposes color neutrality at the nucleon size. In addition, in eq. (3.11) we assume that f (x ⊥ − y ⊥ ) behaves in such a way that its Fourier transformf (k ⊥ ) tends to unity in the infrared limit. This requirement, along with the assumed 'slow' behavior for h(b ⊥ ), result in this factor being approximately unaffected by the differential operators in both eq. (3.10) and eq. (3.11) (see appendix A for more details about these assumptions). Substituting in eq. (3.7), we finally get: where the last factor corresponds to the dipole function in the adjoint representation [40]:

JHEP01(2019)073
Here we introduced the factor Note that in eq. (3.15) we applied the same approximation as in eq. (3.10) in order to obtain the factorization of Now, taking the limit y ⊥ → x ⊥ : we will identify functions integrated in the longitudinal direction from −∞ to ∞ by simply omitting their longitudinal dependence) and substituted the following expression: Here the double derivative ∂ 2 L(0 ⊥ ) is a model-dependent constant (see appendix A for details). We apply this result for both nuclei in eq. (3.4), obtaining: whose dependence on the transverse position is a consequence of our generalized MV model approach, where we assume finite nuclei. Note that we label both factors µ 2 and h according to the corresponding nucleus, which potentially allows for the use of different nuclear profiles for target and projectile. We absorb these quantities in the definition of the following momentum scale: which characterizes each colliding nucleus. Performing this substitution we obtain: The EMT two-point correlator in the classical approximation The next step in our calculation is the computation of T µν We start by expanding the product of energy densities: Here we defined the transverse and color structure tensors respectively as: 3) and adopted a shorthand notation for the gluon fields α i,a (x ⊥ ) ≡ α i,a x . As the average operation is performed independently for both nuclei, the building block of 0 (x ⊥ ) 0 (y ⊥ ) reads: This is an extended and more complicated version of eq. (3.6), with twice as many color sources depending on different longitudinal coordinates. The correlations between its different elements result in the following sum: The details of the above decomposition are explained in appendix B. For simplicity we momentarily adopted a shorthand notation that omits the longitudinal coordinate dependence and the differential operators 1/∇ 2 , ∂ i . In the previous expression a major source of JHEP01(2019)073 difficulty stands out: unlike the average featured in eq. (3.6), the one in eq. (4.4) gets nontrivial contributions from correlators connecting external color sources with those arranged inside Wilson lines. We name these 'connected' correlators, and indicate them as . . . c . Based on the diagrammatic rules derived in [41] we are able to compute these contributions and express them in terms of the following function: . Here we recover the notation used in the previous section for the adjoint dipole correlator and extend it to the case of three Wilson lines [40] as: The second longitudinal coordinate in the dependence of both C (2) adj and C adj stands for the lower limit of the integral contained in their definition: We also introduced the adjoint Wilson line quadrupole tensor: The fully connected term (last term of eq. (4.5)) vanishes and thus we are able to write all connected contributions in terms of C ij;kl ab;cd . The remaining contributions result from the factorization of external and internal color sources (first term after the equal sign in eq. (4.5)):

JHEP01(2019)073
This term can be further expanded by application of Wick's theorem, which tells us that the external source correlator breaks down into the following sum of pairwise contractions: Following this decomposition, eq. (4.10) yields three terms that we can address in terms of the 'disconnected' function, which we derive explicitly in the following lines: where b ⊥ = (x ⊥ + y ⊥ )/2. Note that both here and in the connected function eq. (4.6) we substituted the result of eq. (3.12), which implies that we adopt the same assumptions over h(b ⊥ ) and f (x ⊥ − y ⊥ ) as in the previous section. We also made use of the knowledge that eventually all the transverse positions that enter this expression will be either x ⊥ or y ⊥ , which allows us to neglect the corrections to an expansion of h((u ⊥ + u ⊥ )/2) around h((x ⊥ +y ⊥ )/2) (see appendix A for details). This approximation was also taken in eq. (4.7), allowing us to extract λ(w − , b ⊥ ) as a common factor of the sum. Going back to eq. (4.12), note that we introduced the following function: where, as was also the case in the connected function eq. (4.6), we encounter double derivatives of L(x ⊥ − x ⊥ ). From their symmetries and dimension, we can parameterize them as: (4.14) In appendix A we obtain expressions for the coefficients A(r ⊥ ) and B(r ⊥ ) in terms of f (r ⊥ ) and provide an explicit calculation in the specific case of the MV model (where f (r ⊥ ) = δ (2) (r ⊥ )). However, for now we prefer to stay in the most general case and leave them undetermined. In order to solve the integral present in D ij;kl ab;cd we consider separately the region where z − > w − and its complementary. Assuming a certain ordering in the integration variables JHEP01(2019)073 allows us to factorize the Wilson line correlator by applying the locality in rapidity implied in eq. (2.6). For instance, in the region z − > w − (see figure 2): Summing the contributions from each integration region z − > w − and w − > z − we get: Having defined these functions, we can rewrite our building block eq. (4.4) as: Note that, in addition to the fully connected correlator, also the first two partially connected terms of eq. (4.5) vanish (see appendix B).
Remarkably, in both D ij;kl ab;cd and C ij;kl ab;cd we find different projections of the adjoint Wilson line quadrupole eq. (4.9), which is a quite complex object. However, the fact that in our calculation we only deal with two transverse coordinates x ⊥ and y ⊥ yields great simplification in some instances. For example, the first term after the equal sign in eq. (4.17) corresponds to:

JHEP01(2019)073
In this case, the projection of the adjoint Wilson line quadrupole can be obtained in a straightforward way. Writing it explicitly: ) and expanding the first pair of adjoint Wilson lines in terms of fundamental Wilson lines as U ab = 2 Tr U † t a U t b , we get: Therefore: In the two remaining disconnected terms the Wilson lines that share a color index depend on different transverse coordinates, which prevents the Fierz identity from simplifying the expression. In other words, while in eq. (4.18) we have which corresponds to the trivial propagation of an eigenvector in color space, in the other two particular cases of D ij;kl ab;cd we find instead, whose calculation requires expressing δ AC δ BD in terms of the eigenvectors of Q ABCD acbd . This is a highly non-trivial problem that we analyze in depth in appendix C.

JHEP01(2019)073
Substituting the result of eq. (4.23) and solving the double integrals, we obtain: where: adj (x ⊥ , y ⊥ )) 2 exp g 2 Γλ(b ⊥ ) (4.28) adj (x ⊥ , y ⊥ ) . (4.29) As for the C ij;kl ab;cd function, only one particular case enters eq. (4.17): Remarkably, the previous expression contains the propagation of the color vector f ACe f BDe by the adjoint Wilson line quadrupole. In appendix C we show that it is actually an eigenvector of Q ABCD abcd , yielding the following straightforward result: (4.31)

JHEP01(2019)073
Substituting and solving the double integrals, we get: which concludes the calculation of the building block α i,a x α k,c x α i ,a y α k ,c y . The final step consists in explicitly expanding the color contractions between these objects (one for each nucleus) and the transverse and color structure tensors defined earlier: The product of the seven terms corresponding to each nucleus (eq. (4.17)) yields a total of 49 terms, which, by application of the symmetries of the tensors A ik;i k jl;j l and F ac;a c bd;b d , can be reduced to: i k jl;j l F ac;a c bd;b d +A ik;i k jl;l j F ac;a c bd;d b + D ik;i k 1 ac;a c (x ⊥ , x ⊥ , y ⊥ , y ⊥ )C jj ;ll 2 bb ;dd (x ⊥ , y ⊥ , x ⊥ , y ⊥ ) +2 D ii ;kk 1 aa ;cc (x ⊥ , y ⊥ , x ⊥ , y ⊥ )C jj ;ll 2 bb ;dd (x ⊥ , y ⊥ , x ⊥ , y ⊥ ) +2 C ii ;kk 1 aa ;cc (x ⊥ , y ⊥ , x ⊥ , y ⊥ )C jj ;ll 2 bb ;dd (x ⊥ , y ⊥ , x ⊥ , y ⊥ ) × A ik;i k jl;j l F ac;a c bd;b d +A ik;i k lj;j l F ac;a c db;b d +A ik;i k jl;l j F ac;a c bd;d b +A ik;i k lj;l j F ac;a c db;d b It is worth mentioning that the terms resulting of the first contraction after the equal sign in eq. (4.34) are identical to the product of the separate averages of 0 (x ⊥ ) and 0 (y ⊥ ): where we approximated h(x ⊥ ) and h(y ⊥ ) with h(b ⊥ ), as repeatedly done throughout the calculation. Therefore, the result of Cov[ ](τ = 0 + ; x ⊥ , y ⊥ ) = 0 (x ⊥ ) 0 (y ⊥ ) − 0 (x ⊥ ) 0 (y ⊥ ) corresponds to the remaining terms. We use the Mathematica package FeynCalc [42,43] to perform the contractions featured in eq. (4.34). After doing so we arrive at the main

JHEP01(2019)073
result of this work: Cov[ ](τ = 0 + ; x ⊥ , y ⊥ ) where the dependencies have been omitted for readability. The covariance of the full EMT is simply obtained from the previous expression as The factors A(r ⊥ ) and B(r ⊥ ) were introduced in eq. (4.14). Explicit expressions for them in the general case are given in appendix A and in eqs. (4.42), (4.43) below for the specific case of the original MV model. Also, in order to make our final result more compact we have defined: For simplicity, in the previous expressions we also defined the following momentum scale: which is related to the one introduced in section 3 by:

N c -expansion in the MV model
Our generalization of the classical approach yields a few aspects of the previous calculation that had to be left undetermined. For instance, the function f (r ⊥ ) featured in the twopoint correlator of eq. (2.6) introduces some ambiguity in the computation of the double derivative of L(r ⊥ ), which is left in terms of the unknown coefficients A(r ⊥ ) and B(r ⊥ ) (eq. (4.14)). In the particular case of the MV model, where f (r ⊥ ) is taken as a Dirac delta, we are able to compute them as: where K 0 is a modified Bessel function. The mass m is an infrared scale that we introduce to regularize the divergent Green's function G(r ⊥ ). For simplicity we choose m to be the same mass scale introduced earlier in eq. (3.13). The leading behavior in the m → 0 limit is: and B MV , being a constant, yields a negligible correction to this logarithm. In the same limit, the leading behavior of Γ(r ⊥ ) and the product of its derivatives corresponds to the following expressions: and for the saturation scale: (4.47) (See appendix A for a detailed derivation of these expressions). These factors exhibit logarithmic divergences of different nature. While Γ and ∂ i x Γ∂ i y Γ diverge only in the infrared limit m → 0, A and Q 2 s are divergent in both infrared and ultraviolet r → 0 limits. The latter case enters our solution explicitly through the terms multiplied by the factor ∂ 2 L(0 ⊥ ) ≡ −2 lim r→0 A(r ⊥ ). In the end those terms will be the only ones yielding a divergence, as the logarithms stemming from A, Γ and ∂ i x Γ∂ i y Γ are exactly cancelled in eq. (4.36). Therefore, the overall effect of taking the MV limit on the complete result of the energy density covariance only consists in replacing all r-depending coefficients (except the saturation scales) with constants. As this substitution does not yield a significant simplification to the final formula, instead of showing that result we prefer to display the first orders of its N c -expansion. Note that in these expressions we are not taking the complete, strict MV limit, which would imply h(b ⊥ ) = 1; instead, we are only assuming locality in the JHEP01(2019)073 transversal plane. The leading order of the expansion, of order N 0 c , reads: The next term, of order N −2 c , reads: (4.49) In order to give a general idea of the magnitude and analytical features of our solution, on figure 3 we draw these functions in the GBW model as a function of the dimensionless product rQ s for Q s1 = Q s2 . Note that in this limit, as we are ignoring all logarithmic factors, we also have Q s =Q s . The N −2 c term yields a small but noticeable negative correction (see red dashed curve of figure 3). As the next terms are negligible, the first two orders of the N c -expansion provide a neat approximation to the complete result (see right plot of figure 3). Comparing this curve with the N 0 c -order term we notice that the large-N c limit yields a 12.5% error in the r → 0 limit, which is a reasonable approximation. In the rQ s 1 limit our result vanishes following a power-law behavior. The leading term of this limit results from a combination of terms included in the first two orders of the N c -expansion presented above, eq. (4.48) and eq. (4.49): (4.50) Note that this power-law tail is a non-trivial feature of our general result (shown in eq. (4.36)) that is also displayed in the particular case of the MV model. Normalizing the previous result with a single one-point correlator we obtain the following expression: . (4.51) In the opposite limit, r → 0, the covariance tends to: s2 , (4.52) and the normalized covariance: (4.53) In both expressions we applied eq. (4.41) in the last step.

The Glasma Graph approximation
An alternative approach to this calculation is proposed in [22] under the Glasma Graph approximation, whereby it is assumed that the four-point correlation functions of the WW fields characterizing the single nucleus solution of the classical Yang-Mills equations of motion can be factorized into products of two-point correlation functions such that:  Figure 4. LEFT: comparison of the normalized covariance of energy density 0 against r Q s for Q s1 = Q s2 , N c = 3 in the exact analytical approach (blue full curve) and the Glasma Graph approximation (red dashed curve). As a visual aid we also indicate the asymptotic behavior in the IR limit, which is 16/[(N 2 c − 1)r 2 Q 2 s ] (green dot-dashed curve). RIGHT: ratio of exact analytical result to the Glasma Graph result.
This Wick theorem-like decomposition is equivalent to assuming that the WW fields obey Gaussian statistics. This is not generally correct, since the dynamics relating the Gaussianly distributed color sources and the generated gluon fields (encoded in the Yang-Mills equations) are non-linear. However, as we shall see, the Glasma Graph method remains a good approximation in the limit of small transverse separations r → 0. In that limit the dynamics linearize and reduce to two-gluon exchanges, effectively mapping a Gaussian distribution (for the color sources) onto another one (for the WW fields).
We compare the normalized covariance from our result (in the strict MV model) with the one computed according to the decomposition defined in eq. (4.54). As can be seen in figure 4, although both results agree exactly in the UV limit r → 0, in the rest of the spectrum our computation yields a harder curve. Another remarkable difference is that, while our result for the normalized covariance shows a slowly vanishing behavior in the infrared limit, the Glasma Graph approximation yields a much steeper tail: (4.55) The ∼ 1/r 4 decreasing behavior displayed by eq. (4.55) is in clear contrast with the ∼ 1/r 2 asymptotic behavior of our result. This potentially implies much different results and physical interpretations for any observable built from this quantity.

Discussion and outlook
In this paper we provided an analytical expression for the covariance of the EMT characterizing the Glasma state produced in the early stages of an ultra-relativistic heavy ion collision. We performed this calculation in a classical framework based on the CGC effective theory, which we introduced by outlining the solution to the Yang-Mills equations

JHEP01(2019)073
for two nuclei at τ = 0 + . In our approach we assumed a non-local two-point correlator of color source densities with an explicit impact parameter dependence. These modifications were introduced as small but non-negligible deviations from the original MV model. In this framework we obtained a remarkably lengthy, but still simple, formula whose general behavior we analyzed in the GBW model. We find that the first order of the N c -expansion (of order N 0 c ) yields a good approximation of the complete result (especially at large correlation distances), and that only the first correction to this term (of order N −2 c ) yields a non-negligible contribution. Finally, we compare our result with a recent calculation performed in the Glasma Graph approximation [22]. In this work they assume a Gaussian-like decomposition of correlators at the level of the gluon fields rather than the color source densities, thus neglecting many contributions to the medium average. We find that this approximation quickly becomes unsatisfactory as we move out of the UV (r → 0) limit. In fact, the most striking difference emerges in the IR (rQ s 1) limit, in which our result for the normalized covariance vanishes following a inverse square law whereas the Glasma Graph approximation yields a much more rapidly decreasing inverse fourth power. Arguably, these relatively long-range correlations could have an impact on any observable based on the integration of this quantity in the transverse plane, such as the mean square eccentricity fluctuations.
The results presented here provide a first step towards a first-principles computation of the initial conditions of the hydrodynamical expansion of QGP. The subsequent τevolution up to thermalization time, as well as the calculation of observables relevant to QGP phenomenology, will be analyzed in a forthcoming publication.

A Operations involving the 2-D Laplacian Green's function
Throughout the computation of the covariance of T µν 0 we encounter several non-trivial calculations involving the Green's function for the 2-dimensional Laplace operator G(x ⊥ − y ⊥ ). For instance, when computing the correlator of two gluon fields (eq. (3.6)), we find: This expression includes two undetermined functions, h(b ⊥ ) and f (x ⊥ − y ⊥ ), introduced in the two-point correlator (eq. (2.6)) in order to generalize the MV model. However, we do not take these functions as completely general. For h(b ⊥ ), in addition to overall good analytical properties, we assume a slowly varying behavior over lengths of the order of a length scale 1/m or smaller (as proposed in [32]): where we take m as the infrared regulator. We require that: where R A is the nuclear radius. Thus, the interaction distances of interest in our calculation obey r = |x ⊥ − y ⊥ | m −1 . This requirement, as well as the assumed behavior for h(b ⊥ ),

JHEP01(2019)073
yield a significant simplification to eq. (A.1). To see this, we expand h ((z ⊥ + u ⊥ )/2) around b ⊥ = (x ⊥ + y ⊥ )/2: Cutting the expansion at first order, eq. (A.1) yields the following terms: First, we focus on the leading order term: In order to further transform L(x ⊥ − y ⊥ ) we go to momentum space. The Green's function G(x ⊥ − y ⊥ ) admits a simple Fourier representation: which we substitute in L(x ⊥ − y ⊥ ), yielding: In the last step we introduced the inverse Fourier transform of f , defined as: Now we turn to the linear term of the expansion (second term of eq. (A.5)), which we want to compare with h(b ⊥ )L(r ⊥ ). By performing a simple variable change, it can be written as: where v ⊥ = z ⊥ + u ⊥ and w ⊥ = z ⊥ − u ⊥ . Substituting eq. (A.7) and performing some transformations, we get to: The integration in v ⊥ yields a distribution derivative of the Dirac delta function:

JHEP01(2019)073
Substituting this result in eq. (A.11) and integrating by parts, we finally obtain: and thus, eq. (A.5) yields: Here we applied the fact that r i We will take this expression as a good approximation of eq. (A.1). The next step in the calculation of the two gluon field correlator eq. (3.6) is the computation of the double derivative: The reasoning behind the last approximate equality follows from the dimension of L(r ⊥ ), its IR behavior, and the fact that we imposed an infrared cut-off mass scale m. In order to be able to discuss L(r ⊥ ) in the infrared region we need to assume a certain behavior off (q ⊥ ) in this regime. We assumef IR ∼ 1, just like in the MV model, as we do not expect other possible choices of models to differ in that regime. Then, we can safely assume that L ∝ m −2 , which makes the term (∂ i x ∂ j y h)L suppressed with respect to ∂ i x ∂ j y L (a dimensionless object). Also, this takes us to ∂ i L ∝ m −1 , making the terms of the form (∂ j h)(∂ i L) negligible as well. Thus, we are left with the following double derivative: From its symmetries and dimension, the previous expression can be parameterized as: A priori, this decomposition is not possible when r → 0. However, as it is a symmetric object in i, j, we can make a different parameterization in this limit: that we can relate to: Now, taking the limit r → 0: .20) and contracting with δ ij :

JHEP01(2019)073
we have C = − 1 2 ∂ 2 ⊥ L(0 ⊥ ), which is the notation we use in the body of the article (the same as in [44]). We can express these coefficients in terms off (q ⊥ ) by computing the following projections of eq. (A.16) and eq. (A.18): e iq r cos θ q 2 cos(2θ) (A.23) Note that, as lim r→0 A(r ⊥ ) = C and lim r→0 B(r ⊥ ) = 0, this parameterization of ∂ i x ∂ j y L(r ⊥ ) is continuous in r. We can relate C to the factor Γ, defined as: by taking the limit r → 0: where we assumed thatf (q ⊥ ) =f (|q ⊥ |).
The MV model. In the specific case where f (z ⊥ −w ⊥ ) = δ 2 (z ⊥ −w ⊥ ), i.e. the MV model, we havef (q ⊥ ) = 1 and thus we can explicitly compute our coefficients: Both A(r ⊥ ) MV and C MV yield an infrared logarithmic divergence, which we deal with by introducing a regularizing mass in the Fourier representation of G(r ⊥ ):

JHEP01(2019)073
where K 0 is a modified Bessel function. For simplicity we choose m to be the same mass scale introduced earlier in eq. (A.2) (although it could be an unrelated infrared scale). In our calculation we will keep only the leading behavior in the m → 0 limit, which is: (where γ is the Euler constant) and thus: The coefficient C MV corresponds to the UV limit of the previous expression (r → 0): and thus: which also exhibits a logarithmic divergence. As for Γ, we have: The leading behavior of the previous expression in the m → 0 yields:

B Correlators of n Wilson lines and m external color sources
In this appendix we expand on the calculation of the correlator featured in eq. (4.4): which we perform by application of the techniques derived in [41].
(...) Figure 5. Schematic representation of the connected correlator H m,n (see eq. (B.5)) for the particular case featured in our calculation. The oblong shapes represent the sum of all possible contractions between an external source and n Wilson lines, whose correlator is represented as a dark square (see figure 6).
This aspect can be comprised in a redefinition of H 1,n as: where F n denotes the correlator of n Wilson lines. Note that in these formulas the bracketed indices represent a set of n color indices (not indices from 'missing' sources, as in eq. (B.3) and eq. (B.4)). By application of the approximations outlined in appendix A, the previous expression can be rewritten as: where b ⊥ = (x ⊥ + y ⊥ )/2. In this step we have made use of the knowledge that all the transverse positions that enter our calculation are either x ⊥ or y ⊥ . Thus, when expanding h((z ⊥ + w ⊥ )/2) around h(b ⊥ ) in eq. (B.7), the linear term of the expansion yields a correction proportional to a product of the form ( . Whether x ⊥ = x ⊥ or x ⊥ = y ⊥ , this term is suppressed with respect to h(b ⊥ ) according to the assumptions detailed in appendix A. Another difference between our calculation and the one performed in the aforementioned paper is that in the latter the insertion of external sources is assumed to take place at a longitudinal position y − that satisfies b − < y − < a − . However, in our particular case the longitudinal coordinate on which the external color sourceρ a depends is the same as the one of the Wilson line that it is attached to, yielding the following simplification of the previous expression (see figure 6): Having defined all the basic pieces of the calculation of a correlator with m external sources and n Wilson lines, we can go back to our particular case. According to the notation used in [41], the correlator in eq. (B.1) corresponds to F m,n with m = 4, n = 4. By direct application of eq. (B.2): For simplicity we momentarily adopted a shorthand notation that omits the longitudinal coordinate dependence and the differential operators 1/∇ 2 , ∂ i . However, it should be kept in mind that the external sources and Wilson lines that share an index depend on the same longitudinal coordinate. The first term after the equal sign, which corresponds to a complete factorization of external sources and Wilson lines, can be further expanded by application of the Wick's theorem, which tells us that the factor involving the external sources breaks down into the following sum of pairwise contractions: The three terms resulting from this expansion are addressed on section 4, where we derive the 'disconnected' function:

JHEP01(2019)073
with b ⊥ = (x ⊥ + y ⊥ )/2. Again, to get to this expression we anticipate that in our particular case the general transverse positions u ⊥ , u ⊥ , v ⊥ , and v ⊥ will take the values x ⊥ or y ⊥ . Here we also introduce the following function: The disconnected function is used to express part of the outcome of eq. (4.4), which consists in the integration of the correlator eq. (B.1) in the longitudinal direction. Specifically: All remaining terms of eq. (B.10) contain the connected correlator . . . c , which can be computed by application of formulas eq. (B.5) and eq. (B.9). However, our case of interest is somewhat more general than the one covered in these equations, as in our correlator each Wilson line depends on a different longitudinal coordinate. Even though this may seem a source of extra difficulty, it actually yields great simplification. For example, let us take what seems to be the most dreadful term of our calculation, namely the fully connected version of the correlator (last term in eq. (B.10)): As we have four different longitudinal coordinates, in order to compute eq. (B.15) we need to consider all regions of the integration space. 4 For example, applying eq. (B.5) in the region where z − > w − > z − > w − we have: where, according to eq. (B.9), the first factor reads: As is also the case for a single point in a 1-dimensional integral or a line in a 2-dimensional one, the regions where two or more of the coordinates have the same values (for example z − = z − > w − > w − ) yield a negligible contribution. Therefore, we must always consider a certain ordering in our integration variables. The same logic was applied when splitting the double integral featured in the disconnected function D ij;kl ab;cd . which vanishes due to the antisymmetric property of the SU(N c ) structure constants. As we have the same contribution from every region of the integration space, eq. (B.15) yields 0. In order to address the remaining six terms we define the 'connected' function: where b ⊥ = (x ⊥ + y ⊥ )/2. The only nonvanishing contributions to this integral come from the regions where z − > w − > w − and z − > w − > w − , as in the other cases eq. (B.9) introduces a vanishing H 1,1 factor. For the z − > w − > w − region we have (see figure 7): where we have applied the longitudinal locality of eq. (2.6) to factorize the correlator the same way we do in eq. (4.15). We focus on the first connected correlator, which contains JHEP01(2019)073 one external source and three Wilson lines, thus corresponding to: Substituting the result for the three-point adjoint correlator from [40]: (where, again, we approximated the function h as h((x ⊥ + y ⊥ )/2)), we get to: Here, the color factor (the one that cancels (2N 2 c C F ) −1 ) comes from the trace of the product of two structure constants. The remaining correlator, which contains an external source and four adjoint Wilson lines, yields: Substituting this expression and summing the contribution from the z − > w − > w − region the 'connected' function finally yields: (B.24)

JHEP01(2019)073
The last step is the reexponentiation, where we assume that the neglected terms complete the expansion. As an example, let us use this technique to calculate the well-known dipole function: In the last step we have made use of the locality in rapidity of the MV model to factorize the correlator of the remaining Wilson line slices Tr U (x ⊥ )U † (y ⊥ ) (n−1) . By iterating the process, we arrive at: Lastly, we assume that the neglected higher order terms add up to an exponential expression, which reads:

C.2 Diagonalization method
One important shortcoming of the technique described above lies in the fact that there is no unique way in which we can arrange the terms resulting from expanding the Wilson lines. This step becomes more problematic as we increase the number of Wilson lines in the correlator. Nevertheless, we can reduce the inherent arbitrariness of the reexponentiation process by formulating the method as a diagonalization problem. This allows us to systematically account for all incoming and outgoing states of the interaction embodied in the medium average . . . . In the next subsection we will make use of this technique to obtain the behavior of the following adjoint Wilson line quadrupole: under different color projections. However, to illustrate the method we will first reproduce the more general result obtained in [45] for three different transverse coordinates: First, we need to expand the adjoint Wilson lines in a longitudinal position x − n . For the sake of simplicity in the following calculations we will momentarily adopt a shorthand

JHEP01(2019)073
From the previous projections we can write eq. (C.14) in the following form: U aceg bdf h = (U aceg a c e g ) (n−1) T a c e g bdf h (x − n ) = (U aceg a c e g ) (n−1) (1 + M (x − n )) a c e g bdf h , (C. 24) where the M a c e g bdf h matrix is of order 1 in ∆x − . The next step of the method consists in iterating the expansion of the Wilson lines n − 1 times. Doing this (and neglecting terms of order (∆x − ) 2 or higher), we get: It is worth reminding that we are omitting some of the dependencies ofM for simplicity; this tensor also depends on the transverse coordinates,M (x − ; z ⊥ , x ⊥ , y ⊥ ). In order to reproduce the notation of [45], we introduce the following functions: and thus we obtain the following expression forM (hereby correcting typos in the matrix given in [45]):  which we diagonalize using Mathematica: The final step is the reexponentiation of eq. (C.25), which is straightforward for a diagonal matrix: U aceg bdf h= (1 +M d ) aceg bdf h −→ U aceg bdf h= (eM d ) aceg bdf h . (C.31)

JHEP01(2019)073
Here the dot stresses that in order to use this result we need to work in the basis defined by the eigenvectors ofM , which in the (u 1 , u 2 , u 3 , w 1 , w 2 , w 3 ) basis looks like: Remarkably, we have t 2 = u 2 = δ ca δ ge , t 3 = −f can f gen , and t 4 = w 1 = d ean d gcn .

C.3 Projections of the quadrupole
Let us now go back to our particular case: U Aa (x ⊥ )U Bb (x ⊥ )U Cc (y ⊥ )U Dd (y ⊥ ) , (C. 33) which can be obtained from the quadrupole studied in the previous subsection by setting z ⊥ ≡ x ⊥ and x ⊥ = y ⊥ ≡ y ⊥ . This simplifies the above result, as R b = 0 and R d = −R a . In this limit,M d adopts the following form: As part of the calculation of T µν (x ⊥ )T σρ (y ⊥ ) , we need to calculate the following projections of the adjoint Wilson line quadrupole: δ AC δ BD U Aa (x ⊥ )U Bb (x ⊥ )U Cc (y ⊥ )U Dd (y ⊥ ) .