Tracing the origin of azimuthal gluon correlations in the color glass condensate

We examine the origins of azimuthal correlations observed in high energy proton-nucleus collisions by considering the simple example of the scattering of uncorrelated partons off color fields in a large nucleus. We demonstrate how the physics of fluctuating color fields in the color glass condensate (CGC) effective theory generates these azimuthal multiparticle correlations and compute the corresponding Fourier coefficients v_n within different CGC approximation schemes. We discuss in detail the qualitative and quantitative differences between the different schemes. We will show how a recently introduced color field domain model that captures key features of the observed azimuthal correlations can be understood in the CGC effective theory as a model of non-Gaussian correlations in the target nucleus.


I. INTRODUCTION
Azimuthal anisotropies of multiparticle correlations observed in small systems such as those produced in p+Pb, d+Au, or 3 He+Au collisions have been computed in various theoretical frameworks. Within different calculations these correlations are either dominantly due to initial state parton correlations in the projectile and target [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] or from final state correlations that are generated by the collective flow of matter produced in the collision [16][17][18][19][20][21][22]. Since both of these frameworks are able to describe key features of the data, disentangling the different effects and understanding how correlations are generated in both approaches is essential to obtain novel insight into the QCD dynamics of ultradense parton systems.
We will focus in this work on multiparticle correlations which are generated in the initial state. We will discuss within the color glass condensate (CGC) effective theory of high energy QCD [23] the origin of these correlations and we will critically examine the assumptions underlying different calculations in this framework. Even though some of the models may appear very different, their common features and their differences can be understood systematically as approximations to the underlying QCD dynamics.
We will orient our discussion of initial state correlations within the "dilute-dense" power counting in the CGC, where the incoming parton densities are assumed to be small in the projectile but large in the target nucleus. In this limit, analytical and numerical computations are comparatively simple and thus permit systematic comparisons between different approximation schemes. We note however that the kinematic region where the strongest azimuthal correlations are seen in experiments correspond more to a "dense-dense" situation as the parton densities are also large in the incoming projectile. A systematic power counting then requires one to solve classical Yang-Mills equations in the presence of projectile and target color sources, which has only been achieved numerically [14,24]. While there are important qualitative differences between dense-dense and dilute-dense systems, we believe that the lessons inferred from analytical and numerical studies the dilute-dense case can nevertheless be valuable for the discussion of the phenomenologically more relevant dense-dense collision systems.
We will here compute the two particle correlation function for quarks scattering off a large nucleus in the dilute-dense CGC description. The main ingredient of this computation is the so-called "dipole-dipole" correlator of light-like Wilson lines and we will compute this quantity in the following approximations: i) the Gaussian two gluon exchange ("Glasma graph") approximation [1][2][3][4][5][25][26][27][28][29] and ii) the nonlinear Gaussian approximation [30][31][32][33][34][35][36][37][38][39][40]. We will determine the second, third, and fourth azimuthal Fourier coefficients of the two particle correlation function within these two approximation schemes and compare the results to numerical lattice simulations of the full correlator in the McLerran-Venugopalan (MV) model [41][42][43] and after renormalization group (JIMWLK) evolution [44][45][46][47] of the MV model initial conditions to higher rapidities [13]. We will further discuss how our computations relate to the "color field domain model" introduced in [10,11,48,49] based on ideas developed in [6,7]. Since many features of this model appear similar to those discussed previously [1,2,4,5], it is important to understand the interpretation and justification for this model from first principles. We will demonstrate that the effects of the color field domain model can be reproduced in the CGC effective theory if non-Gaussian correlations are assumed to play an important role.
This paper is organized as follows. In the next section, we shall discuss the physical picture of how initial state multiparticle correlations are generated and derive the leading order expressions for single and double inclusive distributions. These can be expressed in terms of correlators of lightlike Wilson lines, which we will calculate in Sec. 3 within the Glasma graph and nonlinear Gaussian approximations. We then compute the Fourier moments v n of two particle correlations in both approximation schemes in Sec. 4 and compare our results with numerical lattice simulations of the full correlation function. In Section 5 the analytical and numerical results obtained are then compared to the color field domain model. We first cast our analytical results in terms of color electric field correlators and make a direct comparison with expressions for the same correlator in the color field domain model. We will demonstrate that our results provide a clear interpretation of the color field domain model which clarifies the discussion in the recent literature. We end with a summary of the results of the paper and an outlook on further research directions in computations of multiparton correlations in high energy QCD.

