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 ($\langle T^{\mu\nu}(x_{\perp})T^{\sigma\rho}(y_{\perp})\rangle$ and $\langle T^{\mu\nu}(x_{\perp})\rangle$, respectively) at proper time $\tau\!=\!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 ($r^{-2}$ vs $r^{-4}$ power-law decay), indicating that the color screening of the correlations in the transverse plane occurs at distances larger than $1/Q_s$ 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 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 varies 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 is, 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 times 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 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], T µν = diag( , , , − ), where 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 shortrange -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 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 Quark Gluon Plasma 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 is 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 from now on) 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 separation, rQ s 1 where r ≡ |x ⊥ − y ⊥ |: The factors Q s1,2 (r ⊥ ) andQ s1,2 -two definitions of the saturation scales characterizing each nuclei-will be introduced later along with the factor L(0 ⊥ ). Eq. (1.2) is one 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 state 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 two-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 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 bring 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. A few examples where our result could be directly applied are, for instance, the multi-parametric fits based on Bayesian 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 Sec. 2 we introduce a generalization of the MV model with relaxed transversal locality and 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 conditions for the ensuing evolution in the future light-cone. In Sec. 3 we calculate the EMT correlator in the previously presented framework. In Sec. 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.35), 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 nontrivial projections of the correlator of four Wilson lines in the adjoint representation and the decomposition of the correlator of four valence color sources and four Wilson lines. We analyze these problems in depth on appendices B and C. Finally, in Sec. 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: Here 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 the color charge density ρ a (x − , x ⊥ ), we might factorize its x − dependence by assuming that the currents are shaped as a Dirac delta on x − = 0 (last equality). This approximation is motivated by the Lorentz contraction applied to the relativistic nuclei. However, we choose not to make any assumptions about the longitudinal structure of 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 a 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 result of the MV model is the following two-point correlator: which establishes that only interactions with the background field that are local in color, rapidity and transverse position yield a nonvanishing contribution to observables. 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 correlation function of the sources: 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 Fig. 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's 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: 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 exponential: 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, two-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: 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 Eq. (2.12), the separate components of the 2 We work in a specific gauge that acts as a sort of 'mix' of the light-cone gauges of both nuclei: the Fock-Schwinger gauge, defined by the condition (x + A − + x − A + )/τ = 0. Note that, as the separate fields of each nuclei already satisfy this condition, the Fock-Schwinger representation does not introduce any physical assumption that is not already present in the single nucleus characterization.
homogeneous Yang-Mills equations [D µ , F µν ] = 0 take the following form [38]: 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 it to the fields prior to the collision (Eq. (2.11)) 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.12) along with the boundary conditions of Eq. (2.14) we obtain the following expression for the EMT at τ = 0 + : where t µν ≡ diag(1, 1, 1, −1) (in terms of the Cartesian coordinate system) and 0 (x ⊥ ) is the energy density at proper time τ = 0 + in a point x ⊥ of the transverse plane. As mentioned earlier, a remarkable aspect of this tensor is the maximum pressure anisotropy denoted by the negative 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 T µν 0 (x ⊥ ) = 0 (x ⊥ ) ×t µν , with: 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 representa- 2) we get: In the last equality we factorize the average over color source densitiesρ 1 andρ 2 , since the sources fluctuations in each nuclei are independent of each other: Thus, the building block of 0 is the average of two gauge fields evaluated in the same transverse coordinate: α i,a (x ⊥ )α j,b (x ⊥ ) . 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ρ a factors: one external and the rest arranged inside the Wilson lines. Since we are assuming Gaussian statistics of the color sources 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 term of the infinite possibilities available is the one where the external sources are factorized from those inside the Wilson lines (see Appendix B for a general analysis of the decomposition of the correlator of m Wilson lines and n 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 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 m −1 , or smaller. Therefore we take m −1 to be an intermediate scale between the inverse saturation scale and the nuclear radius R A : As discussed below, one can think of m −1 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]: 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 structure models 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 difficulty stands out: unlike the average featured in Eq. (3.6), the one in Eq. (4.4) gets nonvanishing 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: We also introduced the adjoint Wilson line quadrupole tensor: It turns out that the fully connected term, i.e. the last term of Eq. (4.5), vanishes (see Appendix B for the explicit calculation), so 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 equality sign in Eq. (4.5)): 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: (4.10) Following this decomposition, Eq. (4.9) yields three terms that we can address in terms of the 'disconnected' function, which we derive explicitly in the following lines: 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.11), 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: In Appendix A we obtain expressions for the coefficients A(r ⊥ ) and B(r ⊥ ) in terms of the transversal f (r ⊥ ) function 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 allows us to factorize the Wilson line correlator by applying the locality in rapidity of Eq. (2.6). For instance, in the region z − > w − : (4.14) 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.8), 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 equality sign in Eq. (4.16) corresponds to: 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: Now, applying the Fierz identity t a ij t a kl = 1 2 (δ il δ jk − 1 Nc δ ij δ kl ): Therefore: (4.20) 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.17) we have which corresponds to the trivial propagation of an eigenvector in color space, in the other particular cases of D ij;kl ab;cd we find instead, whose calculation requires expressing δ AC δ BD in terms of eigenvectors of Q ABCD acbd . This is a highly non-trivial problem that we analyze in depth in Appendix C. Substituting the result of Eq. (4.22) and solving the double integrals, we obtain: where: adj (x ⊥ , y ⊥ )) . (4.28) As for the C ij;kl ab;cd function, only one particular case enters Eq. (4.16): Remarkably, the previous expression contains the propagation of the color vector f ACe f DBe 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: 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 following 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.16)) 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: +D ii ;kk 1 aa ;cc (x ⊥ , y ⊥ , x ⊥ , y ⊥ )D 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 jl;l j F ac;a c bd;d b It is worth mentioning that the terms resulting of the first contraction after the equality sign in Eq. (4.33) are identical to the product of the separate averages of 0 (x ⊥ ) and 0 (y ⊥ ): where we approximate h(x ⊥ ) and h(y ⊥ ) with h(b ⊥ ), as repeatedly done throughout the calculation. Therefore, the result of Cov[ ](τ = 0 corresponds to the last six lines of Eq. (4.33). We use the Mathematica package FeynCalc [42,43] to perform the contractions appearing in Eq. (4.33). After doing so we arrive at the main result of this work for the energy density covariance: The covariance of the full EMT is simply obtained from the previous expression as The factors A(r ⊥ ) and B(r ⊥ ) (A and B in the above expression for simplicity) were introduced in Eq. (4.13). Explicit expressions for them in the general case are given in Appendix A and in Eqs. (4.41), (4.42) below for the specific case of the original MV model. Also, 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 linked to the one introduced in Sec. 3 by:

N c -expansion in the MV model
Since the full result for the covariance of the EMT provided in Eq. (4.35) is somewhat lengthy, in this and the next sections we discuss a few simplifying limits in the context of the original 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 two-point 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.13)). 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 (see Appendix A for details). The gluon 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: 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 infrared logarithms stemming from A, Γ and ∂ i x Γ∂ i y Γ are exactly cancelled. Therefore, the overall effect of taking the MV limit on the complete result of the energy density covariance (shown in Eq. (4.35)) 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 throughout this section we are not taking the complete, strict MV limit, which would imply h(b ⊥ ) = 1; instead, we are only assuming locality in the transversal plane. The leading order of the expansion, of order N 0 c , reads: The next term, of order N −2 c , reads: (4.48) In order to give a general idea of the magnitude and analytical features of our solution, on Fig. 3 we draw these functions in the GBW limit 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 Fig. 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 Fig. 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.47) and Eq. (4.48): (4.49) Note that this power-law tail is a non-trivial feature of our general result (shown in Eq.(4.35)) 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.50) In the opposite limit, r → 0, the covariance tends to: s2 , (4.51) and the normalized covariance: (4.52) In the last equality of these expressions we apply Eq. (4.40).  Figure 4. LEFT: Comparison of the normalized covariance of energy density 0 against r Q s for Q s1 = Q s2 in the GBW model (blue full curve) and the Glasma Graph-inspired approximation (red dashed curve). As a visual aid we also indicate the asymptotic behavior in the IR limit, which is

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: (4.53) 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.53). As can be seen in Fig. 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-inspired approximation yields a much steeper tail: (4.54) The ∼ 1/r 4 decreasing behavior displayed by Eq. (4.54) 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 a 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 for two nuclei at τ = 0 + . In our approach we assumed an explicit impact parameter dependence in the average of two color source densities, as well as a generalization of the transverse profile of the interaction. 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 (r Q 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 should be understood as a starting point towards a firstprinciples 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. Specifically, we suggest that our results could be straightforwardly applied to the calculation of eccentricity fluctuations through the expressions presented in [25].

A Operations involving the 2-D Laplacian Green's function
Throughout the computation of the covariance of T µν 0 we encounter several nontrivial 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 m −1 or smaller: 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 ⊥ ), yield a significant simplification to Eq. (A.1). To see this, we expand h Cutting the expansion at first order, Eq. (A.1) yields the following terms: where b ⊥ = (z ⊥ + u ⊥ )/2. 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 line we introduced the inverse Fourier transform of f , defined as: Now we turn to the linear term of the expansion (second line in 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: 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 application of the following derivatives: 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 don't 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: and contracting with δ ij : 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.

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 gluon mass in the Fourier representation of G(r ⊥ ): 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): which also exhibits a logarithmic divergence.

