Initial correlations of the Glasma energy-momentum tensor

JHEP01(2019)073 Published for SISSA by Springer Received:October 13, 2018 Revised:December 12, 2018 Accepted:December 19, 2018 Published:January 8, 2019 Initial correlations of the Glasma energy-momentum tensor Javier L. Albacete,aPablo Guerrero-Rodr´ıgueza,b and Cyrille Marquetb aCAFPE & Dpto. de F´ısica Te´orica y del Cosmos, Universidad de Granada, E-18071 Campus de Fuentenueva, Granada, Spain bCentre de Physique Th´eorique, ´ Ecole Polytechnique, CNRS, Universit´e Paris-Saclay, F-91128 Palaiseau, France E-mail: [email protected],[email protected],[email protected] Abstract: 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 (hTµν(x⊥)Tσρ(y⊥)iand hTµν(x⊥)i, 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/r2vs. 1/r4power-law decay), indicating that the color screening of the correlations in the transverse plane occurs at distances larger than 1/Qsby a logarithmic factor sensitive to the infrared. In the Glasma flux tube picture, this implies that the color domains are larger than originally estimated. Keywords: Heavy Ion Phenomenology ArXiv ePrint: 1808.00795 Open Access,c The Authors. Article funded by SCOAP3.https://doi.org/10.1007/JHEP01(2019)073 JHEP01(2019)073 Contents 1 Introduction 1 2 The classical approach to gluon production in heavy ion collisions 5 3 The EMT one-point correlator in the classical approximation 9 4 The EMT two-point correlator in the classical approximation 13 4.1 Nc-expansion in the MV model 21 4.2 The Glasma Graph approximation 23 5 Discussion and outlook 24 A Operations involving the 2-D Laplacian Green’s function 25 B Correlators of nWilson lines and mexternal color sources 29 C The correlator of four Wilson lines in the adjoint representation 37 C.1 Reexponentiation method 37 C.2 Diagonalization method 38 C.3 Projections of the quadrupole 42 1 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 – 1 – 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-xgluons 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–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–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µν 0= [diag(0, 0, 0,−0)]µν , where 0is 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/Qs, 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 – 2 – 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: Cov[Tµν ](τ= 0+;x⊥, y⊥)≡ hTµν 0(x⊥)Tσρ 0(y⊥)i−hTµν 0(x⊥)ihTσρ 0(y⊥)i,(1.1) 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, hρa(x− , x⊥)ρb(y− , y⊥)i, 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, rQs1 with r≡ |x⊥−y⊥|: lim rQs1Cov[T00 MV ](0+;x⊥, y⊥) = 2N2 c−14π ∂2L(0⊥)2¯ Q4 s 1 Q2 s 2 +¯ Q4 s 2 Q2 s 1  g4N2 cr2.(1.2) The factors Qs 1 , 2 (r⊥, b⊥) and ¯ Qs 1 , 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 – 3 – 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/Qsare created. In our result, although the transverse correlation length is parametrically of the order of 1/Qs, the correlations decrease only according to a 1/r2power-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/Qs, which is what happens under the Glasma Graph approximation [22] (that features a 1/r4 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/Qs, as it features a logarithmic enhancement ln(Qs/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 1 hT00 MV (0+, x⊥)iZd2r⊥Cov[T00 MV ](0+;x⊥, x⊥−r⊥)≃¯ Q2 s 1 (x⊥) + ¯ Q2 s 2 (x⊥) αsNc .(1.3) 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 Rd2x⊥T00 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π/(N2 c−1))[ ¯ Q2 s 1 (x⊥)+ ¯ Q2 s 2 (x⊥)]/Rd2x⊥¯ Q2 s 1 (x⊥)¯ Q2 s 2 (x⊥). 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/(N2 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 – 4 – 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 2we 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 3we calculate the EMT correlator in the previously presented framework. In section 4we 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 Nc-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 mcolor sources and nWilson lines. We analyze these problems in depth on appendices Band C. Finally, in section 5we discuss the physical implications and phenomenological applications of our result, as well as its role in future works. 2 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-xgluons carried by each nuclei with gauge fields Aµ 1,2 (x) whose dynamics follow from the classical Yang-Mills equations: [Dµ, Fµν] = Jν=Jν,a ta .(2.1) The source of the fields is a color current Jν,a that represents the flow of large-xvalence 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: Jν,a(x−, x⊥) = δν+ρa(x−, x⊥)≈δν+δ(x−)ρa(x⊥),(2.2) where ρais 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 – 5 – 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 h. . .i. 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: hO[ρ]i=1 NZ[dρ]W[ρ]O[ρ],(2.3) where Nis a normalization constant equal to R[dρ]W[ρ]. The main assumption adopted in the MV model is that in nuclei with large mass numbers the valence partons that enter eq. (2.1) through ρaemerge from a large number of separate nucleons and therefore are uncorrelated. Thus, invoking the central limit theorem, this model approximates W[ρ] with a Gaussian distribution: hO[ρ]i MV =R[dρ] exp n−Rdx−d2x⊥1 2µ2(x−)Tr ρ2(x−, x⊥)oO[ρ] R[dρ] exp n−Rdx−d2x⊥1 2µ2(x−)Tr [ρ2(x−, x⊥)]o.(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: hρa(x−, x⊥)ρb(y−, y⊥)i MV =µ2(x−)δabδ(x−−y−)δ2(x⊥−y⊥).(2.5) 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: hρa(x−, x⊥)ρb(y−, y⊥)i=µ2(x−)h(b⊥)δabδ(x−−y−)f(x⊥−y⊥) ≡λ(x−, b⊥)δabδ(x−−y−)f(x⊥−y⊥),(2.6) 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–36]. These extensions of the original MV model might prove especially useful in subsequent phenomenological applications. In section 3we 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 – 6 – JHEP01(2019)073 t x3 x+ x  1 2 0 3 Figure 1. Space-time diagram of the collision of two ultra-relativistic heavy ion nuclei. The two diagonal lines represent the trajectory of the nuclei. The points below them (quadrant 0) represent a region where the projectiles have not yet arrived. By choosing the gauge fields to vanish in the remote past, in this region we have Aµ= 0. As quadrants 1 and 2 represent regions where only one of the nuclei has arrived, the dynamics of the gauge fields there are described by the Yang-Mills equations with a single source. However, in quadrant 3 we need to take into account both sources. 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: A± 1 = 0 (2.7) Ai 1 =θ(x−) Z∞ −∞ dz−U† 1 (z−, x⊥)∂i˜ρ 1 (z−, x⊥) ∇2U 1 (z−, x⊥)≡θ(x−)αi 1 (x⊥),(2.8) which is a non-abelian Weizs¨acker-Williams (WW) field. Here ˜ρis the color charge density in the covariant gauge1and Uis the Wilson line, an SU(Nc) 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 1Providing that we average gauge invariant observables, the specific gauge in which we work does not affect the result of hOi, as both the Gaussian weight W[ρ] and the functional measure [dρ] are gauge invariant objects. – 7 – JHEP01(2019)073 exponential: U 1 (x−, x⊥)=P−exp (−igZx− x− 0 dz−1 ∇2˜ρ 1 (z−, x⊥)) = P−exp (−igZx− −∞ dz− Zdz2 ⊥G(z⊥−x⊥)˜ρ 1 (z−, z⊥)),(2.9) 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− 0is arbitrary, with different choices giving us solutions Aiconnected 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 direction2(indicated with the label 2), we have: A± 2 = 0 (2.10) Ai 2 =θ(x+) Z∞ −∞ dz+U† 2 (z+, x⊥)∂i˜ρ 2 (z+, x⊥) ∇2U 2 (z+, x⊥)≡θ(x+)αi 2 (x⊥),(2.11) where: U 2 (x+, x⊥)=P+exp (−igZx+ −∞ dz+ Zdz2 ⊥G(z⊥−x⊥)˜ρ 2 (z+, z⊥)).(2.12) Thus, the total gauge field outside the light-cone reads: A±= 0 (2.13) Ai=θ(x−)θ(−x+)αi 1 (x⊥) + θ(x+)θ(−x−)αi 2 (x⊥).(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]: A±=±x±α(τ= 0+ , x⊥) (2.15) Ai=αi(τ= 0+ , x⊥),(2.16) where we adopted the comoving coordinate system, defined by proper time τ=√2x+x− and rapidity η=1 2ln(x+/x−). Substituting the above expressions, the separate components of the homogeneous Yang-Mills equations [Dµ, Fµν] = 0 take the following form [38]: ν=τ−→ igτ [α, ∂τα]−1 τDi, ∂ταi= 0 (2.17) ν=η−→ 1 τ∂τ 1 τ∂τ(τ2α)−Di,Di, α= 0 (2.18) ν=j−→ 1 τ∂τ(τ∂ταj)−igτ2α, Dj, α−Di, Fij= 0.(2.19) 2We 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. – 8 – 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: D˜ρi,e x˜ρk,f x˜ρi0 ,e0 y˜ρk0 ,f0 yE=h˜ρi,e x˜ρk,f xih˜ρi0 ,e0 y˜ρk0 ,f0 yi+h˜ρi,e x˜ρi0 ,e0 yih˜ρk,f x˜ρk0 ,f0 yi+h˜ρi,e x˜ρk0 ,f0 yih˜ρk,f x˜ρi0 ,e0 yi. (4.11) 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: Dij;kl ab;cd(u⊥, u0 ⊥, v⊥, v0 ⊥) =Z∞ −∞ dz−dz−0dw−dw−0*∂i˜ρa0 (z−, u⊥) ∇2 ∂j˜ρb0 (z−0, u0 ⊥) ∇2+*∂k˜ρc0 (w−, v⊥) ∇2 ∂l˜ρd0 (w−0, v0 ⊥) ∇2+ ×DUa0a(z−, u⊥)Ub0b(z−0, u0 ⊥)Uc0c(w−, v⊥)Ud0d(w−0, v0 ⊥) E =Z∞ −∞ dz− dz−0 dw− dw−0δa0b0λ(z−, b⊥)δ(z−−z−0 )∂i u∂j u0L(u⊥−u0 ⊥)δc0d0 λ(w−, b⊥) ×δ(w−−w−0)∂k v∂l v0L(v⊥−v0 ⊥)DUa0a(z−, u⊥)Ub0b(z−0, u0 ⊥)Uc0c(w−, v⊥)Ud0d(w−0, v0 ⊥) E =Tij;kl (u⊥, u0 ⊥, v⊥, v0 ⊥) Z∞ −∞ dz− dw− λ(z−, b⊥)λ(w−, b⊥) ×DUa0a(z−, u⊥)Ua0b(z−, u0 ⊥)Uc0c(w−, v⊥)Uc0d(w−, v0 ⊥)E,(4.12) 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⊥+u0 ⊥)/2) around h((x⊥+y⊥)/2) (see appendix Afor 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: Tij;kl(u⊥, u0 ⊥, v⊥, v0 ⊥)≡∂i u∂j u0L(u⊥−u0 ⊥)∂k v∂l v0L(v⊥−v0 ⊥),(4.13) where, as was also the case in the connected function eq. (4.6), we encounter double derivatives of L(x⊥−x0 ⊥). From their symmetries and dimension, we can parameterize them as: ∂i x∂j yL(r⊥) = A(r⊥)δij +B(r⊥)δij 2−rirj r2.(4.14) In appendix Awe 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 Dij;kl ab;cd we consider separately the region where z−>w−and its complementary. Assuming a certain ordering in the integration variables – 15 – JHEP01(2019)073 ⌦ ↵ 1 w  z  a0 b0 c0 d0 a b c d = ⌦ ↵ a0 b0 A B w  z  ⌦ ↵ w  1 A B c0 d0 a b c d ⌦ ↵ 1 w  z  a0 c0 a b c d = a0 c0 ⌦ ↵ a0 A B w  z  ⌦ ↵ w  1 A B c0 a b c d a0 c0 Figure 2. Schematic representation of the correlator factorization performed in eq. (4.15). 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): DUa0a(z−, u⊥)Ua0b(z−, u0 ⊥)Uc0c(w−, v⊥)Uc0d(w−, v0 ⊥)E =DUa0A(z−, w−;u⊥)Ua0B(z−, w−;u0 ⊥)E ×DUAa(w−, u⊥)UBb(w−, u0 ⊥)Uc0c(w−, v⊥)Uc0d(w−, v0 ⊥)E =C(2) adj(z−, w−;u⊥, u0 ⊥)DUAa(w−, u⊥)UAb(w−, u0 ⊥)Uc0c(w−, v⊥)Uc0d(w−, v0 ⊥)E =C(2) adj(z−, w−;u⊥, u0 ⊥)QAAc0c0 abcd (w−;u⊥, u0 ⊥, v⊥, v0 ⊥).(4.15) Summing the contributions from each integration region z−>w−and w−>z−we get: Dij;kl ab;cd(u⊥, u0 ⊥, v⊥, v0 ⊥) = Tij;kl (u⊥, u0 ⊥, v⊥, v0 ⊥) Z∞ −∞ dz−Zz− −∞ dw− λ(z−, b⊥)λ(w−, b⊥) ×C(2) adj(z−, w−;u⊥, u0 ⊥) + C(2) adj(z−, w−;v⊥, v0 ⊥)δABδCDQABCD abcd (w−;u⊥, u0 ⊥, v⊥, v0 ⊥). (4.16) Having defined these functions, we can rewrite our building block eq. (4.4) as: hαi a(x⊥)αk c(x⊥)αi0a0(y⊥)αk0c0(y⊥)i =Dik;i0k0 ac;a0c0(x⊥, x⊥, y⊥, y⊥) + Dii0;kk0 aa0;cc0(x⊥, y⊥, x⊥, y⊥) +Dik0;ki0 ac0;ca0(x⊥, y⊥, x⊥, y⊥) + Cii0;kk0 aa0;cc0(x⊥, y⊥, x⊥, y⊥) + Cik0;ki0 ac0;ca0(x⊥, y⊥, x⊥, y⊥) +Ckk0;ii0 cc0;aa0(x⊥, y⊥, x⊥, y⊥) + Cki0;ik0 ca0;ac0(x⊥, y⊥, x⊥, y⊥).(4.17) 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 Dij;kl ab;cd and Cij;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: Dik;i0k0 ac;a0c0(x⊥, x⊥, y⊥, y⊥)=2Tik;i0k0 (x⊥, x⊥, y⊥, y⊥) Z∞ −∞ dz−Zz− −∞ dw− λ(z−, b⊥)λ(w−, b⊥) QAACC aca0c0(w−;x⊥, x⊥, y⊥, y⊥).(4.18) – 16 – JHEP01(2019)073 In this case, the projection of the adjoint Wilson line quadrupole can be obtained in a straightforward way. Writing it explicitly: QAACC aca0c0(w−;x⊥, x⊥, y⊥, y⊥) = DUAa(w−, x⊥)UAc(w−, x⊥)UCa0(w−, y⊥)UCc0(w−, y⊥)E (4.19) and expanding the first pair of adjoint Wilson lines in terms of fundamental Wilson lines as Uab = 2 Tr U†taUtb, we get: UAaUAc = 4 U† ijtA jkUklta liU† i0j0tA j0k0Uk0l0tc l0i0. Now, applying the Fierz identity ta ijta kl =1 2(δilδjk −1 Ncδijδkl) = 2 δjk0δkj0−1 Nc δjkδj0k0U† ijUklta liU† i0j0Uk0l0tc l0i0 = 2 U† ijUklU† i0kUjl0−1 Nc U† ijUjlU† i0j0Uj0l0ta litc l0i0 = 2 δil0δi0l−1 Nc δilδi0l0ta litc l0i0= 2 Tr{tatc}− 1 Nc Tr{ta}Tr{tc}=δac.(4.20) Therefore: Dik;i0k0 ac;a0c0(x⊥, x⊥, y⊥, y⊥)=2Tik;i0k0 (x⊥, x⊥, y⊥, y⊥)δacδa0c0 Z∞ −∞ dz− Zz− −∞ dw− λ(z−, b⊥)λ(w−, b⊥) =Tik;i0k0 (x⊥, x⊥, y⊥, y⊥)δacδa0c0¯ λ2(b⊥) =1 4δikδi0k0 ∂2L(0⊥)2δacδa0c0¯ λ2(b⊥).(4.21) 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 δABδCDQABCD abcd (w−;x⊥, x⊥, y⊥, y⊥) = δabδcd,(4.22) which corresponds to the trivial propagation of an eigenvector in color space, in the other two particular cases of Dij;kl ab;cd we find δACδBDQABCD acbd (w−;x⊥, x⊥, y⊥, y⊥) (4.23) instead, whose calculation requires expressing δAC δBD in terms of the eigenvectors of QABCD acbd . This is a highly non-trivial problem that we analyze in depth in appendix C. – 17 – JHEP01(2019)073 Substituting the result of eq. (4.23) and solving the double integrals, we obtain: Dij;kl ab;cd(x⊥, y⊥, x⊥, y⊥) = 2δabδcdN2 c−4 2N2 c f 1 +2 N2 c f 2 +Nc+ 2 4Nc f 3 +Nc−2 4Nc f 4  +δacδbd1 N2 c−1f 5 −Nc+ 2 2Nc(Nc+ 1)f 3 +Nc−2 2Nc(Nc−1)f 4  +δadδbc−N2 c−4 2N2 c f 1 −2 N2 c f 2 +Nc+ 2 4Nc f 3 +Nc−2 4Nc f 4  +dabmdcdm−1 Nc f 1 +1 Nc f 2 +1 4f 3 −1 4f 4 +dadmdcbm1 Nc f 1 −1 Nc f 2 +1 4f 3 −1 4f 4  +dacmdbdmNc N2 c−4f 2 −Nc+ 4 4(Nc+ 2)f 3 +Nc−4 4(Nc−2)f 4 Tij;kl (x⊥, y⊥, x⊥, y⊥),(4.24) where: f 1 =2 (Ncg2Γ)21−C(2) adj(x⊥, y⊥)2(4.25) f 2 =2 Ncg2Γ2 Ncg2Γ1−C(2) adj(x⊥, y⊥)−¯ λ(b⊥)C(2) adj(x⊥, y⊥)(4.26) f 3 =4 Nc(Nc+ 2)g4Γ21−C(2) adj(x⊥, y⊥) −2 (Nc+ 2)(Nc+ 1)g4Γ21−(C(2) adj(x⊥, y⊥))2exp −g2Γ¯ λ(b⊥)(4.27) f 4 =4 Nc(Nc−2)g4Γ21−C(2) adj(x⊥, y⊥) −2 (Nc−2)(Nc−1)g4Γ21−(C(2) adj(x⊥, y⊥))2exp g2Γ¯ λ(b⊥)(4.28) f 5 =2 Ncg2Γ¯ λ(b⊥)−2 Ncg2Γ1−C(2) adj(x⊥, y⊥).(4.29) As for the Cij;kl ab;cd function, only one particular case enters eq. (4.17): Cij;kl ab;cd(x⊥, y⊥, x⊥, y⊥) =g2 2∂i x∂j yL(x⊥−y⊥)h3(b⊥)∂k xΓ(x⊥−y⊥)∂l yΓ(y⊥−x⊥) ×Z∞ −∞ dz− Zz− −∞ dw− Zw− −∞ dw−0µ2(z−)µ2(w−)µ2(w−0)C(2) adj(z−, w−0;x⊥, y⊥) ×fACefBDeQABCD abcd (w−0;x⊥, y⊥, x⊥, y⊥).(4.30) Remarkably, the previous expression contains the propagation of the color vector fACefBDe by the adjoint Wilson line quadrupole. In appendix Cwe show that it is actually an eigenvector of QABCD abcd , yielding the following straightforward result: fACefBDeQABCD abcd (w−0;x⊥, y⊥, x⊥, y⊥) = fABefCDeQABCD acbd (w−0;x⊥, x⊥, y⊥, y⊥) =facefbdeC(2) adj(w−0;x⊥, y⊥).(4.31) – 18 – JHEP01(2019)073 Substituting and solving the double integrals, we get: Cij;kl ab;cd(x⊥, y⊥, x⊥, y⊥) = facefbde∂i x∂j yL(x⊥−y⊥)∂k xΓ(x⊥−y⊥)∂l yΓ(y⊥−x⊥) ×4 Γ3g4N3 c−¯ λ2(b⊥) 2ΓNc +4 Γ3g4N3 c +2¯ λ(b⊥) Γ2g2N2 cC(2) adj(x⊥, y⊥), (4.32) which concludes the calculation of the building block hαi,a xαk,c xαi0 ,a0 yαk0 ,c0 yi. 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: h0(x⊥)0(y⊥)i=Aik;i0k0 jl;j0l0Fac;a0c0 bd;b0d0hαi,a xαk,c xαi0 ,a0 yαk0 ,c0 yi 1 hαj,b xαl,d xαj0 ,b0 yαl0 ,d0 yi 2 .(4.33) 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 Aik;i0k0 jl;j0l0and Fac;a0c0 bd;b0d0, can be reduced to: h0(x⊥)0(y⊥)i =1 2Dik;i0k0 1ac;a0c0(x⊥,x⊥,y⊥,y⊥)Djl;j0l0 2bd;b0d0(x⊥,x⊥,y⊥,y⊥)Aik;i0k0 jl;j0l0Fac;a0c0 bd;b0d0 +Dik;i0k0 1ac;a0c0(x⊥,x⊥,y⊥,y⊥)Djj0;ll0 2bb0;dd0(x⊥,y⊥,x⊥,y⊥) +Dii0;kk0 1aa0;cc0(x⊥,y⊥,x⊥,y⊥)Djj0;ll0 2bb0;dd0(x⊥,y⊥,x⊥,y⊥)hAik;i0k0 jl;j0l0Fac;a0c0 bd;b0d0+Aik;i0k0 jl;l0j0Fac;a0c0 bd;d0b0i +Dik;i0k0 1ac;a0c0(x⊥,x⊥,y⊥,y⊥)Cjj0;ll0 2bb0;dd0(x⊥,y⊥,x⊥,y⊥) +2Dii0;kk0 1aa0;cc0(x⊥,y⊥,x⊥,y⊥)Cjj0;ll0 2bb0;dd0(x⊥,y⊥,x⊥,y⊥) +2Cii0;kk0 1aa0;cc0(x⊥,y⊥,x⊥,y⊥)Cjj0;ll0 2bb0;dd0(x⊥,y⊥,x⊥,y⊥) ×hAik;i0k0 jl;j0l0Fac;a0c0 bd;b0d0+Aik;i0k0 lj;j0l0Fac;a0c0 db;b0d0+Aik;i0k0 jl;l0j0Fac;a0c0 bd;d0b0+Aik;i0k0 lj;l0j0Fac;a0c0 db;d0b0i+[1 ↔2].(4.34) 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⊥): Dik;i0k0 ac;a0c0(x⊥, x⊥, y⊥, y⊥)Djl;j0l0 bd;b0d0(x⊥, x⊥, y⊥, y⊥)Aik;i0k0 jl;j0l0Fac;a0c0 bd;b0d0 =g4(∂2L(0⊥))4 N4 cC2 F ¯ λ2 1 (b⊥)¯ λ2 2 (b⊥) =1 g4α4 s(4π ∂2L(0⊥))4 N4 cC2 F ¯ λ2 1 (b⊥)¯ λ2 2 (b⊥) = C2 F g4¯ Q4 s 1 ¯ Q4 s 2 (4π ∂2L(0⊥))4≈h0(x⊥)ih0(y⊥)i, (4.35) 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⊥) = h0(x⊥)0(y⊥)i−h0(x⊥)ih0(y⊥)i 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 – 19 – JHEP01(2019)073 result of this work: Cov[](τ= 0+;x⊥, y⊥) ≡ h0(x⊥)0(y⊥)i−h0(x⊥)ih0(y⊥)i =∂i xΓ∂i yΓ(N2 c−1)A(4A2−B2) 16N2 cΓ5g4(p 1 q 2 +p 2 q 1 ) +(N2 c−1)(16A4+B4) 2N2 cΓ4g4p 1 p 2 +(∂i xΓ∂i yΓ)2(N2 c−1)A2 64N2 cΓ6g4q 1 q 2 +(N2 c−1)(4A2+B2) 2N2 cΓ2g44π ∂2L(0⊥)2¯ Q4 s 1 (Q2 s 2 r2−4+4e−Q2 s 2 r2 4)+ [1 ↔2] +(4A2+B2)2 g4Γ4N2 c "N6 c+ 2N4 c−19N2 c+ 8 (N2 c−1)2−4N6 c−3N4 c−26N2 c+ 16 (N2 c−1)(N2 c−4) e−Q2 s 1 r2 4 +(Nc−1)(Nc+ 3)N3 c (Nc+ 1)2(Nc+ 2)2Nc 2e−(Nc+1)r2Q2 s 2 2Nc+ (Nc+ 2) −2(Nc+ 1)e−Q2 s 2 r2 4e−(Nc+1)r2Q2 s 1 2Nc +(Nc+ 1)(Nc−3)N3 c (Nc−1)2(Nc−2)2Nc 2e−(Nc−1)r2Q2 s 2 2Nc+ (Nc−2) −2(Nc−1)e−Q2 s 2 r2 4e−(Nc−1)r2Q2 s 1 2Nc +r4 2Q2 s 1 Q2 s 2 −4r2Q2 s 1 1−e−Q2 s 2 r2 4+ 4(N2 c−8)(N2 c−1)(N2 c+ 4) (N2 c−4)2e−(Q2 s 1 +Q2 s 2 )r2 4# + [1 ↔2]!,(4.36) where the dependencies have been omitted for readability. The covariance of the full EMT is simply obtained from the previous expression as Cov[Tµν ](0+;x⊥, y⊥) = Cov[](0+;x⊥, y⊥)×tµνtσρ.(4.37) 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 Aand 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: p 1 , 2 ≡e−Q2 s 1 , 2 r2 4(Q2 s 1 , 2 r2+ 4) −4 (4.38) q 1 , 2 ≡e−Q2 s 1 , 2 r2 4Q4 s 1 , 2 r4+ 8Q2 s 1 , 2 r2+ 32−32.(4.39) For simplicity, in the previous expressions we also defined the following momentum scale: r2Q2 s 4=g2Nc 2Γ(r⊥)¯ λ(b⊥),(4.40) which is related to the one introduced in section 3by: Q2 s(r⊥, b⊥)r→0 =¯ Q2 s(b⊥)−4π ∂2L(0⊥).(4.41) (See appendix A). Since eq. (4.36) is somewhat lengthy, in the next sections we discuss a few simplifying limits in the context of the original MV model. – 20 – JHEP01(2019)073 4.1 Nc-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: A(r⊥) MV =−1 2G(r⊥) = 1 4πK 0 (mr) (4.42) B(r⊥) MV =1 4π,(4.43) where K 0 is a modified Bessel function. The mass mis an infrared scale that we introduce to regularize the divergent Green’s function G(r⊥). For simplicity we choose mto be the same mass scale introduced earlier in eq. (3.13). The leading behavior in the m→0 limit is: A(r⊥) MV ≈1 8πln4 m2r2,(4.44) 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: Γ(r⊥) MV =1 2πm2−r 2πmK 1 (mr)≈r2 8πln 4 m2r2(4.45) ∂i xΓ∂i yΓ MV ≈ − r2 16π2ln m2r2 42 ,(4.46) and for the saturation scale: Q2 s(r, b⊥) MV ≈¯ Q2 s(b⊥) ln4 m2r2.(4.47) (See appendix Afor 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, Aand Q2 sare divergent in both infrared and ultraviolet r→0 limits. The latter case enters our solution explicitly through the terms multiplied by the factor ∂2L(0⊥)≡−2 lim r→0A(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 Nc-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 – 21 – JHEP01(2019)073 transversal plane. The leading order of the expansion, of order N0 c, reads: Cov[ MV ](0+;x⊥,y⊥)N0 c =1 g4r8e−r2 2(Q2 s 1 +Q2 s 2 )16+32e Q2 s 1 r2 2 −64e Q2 s 1 r2 4−4er2 4(2Q2 s 1 +Q2 s 2 )Q4 s 2 r4−24π∂2L(0⊥)2¯ Q4 s 1 r4+8Q2 s 2 r2+48 +1 8er2 4(Q2 s 1 +Q2 s 2 )Q4 s 1 Q4 s 2 r8+(4Q2 s 1 Q2 s 2 r6+128r2)Q2 s 1 +Q2 s 2 +16r4Q2 s 1 +Q2 s 2 2+1024 +2er2 2(Q2 s 1 +Q2 s 2 )¯ Q4 s 1 r4Q2 s 2 r2−44π∂2L(0⊥)2+40+[1 ↔2].(4.48) The next term, of order N−2 c, reads: Cov[ MV ](0+;x⊥,y⊥)N−2 c =1 N2 cg4r8e−r2 2(Q2 s 1 +Q2 s 2 )2Q2 s 1 r2+Q2 s 2 r2+82 +4Q2 s 1 r2(8+Q2 s 1 r2)e Q2 s 2 r2 2−8(8+Q2 s 1 r2)(4+Q2 s 1 r2)e Q2 s 2 r2 4 +4er2 4(2Q2 s 1 +Q2 s 2 )Q4 s 2 r4−2(4π∂2L(0⊥))2¯ Q4 s 1 r4+8Q2 s 2 r2+16Q2 s 1 r2 −1 8er2 4(Q2 s 1 +Q2 s 2 )Q4 s 1 Q4 s 2 r8+(4Q2 s 1 Q2 s 2 r6+128r2)Q2 s 1 +Q2 s 2 +16r4Q2 s 1 +Q2 s 2 2−1024 −2er2 2(Q2 s 1 +Q2 s 2 )¯ Q4 s 1 r4(Q2 s 2 r2−4)(4π∂2L(0⊥))2+32Q2 s 1 r2−4Q2 s 1 Q2 s 2 r4+[1 ↔2]. (4.49) In order to give a general idea of the magnitude and analytical features of our solution, on figure 3we draw these functions in the GBW model as a function of the dimensionless product rQsfor Qs 1 =Qs 2 . Note that in this limit, as we are ignoring all logarithmic factors, we also have Qs=¯ Qs. The N−2 cterm 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 Nc-expansion provide a neat approximation to the complete result (see right plot of figure 3). Comparing this curve with the N0 c-order term we notice that the large-Nclimit yields a 12.5% error in the r→0 limit, which is a reasonable approximation. In the rQs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 Nc-expansion presented above, eq. (4.48) and eq. (4.49): lim rQs1Cov[ MV ](0+;x⊥, y⊥) = 2N2 c−14π ∂2L(0⊥)2¯ Q4 s 1 Q2 s 2 +¯ Q4 s 2 Q2 s 1  g4N2 cr2.(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 – 22 – JHEP01(2019)073 � ��� ��� ��� ��� ��� ��� ��� ��� � � � � � �� �������� � ��� �� ������ rQ s g4Cov/Q8 s O(N0 c) O(N0 c)+O(N2 c) ����� ����� � ����� � � � � � �� �� ����� � ������� Ratio rQ s Figure 3. LEFT: sum of the first two orders of the Nc-expansion of the energy density covariance against rQsfor Qs 1 =Qs 2 and Nc= 3. Blue full curve: N0 c-order term. Red dashed curve: sum of N0 cand N−2 c-order terms. RIGHT: ratio between the complete result and the sum of the first two orders of the Nc-expansion. the previous result with a single one-point correlator we obtain the following expression: lim rQs1Cov[](0+;x⊥, y⊥) h0(x⊥)i MV =4 g2Ncr2¯ Q2 s 1 Q2 s 2 ¯ Q2 s 2 +¯ Q2 s 2 Q2 s 1 ¯ Q2 s 1 .(4.51) In the opposite limit, r→0, the covariance tends to: lim r→0Cov[ MV ](0+;x⊥, y⊥) = C F 2Ncg4Q4 s 1 Q4 s 2 +4π∂2L(0⊥)2¯ Q4 s 1 Q4 s 2 +¯ Q4 s 2 Q4 s 1  =3C F 2Ncg44π∂2L(0⊥)4¯ Q4 s 1 ¯ Q4 s 2 ,(4.52) and the normalized covariance: lim r→0Cov[](0+;x⊥, y⊥) h0(x⊥)ih0(y⊥)i MV =Q4 s 1 Q4 s 2 +4π∂2L(0⊥)2¯ Q4 s 1 Q4 s 2 +¯ Q4 s 2 Q4 s 1  (4π∂2L(0⊥))4¯ Q4 s 1 ¯ Q4 s 2 (N2 c−1) =3 N2 c−1. (4.53) In both expressions we applied eq. (4.41) in the last step. 4.2 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: hαi,a(x⊥)αk,c(x⊥)αi0 ,a0(y⊥)αk0 ,c0(y⊥)i GG =hαi,a(x⊥)αk,c(x⊥)ihαi0 ,a0(y⊥)αk0 ,c0(y⊥)i +hαi,a(x⊥)αi0 ,a0(y⊥)ihαk,c(x⊥)αk0 ,c0(y⊥)i +hαi,a(x⊥)αk0 ,c0(y⊥)ihαk,c(x⊥)αi0 ,a0(y⊥)i.(4.54) – 23 – JHEP01(2019)073 ���� ���� ���� ����� ��� � � � � � �� �� ���� �� � ������� ���� ������ ������ ����� ������������� ���������� �� �������� Cov/(h✏0(x?)ih✏0(y?)i) rQ s � ��� � �� �� �� ��� � � � � � �� �� ����� � ������� rQ s Ratio Figure 4. LEFT: comparison of the normalized covariance of energy density  0 against r Qsfor Qs 1 =Qs 2 ,Nc= 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/[(N2 c−1)r2Q2 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: lim rQs1Cov[](0+;x⊥, y⊥) h0(x⊥)ih0(y⊥)i GG =16( ¯ Q4 s 1 +¯ Q4 s 2 ) r4(N2 c−1) ¯ Q4 s 1 ¯ Q4 s 2 (4π∂2L(0⊥))2.(4.55) The ∼1/r4decreasing behavior displayed by eq. (4.55) is in clear contrast with the ∼1/r2 asymptotic behavior of our result. This potentially implies much different results and physical interpretations for any observable built from this quantity. 5 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 – 24 – JHEP01(2019)073 b1 b2 bn (↵1)1 (↵1)2 (↵1)n b c  1 (↵p)1 (↵p)2 (↵p)n (↵p+1)1 (↵p+1)2 (↵p+1)n (...) (↵m1)1 (↵m1)2 (↵m1)n a1 a2 an c  p c  p+1 c  m1 a  ⇢c1 ⇢c2 ⇢cp+1 ⇢cp+2 ⇢cm (...) Figure 5. Schematic representation of the connected correlator Hm,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 nWilson lines, whose correlator is represented as a dark square (see figure 6). This aspect can be comprised in a redefinition of H1,n as: H1,n(b−,a−|{b},{a})i≡g n X j=1 µ2(y−)Fn(b−,y−|{b}{β})|βj=dFn(y−,a−|{β}{a})|βj=d0 ×∂i yZdz⊥dw⊥G(z⊥−xj⊥)G(w⊥−y⊥)f(z⊥−w⊥)hz⊥+w⊥ 2fcdd0, (B.7) where Fndenotes the correlator of nWilson lines. Note that in these formulas the bracketed indices represent a set of ncolor 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: H1,n(b−, a−|{b},{a})i ≈g n X j=1 µ2(y−)Fn(b−, y−|{b}{β})|βj=dFn(y−, a−|{β}{a})|βj=d0h(b⊥)∂i yL(xj⊥−y⊥)fcdd0 =gλ(y−, b⊥) n X j=1 ∂i yL(xj⊥−y⊥)fcdd0Fn(b−, y−|{b}{β})|βj=dFn(y−, a−|{β}{a})|βj=d0, (B.8) 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 (x0 ⊥ −b⊥)i∂ih(b⊥). Whether x0 ⊥=x⊥or x0 ⊥=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 ˜ρadepends is the same as the one of the Wilson line that it is attached to, yielding the following simplification of the – 31 – JHEP01(2019)073 b1 b2 bn a1 a2 an = b a  b1 b2 bn a1 a2 an b a  b1 b2 bn a1 a2 an b a  + + (...) + b1 b2 bn a1 a2 an b a  ⇢c ⇢c ⇢c ⇢c b0 b0 b0 Figure 6. Schematic representation of eq. (B.9). The circles in the right hand of the equation represent couplings of the external source ρcto Wilson lines inside the correlator (dark square). Each of these couplings multiplies the correlator by a gλ(b−, b⊥)fcbjb0∂yL(xj⊥−y⊥) factor. previous expression (see figure 6): H1,n(b−, a−|{b},{a})i=gλ(b−, b⊥) n X j=1 ∂i yL(xj⊥−y⊥)fc bjb0Fn(b−, a−|{β}{a})|βj=b0.(B.9) Having defined all the basic pieces of the calculation of a correlator with mexternal sources and nWilson lines, we can go back to our particular case. According to the notation used in [41], the correlator in eq. (B.1) corresponds to Fm,n with m= 4, n= 4. By direct application of eq. (B.2): D˜ρi,e xUea x˜ρk,f xUfc x˜ρi0 ,e0 yUe0 a0 y˜ρk0 ,f0 yUf0 c0 yE=D˜ρi,e x˜ρk,f x˜ρi0 ,e0 y˜ρk0 ,f0 yEDUea xUfc xUe0 a0 yUf0 c0 yE +D˜ρi,e x˜ρk,f xED˜ρi0 ,e0 y˜ρk0 ,f0 yUea xUfc xUe0 a0 yUf0 c0 yEc +D˜ρi0 ,e0 y˜ρk0 ,f0 yED˜ρi,e x˜ρk,f xUea xUfc xUe0 a0 yUf0 c0 yEc +D˜ρi,e x˜ρi0 ,e0 yED˜ρk,f x˜ρk0 ,f0 yUea xUfc xUe0 a0 yUf0 c0 yEc +D˜ρi,e x˜ρk0 ,f0 yED˜ρk,f x˜ρi0 ,e0 yUea xUfc xUe0 a0 yUf0 c0 yEc +D˜ρk,f x˜ρi0 ,e0 yED˜ρi,e x˜ρk0 ,f0 yUea xUfc xUe0 a0 yUf0 c0 yEc +D˜ρk,f x˜ρk0 ,f0 yED˜ρi,e x˜ρi0 ,e0 yUea xUfc xUe0 a0 yUf0 c0 yEc +D˜ρi,e xUea x˜ρk,f xUfc x˜ρi0 ,e0 yUe0 a0 y˜ρk0 ,f0 yUf0 c0 yEc .(B.10) 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: D˜ρi,e x˜ρk,f x˜ρi0 ,e0 y˜ρk0 ,f0 yE=h˜ρi,e x˜ρk,f xih˜ρi0 ,e0 y˜ρk0 ,f0 yi+h˜ρi,e x˜ρi0 ,e0 yih˜ρk,f x˜ρk0 ,f0 yi+h˜ρi,e x˜ρk0 ,f0 yih˜ρk,f x˜ρi0 ,e0 yi. (B.11) The three terms resulting from this expansion are addressed on section 4, where we derive the ‘disconnected’ function: Dij;kl ab;cd(u⊥, u0 ⊥, v⊥, v0 ⊥) = Tij;kl (u⊥, u0 ⊥, v⊥, v0 ⊥) Z∞ −∞ dz−Zz− −∞ dw− λ(z−, b⊥)λ(w−, b⊥) ×C(2) adj(z−, w−;u⊥, u0 ⊥) + C(2) adj(z−, w−;v⊥, v0 ⊥)δABδCDQABCD abcd (w−;u⊥, u0 ⊥, v⊥, v0 ⊥), (B.12) – 32 – 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⊥,u0 ⊥,v⊥, and v0 ⊥will take the values x⊥or y⊥. Here we also introduce the following function: Tij;kl(x⊥, x0 ⊥, y⊥, y0 ⊥)≡∂i x∂j x0L(x⊥−x0 ⊥)∂k y∂l y0L(y⊥−y0 ⊥).(B.13) 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: Z∞ −∞ dz−dw−dz−0dw−0 *∂i˜ρe(z−, x⊥) ∇2 ∂k˜ρf(w−, x⊥) ∇2 ∂i0˜ρe0(z−0, y⊥) ∇2 ∂k0˜ρf0(w−0, y⊥) ∇2+ ×DUea(z−, x⊥)Ufc(w−, x⊥)Ue0a0(z−0, y⊥)Uf0 c0(w−0, y⊥)E =Dik;i0k0 ac;a0c0(x⊥, x⊥, y⊥, y⊥) + Dii0;kk0 aa0;cc0(x⊥, y⊥, x⊥, y⊥) + Dik0;ki0 ac0;ca0(x⊥, y⊥, x⊥, y⊥).(B.14) All remaining terms of eq. (B.10) contain the connected correlator h. . .ic, 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)): Z∞ −∞ dz−dw−dz−0dw−0 *∂i˜ρe(z−, x⊥) ∇2 ∂k˜ρf(w−, x⊥) ∇2 ∂i0˜ρe0 (z−0, y⊥) ∇2 ∂k0˜ρf0 (w−0, y⊥) ∇2 ×Uea(z−, x⊥)Ufc(w−, x⊥)Ue0a0(z−0, y⊥)Uf0 c0(w−0, y⊥)Ec.(B.15) As we have four different longitudinal coordinates, in order to compute eq. (B.15) we need to consider all regions of the integration space.4For example, applying eq. (B.5) in the region where z−>w−>z−0 >w−0 we have: H4,4(z−,−∞|e, f, e0, f0;a, c, a0, c0) = H1,1(z−, w−|e;α 1 )iH1,2(w−, z−0|α 1 , f ;α 2 , β 1 )k ×H1,3(z−0, w−0|α 2 , β 1 , e0;α 3 , β 2 , γ 1 )i0 ×H1,4(w−0,−∞|α 3 , β 2 , γ 1 , f0;a, c, a0, c0)k0,(B.16) where, according to eq. (B.9), the first factor reads: H1,1(z−, w−|e;α 1 )i=g λ(z−, b⊥)∂i xL(0⊥)feeα Uαα1(z−, w−;x⊥)= 0,(B.17) 4Namely the regions where z−> z−0 > w−> w−0,z−> z−0 > w−0 > w−, etcetera. 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−0 >w−>w−0) 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 Dij;kl ab;cd. – 33 – JHEP01(2019)073 a0 b0 c0 d0 a b c d 1 A B C ⇢c 0 ⇢d0 z  w  a0 c0 d0 a b c d 1 A B C ⇢c 0 ⇢d0 a0 w  w0 Figure 7. Schematic representation of the connected correlator H2,4factorized in eq. (B.19). which vanishes due to the antisymmetric property of the SU(Nc) 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: Cij;kl ab;cd(u⊥, u0 ⊥, v⊥, v0 ⊥) =Z∞ −∞ dz−dz−0dw−dw−0 *∂i˜ρa0(z−, u⊥) ∇2 ∂j˜ρb0(z−0, u0 ⊥) ∇2+ ×*∂k˜ρc0(w−, v⊥) ∇2 ∂l˜ρd0(w−0, v0 ⊥) ∇2Ua0a(z−, u⊥)Ub0b(z−0, u0 ⊥)Uc0c(w−, v⊥)Ud0d(w−0, v0 ⊥)+c =∂i u∂j u0L(u⊥−u0 ⊥)Z∞ −∞ dz−dw−dw−0λ(z−, b⊥) ×*∂k˜ρc0(w−, v⊥) ∇2 ∂l˜ρd0(w−0, v0 ⊥) ∇2Ua0a(z−, u⊥)Ua0b(z−, u0 ⊥)Uc0c(w−, v⊥)Ud0d(w−0, v0 ⊥)+c , (B.18) where b⊥= (x⊥+y⊥)/2. The only nonvanishing contributions to this integral come from the regions where z−> w−> w−0 and z−> w−0 > w−, as in the other cases eq. (B.9) introduces a vanishing H1,1factor. For the z−>w−>w−0 region we have (see figure 7): ∂i u∂j u0L(u⊥−u0 ⊥)Z∞ −∞ dz−Zz− −∞ dw−Zw− −∞ dw−0λ(z−, b⊥)C(2) adj(z−, w−;u⊥, u0 ⊥) ×*∂k˜ρc0(w−, v⊥) ∇2Ua0A(w−, w−0;u⊥)Ua0B(w−, w−0;u0 ⊥)Uc0C(w−, w−0;v⊥)+c ×*∂l˜ρd0(w−0, v0 ⊥) ∇2UAa(w−0, u⊥)UBb(w−0, u0 ⊥)UCc(w−0, v⊥)Ud0d(w−0, v0 ⊥)+c ,(B.19) 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 – 34 – JHEP01(2019)073 one external source and three Wilson lines, thus corresponding to: H1,3(w−, w−0|a0, a0, c0;A, B, C)k/(g λ(w−, b⊥)) =∂k vL(v⊥−u⊥)fc0a0αDUαA(w−, w−0;u⊥)Ua0B(w−, w−0;u0 ⊥)Uc0C(w−, w−0;v⊥)E +∂k vL(v⊥−u0 ⊥)fc0a0αDUa0A(w−, w−0;u⊥)UαB(w−, w−0;u0 ⊥)Uc0C(w−, w−0;v⊥)E+ 0. (B.20) Substituting the result for the three-point adjoint correlator from [40]: hUaa0(w−, x 1 ⊥)Ubb0(w−, x 2 ⊥)Ucc0(w−, x 3 ⊥)i =1 2N2 cC F fabcfa0b0c0+N2 c N2 c−4dabcda0b0c0exp   −g2Nc 4¯ λ(w−, b⊥)X i>j Γ(xi ⊥−xj ⊥)   ≡1 2N2 cC F fabcfa0b0c0+N2 c N2 c−4dabcda0b0c0C(3) adj(w−;x 1 ⊥, x 2 ⊥, x 3 ⊥),(B.21) (where, again, we approximated the function has h((x⊥+y⊥)/2)), we get to: H1,3(w−, w−0|a0, a0, c0;A, B, C)k=gfABCC(3) adj(w−, w−0;u⊥, u0 ⊥, v⊥) ×λ(w−, b⊥)∂k vL(v⊥−u0 ⊥)−L(v⊥−u⊥).(B.22) Here, the color factor (the one that cancels (2N2 cC 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: H1,4(w−0|{A, B, C, d0},{a, b, c, d})l/(g λ(w−0, b⊥)) =∂l v0L(v0 ⊥−u⊥)fd0AαDUαa(w−0, u⊥)UBb(w−0, u0 ⊥)UCc(w−0, v⊥)Ud0d(w−0, v0 ⊥)E +∂l v0L(v0 ⊥−u0 ⊥)fd0BαDUAa(w−0, u⊥)Uαb(w−0, u0 ⊥)UCc(w−0, v⊥)Ud0d(w−0, v0 ⊥)E +∂l v0L(v0 ⊥−v⊥)fd0CαDUAa(w−0, u⊥)UBb(w−0, u0 ⊥)Uαc(w−0, v⊥)Ud0d(w−0, v0 ⊥)E+ 0. (B.23) Substituting this expression and summing the contribution from the z−>w−0 >w−region the ‘connected’ function finally yields: Cij;kl ab;cd(u⊥, u0 ⊥, v⊥, v0 ⊥) =g2h3(b⊥)∂i u∂j u0L(u⊥−u0 ⊥)Z∞ −∞ dz− Zz− −∞ dw− Zw− −∞ dw−0µ2(z−)µ2(w−)µ2(w−0) ×C(2) adj(z−, w−;u⊥, u0 ⊥) h∂k vL(v⊥−u0 ⊥)−L(v⊥−u⊥)C(3) adj(w−, w−0;u⊥, u0 ⊥, v⊥) ×∂l v0fAeDfCBeL(v0 ⊥−u⊥)+fACefDBeL(v0 ⊥−u0 ⊥) + fABefeCDL(v0 ⊥−v⊥) ×QABCD abcd (w−0;u⊥, u0 ⊥, v⊥, v0 ⊥)i+    l←→ k c←→ d v⊥ ←→ v0 ⊥    !.(B.24) – 35 – JHEP01(2019)073 Now we can rewrite eq. (4.4) in terms of Dij;kl ab;cd and Cij;kl ab;cd as: hαi,a(x⊥)αk,c(x⊥)αi0 ,a0(y⊥)αk0 ,c0(y⊥)i =Dik;i0k0 ac;a0c0(x⊥, x⊥, y⊥, y⊥) + Dii0;kk0 aa0;cc0(x⊥, y⊥, x⊥, y⊥) +Dik0;ki0 ac0;ca0(x⊥, y⊥, x⊥, y⊥) + Cik;i0k0 ac;a0c0(x⊥, x⊥, y⊥, y⊥) + Ci0k0;ik a0c0;ac(y⊥, y⊥, x⊥, x⊥) +Cii0;kk0 aa0;cc0(x⊥, y⊥, x⊥, y⊥) + Cik0;ki0 ac0;ca0(x⊥, y⊥, x⊥, y⊥) + Cki0;ik0 ca0;ac0(x⊥, y⊥, x⊥, y⊥) +Ckk0;ii0 cc0;aa0(x⊥, y⊥, x⊥, y⊥).(B.25) The fact that in our particular case the Wilson lines depend on only two transverse coordinates (x⊥and y⊥) yields a significant simplification in the final expression. For example, by taking u⊥=u0 ⊥in eq. (B.24) we can see that the first two connected terms of eq. (B.25) yield 0. The next four terms take the following form: Cij;kl ab;cd(x⊥, y⊥, x⊥, y⊥) = 2g2h3(b⊥)∂i x∂j yL(x⊥−y⊥)Z∞ −∞ dz− Zz− −∞ dw− Zw− −∞ dw−0µ2(z−)µ2(w−)µ2(w−0) ×C(2) adj(z−, w−;x⊥, y⊥)∂k x(L(x⊥−x⊥)−L(x⊥−y⊥))C(3) adj(w−, w−0;x⊥, y⊥, x⊥) ×∂l y(L(y⊥−y⊥)−L(y⊥−x⊥)) fACefBDeQABCD abcd (w−0;x⊥, y⊥, x⊥, y⊥),(B.26) where we applied the Jacobi identity of SU(Nc). The previous expression contains a trivial projection of the adjoint Wilson line quadrupole: fACefBDeQABCD abcd (w−0;x⊥, y⊥, x⊥, y⊥) = fABefCDeQABCD acbd (w−0;x⊥, x⊥, y⊥, y⊥) =facefbdeC(2) adj(w−0;x⊥, y⊥).(B.27) (See appendix Cfor the detailed computation). Also, in the chosen limit the adjoint Wilson line tripole yields: C(3) adj(w−, w−0;x⊥, y⊥, x⊥) = C(2) adj(w−, w−0;x⊥, y⊥).(B.28) Both correlators tend to the dipole function in the case of only two different transverse coordinates. We are left with a product of three dipole functions that combine as: C(2) adj(z−, w−;x⊥, y⊥)C(2) adj(w−, w−0;x⊥, y⊥)C(2) adj(w−0;x⊥, y⊥) = exp −g2Nc 2Γ(x⊥−y⊥)h(b⊥)¯µ2(z−, w−) + ¯µ2(w−, w−0) + ¯µ2(w−0) = exp (−g2Nc 2Γ(x⊥−y⊥)h(b⊥) Zz− w− du−µ2(u−)+Zw− w−0 du−µ2(u−)+Zw−0 −∞ du−µ2(u−)!) =C(2) adj(z−;x⊥, y⊥),(B.29) – 36 – JHEP01(2019)073 and therefore: Cij;kl ab;cd(x⊥,y⊥,x⊥,y⊥) = g2 2facefbdeh3(b⊥)∂i x∂j yL(x⊥−y⊥)∂k xΓ(x⊥−y⊥)∂l yΓ(y⊥−x⊥) ×Z∞ −∞ dz−Zz− −∞ dw−Zw− −∞ dw−0µ2(z−)µ2(w−)µ2(w0−)C(2) adj(z−;x⊥,y⊥). (B.30) Solving the double integral, we obtain: Cij;kl ab;cd(x⊥, y⊥, x⊥, y⊥) = facefbde∂i x∂j yL(x⊥−y⊥)∂k xΓ(x⊥−y⊥)∂l yΓ(x⊥−y⊥) ×4 Γ3g4N3 c−¯ λ2(b⊥) 2ΓNc +4 Γ3g4N3 c +2¯ λ(b⊥) Γ2g2N2 cC(2) adj(x⊥, y⊥). (B.31) C The correlator of four Wilson lines in the adjoint representation C.1 Reexponentiation method 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 nlayers of length ∆x−. Due to the properties of path ordered exponentials, this leads to the factorization of the Wilson line into a product of nindependent contributions from each zone: U(x−, x⊥)ij ≡U(n) ij ≡(Un(x− n, x⊥)Un−1(x− n−1, x⊥). . . U1(x− 1, x⊥))ij,(C.1) 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 nfactors to order g2: U(x−, x⊥)ij ≈δik +ig ˜ A+a(x− n, x⊥)ta ik∆x−−g2C F 2λ(x− n, b⊥)L(0⊥)∆x−δikU(n−1) kj ,(C.2) where ˜ A+a(x−, x⊥) = −˜ρa(x−, x⊥) ∇2(C.3) is the only non-trivial component of the gluon field expressed in the covariant gauge (defined by the condition ∂µ˜ Aµ= 0). Note that in the g2-order term of eq. (C.2) we already applied the two-point correlator, whose discretized version reads: h˜ A+a(x−, x⊥)˜ A+b(y−, y⊥)i=λ(x−, b⊥)δabL(x⊥−y⊥)δx−y− ∆x−.(C.4) We iterate this process n−1 more times neglecting terms of order (∆x−)2or higher. Then, we rearrange the resulting terms in the form of the first orders of an expanded exponential. – 37 – 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: DTr nU(x⊥)U†(y⊥)oE ≈U(n−1) kl (x⊥)U(n−1)† lj (y⊥)δik +ig ˜ A+a(x− n, x⊥)ta ik∆x−−g2C F 2λ(x− n, b⊥)L(0⊥)∆x−δik ×δji −ig ˜ A+b(x− n, y⊥)tb ji∆x−−g2C F 2λ(x− n, b⊥)L(0⊥)∆x−δji =DTr nU(x⊥)U†(y⊥)oE(n−1) 1−g2 2C F ∆x−λ(x− n, b⊥)Γ(x⊥−y⊥).(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: DTr nU(x⊥)U†(y⊥)oE≈ 1−g2 2C F Γ(x⊥−y⊥)h(b⊥) n X i=1 ∆x−µ2(x− i)! =1−g2 2C F Γ(x⊥−y⊥)¯ λ(x−, b⊥).(C.6) Lastly, we assume that the neglected higher order terms add up to an exponential expression, which reads: DTr nU(x⊥)U†(y⊥)oE= exp −g2 2C F Γ(x⊥−y⊥)¯ λ(x−, b⊥).(C.7) 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 h. . .i. In the next subsection we will make use of this technique to obtain the behavior of the following adjoint Wilson line quadrupole: DUAa(x⊥)UBb(x⊥)UCc(y⊥)UDd(y⊥)E,(C.8) 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: DUab(z⊥)Ucd(z⊥)Uef (x⊥)Ugh(y⊥)E.(C.9) 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 – 38 – JHEP01(2019)073 notation similar to the one used in [45]. We absorb the g∆x−factor in the definition of our fields: g˜ A+a(x−, x⊥)∆x−≡˜ A+a(x−, x⊥),(C.10) which yields the following two-point function: h˜ A+a(x−, x⊥)˜ A+b(y−, y⊥)i=δx−y−δabBxy(x−, b⊥),(C.11) where, due to the discretization of the rapidity range, the Kronecker delta δx−y−now takes the place of the Dirac delta. For simplicity we also introduced: Bxy(x−, b⊥)≡g2∆x−λ(x−, b⊥)L(x⊥−y⊥).(C.12) Using this notation the expansion to order g2of the adjoint Wilson line looks like: Uab(x−, x⊥)=(Uab1)(n−1) δb1b1−Nc 2Bx(x− n, b⊥)−˜ Ag(x− n, x⊥)fb1gb.(C.13) Performing this expansion for every Wilson line in eq. (C.9) and neglecting terms of order (∆x−)2or higher we get: DUab(z⊥)Ucd(z⊥)Uef (x⊥)Ugh(y⊥)E =DUaa0(z⊥)Ucc0(z⊥)Uee0(x⊥)Ugg0(y⊥)E(n−1) ×δa0bδc0dδe0fδg0h1−Nc 2(2Bz+Bx+By)+δa0bδc0dfe0mf fg0mhBxy +δa0bδe0ffc0mdfg0mhBzy +δa0bδg0hfe0mf fc0mdBzx +δe0fδc0dfa0mbfg0mhBzy +δg0hδc0dfe0mf fa0mbBzx +δe0fδg0hfa0mbfc0mdBz.(C.14) We express the previous lines as a matrix equation: Uaceg bdfh = (Uaceg a0c0e0g0)(n−1)Ta0c0e0g0 bdfh , for which we introduce the following color vector basis: u1=δeaδgc u2=δcaδge u3=δgaδec w1=deamdgcm w2=dcamdgem w3=dgamdecm z1=deamfgcm z2=dcamfgem z3=dgamfecm.(C.15) It can be shown via color algebra arguments that this ensemble covers the entirety of possible interactions embodied in Ta0c0e0g0 bdfh (see [46]).5The last three (z1,z2,z3) form a basis that does not mix with the rest of the vectors in the Gaussian model we are considering. Thus, if we expressed Ta0c0e0h0 bdfh in this 9-dimensional space it would look like a block diagonal matrix with a 6 ×6 part corresponding to the vectors ui,wiand a 3 ×3 sector corresponding to the ziset. In our specific calculation, the vectors that we are 5In [46], the author mentions only 8 such tensors, but that is because he is dealing with SU(3), and there exists a relation between the SU(Nc) generators, valid only for Nc= 3, which reduces the number of independent rank 4 tensors from 9 to 8 in that case. – 39 – JHEP01(2019)073 interested in propagating live in the 6-dimensional space defined by the first two sets, so it will be enough to consider Ta0c0e0h0 bdfh in the basis formed by uiand wi. To build this matrix we propagate these six vectors using eq. (C.14): Ta0c0e0h0 bdfh δe0a0δg0c0=δfbδhd 1−Nc 2(2Bz+Bx+By−2Bzx −2Bzy) +fbfmfdhm (Bz+Bxy −Bzx −Bzy) =δfbδhd 1−g2Nc 2∆x−λ(x− n, b⊥)(Γ(z⊥−x⊥) + Γ(z⊥−y⊥)) +fbfmfdhm g2 2∆x−λ(x− n, b⊥) (Γ(z⊥−x⊥) + Γ(z⊥−y⊥)−Γ(x⊥−y⊥)). (C.16) The SU(Nc) factor fbfmfdhm, as well as the ones resulting from permutations of its indices, can be expressed in terms of our basis vectors by means of the following identity: fabmfcdm =2 Nc (δacδbd −δadδbc) + dacedbde −dadedbce,(C.17) and thus the propagation of u1reads: Tu1=u11−g2Nc 2∆x−λ(x− n,b⊥)(Γ(z⊥−x⊥)+Γ(z⊥−y⊥)) +g2 2∆x−λ(x− n,b⊥)2 Nc (u2−u3)+w2−w3(Γ(z⊥−x⊥)+Γ(z⊥−y⊥)−Γ(x⊥−y⊥)). (C.18) Repeating this process for the remaining vectors, we obtain: Tu2=u21−g2Nc 2∆x−λ(x− n,b⊥)Γ(x⊥−y⊥)(C.19) Tu3=u31−g2Nc 2∆x−λ(x− n,b⊥)(Γ(z⊥−x⊥)+Γ(z⊥−y⊥)) +g2 2∆x−λ(x− n,b⊥)2 Nc (u2−u1)+w2−w1(Γ(z⊥−x⊥)+Γ(z⊥−y⊥)−Γ(x⊥−y⊥)) (C.20) Tw1=w11−g2Nc 8∆x−λ(x− n,b⊥)(Γ(x⊥−y⊥)+3Γ(z⊥−x⊥)+3Γ(z⊥−y⊥)) +g2 2∆x−λ(x− n,b⊥)2 Nc−Nc 4(w2−w3)+4 N2 c−1(u2−u3) ×(Γ(x⊥−y⊥)−Γ(z⊥−x⊥)−Γ(z⊥−y⊥)) (C.21) Tw2=w21−g2Nc 4∆x−λ(x− n,b⊥)(Γ(x⊥−y⊥)+Γ(z⊥−x⊥)+Γ(z⊥−y⊥))(C.22) Tw3=w31−g2Nc 8∆x−λ(x− n,b⊥)(Γ(x⊥−y⊥)+3Γ(z⊥−x⊥)+3Γ(z⊥−y⊥)) +g2 2∆x−λ(x− n,b⊥)2 Nc−Nc 4(w2−w1)+4 N2 c−1(u2−u1) ×(Γ(x⊥−y⊥)−Γ(z⊥−x⊥)−Γ(z⊥−y⊥)).(C.23) – 40 – Reproduced with permission of copyright owner. Further reproduction prohibited without permission.