II. MULTIPARTICLE CORRELATIONS FROM FLUCTUATING COLOR FIELDS
We begin our discussion of the physics of initial state correlations with the simplest possible example of the high energy scattering of individual (uncorrelated) quarks off a large nucleus. Our general picture is that each parton scatters independently off the color field of the nucleus receiving a transverse momentum kick in the process. As noted previously [6,7,25], the color fields fluctuate from event to event and are locally organized in domains of size ∼ 1/Q s as illustrated in Fig. 1. When two (or more) quarks scatter off the same domain, they will receive a similar kick whenever they are in the same color state. This leads to a correlation which is suppressed by 1/N c 2 (in the limit of large N c ) and the number of domains Q 2 s S ⊥ , where S ⊥ denotes the transverse area probed by the projectile. We will now discuss this physical picture in in more detail and further develop its quantitative implementation along the lines of the discussion in Ref. [13].

A. Single quark scattering
Within the CGC formalism, the color fields inside the target nucleus are determined by the solution of the classical Yang-Mills equations where the eikonal current J µ is given in terms of the density of color charges ρ inside the target nucleus as The solution to the classical Yang-Mills equations takes the well known form [50] A where The scattering of an incoming quark inside the projectile can be described to leading order accuracy in α s by the solution of the Dirac equation in the presence of the background field of the target in Eq. (3). One finds that the forward scattering amplitude of a quark with momentum p to scatter off the color fields in the target is given by 1 where denotes the Wilson line at a spatial position x in the fundamental representation. Within the leading order dilute-dense framework, it is then straightforward to compute the single inclusive distributions of quarks in a high energy projectile after scattering from a nuclear target whereρ is the reduced one particle density matrix in the probe. This general expression can be rewritten explicitly as 2 1 Since we are primarily interested in the transverse coordinate dependence, we have omitted a delta function for longitudinal momentum conservation as well as the spin structure to lighten the notation. We refer to [51] for the complete expression. 2 Our expression generalizes the one given in [51] by replacing the collinear quark distribution with a Wigner function Wq(b, k) that is a function of both the k of quarks in the projectile and their impact parameter b. Equivalent expression for gluons, differing only by the representation of the Wilson lines, have explicitly been derived in [31,52].
The Wigner function W q (b, k) characterizes the transverse momentum and position distribution of incoming quarks inside the projectile and is defined to be (9) For the illustrative purpose of this paper, it is sufficient to choose a Gaussian form with a dimensionful constant B characterizing the transverse area of the projectile. The dynamics of interest to us is given by the unintegrated gluon distribution of the target nucleus which represents the distribution of momentum transfers from the target that contribute to the corresponding momentum distribution of the scattered quark. Here is the expectation value of the dipole operator which results from the product of the forward scattering amplitude in Eq. (5) and its complex conjugate equivalent, to obtain the single inclusive probability [53][54][55]. Written fully in coordinate space our expression for the single inclusive multiplicity (8) has the form which is recognizable as the one used in [13] up to a constant factor which cancels in the anisotropy coefficients v n . Since we are interested in the scattering of a small probe -such as a quark inside a proton -off a large nucleus, it is reasonable to neglect the impact parameter b = (x + y)/2 dependence in the target. Because on average there is no preferred direction in the transverse plane, the expectation value D(x, y) then depends only on the magnitude of the transverse separation r = x − y of the dipole 3 , 3 The dipole expectation value can in general depend on the rota-Since the dipole coordinate r is the conjugate variable to the momentum transfer p − k for a single quark scattering, the symmetry in Eq. (15) ensures that the momentum transfer to each individual quark is on average symmetric with respect to the azimuthal angle.

B. Double inclusive spectrum and multiparticle correlations
We will now consider the case where two quarks scatter independently off the same nucleus and shall study the correlations between the two scattered quarks. We will make the simplifying assumption that the momenta of the two incoming quarks are initially uncorrelated 4 -the two particle distribution W qq (b 1 , k 1 , b 2 , k 2 ) of incoming quarks factorizes into the product of single quark distributions, The double inclusive distribution of scattered quarks then takes the form While we assumed the transverse momenta k 1 and k 2 of the two incoming quarks to be uncorrelated, one observes from Eq. (17) that this is no longer the case for the momenta p 1 and p 2 of the scattered quarks.
Since both quarks scatter off the same nucleus, the momentum transfers p 1 − k 1 and p 2 − k 2 are correlated tional invariants |b|,|r| and b · r. Once one performs the averaging over the b distribution of the projectile, only the dependence on r remains. We note however, that such an averaging can effectively introduce non-Gaussian correlations between several dipoles. These can affect multiparticle correlations at a characteristic momentum scale given by the inverse impact parameter dependence. 4 We are mostly concentrating on the near-side "ridge" correlation for semihard momenta ∼ Qs. Thus we are neglecting back-toback momentum correlations [56,57] that are particularly important for the away-side "jet" peak [3,58,59] and contributions at the nonperturbative small intrinsic transverse momentum scale of the probe. Such correlations should be taken into account in a full comparison with data.
with each other, giving rise to azimuthal correlations of the scattered quarks.
We note that the above model can be generalized in a straightforward way to study correlation functions involving more than two particles. Such higher order correlation functions can be related to higher order correlation functions of dipole correlators. For example, four quark correlations in this model involve expectation values of products of four dipoles in the fundamental representation.