B Correlators of m Wilson lines and n 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]. In one of the appendices of this work they analyze the general case of the correlator of m Wilson lines and n color charge densities 3 : is the correlator of m−1 color charge densities. In the notation adopted here, the indices corresponding to sources that are 'missing' from the correlators are indicated by brackets {...}. We also have: which is the 'connected' correlator of n Wilson lines with j insertions of external sources at the positions J 1 , J 2 , ... J j (with J 1 < J 2 < ... < J j ). This is a special kind of correlator that does not include contractions between color sources outside the Wilson lines. Therefore, when computing it, any of these external sources can only be linked to those arranged inside Wilson lines. An iterative formula is derived for this object: where H 1,n is the basic building block of the 'connected' correlators, having only one external source being linked to those inside n Wilson lines. Applying our generalized version of the MV model (embodied in the two-point correlator Eq. (2.6)), it yields the following expression: 2) is derived for the case where m is even. In [41] the odd m formula is also provided.
However, a fundamental difference between our calculation and the one featured in [41] is that in our case the external sources are affected by the differential operators 1/∇ 2 and ∂ i . 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 this formula 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: 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 Fig. 5): (B.9) 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. 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 equality 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: 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: (B.14) 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's 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)): region where z − > w − > z − > w − we have: where, according to Eq. (B.9), the first factor reads: 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: where we have applied the longitudinal locality of Eq. (2.6) to factorize the correlator the same way we do in Eq. (4.14). We focus on the first connected correlator, which contains 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. For the remaining correlator, which contains an external source and four adjoint Wilson lines: Unsurprisingly, both correlators tend to the dipole in the limit of only two transverse coordinates. We are left with a product of three dipole functions that combine as: and therefore: Solving the double integral, we obtain: adj (x ⊥ , y ⊥ ) . Before addressing the problem of the adjoint Wilson line quadrupole we will briefly describe and apply a general method for the computation of Wilson line correlators. This technique, first applied in [45], is based in the discretization of the x − direction into n layers of length ∆x − . Due to the properties of path ordered exponentials, this leads to the factorization of the Wilson line into a product of n independent contributions from each zone: assuming that ∆x − is equal to or shorter than the correlation length of the gluon field fluctuations. This assumption is trivially satisfied in the MV model (and also in our generalized version), where interactions are local in rapidity, allowing us to take the limit ∆x − → 0. As a first step we expand one of these n factors to order g 2 : is the only nontrivial component of the gluon field expressed in the covariant gauge (defined by the condition ∂ µÃ µ = 0). Note that in the g 2 -order term of Eq. (C.2) we already applied the two-point correlator, whose discretized version reads: We iterate this process n times neglecting terms of order (∆x − ) 2 or higher. Then, we rearrange the resulting terms in the form of the first orders of an expanded exponential. The last step is the reexponentiation, where we assume that the neglected terms complete the expansion. As an example, let's use this technique to calculate the well-known dipole function: (C.5) 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 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 ... . We make use of this technique to obtain the behavior of the following adjoint Wilson line quadrupole: under different color projections. However, 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 notation (similar to the one used in [45]). We absorb the g∆x − factor in the definition of our fields: which yields the following two-point function: where, due to the discretization of the rapidity range, the Kronecker delta δ x − y − now occupies the place of the Dirac delta. For simplicity we also introduced: (C.12) Using this notation the expansion to second order in g of the adjoint Wilson line looks like: Performing this expansion for every Wilson line in Eq. (C.9) and neglecting terms of order (∆x − ) 2 or higher we get: U ab (z ⊥ )U cd (z ⊥ )U ef (x ⊥ )U gh (y ⊥ ) = U aa (z ⊥ )U cc (z ⊥ )U ee (x ⊥ )U gg (y ⊥ ) (n−1) × δ a b δ c d δ e f δ g h 1 − N c 2 (2B z + B x + B y ) + δ a b δ c d f e mf f g mh B xy +δ a b δ e f f c md f g mh B zy + δ a b δ g h f e mf f c md B zx + δ e f δ c d f a mb f g mh B zy +δ g h δ c d f e mf f a mb B zx + δ e f δ g h f a mb f c md B z . (C.14) We express the previous lines as a matrix equation: U aceg bdf h = (U aceg a c e g ) (n−1) T a c e g bdf h , for which we introduce the following color vector basis: u 1 = δ ea δ gc , u 2 = δ ca δ ge , u 3 = δ ga δ ec , w 1 = d eam d gcm , w 2 = d cam d gem , w 3 = d gam d ecm , z 1 = d eam f gcm , z 2 = d cam f gem , z 3 = d gam f ecm .
(C.15) T w 1 = w 1 1 − g 2 N c which we diagonalize: The final step is the reexponentiation of Eq. (C.25), which is straightforward for a diagonal matrix: 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's now go back to our particular case: which can be obtained from the quadrupole calculated in the previous section 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: The first of them corresponds to the propagation of the eigenvector t 3 , which is straightforward to compute: The case of Eq. (C.36) corresponds to the propagation of u 1 , which is not an eigenvector and thus requires that we express it in terms of the t i set: and then: where for simplicity we omitted the dependencies ofλ. Expanding the eigenvectors in terms of our original basis we obtain: = δ ac δ bd N 2 c − 4 2N 2 c e −g 2 NcΓλ + 2 N 2 c e −g 2 Nc 2 Γλ + N c + 2 4N c e −g 2 (Nc+1)Γλ + N c − 2 4N c e −g 2 (Nc−1)Γλ Note that in order to use these results in the calculation of D ij;kl ab;cd (x ⊥ , y ⊥ , x ⊥ , y ⊥ ) and C ij;kl ab;cd (x ⊥ , y ⊥ , x ⊥ , y ⊥ ) we need to permute the indices b and c (see Eq. (4.23)).