III. DIPOLE-DIPOLE CORRELATOR
The discussion in the previous section shows that all the features of two-particle correlations are encoded in the expectation value of the dipole-dipole correlator D(x, y)D(u, v) . We will now study the properties of this correlator and discuss approximation schemes that have been frequently employed in the literature.

A. Glasma graph approximation
The basic properties of the two particle correlation can be understood by studying the interaction between the incoming quark and the target nucleus in terms of multiple gluon exchanges. Clearly the dynamics of these gluon exchanges depends on the nature of color charge correlations in the target and needs to be specified. A simple model of such correlations is the MV model [41][42][43]. In this model, the underlying distribution of color charges in the nucleus is assumed to be a Gaussian distribution such that all multigluon correlations are uniquely determined by the two gluon correlation function.
A further approximation that simplifies the computation considerably is to assume that each of the quarks in the projectile exchanges only two gluons with the target nucleus. We will refer to this combination of the two gluon exchange approximation and Gaussian statistics as the Glasma graph approximation, a term first introduced in Ref. [2]. This approximation has been used in a number of phenomenological studies of ridge correlations in high energy collisions [1,2,4,5,[25][26][27][28][29].
More specifically, the Glasma graph approximation can be understood by examining the expression for the dipole-dipole correlator in terms of the color field. When the density of color charges in the target nucleus gρ is small -corresponding to a dilute-dilute situation -one can perform an expansion of the path ordered Wilson line around the identity matrix, representing an expansion in terms of the number of gluons exchanged between the projectile and the target. In order to keep the notation as light as possible, we will denote the Wilson lines as V (x) = exp(−iΛ(x)) in the following 5 , so that the expansion takes the form Evaluating the color traces for Λ( the dipole operator to O(Λ 3 ) within this approximation takes the form In the Gaussian approximation for the correlations of Λ the C and P odd O(Λ 3 ) term vanishes in the expectation value of the dipole operator and only the O(Λ 2 ) remains. We will denote the ΛΛ correlator as which defines the correlation function γ(r). With this definition, we obtain [53] D Following the same logic as previously for the dipole expectation value, we will now compute the dipole-dipole correlator by expanding up to O(Λ 4 ) in the coupling constant. The assumption of Gaussian statistics means that the four point correlation function can be expressed in terms of the two point function in Eq. (21) as

Thus one obtains
5 Careful path ordering will not modify the results obtained as long as the correlations are local in x + . which, to this order of the approximation, is consistent with One observes from Eq. (24) that the leading N c contribution to the dipole-dipole correlator factorizes into the product of single inclusive averages -corresponding to the independent scattering of two quarks. Diagrammatically this corresponds to the disconnected contribution in Fig. 2. Genuine correlations are contained in the second term of Eq. (24) and suppressed by a factor of 1/(N c 2 −1) as pointed out previously in [25,27,28]. These correspond diagrammatically to the connected graphs in Fig. 2 and feature all possible contractions between the x, y and u, v dipoles.
When calculating two particle correlations using Eq. (17), one needs a Fourier transform of the dipoledipole correlator in Eq. (24) to obtain the momentum transfer p − k from the interaction with the target.
One finds that the terms in Eq.(24) that do not depend on all four coordinates x, y, u, v only contribute for very soft momenta p, k on the order of the incoming momentum of the projectile p T , q T ∼ 1/ √ B. Since we are focusing on the dominant semihard gluons with p T , q T Q s 1/ √ B we will neglect these contributions in the following and replace Eq. (24) with Equation (25) is in fact the form used in the Glasma graph papers [1,2,4,5,[25][26][27][28][29]. In particular, in the comparisons to data performed in [2,4,5], the expectation value of the dipole correlator is taken to satisfy the Balitsky-Kovchegov (BK) equation [60,61]. Thus the Glasma graph approximation, as employed in phenomenological computations, corresponds to the merging of gluon ladders from multiple sources (all localized on a transverse scale ∼ 1/Q s ) into a single ladder exchange represented by the dipole correlator.

B. Nonlinear Gaussian approximation
If we restrict ourselves to Gaussian correlations of gluon fields in the target, it is possible to resum the multiple gluon exchange contributions to the dipole amplitude and the dipole-dipole correlator analytically to all orders 6 . Generalizing the notation from the glasma graph case, one obtains [53] This expression is a part of the nonlinear Gaussian approximation because it includes all orders in Λ, evaluated with a the Gaussian ΛΛ correlator in Eq. (21). This is to be contrasted with the two-gluon exchange approximation Eq. (22) which was expanded to the lowest order. Using a well known algorithm [31, 32, 34-36, 38, 38-40, 55] for computing higher point correlators one can obtain the dipole-dipole operator expectation value to all orders in Λ assuming a Gaussian ΛΛ correlator as [62] D where Here ∆ is short for ∆(x, y; u, v) and is given by We will refer to Eq. (27) as the "nonlinear Gaussian approximation" for the dipole-dipole correlator. It is exact in the McLerran-Venugopalan (MV) model where dipole and multipole correlators depend nonlinearly on the Λ fields, but the correlators of Λ's are Gaussian. One can easily check that the two-gluon exchange limit of the nonlinear Gaussian (27) is the same as the two-gluon exchange approximation (23).
One can obtain sight into this relation by taking the large N c limit for a constant D(x − y): which shows again that the leading N c contribution corresponds to the independent scattering of two quarks, while genuine correlations are N c suppressed. While the nonlinear Gaussian approximation reproduces the double gluon exchange (Glasma graph) approximation at O(Λ 4 ) in the dilute limit, it also contains a series of higher order terms. An important subset of higher order contributions corresponds to the diagrams which separately break the x − y → y − x and u − v → v − u symmetries. As we will discuss shortly these contributions are responsible for generating the odd moments (v 3 , v 5 , ...) in the Fourier expansion of the correlation function. One finds that the leading contribution to the x ↔ y antisymmetric part can be associated with the square of the C and P odd contribution to the dipole operator in Eq. (20). Note that while the expectation value of the odd term is zero in the Gaussian approximation, the expectation value of its square is not, but is The contributions from the odd terms correspond diagrammatically to the processes depicted in Fig. 3 and take the form While for N c = 2 the antisymmetric contribution vanishes identically to all orders, it is nonvanishing for N c ≥ 3 and suppressed by a factor of 1/N c 2 relative to the disconnected contribution in the large N c limit.

IV. AZIMUTHAL CORRELATIONS IN QUARK NUCLEUS SCATTERING
We will now discuss the azimuthal correlations of quarks scattering off a large nucleus and present comparisons of results for these correlations within different approximation schemes introduced in the previous sections. To further quantify the correlations introduced by the scattering, we will decompose the double inclusive distribution in Eq. (17) into Fourier modes in the relative azimuthal angle ∆φ between the two scattered quarks, The familiar coefficients v n {2}(p T ) can be obtained from the V n∆ using [63] v where p Ref T represents a range in transverse momentum corresponding to an experimental reference bin.

A. Analytic estimates
Before we turn to a discussion of numerical results, it is useful to obtain further analytic insight into the correlation functions themselves by considering the limit of a large probe with small intrinsic transverse momentum (Q 2 s B 1). With the Gaussian Wigner distribution introduced in Sec. II it is convenient to absorb a part of the Wigner distribution into the definition of a modified dipole distributioñ such that the single inclusive distribution (8) becomes Within the Glasma graph approximation the double inclusive distribution can then be evaluated by combining Eqs. (17) and (25) as which further reduces to in the limit of a large probe or at high momenta, where the intrinsic transverse momentum of the probe can be neglected and we can approximate 2πBe − B 2 k 2 → (2π) 2 δ (2) (k). While the first term inside the square bracket corresponds to the disconnected contribution and does not contain any correlations, the connected term gives rise to a correlation which is suppressed by 1/(N c 2 − 1) and 1/(Q 2 s B). One also observes from Eqs. (36) and (37) that the two particle correlation function has a similar structure as the collinear limit of the Glasma graph computation [1,2,4,5,[25][26][27][28][29]. It features a near side (p ≈ q) contribution as well as one on the away side (p ≈ −q), which in terms of the Fourier decomposition in Eq. (32) give rise to even harmonics v 2 , v 4 , .... Odd harmonics v 3 , v 5 , ... are not present in Eqs. (36) and (37) as the above expressions are manifestly symmetric under p → −p.

B. Numerical results
In order to establish the quality of the different approximations, we will now compare the results for the Fourier coefficients v n to numerical lattice computations that fully evaluate the correlation functions of Wilson lines. We begin with a comparison in the McLerran-Venugopalan (MV) model where the Wilson lines are generated from a Gaussian ensemble of fluctuating color charges. Following the numerical procedure of Ref. [64], a set of Wilson line configurations is generated according to Eq. (6); these are then employed to extract numerically the expectation value of the dipole operator D(r).
With the expectation value of the dipole operator we can then compute the double inclusive spectrum in Eq. (17) using the dipole-dipole correlator in the Glasma graph approximation of Eq. (25) and in the nonlinear Gaussian approximation of Eq. (27). We also compute the double inclusive spectrum directly from the lattice Wilson lines using the procedure described in Ref. [13]. For each of these three different double inclusive spectra, we determine the Fourier coefficients v n using Eqs. (32) and (33). We choose the reference momentum to be p Ref Our results for the Fourier coefficients v n (p T ) in the MV model are shown in the left panel of Fig. 4.
We also perform, as discussed in [13], the JIMWLK rapidity evolution of the Wilson lines for y = 7.6 units in rapidity for SU(3) and for y = 12.4 units for SU(2) with the same running coupling formula and an initial saturation scale Q s /Λ QCD = 3.7 in both cases. We use the running coupling prescription for the JIMWLK equation proposed in Ref. [65]. We then compute again the double inclusive spectrum and the Fourier harmonics directly from the lattice Wilson lines, as well as from the lattice result for the dipole. The results including JIMWLK evolution are shown in the right panel of Fig. 4. For the MV model the probe size is BQ 2 s = 3.7 and for the JIMWLK simulations BQ 2 s = 2.5; for a discussion of the B-dependence see [13].
We note that for the MV model case the nonlinear Gaussian approximation agrees perfectly with the direct numerical calculation for all v n up to the numerically accessible values of p T . This is of course a simple numerical check of the analytical expressions and the agreement should be exact because non-Gaussianities are absent by definition in the MV model. In contrast, the Glasma graph approximation deviates significantly from the numerical result. As discussed previously, the Glasma graph result does not have any odd harmonics. However the even Fourier coefficients v 2 and v 4 too show significant deviations from the exact numerical result especially around p T ∼ 2Q s , demonstrating that nonlinearities are significant even at larger p T . The Glasma graph approximation should be exact in the high momentum limit when |p|, |q|, |p − q| and |p + q| are all large. However, even for a large |p| = |q| the azimuthal harmonics receive contributions from |p ± q| Q s where the two- gluon exchange approximation is not very accurate.
With JIMWLK rapidity evolution, the distribution of color charges is no longer explicitly Gaussian. Therefore, although non-Gaussian contributions were not seen in the operators studied previously in [37], it is possible that JIMWLK evolution of the azimuthal anisotropies will introduce non-Gaussian contributions. Indeed this is seen in Fig. 4 (right) where we observe a deviation of the nonlinear Gaussian approximation from the numerical result from solving the JIMWLK equations. Furthermore we see that both the numerical and the nonlinear Gaussian results for all v n are reduced by the JIMWLK evolution. On the contrary, the Glasma graph results, after JIMWLK rapidity evolution of the Wilson lines, are roughly the same as those for the MV model. Thus while better agreement of the Glasma graph approximation with the nonlinear Gaussian and full JIMWLK results is seen after rapidity evolution, this agreement is a fortuitous numerical coincidence.
We have also analyzed the N c dependence of our results for v n employing the full numerical calculation of the dipole-dipole correlator. Our results for the SU(2) and SU(3) gauge theory are shown in Fig. 5, where we scale the Fourier coefficients v 2 and v 4 by the color factor N c 2 − 1. This is because azimuthal correlations in the double inclusive spectrum contain an overall factor of 1/(N c 2 − 1) (see e.g. Eq. (37)) and the v n 's are related via a square root to the Fourier coefficients V n∆ in the expansion of the double inclusive spectrum. In fact, we find that this scaling works nearly perfectly both in the MV model case (Fig. 5 (left)) and the JIMWLK evolved case (Fig. 5 (right)); in the former, small differences are seen in v 2 and v 4 for only for large p T , where lattice cutoff effects can already have an effect.
Finally, we demonstrate the dependence of our results on different choices for the reference transverse momen- tum in Fig. 6. This is an interesting exercise because similar studies can be performed with the experimental data and will help to distinguish between different models. We find a clear suppression of the signal to the previously regarded case p Ref T = p T when employing a fixed reference bin 0.5 Q s < p Ref T < 3 Q s . One observes that for the Gaussian correlations in the MV model, this effect is particularly strong at large p T in the Glasma graph approximation. In case of the JIMWLK-evolved results, all approximations show a similarly strong suppression of the signal at large p T when using the stated fixed reference momentum bin.
The decorrelation in p T observed in Fig. 6 is fairly fast and appears incompatible with the experimental observations. Experimentally only a slow decorrelation can be seen in the data when comparing experimental results for v n in p+Pb collisions using different methods [63,66]. However a number of caveats are in order with regard to this comparison. We emphasize that our results are for quarks or more generally on the parton level. While hadronization effects will weaken the strong dependence on the choice of reference momentum observed on the parton level, a quantitative description of the experimental data in initial state frameworks will also be quite sensitive to the choice of the fragmentation scheme. The role of fragmentation in such correlations deserves a more detailed study in the future (see also [67][68][69]).

V. THE CGC AND THE COLOR FIELD DOMAIN MODEL
We will now discuss the relation of the azimuthal correlations derived in the CGC framework to those computed recently in the color field domain model introduced in [10,11,48,49]. Since the latter qualitatively describes some key features of the ridge data in proton-lead colli-sions at the LHC, it is interesting to compare and contrast this model with the CGC based calculations we discussed thus far.

A. Electric fields in the Glasma graph approximation
The color field domain model is usually formulated in terms of transverse color electric fields and their correlators. To achieve an "apples-to-apples" comparison, we will first show how our previous discussion in terms of dipole-dipole correlators can be formulated in terms of color electric fields.
The classical color electric fields in the target, as shown in Fig. 1, can be expressed in terms of the Wilson line correlators as Performing a short distance expansion of the dipole operator, one obtains where we denote r = x − y and b = (x + y)/2. Within the weak field limit of the Glasma graph approximation outlined previously, one can evaluate the correlator of color electric fields as where we have used Eq. (21). We emphasize however that this equivalence is valid only in the combined limit of weak fields and short distances. Using the electric field correlator then yields the following expression for the dipole operator One can similarly use Eq. (39) and express the dipoledipole correlator in the short distance limit as This expression shows that the expectation value of the dipole-dipole correlator is sensitive to fluctuations of the color electric fields as characterized by the four point correlator.
In the Glasma graph approximation, the electric field is linearly proportional to the charge density and thus has Gaussian correlations The four point correlation function can be expressed in terms of the two point function as Combining the above expressions, the dipole-dipole correlator takes the form which agrees precisely with the expansion of the result in the double gluon exchange approximation in Eq. (24) in the short distance limit |r 1 | ∼ |r 2 | |b 1 | ∼ |b 2 |.
This simple calculation shows that the physics of fluctuating color electric field domains is implicitly contained in the conventional Glasma graph picture. While in the short distance (large momentum) limit the dipole-dipole correlator is expressed in terms of two and four point correlators of electric fields, there is a one-to-one mapping between the statistical properties of these electric fields and those of the Λ's in the Glasma graph calculation.

B. The color field domain model and non-Gaussian correlations
We focused thus far on conventional models based on Gaussian correlations of color fields inside a large nucleus. It is now interesting to understand how these relate to the color field domain model [10,11,48,49]. In our language, the color field domain model is obtained by replacing the electric field correlator in Eq. (40) by This correlator in the color field domain model depends explicitly on the effective degree of polarization A and the unit vectorâ characterizing the direction of the color electric field. Expectation values of operators within the color field domain model are computed by a two step averaging procedure. In the first step, one performs a Gaussian average with the modified two point correlation function in Eq. (46). The second step consists of an average over all possible directions of the chromoelectric fields,â, such that and â iâjâkâl â = δ ij δ kl + δ ik δ jl + δ il δ jk 8 .
Implicit in this two step procedure is a physical assumption about time scales. One assumes that partons within a color domain of size ∼ 1/Q s generate a color electric field oriented in a particular directionâ with a likelihood A ranging from 0 − 100% , which is long lived on the time scale of the interaction such that partons in the projectile are collimated relative to the directionâ of this color electric field. A similar picture is implicit in the work of [12] where the net momentum transfer from the projectile partons to the target takes the place of the color electric field.
When the life time of the degrees of freedom responsible for the breaking of rotational symmetry within each color domain is much larger than the time scale over which one performs the average over the color electric fields inside the target in Eq. (46), the orientation of the domain appears frozen on that time scale and a separate averaging is justifiable. However it is not a priori evident that such a separation of time scales exists. In particular, it is not clear what would be the intrinsic or dynamical scale that separates the time scale over which the color electric fields align themselves from the time scale over which one performs the average over the different orientations of the field and why such an average should be a Gaussian average.
Let us now discuss explicitly the calculation of the four point correlator of color electric fields in the color field domain model. The only correlators that are related to physical observables are the ones averaged over all unobservable degrees of freedom, including the direction ofâ. Therefore, to understand the correlation structure of the model we must compare the two-and four-point functions of the electric field after the full two step average.
Carrying out the two step averaging procedure for two point correlators by averaging Eq. (46) overâ using Eq. (47) one obtains The result is independent of the effective polarization A and the normalization has been chosen in such a way that the two-point function Eq. (48) agrees with the previous result in Eq. (40) for a rotationally invariant system in the short distance limit 7 .
The short distance expansion of the dipole-dipole correlator in Eq. (42) involves the (double) average of the four point correlation function which, after the first step, can be expressed as Performing the second average with Eq. (48), one obtains . 7 Note that for a rotationally invariant correlation function γ(r) in the short distance limit r → 0 we can replace Comparing the two point and four point correlation functions of the color electric fields in Eqs. (48) and (51), one sees explicitly that the four point correlation function can not be expressed in terms of the two point function, i.e. the color field domain model is non-Gaussian. Specifically one finds that the dipole-dipole correlator is a sum of a Gaussian piece (present already in Eq. (45)) and non-Gaussian terms proportional to A 2 induced by the two step averaging procedure. The non-Gaussian terms are referred to as "disconnected contributions" in Ref. [48]. In addition to the small dipole limit r 2 1 , r 2 2 1/Q 2 s that is assumed in this discussion, one can additionally take the limit b 2 1 , b 2 2 1/Q 2 s . In this case the two-gluon exchange approximation becomes exact and thus the Glasma graphs and the nonlinear Gaussian are equivalent to each other. Even in this limit the A-terms are not suppressed in any way and the color field domain model remains different from the other approaches considered in this paper.
It is obvious that the (r 1 · r 2 ) 2 terms introduce an additional A-dependent cos 2φ correlation between the two dipoles. Upon Fourier transformation, this modifies the angular structure of the correlation between the two produced particles. What remains unclear at this stage is the physical origin of this particular form of non-Gaussian correlators as well as the magnitude of the non Gaussianity characterized by the additional A parameter in this model.

C. Interpretations of the color field domain model
We noted previously but wish to emphasize again that fluctuating domains of color electric field are present also in the Glasma graph or non-linear Gaussian approximation and they are not a new physical feature added by the color field domain model. What is different in the color field domain model is that the direction of the chromoelectric field is treated explicitly as a long lived degree of freedom. By modifying the correlation function of electric fields according to Eq. (46) and performing a separate average over the orientation of the electric fields the statistics of these domains is altered significantly. Most importantly, this can lead to sizable non-Gaussian correlations depending on the magnitude of the parameter A in this model. Since the single inclusive distribution is not sensitive to A, its value can only be determined from correlation measurements. An observable that would be particularly sensitive to the presence of intrinsic non-Gaussianities would be the four-particle cumulant flow coefficient, as discussed in [48]. We will now discuss three possible interpretations of the non-Gaussian correlations represented by the A-term in the color domain model. 1. The color electric field is a nonlinear function of the color charge density. Thus even if the color charges (or, equivalently the Λ's) have a Gaussian distribution, the color electric fields and thus the dipole operators can have non-Gaussian correlations. One possible interpretation of the color field domain model is as an effective way to account for the non-linear relation between electric fields and color charge densities in the target in an otherwise linearized calculation. The Glasma graph approximation of Eq. (38) assumes that the electric field is linearly proportional to the color charge and thus the electric fields have Gaussian correlations. The nonlinear Gaussian approximation, on the other hand, has Gaussian correlations for color charges but not for the electric fields. In this interpretation of the color field domain model, the corrections encoded in A should be proportional to the difference between the Glasma graph and the nonlinear Gaussian computations in this paper. The total anisotropy of the azimuthal two-particle distribution calculated from the MV model (see e.g. [11]), is then the sum of the Glasma graphs and the Aterm. In this interpretation one has parametrically A 2 ∼ 1/(N c 2 − 1), since correlations are N csuppressed (relative to the uncorrelated term) in both the nonlinear Gaussian and the Glasma graph approximations. Our numerical results in Fig. 4 show that such non-linear corrections can indeed be sizeable and should be taken into account as a correction to the Glasma graph result. However, if this is the interpretation it would seem more natural to directly use the non-linear Gaussian approximation rather than introducing an additional parameter. In particular, it is not obvious whether a constant A could have a similar momentum dependence as the nonlinear Gaussian approximation, given that the color field domain model differes from the nonlinear Gaussian even in the small distance limit.

A second possible interpretation of the A-term is
that it represents non-Gaussian correlations that can emerge from JIMWLK rapidity evolution even when starting from a Gaussian initial condition. This contribution is, in our present calculation, represented by the difference between the full JIMWLK result and the nonlinear Gaussian. Indeed we see signs of a ∼ 10% deviation between the two for p T Q s . However, this difference is relatively small in practical terms and might not have a significant influence on phenomenology. From a theory perspective, it is to our knowledge the first instance observed in the literature of a meaningful breaking of the Gaussian approximation to JIMWLK. We see no obvious reason why the deviation from Gaussianity seems larger here than in the observables studied in Ref. [37]. This issue might call for additional studies in the future including a more systematical check of discretization effects in the lattice calculations.
3. Finally the most intriguing possibility is that the A-term represents an intrinsic non-Gaussian correlation that is present in the initial condition for JIMWLK evolution and survives substantially after evolution. This is the interpretation suggested in [10,48]. The possible existence of such non-Gaussian correlators was previously suggested in [6,7] and later studied in [70,71]. While the Gaussian MV model can be justified on quite general grounds [72] as arising, due to the central limit theorem, from a superposition of a large amount of uncorrelated color charges in a heavy nucleus, deviations from Gaussian statistics are naturally expected for a small number of large x degrees of freedom. The existence and persistence of such a non-Gaussianity at small x would thus be a signal of remarkably strong long range rapidity correlations inside the gluon cascade building up the strong color fields at small x. The computations in this paper do not address this possibility of an intrinsic non-Gaussian four-particle correlation because we have been working in the MV model+JIMWLK evolution setup where such correlations are absent in the initial conditions.

VI. SUMMARY AND CONCLUSIONS
We explored a number of calculational schemes to compute two particle correlations of quarks scattering off a highly energetic nucleus. All cases correspond to different approximations within the dilute-dense limit of the color glass condensate framework. The two-particle correlations are quantified in terms of Fourier coefficients in an expansion in relative azimuthal angle of the double inclusive distribution of scattered quarks. This distribution is proportional to the dipole-dipole correlator. The study of the properties of this correlator in the various approximation schemes was the primary objective of this work.
The simplest approximation scheme considered was the glasma graph approximation. In this case, the lightlike Wilson lines in the dipole-dipole correlator are expanded to lowest order, restricting the interaction with the target to two gluon exchange. One further assumes that the gluon correlations in the target are Gaussian correlations. This approximation scheme has been used previously in the literature to study azimuthally collimated double in-clusive gluon production in p+p and p+Pb collisions.
Another approximation scheme, of greater complexity, is the nonlinear Gaussian approximation. In this case, all multigluon exchanges are resummed to all orders to obtain a complicated analytical expression for the dipoledipole correlator. This expression is exact as long as there are only Gaussian correlations in the target. This is for instance the case for the MV model.
These analytical results are a good benchmark for numerical studies wherein the Wilson lines are computed on 2+1-dimensional lattices for Gaussian distributed sources -good agreement is expected and achieved. With the MV initial conditions for the Wilson lines at a given rapidity, the JIMWLK equations are solved on the lattice to determine the Wilson lines at larger rapidities. These then allow one to determine in principle expectation values of n-dipole correlators as a function of rapidity.
We studied the Fourier harmonics v 2 , v 3 , and v 4 that are extracted from azimuthal two particle correlations in the various approximation schemes. The Glasma graph approximation and the nonlinear Gaussian approximations differ appreciably for p T ∼ 2Q s in the MV model, indicating the importance of coherent multiple scattering effects. Since the Glasma graph approximation is at the heart of most comparisons to experimental data, this calls for a more detailed study to further quantify its theoretical uncertainties. However, we believe that in terms of the phenomenological consequences most of this difference can be accomodated within the uncertainties in the overall normalization of the Glasma graph calculations [1][2][3][4][5]. A significant difference between the two approximation schemes is that the symmetry constraints inherent in the Glasma graph approximation do not produce any odd harmonics in contrast to the nonlinear Gaussian approximation which generates all odd harmonics of the azimuthal double inclusive distribution. With JIMWLK rapidity evolution, the differences in the v 2,4 coefficients computed in the Glasma graph and nonlinear Gaussian schemes decrease significantly. Furthermore, we find that the coefficients in the two schemes are quite close to those computed by solving JIMWLK numerically without any approximation to the dipole-dipole correlator. Since we do not currently have a good interpretation for this better agreement, we believe that it may to some extent be accidental.
We analyzed the dependence of our results on the number of colors by comparing computations for SU(3) and SU(2) gauge fields and found precisely the expected scaling of V n∆ (p T ) with 1/(N c 2 − 1). This result confirmed that the azimuthal angle dependent correlations are suppressed parametrically by 1/N c 2 . We studied the dependence of the Fourier coefficients on the reference transverse momentum p Ref Q s , with the Glasma graph approximation showing the largest differences. JIMWLK rapidity evolution seems to increase the difference between the different p Ref T choices. We note however, that the choice of the fragmentation scheme can qualitatively influence the comparison of model computations of gluon correlations to the hadron correlation data. This topic deserves a more detailed study in the future.
Finally we analyzed in detail the relation of the color glass condensate based computations to a color domain model which captures qualitative features of the multiparticle azimuthal correlations observed in protonnucleus collisions. In this framework, the dipole-dipole correlator is modified to include an additional term that models the polarization of gluon fields in individual domains of color charge within the target. We conclude that this term, proportional to the polarization parameter A, introduces non-Gaussian correlations amongst the color electric fields inside the target nucleus. On the other hand, if such non-Gaussianities are not explicitly introduced, the color domain model reduces to the MV model in the Glasma graph approximation.
We also discussed possible origins of non-Gaussian correlations of the color fields of a large nucleus. One possibility is that these correspond to non-Gaussian correlations induced by JIMWLK rapidity evolution of Gaussian correlations at the initial rapidity. However our estimates of this effect suggest that such non-Gaussian correlations are too small to be relevant phenomenologically. A more interesting possibility is that the A polarization term introduced in this model arises from intrinsic four point correlations that are significant in the initial condition and whose magnitude is preserved with rapidity evolution. Such correlations would be interesting to study in the future.