JIMWLK evolution in the Gaussian approximation

We demonstrate that the Balitsky-JIMWLK equations describing the high-energy evolution of the n-point functions of the Wilson lines (the QCD scattering amplitudes in the eikonal approximation) admit a controlled mean field approximation of the Gaussian type, for any value of the number of colors Nc. This approximation is strictly correct in the weak scattering regime at relatively large transverse momenta, where it reproduces the BFKL dynamics, and in the strong scattering regime deeply at saturation, where it properly describes the evolution of the scattering amplitudes towards the respective black disk limits. The approximation scheme is fully specified by giving the 2-point function (the S-matrix for a color dipole), which in turn can be related to the solution to the Balitsky-Kovchegov equation, including at finite Nc. Any higher n-point function with n greater than or equal to 4 can be computed in terms of the dipole S-matrix by solving a closed system of evolution equations (a simplified version of the respective Balitsky-JIMWLK equations) which are local in the transverse coordinates. For simple configurations of the projectile in the transverse plane, our new results for the 4-point and the 6-point functions coincide with the high-energy extrapolations of the respective results in the McLerran-Venugopalan model. One cornerstone of our construction is a symmetry property of the JIMWLK evolution, that we notice here for the first time: the fact that, with increasing energy, a hadron is expanding its longitudinal support symmetrically around the light-cone. This corresponds to invariance under time reversal for the scattering amplitudes.


Introduction
The final state of an ultrarelativistic hadron-hadron collision, as currently explored at RHIC and the LHC, is characterized by an extreme complexity in terms of the number and the distribution of the produced particles. The study of multiparticle correlations represents an essential tool for organizing this complexity and extracting physical information out of it. In particular, a recent measurement at RHIC of di-hadron correlations in deuteron-gold collisions [1] revealed an interesting phenomenon -the azimuthal correlations are rapidly suppressed when increasing the rapidity towards the fragmentation region of the deuteron -, which is qualitatively [2][3][4][5] and even semi-quantitatively [6,7] consistent with the physical picture of gluon saturation in the nuclear wavefunction. For this interpretation to be firmly established, one needs a more precise understanding of the multi-particle correlations in the high-energy scattering and, in particular, of their evolution with increasing rapidity. This triggered new theoretical studies [8][9][10][11][12] of many-body correlations in the color glass condensate (CGC), which is the QCD effective theory for high-energy evolution and gluon saturation, to leading logarithmic accuracy at least.
The central ingredient in the CGC theory is the JIMWLK (Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov, Kovner) equation [13][14][15][16][17][18][19][20], a functional renormalization group equation of the Fokker-Planck type which describes the non-linear evolution of the gluon distribution in the hadron wavefunction with increasing rapidity, or decreasing the gluon longitudinal momentum fraction x. When applied to an asymmetric, 'dilute-dense', scattering (like a pA collision), all orders. But within the context of the JIMWLK evolution, such Gaussian approximations seem to be prohibited by the highly non-linear structure of the evolution equation, which is the mathematical expression of gluon saturation.
In spite of this theoretical prejudice, the numerical results in Ref. [11] suggest that a Gaussian approximation to the JIMWLK evolution may nevertheless work. Another piece of evidence in that sense emerges from the recent analytic study in Ref. [12]. There, we have constructed an approximate version of the Balitsky-JIMWLK hierarchy which is simple enough to allow for explicit solutions. Then we have showed that, for the special configurations of the quadrupole considered in Ref. [11], these approximate solutions coincide with the respective predictions of the MV model extrapolated to high energy. But in that context too, the similarity with the MV model appears as merely an 'accident', with no deep motivation: the simplified hierarchy proposed in Ref. [12] is generated by the 'virtual' piece of the JIMWLK Hamiltonian, which is non-linear in the target field and therefore seems incompatible with a Gaussian approximation. Moreover, the approximations in Ref. [12] have been justified only in the limit where the number of colors N c is large (formally, N c → ∞). This does not explain the observation in Ref. [11] that the numerical solutions to the JIMWLK evolution for N c = 3 are better reproduced by the finite-N c version of the MV model (with N c = 3, of course) than by its large-N c limit.
Our purpose in the present analysis is to clarify such 'coincidences' and 'apparent contradictions' by resolving the aforementioned tensions between the simplified hierarchy proposed in Ref. [12], the Gaussian approximation, and the large-N c limit. The results that we shall obtain can be summarized as follows. We shall demonstrate that the JIMWLK equation admits indeed an approximate Gaussian solution for the CGC weight function, that this solution is unique within the limits of its accuracy, and that it is tantamount to a simplified system of evolution equations, which are linear (while being consistent with unitarity) and local in the transverse coordinates. In the limit where N c → ∞, these new equations reduce to those previously proposed in Ref. [12]. The ultimate outcome of our analysis is a global approximation to the Balitsky-JIMWLK hierarchy, which is valid for any N c and allows one to construct explicit, analytic, solutions for all the n-point functions of the Wilson lines. These approximate solutions are strictly correct in the limiting regimes at very large (k ⊥ ≫ Q s (Y )) and, respectively, very small (k ⊥ ≪ Q s (Y )) transverse momenta, and provide a smooth (infinitely differentiable) interpolation between these limits. Here, Q s (Y ) denotes the saturation momentum in the target at a rapidity Y equal with the rapidity separation between the target and the projectile.
To describe our results in more detail, let us first explain the distinction between 'real' and 'virtual' terms in the Balitsky-JIMWLK equations. The 'real' terms describe the evolution of the projectile via the emission of small-x gluons, whereas the 'virtual' terms express the probability for the projectile not to evolve, i.e. not to radiate such (small-x) gluons. The 'virtual' terms dominate the evolution in the approach towards the unitarity (or 'black disk') limit, since in that regime the scattering is strong and the projectile has more chances to survive unscattered if it remains 'simple' -i.e., if it does not evolve by emitting more gluons. By the same token, the 'virtual' terms control the evolution of the many-body correlations which, within the context of JIMWLK, are built exclusively via non-linear effects (multiple scattering and gluon recombination) in the regime of strong scattering. More precisely, the 'real' terms are important for that process too -they include the non-linear effects responsible for unitarity and saturation -, but deeply at saturation their role becomes very simple: they merely prohibit the emission of new gluons with low transverse momenta k ⊥ Q s (Y ). Thus, one can follow the evolution of correlations at saturation by keeping only the 'virtual' terms in the Balitsky-JIMWLK equations, but supplementing them with a phase-space cutoff which expresses the effect of the 'real' terms. (This is strictly correct in a 'leading-logarithmic approximation' to be detailed in Sect. 3.4.) Moreover, since the simplified equations thus obtained are linear, they can be extended to also cover the BFKL evolution in the weak scattering regime at k ⊥ ≫ Q s (Y ). Indeed, in that regime and to the accuracy of interest, the n-point functions of the Wilson lines reduce to linear combinations of the dipole scattering amplitude, with the latter obeying the BFKL equation. The BFKL dynamics involves both 'real' and 'virtual' terms, but it can be effectively taken into account by tuning the kernel in the 'virtual' terms -namely, by requiring this kernel to approach the solution to the BFKL equation at large k ⊥ .
The above considerations, to be substantiated by the subsequent analysis, explain why it is possible to approximate the Balitsky-JIMWLK equations by simpler equations which are linear and whose overall structure is inherited from the 'virtual' terms in the original equations. Similar considerations have underlined our previous construction in Ref. [12], but their generalization to finite N c (that we shall provide in this paper) turns out to be highly non-trivial.
Another subtle aspect of our present analysis is the recognition of the fact that the simplified equations that we shall propose (for either finite or infinite N c ) correspond to a Gaussian approximation for the CGC weight function. A priori, the association of a linear system of equations with a Gaussian approximation may look natural, but in the present case this is complicated by the fact that, as alluded to before, the 'virtual' piece of the JIMWLK Hamiltonian is non-linear in the target field to all orders. Such a non-linear structure seems to preclude any Gaussian solution. The resolution of this mathematical puzzle turns out to be interesting on physical grounds, as it sheds new light on the physical picture of the JIMWLK evolution. Namely, we shall show that the Wilson lines within the 'virtual' Hamiltonian do not represent genuine non-linear effects associated with saturation, rather they express the physical fact that, with increasing energy, the longitudinal support of the target expands symmetrically around the light-cone. That is, in contrast to a widespread opinion in the literature (see e.g. [28,[36][37][38]40]), which was based on a misinterpretation of the mathematical structure of the JIMWLK equation, the gluon distribution in the target expands simultaneously towards larger and respectively smaller values of x − , in such a way to remain symmetric around x − = 0. (We assume the target to propagate along the positive x 3 axis at nearly the speed of light and we define ) In turn, this symmetry has physical consequences for the multi-partonic scattering amplitudes: it implies that the n-point functions of the Wilson lines with n ≥ 4 obey a special permutation symmetry -the mirror symmetry -which expresses their invariance under time reversal.
To summarize, a Gaussian weight function which is symmetric in x − and whose kernel is energy-dependent and interpolates between the solution to the BFKL equation at high transverse momenta k ⊥ ≫ Q s and the JIMWLK (or 'dipole') kernel at low momenta k ⊥ ≪ Q s , provides a reasonable approximation to the JIMWLK equation, which is strictly correct in the limiting regimes alluded to above (for any value of N c ). Within its limits of validity, this approximation is essentially unique: different constructions for the kernel can differ from each other only in the transition region around saturation, which is anyway not under control within the present approximation.
In practice, it is convenient to trade this kernel for the dipole S-matrix, which in turn can be obtained either as the solution to the Balitsky-Kovchegov (BK) equation [21,41], or by solving a self-consistency condition similar to that in Ref. [36]. (The differences between these two expressions for the kernel should be viewed as an indicator of the stability of the approximation scheme.) Then the n-point functions with n ≥ 4 (quadrupole, sextupole etc) can be determined in terms of the 2-point function (the dipole S-matrix) by solving the evolution equations associated with the Gaussian weight function. These equations become particularly simple at large N c , where they reduce to the equations proposed in Ref. [12] and can be explicitly solved for arbitrary configurations of the n external points in the transverse plane.
For finite N c and for generic configurations, the equations are more complicated, as they couple the evolution of the various n-point functions with the same value of n. (For instance, the quadrupole mixes under the evolution with a system of two dipoles.) Yet, explicit solutions can be obtained under the simplifying assumption that the kernel of the Gaussian is a separable function of the rapidity and the transverse coordinates (plus an arbitrary function of Y ; see Sect. 4.2 for details). This is certainly not the case for the actual kernel (say, as given by the solution to the BK equation), but it is a good piecewise approximation to it and it is furthermore true for the MV model 2 , that we shall take as our initial condition at low energy. So, not surprisingly, the expressions for the n-point functions that we shall obtain within this scenario are formally similar to the respective predictions of the MV model. One can reverse this last argument as follows: given that the Gaussian weight function is a good, piecewise approximation to the JIMWLK evolution, as we shall demonstrate, and that the kernel of this Gaussian can be taken to be separable within the relevant kinematical regimes, we expect the predictions of this approximation to be very close to those of the MV model extrapolated to high energy.
For special configurations which are highly symmetric, exact solutions can be obtained at finite N c even without assuming separability. We shall study various examples of this type for the 4-point function and the 6-point function, and thus find some surprising factorization properties, that would be interesting to test against numerical solutions to the JIMWLK equation. For one particular configuration of the 6-point function, the exact numerical result is already known [11] and our respective analytic solution appears to agree with it quite well.

The Balitsky-JIMWLK evolution equations
In this section, we shall briefly review the general formulation of the JIMWLK evolution and then use the evolution equations satisfied by the dipole and the quadrupole S-matrices in order to illustrate various properties of the evolution, which are important for what follows: the role and origin of the 'real' and 'virtual' terms, the factorization of multi-trace observables at large N c , and, especially, the symmetric expansion of the longitudinal support of the target and the ensuing, 'mirror', symmetry of the n-point functions with n ≥ 4.

JIMWLK evolution: a brief reminder
The color glass condensate is an effective theory for the small-x part of the wavefunction of an energetic hadron: the gluons carrying a small fraction x ≪ 1 of the hadron's longitudinal momentum are described as a random distribution of classical color fields generated by sources with larger momentum fractions x ′ ≫ x. Given the high-energy kinematics, in particular the fact that the distribution of the color sources is 'frozen' by Lorentz time dilation, this color field can be chosen (in a suitable gauge) to have a single non-zero component, namely A µ a (x) = δ µ+ α a (x − , x) for a hadron moving along the positive z axis 3 (a 'right mover'). All the correlations of this field are encoded into a functional probability distribution, the 'CGC weight function' W Y [α], which contains information about the evolution of the color sources with increasing 'rapidity' Y ≡ ln(1/x), from some initial value Y 0 up to the value Y of interest. In the high energy regime where α s (Y − Y 0 ) 1 and to leading logarithmic accuracy with respect to the large logarithm Y − Y 0 = ln(x 0 /x), this evolution is described by a functional renormalization group equation for W Y [α], known as the JIMWLK equation. The latter can be given a Hamiltonian form, where H is the JIMWLK Hamiltonian -a second-order, functional differential operator, whose most convenient form for the present purposes is that given in [42] and reads where we use the notation u... ≡ d 2 u . . . to simplify writing, M is the 'dipole kernel', and V † and V are Wilson lines in the adjoint representation: with P denoting path-ordering in x − . The above form of the Hamiltonian is valid only when acting on gauge-invariant functionals of α a , which will always be the case throughout our analysis. In fact, the observables of interest are gauge-invariant products of Wilson lines (see below). The functional derivatives in Eq. (2.2) are understood to act at the largest value of x − , that is, at the upper end point of path-ordered exponentials like that in Eq. (2.4) (see e.g. Eq. (2.11)). These derivatives do not commute with each other, but their commutator is proportional to δ uv (cf. Eq. (2.27)) and thus vanishes when multiplied by M uvz ; hence, there is no ambiguity concerning the ordering of the functional derivatives in Eq. (2.2). One can also notice that the last two terms in the JIMWLK Hamiltonian, i.e. those proportional to V † u V z and to V † z V v respectively, are in fact identical to each other, as it can be checked by exchanging u ↔ v and a ↔ b and by using the property V † ac z = V ca z for color matrices in the adjoint representation.
To fully specify the problem, one also needs an initial condition for Eq. (2.1) at Y = Y 0 ; at least for a sufficiently large nucleus, this initial condition is provided by the McLerran-Venugopalan (MV) model [29,30] (see Sect. 3.2 below). Physical observables, like scattering amplitudes for external projectiles, are represented by gauge invariant operatorsÔ[α] built with the field α a , whose target expectation values are computed via functional averaging with the CGC weight function: (2.5) By taking a derivative in this equation with respect to Y , using Eq. (2.1), and integrating by parts within the functional integral over α, one obtains an evolution equation for the observable, in which the JIMWLK Hamiltonian acts on the operatorÔ[α] : Unlike the JIMWLK equation (2.1), this is not a functional equation anymore, but an integrodifferential equation. However, due to the non-linear structure of the Hamiltonian (2.2) with respect to the field α a , Eq. (2.6) is generally not a closed equation -the action of H onÔ generates additional operators in the right hand side -, but just a member of an infinite hierarchy of coupled equations -the Balitsky-JIMWLK equations. Although mathematically equivalent, the functional equation (2.1) and the Balitsky-JIMWLK hierarchy offer complementary perspectives over the high-energy evolution. Eq. (2.1) depicts the evolution of the target via the emission of an additional gluon with rapidity between Y and Y + dY , in the background of the color field α generated via previous emissions, at rapidities Y ′ ≤ Y . The Wilson lines within the structure of the Hamiltonian (2.2) describe the scattering between this new gluon and the background field, in the eikonal approximation. The Balitsky-JIMWLK hierarchy rather refers to the evolution of the projectile, more precisely, of the operator which describes its scattering off the target. This scattering is again computed in the eikonal approximation, so the operator O is naturally built with Wilson lines -one such a line for each parton within the projectile.

Evolution equations for the dipole and the quadrupole
To be more explicit, we shall consider two specific projectiles: a 'color dipole' made with a quark-antiquark (qq) pair and a 'color quadrupole' made with two qq pairs. In both cases, the overall color state of the partonic system is a color singlet. The S-matrix operators describing the forward scattering of these projectiles off the CGC target read for the color dipole and, respectively, for the color quadrupole. In these equations, V † and V are Wilson lines similar to those in Eq. (2.4), but in the fundamental representation. The results that we shall obtain for these two partonic systems will be easy to extend to projectiles made with n qq pairs, for whicĥ S (2n) As we shall see, within the high-energy evolution, such single-trace operators mix with the multi-trace operators, of the form In order to construct evolution equations according to Eq. (2.6), we need the action of the functional derivatives w.r.t. α a on the Wilson lines. This reads (with δ xu = δ (2) (  [36,40]. By acting with these Hamiltonian pieces on the dipole S-matrix, one finds (withᾱ ≡ α s N c /π).
and respectively (recall that the last two terms in Eq. (2.2), which define H real , are actually identical with each other) where the second line follows after reexpressing the adjoint Wilson line in terms of fundamental ones, according to (2.14) and then using the Fierz identity tr t a A t a B = 1 2 trA trB − 1 2N c tr(AB). (2.15) By adding together the above results, one sees that the terms proportional to 1/N 2 c , that would be suppressed at large N c , exactly cancel between 'real' and 'virtual' contributions, and we are left with This equation has the following physical interpretation: the first term in the right hand side, which is quadratic inŜ and has been generated by the 'real' piece of the Hamiltonian, cf. Eq. (2.13), describes the splitting of the original dipole (x 1 , x 2 ) into two new dipoles (x 1 , z) and (z, x 2 ), which then scatter off the target. More precisely, the evolution step consists in the emission of a soft gluon, so the original dipole gets replaced by a quark-antiquark-gluon system which is manifest in the first line of Eq. (2.13), but in large-N c limit (to which refers the first term in the second line of Eq. (2.13)), this emission is equivalent to the dipole splitting alluded to above. As for the second term in Eq. (2.16), i.e. the negative term linear inŜ which has been produced by H virt , it describes the reduction in the probability that the dipole survive in its original state -that is, the probability for the dipole not to emit. In what follows, we shall often refer to the terms produced by H virt (H real ) as the 'virtual' ('real') terms, but one should keep in mind that not all such terms are actually visible in the evolution equation in their standard form in the literature (to be also used in this paper): some of these terms may have canceled between 'real' and 'virtual' contributions. A similar discussion applies to the evolution equation for the quadrupole, which reads Namely, the terms involving ŜQ Y in the right hand side are 'real' terms describing the splitting of the original quadrupole into a new quadrupole plus a dipole, and have been all generated by the action of the last two terms in the Hamiltonian (2.2). The 'virtual' terms involving Q Y and ŜŜ Y are necessary for probability conservation, and have been generated by the first two terms in the Hamiltonian. Once again, all the terms subleading at large N c (as separately generated by H virt and H real ) have canceled in the final equation. The above features are generic: they apply to the evolution equations obeyed by all the single-trace observables like Eq. (2.9). As visible on Eqs. (2.16) and (2.17), these equations are generally not closed : they couple single-trace observables with the multi-trace ones. E.g., the equation for the quadrupole also involves the 4-point function ŜŜ Y and the 6-point function ŜQ Y , which in turn are coupled (via the respective evolution equations) to even higher-point correlators. The equations obeyed by the multi-trace observables exhibit an interesting new feature: they involve genuine 1/N 2 c corrections, as generated when the two functional derivatives in Eq. (2.2) act on Wilson lines which belong to different traces (see e.g. Appendix F in [43] for an example). At large N c , these corrections can be neglected and then it is easy to check that the hierarchy admits the factorized solution .. , (2.18) provided this factorization is already satisfied by the initial conditions. Then the hierarchy drastically simplifies: it breaks into a set of equations which can be solved one after the other (at least in principle). Namely, Eq. (2.16) becomes a closed equation for Ŝ Y (the BK equation [21,41]), Eq. (2.17) becomes an inhomogeneous equation for Q Y with coefficients which depend upon Ŝ Y [3], and so on. In practice, however, the resolution of these equations is hindered by their strong non-locality in the transverse coordinates. So far, only the (numerical) solution to the BK equation has been explicitly constructed.

The mirror symmetry
In this subsection, we shall discuss a symmetry property of the JIMWLK equation, which has not been noticed in the previous literature and which has far-reaching consequences: the symmetry of the target field distribution (the CGC) under reflection in x − . To start with, we shall identify a mirror symmetry in the evolution equation (2.17) for the quadrupole, that can be easily demonstrated in the large-N c limit, but is likely to hold for any N c . (It does so, at least, in the Gaussian approximation that we shall later construct.) Specifically, if the quadrupole S-matrix Q x 1 x 2 x 3 x 4 Y is symmetric under the exchange of the two antiquark Wilson lines (that is, the Wilson lines at x 2 and x 4 ) at the initial rapidity Y 0 -a condition which is indeed satisfied within the MV model [8] -, then this symmetry is preserved by the evolution. That is, for any Y ≥ Y 0 , one has A similar property holds for the exchange of the quark Wilson lines at x 1 and x 3 , but this is not independent of Eq. (2.19), sinceQ x 3 x 2 x 1 x 4 =Q x 1 x 4 x 3 x 2 by the cyclic symmetry of the trace. To demonstrate Eq. (2.19), let us consider the respective anti-symmetric piece: By using Eq. (2.17), it is easy to see that the associated expectation value obeys the following evolution equation: At large N c , where one can factorize ŜQ asym Y ≃ Ŝ Y Q asym Y , Eq. (2.21) becomes a homogeneous equation which implies that Q asym Y = 0 at any Y provided this condition was satisfied at Y 0 . In turn, this implies the symmetry property (2.19). By inspection of the higher equations in the Balitsky-JIMWLK equations, one can check that a similar symmetry holds for all the n-point functions of the Wilson lines. For instance, the equation obeyed by the sextupole S-matrix Ŝ (6) Y is explicitly shown in Appendix B of Ref. [12]. From this equation, one can read the following symmetry property: The generalization of this property to the 2n-point function shown in Eq. (2.9) reads To better understand the content of this symmetry, it is useful to give a pictorial representation for it. To that aim, consider a generic configuration of the quadrupole operatorQ x 1 x 2 x 3 x 4 in the transverse plane, as illustrated in Fig. 1, and join the four points by oriented lines, which follow the direction of color multiplication. In this way, one constructs a closed, oriented, contour, whose orientation indicates the flow of color within the operator. By repeating this procedure for the 'permuted' operatorQ x 1 x 4 x 3 x 2 , one obtains a similar contour, where however the orientation of the color flow is reversed. One can similarly check that, for a general n-point function, the symmetry property (2.23) refers to changing the contour orientation, say from clockwise to counterclockwise. Such a change would also result from the reflection in a mirror, so we shall refer to the symmetry property (2.23) as the 'mirror symmetry'. Additional arguments in the favor of this name will be given below.
There are several reasons why this this symmetry is so important for us here. First, as we shall shortly argue, this corresponds to an important symmetry property of the scattering amplitudes: their invariance under time-reversal. Second, the way how this symmetry is actually preserved by the JIMWLK evolution is very interesting as it sheds light on the physical picture of the target field distribution: with increasing Y , the color glass condensate expands symmetrically around x − = 0. Third, this symmetry will later guide us in the construction of a mean field approximation to the Balitsky-JIMWLK hierarchy.
To understand the relation to time-reversal, let us present another pictorial representation for the two quadrupole operators which enter Eq. (2.19): this is shown in Fig. 2, where the transverse space is schematically represented by the vertical axis, whereas the horizontal axis refers to x − -the light-cone time for the projectile. The Wilson lines are now explicitly shown, as the oriented horizontal lines extending along the x − axis and connected with each other, via matrix multiplication, at x − → ±∞. Once again, the orientation of these lines corresponds to the direction of the color flow. Clearly, these two figures get exchanged with each other when inverting the arrow of time. In the left figure, corresponding toQ x 1 x 2 x 3 x 4 , the 4-body system starts at x − → −∞ as a set of 2 dipoles, (x 1 , x 2 ) and (x 3 , x 4 ), which then exchange color with each other at x − → ∞ and thus reconnect into the new dipoles (x 1 , x 4 ) and (x 3 , x 2 ). In the right figure, the opposite process happens: the system starts with the dipoles (x 1 , x 4 ) and (x 3 , x 2 ), which then reconnect at x − → ∞ into the dipoles (x 1 , x 2 ) and (x 3 , x 4 ), thus yielding the quadrupoleQ x 1 x 4 x 3 x 2 . Hence, the symmetry property (2.19) corresponds indeed to invariance under time-reversal, as anticipated.
We now turn to the physical interpretation of the mirror symmetry in the context of the JIMWLK evolution. One can check that the symmetric structure of the virtual terms in Eq. (2.17) stems from the combined action of the first two terms in Eq. (2.2). Half of the 'virtual' terms are generated by the first term, proportional to the color unity matrix, but by themselves these terms do not show the mirror symmetry; this symmetry is recovered only after adding the other half of the 'virtual' terms, as generated by the second term in Eq. (2.2), proportional to V † u V v . As an example, consider two of the 'virtual' terms in the r.h.s. of Eq. (2.17) whose coefficients get exchanged with each other under the exchange The first of them is generated when acting with the first term in the Hamiltonian on the pair (x 1 , x 4 ) of the quadrupole, whereas the second one emerges from the action of the second term in H on the pair (x 1 , x 2 ).
Hence, to elucidate this symmetry, one needs to better understand the action of the JIMWLK Hamiltonian. As manifest from Eq. (2.11), the functional derivatives within the Hamiltonian act as generators of infinitesimal color rotations of the Wilson lines at their upper end point in x − ; that is, they act as Lie derivatives for the color group SU(N c ). These color rotations express the evolution of the target color field α a (x − , x) with increasing rapidity: performing one infinitesimal step in the evolution, from Y to Y + dY , amounts to 'integrating out' one layer of quantum fluctuations within the target wavefunction -the gluons with longitudinal momentum fractions between x = e −Y and x ′ = e −(Y +dY ) -and results in adding one additional layer to the classical color field α a (x − , x). The fact that the JIMWLK Hamiltonian acts on the Wilson lines via color rotations at the largest value in x − means that the new layer of color field is located at larger x − as compared to the previous layers.
This argument makes it tempting to conclude that, with increasing Y , the support of the target field α a (x − , x) extends only towards increasing x − , thus yielding a field distribution which is asymmetric in x − . This was indeed the prevailing viewpoint in the original literature on the JIMWLK evolution (see e.g. [36,40]), but now we shall argue that this is actually not quite right: although the functional derivatives in Eq. (2.2) have a one-sided action which amounts to color rotations at the largest value of x − alone, the overall structure of the Hamiltonian is such that the target field is nevertheless built symmetrically in x − . In fact, it is precisely this symmetry of the target field distribution under reflection in x − which is responsible for the mirror symmetry in the evolution equations.
To see that, it is useful to notice that Eq. (2.2) can be alternatively rewritten as [44,45] where J a Lu and J a Rv are functional differential operators acting as 'left' and 'right' Lie derivatives -that is, the generators of infinitesimal color rotations at the largest and, respectively, smallest value of x − . They are defined as and satisfy where the second equation follows from the first one after using Eq. (2.14). These equations imply the following commutation relations showing that the two sets of generators satisfy two independent SU(N c ) Lie algebras. The physical interpretation of the various terms in Eq. (2.24) is quite transparent: the action of J a Lu on the quark Wilson line V † x (the first equation in Eq. (2.26)) amounts to the addition of an infinitesimal layer of color field at the largest values of x − , whereas the action of J a Ru is tantamount to a corresponding addition at the smallest values of x − . Hence, the manifest symmetry of Eq. (2.24) under the exchange L ↔ R implies that, during the highenergy evolution, the distribution of the target color field α a (x − , x) -by which we mean its support and correlations -is built symmetrically in x − around x − = 0. Moreover, this is also the origin of the mirror symmetry since, as previously noticed, the latter follows from the combined action of the first two terms in the Hamiltonian (2.2), or (2.24) -those which get interchanged with each other under the permutation L ↔ R of the Lie derivatives.
To better appreciate the differences between an evolution which is symmetric in x − and one which is not, it is instructive to consider the evolution of a Wilson line, say for a quark projectile. When computing the target average (2.5) with the CGC weight function at rapidity Y , the support of the color field α a ( . (This follows from the uncertainty principle: the softest gluon modes that have been integrated over have longitudinal momentum p + = xP with x = e −Y , with P = the total hadron momentum; hence, they are delocalized in x − over a distance ∼ 1/p + ∝ e Y .) Thus, the quark Wilson line can be equivalently rewritten as This makes it manifest that, with increasing Y , the Wilson line 'grows' simultaneously at its both endpoints. To visualize the effect of one step in the evolution (Y → Y + dY ), it is useful to discretize rapidity by writing Y n = nǫ, with ǫ an infinitesimal rapidity interval. Then under one additional step n → n + 1, the upper bound of the support extends as represent the additional fields generated in this evolution step per unit of space-time rapidity. The infinitesimal gauge rotations associated with these new fields can be expanded in powers of ǫ. Strictly speaking, this expansion must be pushed to quadratic order in ǫ, to match with the fact that the evolution Hamiltonian (2.24) involves second order functional derivatives. However, the quadratic terms arising from the expansion of a given Wilson line do not contribute to the evolution of gauge-invariant observables: they would yield 'tadpole' contributions ∝ δ uv , but the dipole kernel in the Hamiltonian vanishes when u = v. In other terms, the two functional derivatives within H must act on different Wilson lines within the observable to give a non-zero result. So, we can restrict the expansion of (2.29) to linear order in ǫ, which yields In Ref. [28], the JIMWLK evolution has been reformulated as a random walk in the space of Wilson lines, which is formally such that one additional step corresponds to an infinitesimal rotation of V † (x) on the 'left' alone. However, by inspection of the manipulations there, one can check that the additional contribution α n+1 (x) to the target field in the (n + 1)th step is such that, in reality, that step simultaneously generate a color precession on the 'left' and on the 'right'. That is, the Langevin process introduced in Ref. [28] does in fact describe a symmetric evolution for the Wilson lines (or for the target field distribution), although this has not been recognized there.

The Gaussian approximation
In this section we shall demonstrate that the JIMWLK equation for the CGC weight function admits an approximate Gaussian solution which properly captures both the BFKL dynamics in the dilute regime at k ⊥ ≫ Q s (Y ) and the approach towards the black disk limit in the saturation regime at k ⊥ ≪ Q s (Y ). Our analysis improves over previous, related, constructions in the literature [36][37][38] at two important levels: (i) we actually justify the Gaussian approximation -including for the description of the higher-point correlation functions and for finite N c -, on the basis of the Balitsky-JIMWLK equations; (ii) we implement the 'mirror' symmetry discussed in Sect. 2.3, that is, we construct a Gaussian distribution which is symmetric in x − at any Y . As we shall see, this last condition is in fact compulsory to achieve a faithful description of the JIMWLK dynamics deeply at saturation.
The material of this section is organized as follows: the Gaussian weight function is introduced in Sect. 3.1, compared to the MV model in Sect. 3.2, and justified in Sects. 3.3 and 3.4 by comparison with piecewise approximations to the JIMWLK equations in the limiting regimes alluded to above.

The Gaussian weight function
The most general Gaussian weight function which is consistent with gauge symmetry 4 and describes a target field distribution which is symmetric in is an even function of x − , assumed to be invertible. The functional δ-function δ Y [α] ensures that the target field vanishes at larger Here, δ(α a (x − , x)) is the usual δ-function and a discretization of the space-time is understood. Finally, the overall normalization factor N Y in Eq.
Eq. (3.1) implies that the only non-trivial correlation of the target fields is their 2-point function, which is moreover local in x − : Within this Gaussian approximation, the locality in x − is required by gauge symmetry: to preserve the latter, any non-locality in x − should be accompanied by gauge links (Wilson lines) built with the field α a , which would spoil Gaussianity. The Gaussian distribution (3.1) is manifestly symmetric in x − around x − = 0 and this symmetry is preserved by the high energy evolution. In fact, Eq. (3.1) depends upon Y only via the two endpoints, , of the support in x − , meaning that the high-energy evolution proceeds via the symmetric expansion of the color field distribution towards both larger and smaller values of x − . Specifically, by using the methods in Refs. [28,36], one can check that the Gaussian weight function (3.1) obeys the following evolution equation: where the 'left' ('right') functional derivatives act on the target field at the largest (smallest) value of x − , that is, at x 2 ) denotes the field correlator per unit space-time rapidity as produced in the last step of the evolution, Eq. To justify the Gaussian Ansatz (3.1) for the CGC weight function, we shall shortly compare the associated evolution equation (3.4) to the actual JIMWLK equation, in different kinematical regimes. In this process, we shall deduce piecewise approximations for the kernelγ Y (x 1 , x 2 ), valid at high (k ⊥ ≫ Q s (Y )) and low (k ⊥ ≪ Q s (Y )) momenta, respectively.

The McLerran-Venugopalan model
Before we turn to the JIMWLK evolution, let us briefly discuss the McLerran-Venugopalan (MV) model [29,30] that we shall take as our initial condition at rapidity Y 0 . Besides providing the initial conditions, this model (and its ad-hoc extrapolation towards high-energy) will serve as a baseline of comparison for the mean-field results that we shall later obtain. Its discussion will also give us an opportunity to clarify some subtle aspects of the Gaussian approximation, like the gauge artifacts in the 'α-representation', cf. Eq. (3.1).
In the MV model one assumes that the color charges in the nucleus are uncorrelated valence quarks. Accordingly the distribution of the color charge density ρ a (x − , x) is a Gaussian with a kernel which is local in the transverse plane: The quantity λ(x − , x) has the meaning of color charge squared per unit transverse area per unit longitudinal distance. In general the nucleus is assumed to be homogeneous in the transverse plane, i.e. the kernel in Eq. (3.7) is taken to be independent of x. Under that assumption, the calculation of expectation values in the MV model is not sensitive to the detailed dependence of the kernel upon x − , but only to its integral which physically represents the color charge squared per unit area. In the above equation, the quantity λ y ≡ |x − |λ(x − ) (the strength of the charge correlator per unit space-time rapidity) has been defined by analogy with Eq. (3.6). The last equality follows after counting the color charges of the valence quarks within a nucleus with atomic number A and transverse area πR 2 A (see e.g. [40]). The fact that it is only the integrated quantity (3.9) which matters arises from the fact that, under the present assumptions, the charge correlator (3.8) is separable as a function of x − and the transverse coordinates. We shall return to this issue in Sects. 4.2 and 4.3.
Eq. (3.7) is gauge invariant, but in order to make contact with the α-representation that we use throughout this paper, we shall henceforward consider it within the class of gauges where the target field is of the form So, for a homogeneous target, Eq. (3.8) implies the following expression for the 2-point function for the color field, in transverse momentum space (we denote r = x 1 − x 2 and k ⊥ = |k|): Here and from now on, we prefer to work with the expressions of the various correlators per unit space-time rapidity, cf. Eq. (3.6), since these are the expressions which most directly enter the mean-field evolution equations like (3.4). Eq. (3.10) raises a potential problem: the Fourier transform of this expression back to the transverse coordinate space is not well defined, as it involves a (quadratic) infrared divergence. This problem reflects the fact that, by itself, the field α a (x − , x) is not invariant under the residual gauge transformations which preserve the structure A µ a = δ µ+ α a for the target field. The infinitesimal version of such a transformation reads with ω a (x − ) an arbitrary function [42]; so, clearly, the color charge density ρ a (x − , x) is invariant under this transformation. Strictly speaking, the general weight function W Y (and, in particular, its Gaussian approximation, Eq. (3.1)) should be written as a functional of ρ, to make gauge symmetry manifest. On the other hand, observables like scattering amplitudes are built with Wilson lines, which are path-ordered exponentials of α. Taken separately, one Wilson line is not gauge invariant (rather, it transforms via color rotations [42]), but the physically relevant operators, which involve a product of such lines, cf. Eq. (2.9), are invariant. Whenever computing the expectation value of such a gauge-invariant operator, there is no problem with using the weight function in the α-representation, as given in Eq. (3.1): all the gauge artifacts cancel out in the final result.
As an example, consider the calculation of the dipole S-matrix within the MV model. The corresponding result is well known and reads (see also Sect. 4.2) where we have assumed the MV model to apply at all the rapidities Y ≤ Y 0 and we defined involves only the following linear combination of the target field correlators which is gauge-invariant, since under a residual gauge transformation the target field α a (x − , x) changes by a x-independent quantity. The sign in the r.h.s. of Eq. (3.13) is such that γ y (x 1 , x 2 ) be positive-semidefinite. The last equality in Eq. (3.13), which involves the color charge correlator λ y , illustrates the fact that the infrared divergences due to gauge artifacts cancel out in the linear combination (3.13). Strictly speaking, the above integral over k still has a logarithmic infrared divergence, but this is milder than the quadratic divergence appearing in the Fourier transform ofγ y (k) in Eq. (3.10). The remaining divergence is not a gauge artifact anymore, but a 'physical' singularity of this model: it reflects the lack of correlations among the color sources. After taking into account the high-energy evolution, transverse correlations get built which screen out this divergence, as we shall shortly see. For completeness, let us estimate the final integral in Eq. (3.12): introducing an infrared cutoff Λ to regularize the remaining infrared divergence and writing r = |x 1 − x 2 |, one finds where Q 2 0 ≡ 2α 2 s C F A/R 2 A is essentially the nuclear saturation scale 5 (as probed by a quarkantiquark dipole) at the initial rapidity Y 0 . Although obtained within the MV model, the above results are generic in the following sense: all the gauge-invariant observables computed in the Gaussian approximation involve the kernelγ y (x 1 , x 1 ) of the Gaussian (the correlator of the target color field) only via the linear combination shown in Eq. (3.13). So, in practice, there is no problem with using the α-representation, as shown in Eq.
In this equation, ρ a L and ρ a R refer to the color charge densities at x − = X − and x − = −X − , respectively. When applied to the evolution of the Wilson-line correlations, Eq. (3.15) amounts to constructing the Wilson lines via symmetric iterations, i.e. via infinitesimal color precessions which proceed simultaneously 'on the left' and 'on the right', as shown in Eq. (2.29). However, within the context of the MV model, this symmetric iteration is merely a choice of a discretization prescription and any other choice is equally good. As a matter of fact, the common choice in the literature in this context (see e.g. [3,8,31,32]) is to perform asymmetric iterations 'on the left' : where this time n refers to a discretization of the x − axis. This procedure is tantamount to solving the following evolution equation In practice, one often takes x − 0 → ∞, since the results are anyway insensitive to the actual value of x − 0 , but only depend upon the integral dX − λ(X − ). The above discussion sheds more light on the role of the 'left-right' symmetry in the evolution equations. So long as the CGC weight function is given (like in the MV model) and the associated evolution equations are merely used as a convenient device to compute expectation values, the symmetric discretization in Eq. (2.29) is not compulsory and it might not even be the most convenient one in practice. However, for the JIMWLK equation and any (mean field) approximation to it, the symmetric iteration is the only one to be correct, since this is how the target field distribution gets actually built via quantum evolution: the 'outer' layers (those located at larger values of |x − |) are constructed after the 'inner' ones (those at smaller |x − |), and the new correlations built in one step depend upon the color field produced in all the previous steps. Hence, it would make no sense to consider an asymmetric evolution, like Eq. (3.17), since this would violate causality in the domain of negative x − .

Weak-scattering regime: the BFKL dynamics
In what follows, we shall study the JIMWLK evolution in two limiting regimes -large transverse momenta k ⊥ ≫ Q s (Y ) in this section and relatively small momenta k ⊥ ≪ Q s (Y ) in the next subsection -with the purpose of showing that, in both regimes, the evolution is consistent with a mean field approximation of the type shown in Eq. (3.4). We recall that Q s (Y ) is the saturation momentum in the target (in a frame in which the target carries most of the total rapidity separation Y ) and it increases with Y very fast. For a multi-point correlation function like the quadrupole (2.8), the statement that the 'transverse momenta are much larger than Q s ' means that all the transverse separations r ij ≡ |x i − x j | between the external points are much smaller than 1/Q s (Y ). Similarly, by 'momenta much smaller than Q s ', we mean that r ij ≫ 1/Q s (Y ) for any pair (x i , x j ) of external points. Very asymmetric configurations, where some of the distances r ij are much larger than 1/Q s (Y ) while the others are much smaller, are strictly speaking not covered by the present analysis and must be separately studied. We shall discuss some examples of that kind in Sect. 4.4 below.
For high transverse momenta k ⊥ ≫ Q s (Y ), the gluon density in the target is low, meaning that the corresponding color field is weak: g dx − α ≪ 1. It is then possible to expand the Wilson lines to lowest non-trivial order in the field in their exponent, within both the JIMWLK Hamiltonian and the operators defining the observables. For an operator like the dipole S-matrix Eq. (2.7), we need to push the expansion in gα up to the second order, since the linear terms vanish after averaging. Introducing the dipole T -matrix operatorT x 1 x 2 ≡ 1 −Ŝ x 1 x 2 , whose expectation value represents the corresponding scattering amplitude, this expansion yields The weak scattering regime corresponds to T Y ≪ 1. Note that Eq.
where it is understood that T Y is evaluated according to Eq. (3.18). More generally, in this dilute regime, all the n-point functions of the type shown in Eqs. (2.9) or (2.10) reduce to linear combinations of dipole amplitudes. This already shows that a Gaussian approximation for the CGC weight function should be indeed possible, to the accuracy of interest. To identify this approximation, let us also consider the weak-field limit of the JIMWLK Hamiltonian. Its obtention is facilitated by observing that Eq. (2.2) can be rewritten as with α a u as defined in Eq. (3.18). After also using (T a ) bc = if abc , one finds H ≃ H BFKL with (

3.22)
This Hamiltonian is supposed to act on operators which are themselves evaluated in the weakscattering regime and hence are quadratic functions of the field α a x , as illustrated in (3.18) and (3.19). Clearly, the only evolution equation of interest for us here is that obeyed by the dipole scattering amplitude (3.18). This is readily obtained as and is recognized as the BFKL equation [46][47][48], that is, the equation obtained after linearizing Eq. (2.16) with respect to T Y . By using its solution, one can compute any other n-point function of the Wilson lines, like Eq. (3.19), in this dilute regime. We now construct the Gaussian approximation which reproduces the BFKL equation. To that aim, we shall compare the mean-field equation for T x 1 x 2 Y generated by Eq. (3.4) with Eq. (3.23) and thus deduce an approximate expression forγ Y (u, v) valid in this linear regime. Notice that the left and right functional derivatives yield identical results when acting on the field α a x which is integrated over x − . Hence, Eqs. (3.4) and (3.18) imply with γ Y (x 1 , x 2 ) defined as in Eq. (3.13). The last equation can be integrated to yield It is easy to check that the same expression for T Y would be obtained by directly evaluating the expectation value in Eq. Returning to the mean-field expression (3.25) for T Y , this must be consistent with the BFKL equation (3.23). This is clearly the case provided the function f Y (x 1 , x 2 ) itself satisfies the BFKL equation: The initial conditions for the above equations can be taken from the MV model, which yields The solution to Eq. (3.26) is by now well understood. Here we will just remind that the BFKL evolution introduces transverse correlations between the 'color sources' (radiated gluons) which ensure that the solution f Y (x 1 , x 2 ) becomes infrared finite after a rapidity evolution Y − Y 0 ∼ 1/ᾱ. In particular, in the window for 'extended geometric scaling' [49][50][51], which holds for transverse momenta relatively close to (but still larger than) the saturation momentum Q s (Y ), one has 6 with γ s ≈ 0.63 (the 'BFKL anomalous dimension 6 Notice that k 4 ⊥ γY (k) = k 4 ⊥γY (k) since in momentum space the difference between γY (k) andγY (k) is proportional to δ (2) (k). at saturation'). Then, clearly, the integral over k in Eq. (3.13) is well defined when computed within the BFKL approximation.
To summarize, the mean-field equation (3.4), where it is understood that the kernel can be replaced asγ Y → −γ Y with the function γ Y (u, v) determined by Eq. (3.26), properly encodes the BFKL evolution of the dipole amplitude in the weak scattering regime. This conclusion holds for any value of the number of colors N c and it extends to all the n-point functions like (2.9) and (2.10) which, in this regime, reduce to linear combinations of dipole amplitudes.
Before concluding this section, let us recall that there are also other aspects of the BFKL dynamics, which cannot be encoded into a Gaussian weight function. They refer to operators more complicated than those in Eq. (2.9), which already at weak scattering involve more than two gluon exchanges; that is, to lowest order in the weak field expansion, they involve polynomials in α of a degree higher than two. (Such operators can be obtained e.g. by subtracting the dipolar contributions to the Wilson-line operators in Eqs. (2.9)-(2.10).) An example of that type is the 'odderon' operator, which describes C-odd exchanges and which in perturbation theory starts with three gluon exchanges. The corresponding evolution equation is correctly encoded (to leading logarithmic accuracy) in the JIMWLK equation [42] -in particular, its low-density limit, known as the 'BKP equation' [52][53][54], is generated by the weak-field limit (3.22) of the JIMWLK Hamiltonian [42] -but this description goes beyond the purpose of a Gaussian approximation, which by construction can encode only the 2-point function of the α field. For instance, to describe odderon effects in the initial conditions, one needs an extension of the MV model allowing for a non-trivial 3-point function [55].
What is however remarkable about the Gaussian approximation that we pursue here is its capacity to encode non-trivial correlations among n Wilson lines with arbitrary n in the strong scattering regime, where the linear relation between the n-point functions and the 2-point function does not hold anymore. This will be discussed in the next subsection.

Strong-scattering regime: the dominance of the 'virtual' terms
For relatively low transverse momenta k ⊥ Q s (Y ), the gluon occupation numbers in the target wavefunction saturate at a large value of order 1/α s , meaning that g dx − α ∼ O(1). This in turn implies that the scattering is strong for projectiles with transverse sizes r 1/Q s . For instance, the dipole scattering amplitude T x 1 x 2 Y becomes of order one when |x 1 − x 2 | 1/Q s . Then Eq. (3.19) implies that, for generic configurations at least, the quadrupole scattering becomes strong when at least one (which necessarily means at least three) of the six transverse distances r ij = |x i − x j | is of order 1/Q s , or larger. Similar considerations apply to the higherpoint correlations. In this regime, the Wilson lines cannot be expanded out anymore. Rather, they resum multiple scattering to all orders in the eikonal approximation.
To correctly describe the high-energy evolution in the presence of gluon saturation and multiple scattering, it is of course essential to keep the non-linear terms in the Balitsky-JIMWLK equations, so like ŜŜ Y in the equation (2.16) for the dipole S-matrix and ŜQ Y in the r.h.s. of Eq. (2.17) for the quadrupole. In fact, these are precisely the terms responsible for the approach towards saturation in the gluon distribution and towards unitarity in the scattering of the projectile. Accordingly, in the transition regime towards saturation/unitarity (i.e. for k ⊥ ∼ Q s (Y )), one has to deal with the whole, infinite, hierarchy of coupled evolution equations: no simple mean-field approximation (like a Gaussian) is possible in that regime. However, the situation drastically simplifies deeply at saturation (k ⊥ ≪ Q s (Y )), where the only role of the non-linear terms in the equation is to forbid further evolution -or, more correctly, to limit the transverse phase-space for the high-energy evolution: gluons with soft momenta k ⊥ ≪ Q s (Y ) can (almost) not be emitted anymore, meaning that domains separated by transverse distances r ≫ 1/Q s (Y ) evolve independently from each other. This leads to considerable simplifications in the Balitsky-JIMWLK equations, which can be most directly recognized by inspection of the projectile evolution.
For multi-partonic projectiles which are such that all the interparticle separations r ij are much larger than 1/Q s (Y ), the associated S-matrices are very small (close to zero) -the more so the larger the number of partons. Roughly speaking, and up to subtleties related to the 1/N 2 c corrections to which we shall later return, a 2-dipole projectile scatters more strongly than a single-dipole one, ŜŜ Y ≪ Ŝ Y , a projectile made with a dipole plus a quadrupole scatters more strongly than the quadrupole alone, ŜQ Y ≪ Q Y , etc. When this happens, the 'virtual' terms dominate the evolution, whereas the 'real' terms can be simply replaced with a lower cutoff ∼ 1/Q s (Y ) on the transverse separation |z − x i | between the newly emitted gluon at z and any of the preexisting partons at x i . Once this is done, the resulting evolution equations are linear and hence admit a Gaussian solution. This is of course related to our previous observation in Sect. 2.3 that the only effect of the 'non-linear terms' (Wilson lines) within H virt is to transform 'left' color precessions into 'right' ones and thus ensure the symmetric expansion of the target field distribution in x − . This also shows that, in this high density regime, where the Wilson lines cannot be expanded anymore and 'left' and 'right' functional derivatives have different mathematical consequences, it is essential to keep trace of the 'mirror' symmetry of the evolution, by using a symmetric Gaussian, as shown in (3.1).
To render these considerations more precise and construct the corresponding Gaussian approximation, we shall develop our mathematical arguments in two steps: (i) at large N c , and (ii) at finite N c .
(i) Large N c : Within the context of the large-N c approximation, the prominence of the 'virtual' terms in the approach towards the black disk limit is quite obvious and has been pointed out at several places in the literature [12,39,56,57]. Specifically, the 'real' terms which survive at large N c involve double-trace operators, which can be factorized to the accuracy of interest: Then we can write e.g.
Now, in equations like (2.16) or (2.17), the transverse position z of the emitted gluon is integrated over, so it can become close to one of the external points x i , in which case Eq. (3.27) does not hold anymore. However, in the high density regime under consideration, such special configurations are disfavoured by the phase-space for the transverse integration. Namely, assuming |x i − x j | ≫ 1/Q s (Y ) for all the pairs (i, j), one can check that the integrals over z receive their dominant contributions from points relatively far apart from all the external points, which satisfy Indeed, the contribution of such a range is enhanced by the large transverse logarithm Hence, to leading logarithmic accuracy in the sense of Eq. (3.29), one can indeed neglect the 'real' terms in the Balitsky-JIMWLK equations at large N c , as anticipated.
(ii) Finite N c : The physical argument at finite N c is the same as at large N c except that, now, one has to take into account the fact that the evolution described by the 'real' terms truly corresponds to the emission of a gluon, and not just to the splitting of, say, one dipole into two dipoles. When this new gluon is sufficiently soft, in the sense that |z − x i | 1/Q s (Y ) for any i, its emission leads to a partonic system with a wider distribution of color charge in the transverse plane, which therefore interacts stronger with the target than the original projectile. But in order to rigorously justify this, one needs to actually estimate the S-matrix for, say, a quark-antiquark-gluon (qqg) system deeply at saturation and show that this is indeed much smaller than the S-matrix of the dipole (qq). To appreciate how subtle this is, let us recall that, when rewriting the 'real' terms in terms of Wilson lines in the fundamental representation (as customary in the Balitsky-JIMWLK equations), one generates single-trace pieces proportional to 1/N 2 c , which by themselves count on the same footing as the 'virtual' terms near the unitarity limit. For instance, the contribution of the 'real' terms to the dipole equation involves the following expectation value (cf. the second line in Eq. (2.13)) where one may naively think that the second, single-trace, term dominates over the first one when all the transverse separations are much larger than 1/Q s (Y ). As another example, we show here some 'real' terms from the evolution equation (2.17) for the quadrupole, namely those arising when acting with H real on the two quarks at x 1 and x 3 : where the second line follows from the first one after using the Fierz identity (2.15). Once again, one may think that the last term in Eq. (3.31), proportional to (1/N 2 c ) Q Y , is the dominant term for large transverse separations ≫ 1/Q s (Y ) (and for finite N c ). If that was indeed the case, there would be a mixing between 'real' and 'virtual' terms deeply at saturation, which would prevent a Gaussian approximation (since the latter could not accommodate the 'real' terms beyond the BFKL approximation).
The situation becomes even more confusing if one recalls that, in the equations obeyed by the single-trace observables, the terms subleading at large N c precisely cancel between 'real' and 'virtual' contributions. In view of this, one may be tempted to argue that the finite-N c corrections are totally irrelevant. But that would be wrong, since there is no similar cancelation in the equations obeyed by the multi-trace operators, like ŜŜ Y or ŜQ Y .
What 'saves' the Gaussian approximation, is the fact that, in spite of appearance, the single-trace components in equations like (3.30) or (3.31) do not dominate over the respective double-trace ones, but merely subtract fake 'single-trace contributions' from the latter, that have been artificially introduced via the Fierz identity. That is, the expression in the first line of Eq. (2.13), which involves an adjoint Wilson line and describes a qqg system, vanishes very fast in the approach towards the black disk limit, where it is suppressed with respect to the corresponding 'virtual' term Ŝ x 1 x 2 Y . But this is not the case for the 2-dipole S-matrix in the second line of Eq. (2.13), which in that regime approaches to (1/N 2 c ) Ŝ x 1 x 2 Y . A similar discussion refers to Eq. (3.31): deeply at saturation, the observable in the first line, which describes a qqqqg partonic system, is suppressed compared to the respective 'virtual' terms, that is, the quadrupole and the pair of dipoles.
In order to demonstrate this while dealing with an infinite hierarchy, we shall provide a self-consistent argument. That is, we start by assuming that the JIMWLK evolution deeply at saturation is controlled by H virt alone and we prove that, under this assumption, the 'real' terms in Eq. (3.30) and Eq. (3.31) vanish exponentially faster than the respective 'virtual' terms in the vicinity of the unitarity limit. We shall give the details of the proof for the dipole evolution, i.e. for the operator in Eq. (3.30), and then briefly discuss its generalization to the quadrupole and higher n-point functions. In this context, by 'H virt ' we mean, of course, the first two terms in the JIMWLK Hamiltonian (2.2) together with the phase-space restriction |z − x i | ≫ 1/Q s (Y ) as introduced by the 'real' terms. That is, we work in the leading-logarithmic approximation in Eqs. (3.28)-(3.29), which enables us to write (3.32) to the accuracy of interest. So, let us calculate the action of H virt on the combination of the operators appearing in Eq. (3.30). This action on the second term has been already computed in Eq. (2.12), that we here rewrite for convenience as with the integral over w understood in the sense of Eq. (3.29). Now, when both derivatives act on the same (either the first or the second) dipole of the first term in Eq. (3.30), we get the following, 'diagonal', contribution 34) and when they act on different dipoles we find the cross term Putting everything together we arrive at where it is crucial to notice that the operator of interest has been reconstructed in the r.h.s. of the equation. It should be clear from the above derivation that this would have not happened without the subtraction of the 1/N 2 c -suppressed dipole. By assumption, the above equation describes the approach towards unitarity of the 'real' piece in the evolution equation (2.16) for the dipole S-matrix. This should be compared to Eq. (2.12), which describes the corresponding approach for the 'virtual' piece (the dipole itself). Clearly, the kernel in Eq. (3.36) is 'twice as large' than that in Eq. (2.12), showing that, deeply at saturation, the expectation value of the 'real' operator in Eq. (3.30) vanish exponentially faster than the 'virtual' term ∝ Ŝ x 1 x 2 Y . Hence, the latter dominates in the evolution equation and in this regime, as anticipated.
This self-consistent argument can be generalized to higher-point correlators, as we now show for the operator which appears in Eq. (3.31) and counts for the evolution of the quadrupole. Acting with H virt , we see that the only new element appearing, when comparing to Eq. (3.36), is operator mixing. Indeed, one finds that we also need to consider the operatorŝ plus permutations of all the operators appearing in Eqs. (3.37), (3.38) and (3.39). Without going into too much detail, one understands that the action of H virt on the above operators leads to where the elements of the 3 × 3 matrix M are proportional to −ᾱ/2π times an integral over w of linear combinations of the dipole kernel. The counting is such that the integrand in the diagonal elements is the sum of three dipole kernels which enter all with a plus sign, plus terms proportional to 1/N 2 c (in analogy with Eq. (3.36)). Furthermore, the integrand in the nondiagonal elements is the sum of dipole kernels with equal number of plus and minus signs, plus again terms proportional to 1/N 2 c . Clearly, the diagonal components are those which control the approach towards the black disk limit and they are larger than those which control the corresponding evolution for the 'virtual' terms in Eq. (2.17), that is Q Y and ŜŜ Y .
Incidentally, the above argument also shows that the two operators in Eqs. (3.38) and (3.39) vanish faster than the quadrupole and the 2-dipole system in the approach towards unitarity. This is interesting since these are precisely the 'real' terms in the evolution equation for Ŝ x 1 x 2Ŝ x 3 x 4 Y , whereas Q Y and ŜŜ Y are the corresponding 'virtual' terms. So, we have also demonstrated the property of interest (the dominance of the 'virtual' terms deeply at saturation) for the evolution of a system of two dipoles with arbitrary coordinates. We are confident that a similar proof applies to the higher-point (single-trace or multi-trace) correlation functions.
It is furthermore instructive to check these arguments via explicit calculations within the Gaussian approximation (3.1). Via methods to be described later, this yields e.g. [37] where we have assumed that Ŝ x i x j and Ŝ x 1 x 3Ŝ x 3 x 2 are equal to 1 as an initial condition, to simplify writing. This formula makes it clear that the operator in the l.h.s. vanishes, roughly, as a 'dipole squared' in the approach towards the unitarity limit. A corresponding argument for the operator (3.37) which enters the evolution of the quadrupole will be given in Sect. 4.4. We thus conclude that the JIMWLK evolution deeply at saturation is indeed correctly described by the 'virtual' Hamiltonian in Eq. (3.42) This applies for k ⊥ ≪ Q s (Y ) and is recognized as the 2-dimensional Coulomb propagator. In turn this implies that the charge-charge correlator λ Y (k) = k 4 ⊥ γ Y (k) vanishes like k 2 ⊥ when k ⊥ → 0, which is the expression of color shielding due to gluon saturation [38,58]: the average color charge squared vanishes when integrated over a transverse area ≫ 1/Q 2 s (Y ). Notice that in some previous versions of the mean field approximation [36,38,39], one has assumed that the JIMWLK Hamiltonian takes an even simpler form in the vicinity of the black disk limit, namely it reduces to the first term in Eq. (3.32), which involves the 'left' derivatives alone. That simplification was motivated [39] by a 'random phase approximation', which assumed that, in the strong field regime deeply at saturation, all the Wilson lines within the Hamiltonian are rapidly oscillating and thus average out to zero. As shown by our present manipulations, this argument is qualitatively correct, but only for the 'real' terms (the last 2 terms) in the JIMWLK Hamiltonian.
To summarize the arguments in this section, the JIMWLK evolution in the two limiting regimes -the weak-scattering regime at low gluon density and the approach towards the black-disk limit deeply at saturation -can be properly encoded, for any value of N c , into a symmetric Gaussian weight function of the type (3.1). In turn, this Gaussian is tantamount to the functional evolution equation shown in Eq. (3.4), or to the following, mean field, Hamiltonian: This has the same operator structure as the 'virtual' piece of the JIMWLK Hamiltonian, but with a different, Y -dependent, kernel, which is essentially the 2-point function of the target color field. This kernel interpolates between the solution to the BFKL equation ( In practice, however, a proper choice for the kernel is probably essential in order to achieve a good global accuracy. In the following section, we shall propose two such choices.

Evolution equations in the Gaussian approximation
In this section we shall describe two methods for constructing smooth, global, expressions for the kernel γ Y (u, v), which are equivalent with each other to the accuracy of interest. Then we shall derive the evolution equations associated with the mean field Hamiltonian (3.43), first for generic N c (in Sect. 4.2), then at large N c (in Sect. 4.3). As before, we shall mostly focus on the evolution of the dipole and of the quadrupole. In the large-N c limit we shall recover the equations previously proposed in Ref. [12]. At finite N c , the general equations are more complicated, but explicit solutions will be presented for special configurations in Sect. 4.4.

Self-consistent constructions of the kernel
Within the Gaussian approximation, it is always possible to trade the kernel γ Y (x 1 , x 2 ) for the dipole S-matrix Ŝ x 1 x 2 Y , for which it is easier to construct global, smooth, approximations in practice. The expression of Ŝ x 1 x 2 Y in the Gaussian approximation is well known in the literature and will be rederived, for completeness, in the next subsection. Here it is preferable to work with the corresponding evolution equation, which reads for a dipole in an arbitrary representation R of the color group. Hence, if one disposes of a numerical solution to the JIMWLK equation, like in Refs. [11,59,60], then one can use the respective estimate for the dipole S-matrix, say, in the fundamental representation, together with Eq. (4.1) to deduce a corresponding estimate for the kernel. In practice, solving the full JIMWLK evolution is quite tedious, so it is customary to rely on its large-N c approximation (for the dipole evolution), namely the BK equation. This is equally good for the present purposes (including at finite N c ) since, as noticed at the end of the previous section, the kernel γ Y (x 1 , x 2 ) is independent of N c . Hence, its limiting behaviors are correctly reproduced by the large-N c version of Eq. (4.1). The latter implies with Ŝ BK x 1 x 2 Y denoting the solution to the BK equation with an initial condition itself evaluated at large N c . (This is important in order to preserve finite-N c accuracy in the limiting kinematical domains where Eq. (4.2) is correct as it stands for any value of N c .) For instance, within the MV model, this is provided by Eqs. (3.11)-(3.12) with C F ≃ N c /2. The function Ŝ BK can be obtained either as an exact, numerical, solution to the BK equation, or as an analytic approximation to it with the correct limiting behaviors.
Although correct to the accuracy of interest, the construction in Eq. (4.2) may look a bit unesthetic at a conceptual level, as it requires an input -the solution of the BK equationwhich seems external to the Gaussian approximation. However, as we now explain, Eq. (4.2) can be also understood as a self-consistency condition internal to the mean field approximation [37]. Specifically, let us start with the first equation in the Balitsky-JIMWLK hierarchy, that is Eq. (2.16) for the dipole, and evaluate all the expectation values there with the Gaussian weight function (3.1), that is, by using Eq. (4.1) with C R = C F together with Eq. (3.41). This leads to an equation for the kernel γ Y (x 1 , x 2 ) which is precisely equivalent to solving the BK equation and then computing the kernel according to Eq. (4.2) (see Ref. [37] for details).
This procedure, which in Ref. [37] has been dubbed 'the Gaussian truncation', is of course not unique: one can similarly start with any equation in the Balitsky-JIMWLK hierarchy, compute all the expectation values there with the Gaussian weight function (3.1), and thereby transform the original equation into an equation for γ Y (x 1 , x 2 ). A different self-consistency condition, which in practice is not more difficult to use than Eq. (4.2), has been originally proposed in Ref. [36]. It amounts to requiring the mean-field Hamiltonian (3.43) to coincide with the JIMWLK Hamiltonian (2.2) on the average : The average here refers, of course, to the Gaussian weight function, which implies that both the l.h.s. and the r.h.s. in the above equation are proportional to δ ab . Introducing the S-matrix for the gluonic dipole operator which it related to the respective fermionic operator as shown in the second equality above, and multiplying Eq. (4.3) by δ ab , we can rewrite the latter as This relation immediately implies γ Y (u, u) = 0 because of the corresponding property of the dipole kernel M uvz . It furthermore implies that γ Y (u, v) is symmetric under u ↔ v, because so is the gluonic dipole S-matrix, as obvious from its definition (4.4).
The self-consistency condition (4.5) is most conveniently written as an equation for Ŝ A Y : by using Eq. (4.1) in the adjoint representation (C R = C A ≡ N c ), one finds The initial condition should be taken from the MV model applied to a gluonic dipole, that is, To summarize, by using either Eq. (4.2) together with the solution to the BK equation, or Eq. (4.1) with R = A together with the solution to Eq. (4.6), one obtains two global approximations for the kernel γ Y (x 1 , x 2 ), which can differ from each other only in the transition region around saturation. In principle, these two approximations are equivalent with each other to the accuracy of interest. In practice, any (numerical) difference between them will have an impact on the calculation of the higher n-point functions, to be described in the remaining part of this section. It is likely that one can optimize the reliability of this whole scheme by choosing a 'good' approximation for the dipole S-matrix used to compute the kernel, like the exact, numerical, solution to the JIMWLK equation, or to the BK equation at least.

Mean-field equations for the dipole and the quadrupole
The evolution equations associated with the mean-field Hamiltonian (3.43) are straightforward to obtain, by using the same techniques as for the JIMWLK Hamiltonian (2.2). Alternatively, given that H MFA has the same operator structure as H virt , the mean-field equations can be directly inferred from the corresponding Balitsky-JIMWLK equations, by keeping only the 'virtual' terms in the latter and replacing everywhere the dipole kernel according to In doing that, one should be careful to restore all the 'virtual' terms in the Balitsky-JIMWLK equations, including those which may have canceled against similar 'real' contributions and hence are not manifest in the final equations (cf. the discussion in Sect. 2.2). Clearly, the resulting equations will inherit the relatively simple structure characteristic of the 'virtual' terms. They form closed systems of equations, which connect only correlation functions with the same number of Wilson lines (or external points) and are local in the transverse coordinates, meaning that they do not mix different transverse configurations. Note also that, although linear, these new equations respect unitarity by construction: the tame of the BFKL growth by the non-linear physics of gluon saturation is already encoded in the kernel of the Gaussian. The fact that the S-matrices for the various projectiles approach the right limit at strong scattering is ensured by the unitarity of the Wilson lines which appear within the respective operators. Consider first the dipole equation. Using the substitution rule (4.7) and the respective 'virtual' term in Eq. (2.12), one finds This is easily solved to give where it is understood that Γ Y 0 is given by the MV model, cf. Eqs. (3.11)- (3.12). The corresponding expressions for a color dipole Ŝ R x 1 x 2 Y in an arbitrary representation R are obtained by replacing C F → C R in the equations above (cf. Eq. (4.1)). In particular, in the 'BKrepresentation' in which the Gaussian kernel is computed according to Eq. (4.2), the dipole S-matrix in the Gaussian approximation and in a generic representation R of the color group is related to the solution Ŝ BK From the discussion in Sect. 3, one expects Eq. (4.9) to have the correct limits at both weak and strong scattering, and it is instructive to explicitly check that. At weak scattering, one has where Y s (r) is the rapidity at which saturation is reached over a transverse size r (that is, Q s (Y s (r)) = 1/r), ωᾱ is the logarithmic derivative of the saturation momentum, and we used ln r 2 Q 2 s (y) = ωᾱ(y − Y s (r)) for y ≥ Y s (r). Eq. (4.11) holds in the leading logarithmic approximation in the sense of Eq. (3.29). The large-N c version of this result has already appeared in the literature [56,57].
We now turn to the quadrupole. The corresponding evolution equation in the MFA is obtained as (4.12) Once again the BFKL limit is easy to check: when all the separations |x i − x j | are much smaller than 1/Q s (Y ), one can replace Q Y ≈ 1 and ŜŜ Y ≈ 1 in the r.h.s. of the above equation, which then reduces to The second line, which follows from the first one after using Eq. (3.24), is the expected result for Q Y at weak scattering, cf. (3.19).
Let us also notice that, within the Gaussian approximation, the dipole and quadrupole Smatrices are invariant under charge conjugation, that is, under the exchange of the quarks with the antiquarks. More precisely, from Eqs. (4.8) and (4.12), and using the fact that the kernel γ Y (x i , x j ) is symmetric, we easily deduce that so long as the above conditions already hold at Y 0 (as is the case within the MV model).
Conversely, C-odd scattering amplitudes, like the odderon, cannot be accounted by the Gaussian approximation, as they would require non-Gaussian effects already at Y 0 (cf. the discussion at the end of Sect. 3.3).
Returning to the full equation (4.12), we notice that this is consistent with mirror symmetry, as it should. (E.g. the r.h.s. is symmetric under the exchange x 1 ↔ x 3 .) That would not have been the case 7 , had we considered a Gaussian Hamiltonian with only left derivatives, as proposed in Refs. [36][37][38]. In fact, without any further assumption, this equation implies that the evolution for the antisymmetric part of the quadrupole defined in Eq. (2.20) is closed and homogeneous. In turn, this means that Q asym Y = 0 in the MFA provided this condition was satisfied at Y 0 . It is furthermore clear that the quadrupole couples to the product of two dipoles, whose evolution in the MFA is in turn determined by Since the equations above involve Ŝ x 1 x 4Ŝ x 3 x 2 Y and Q x 1 x 4 x 3 x 2 Y we need also to consider them with the their indices x 2 and x 4 interchanged. (Actually, Q x 1 x 4 x 3 x 2 Y coincides with Q x 1 x 2 x 3 x 4 if one assumes mirror symmetry at Y 0 , but here we prefer to keep the discussion general.) Thus we arrive at a homogeneous system of first order differential equations To see this, let us assume for the sake of the argument the large-Nc limit where the second term is absent.
Then we find that in the first term only γY (x3, x2) and γY (x1, x4) are present and that the fourth term is absent. Clearly the aforementioned symmetry of the evolution equation is lost.
with M Y a 4 × 4 matrix. Its elements are proportional to g 2 γ Y (x i , x j ) accompanied by color factors and can be read from Eqs. (4.12) and (4.15). The difficulty that appears now is that one cannot solve Eq. (4.16) for a generic dependence of γ Y (x i , x j ) on Y , since in general the matrices M Y at different rapidities Y do not commute with each other, that is [M Y 1 , M Y 2 ] = 0.
(More precisely, one could write down a formal solution which involves the rapidity-ordered exponential of the mixing matrix, but we do not find that very useful in practice.) There are special cases where the rapidity integration in the equation above can be explicitly performed, leading to a simpler expression. The large-N c limit to be discussed in the next subsection is one such a special case. Another one is when the kernel γ Y (x i , x j ), is a separable function of Y and the transverse coordinates, plus an arbitrary function of Y : This property is manifestly satisfied within the (homogeneous) MV model, as noticed after Eq. (3.9), and it is also approximately satisfied by the solution to the BK equation, at least in particular kinematical regimes: in the window for extended geometric scaling, where γ Y (r) ∝ (r 2 Q 2 s (Y )) γs with γ s ≈ 0.63 [49][50][51], and also deeply at saturation where γ Y (r) ∝ ln[r 2 Q 2 s (Y )], cf. Eq. (4.11). Presumably, this is a reasonable approximation for all the dipole sizes. Furthermore, this is also fulfilled in some widely used dipole models, like the GBW model [61,62], where Γ Y (r) = r 2 Q 2 s (Y )/4. The role of separability in simplifying the results of the high-energy evolution in a similar context has been previously recognized in Ref. [35].
Within such a simplified scenario, the Y -dependence factorizes out from the mixing matrix and the resulting, Y -independent, matrix M can be diagonalized. Then one can explicitly solve the system of equations. Within the context of the MV model this was done in Refs. [8,31,35]. In that case, one had to deal with a 2×2 system only, because the Wilson lines were constructed via 'left' iterations alone, cf. Eqs. (3.16)- (3.17). That is, the associated, effective, evolution equations were not explicitly 'mirror-symmetric', unlike the above equations (4.12) and (4.15). However, this symmetry is recovered in the final results, because the color charge distribution in the MV model in symmetric in x − (cf. the discussion after Eq. (3.17)). Moreover, for a separable kernel, the final results for all the n-point functions depend only upon the integral Y dy γ y (a property which looks natural in the case of the dipole, cf. Eq. (4.9), but which in general requires separability). Hence, the results obtained in [8,31,35] within the MV model can be transposed to a more general Gaussian which is still 'separable', by simply replacing the kernel in the final formulae according to Eq. (4.17). We shall not write here the respective general results, but refer to [8] for Q x 1 x 2 x 3 x 4 Y and to [35] for But separability is not always needed in order to obtain explicit solutions at finite N c : for special configurations of the 4 external points x i in the transverse plane, the matrix M Y in Eq. (4.16) may happen to simplify independently of the structure of the kernel γ Y (x i , x j ). As an example consider the evolution of the 2-dipole S-matrix Ŝ x 1 x 3Ŝ x 3 x 2 Y , in which the quark in one dipole and the antiquark in the other dipole are located at the same point x 3 . The expression of this correlator in the Gaussian approximation has been shown in Eq. (3.41) and we would like to check that here. By identifying x 2 with x 3 in Eq. (4.15) and then relabeling x 4 as x 2 for convenience, we see that the two quadrupole operators in the r.h.s. there reduce to single dipoles and then this equation decouples from the evolution of the quadrupole : The last, inhomogeneous, term in the r.h.s. is particularly interesting, as it describes a process in which the two dipolesŜ x 1 x 3 andŜ x 3 x 2 having one common leg merge with each other into a single dipoleŜ x 1 x 2 . This process is suppressed at large N c (as expected, since dipoles cannot merge with each other in the large N c limit [27]) but for finite N c it controls the approach towards the unitarity limit, since a single dipole scatter less than a system of two dipoles. Using Eq. (4.8) to express γ Y (x i , x j ) in terms of the logarithmic derivative of Ŝ x i x j Y we can rewrite the above equation as where we temporarily defined ε = 1/(N 2 c − 1). This is a first order inhomogeneous linear differential equation which can be readily solved, with the result shown in Eq. (3.41). Some other special configurations, which in particular involve the quadrupole, will be studied in Sect. 4.4.
The generalization of the above considerations to an arbitrary n-point function like Eq. (2.9) is straightforward. For instance, the mean-field version of the equation obeyed by the sextupole S-matrix Ŝ (6) Y can be inferred from the results in Appendix B of Ref. [12].

Quadrupole evolution at large N c
In the large-N c limit, the hierarchy generated by the Gaussian approximation drastically simplifies: it reduces to a triangular hierarchy, in which the equations can be successively decoupled from each other and explicitly solved. Let us illustrate that on the example of the 4-point functions -the quadrupole Q x 1 x 2 x 3 x 4 Y and the 2-dipole system Ŝ x 1 x 2Ŝ x 3 x 4 Y -which in general mix under evolution, as shown in Eq. (4.16). At large N c , one can ignore the last two lines in the r.h.s. of Eq. (4.15), meaning that the 2-dipole system evolves independently of the quadrupole. The corresponding entries in the mixing matrix in Eq. (4.16) are now equal to zero, so that this matrix becomes triangular, as anticipated. Specifically, Eq. (4.15) reduces to which is immediately solved as and where we have defined As expected, this implies that the 2-dipole S-matrix factorizes at large N c provided it did so in the initial conditions at Y 0 : Consider now the evolution of the quadrupole: by keeping in (4.12) only the leading terms at large N c (which in particular means using the factorization property (4.23)), one finds which is an ordinary, first order, inhomogeneous differential equation. The dipole S-matrix is given by (4.9) with C F ≃ N c /2 (as appropriate at large N c ) and acts as a source for the evolution of the quadrupole. As explained in Sect. 4.1, in practice it is preferable to view Eq. (4.9) as a definition of the Gaussian kernel in terms of the dipole S-matrix, since the latter is a more directly relevant physical quantity. By using Eq. (4.2) to express γ Y (x 1 , x 2 ) in terms of the logarithmic derivative of Ŝ x 1 x 2 Y we can rewrite Eq. (4.24) as The general solution of this equation is easily found as [12] This solution is already explicit, but it takes an even simpler form if one assumes separability, in the sense of the discussion after Eq. (4.16). In that case, the integral over y in Eq. (4.26) can be exactly performed, to yield where we have denoted Notice that the function L also depends upon Y and Y 0 , but because of separability the ratio between two L's is a function of the transverse coordinates alone. Eq. (4.27) depends upon the kernel γ y (x i , x j ) only via its integral over y. This is a consequence of separability, as already noticed in Sect. 3.2 in the context of the MV model. As a matter of facts, Eq. (4.27) is quite similar to the corresponding expression in the MV model [3,8] and it becomes formally identical to it once we assume an initial condition of the MV type. Specifically, Eq. (4.27) reduces to 29) provided this functional relation is already satisfied at Y 0 , as is indeed the case in the MV model and for large N c [3,8]. Note that there is an alternative way to deduce Eq. (4.29) from Eq. (4.27), which makes no reference to the MV model. Namely, if one assumes Eq. (4.27) to capture the whole evolution from Y → −∞ (where the Wilson lines reduce to the unit matrix) up to the rapidity Y of interest, then one can use Q Y 0 → 1 and Ŝ Y 0 → 1 for Y 0 → −∞ to check that the expression within the square brackets in Eq. (4.27) vanishes for that particular initial condition. In general, i.e. without assuming separability, one expects the ratio of two L's to depend very weakly on Y . If so, it might be still a good approximation to use the simpler formula (4.29) for the quadrupole rather than the general one in Eq. (4.26), which is more involved. For that purpose, the function L in Eq. (4.29) should be defined by Eq. (4.28) with Y 0 → −∞, and hence Given a smooth approximation for Ŝ Y , such as the numerical solution to the BK equation, Eq. (4.26) (or (4.29)) provides a correspondingly smooth approximation for Q Y , which is guaranteed to be correct whenever all the transverse separations r ij ≡ |x i − x j | are either much smaller, or much larger, than 1/Q s (Y ), and for large N c . On the other hand, the present approximations are strictly speaking not under control in the transition region around saturation (r ij ∼ 1/Q s (Y )), nor for very asymmetric configurations, such that some distances r ij are much larger than 1/Q s while the other ones are much smaller. Some very asymmetric but relatively simple configurations will be discussed in the next subsection, directly for finite N c .

Special configurations at finite N c
In this subsection, we shall study some special configurations of the 4-point function and the 6-point function in the transverse plane, which because of their degree of symmetry allow for explicit, and relatively simple, solutions without additional assumptions like separability or large-N c . First we shall consider the class of configurations introduced in [12] for which the only constraints are r 13 = r 14 and r 23 = r 24 . For example, three such configurations are shown in Figs. 3.a, 3.b and 3.c. From these figures, it should be clear that there is a high degree of variability (concerning both shapes and sizes) within this particular class. By using the constraints aforementioned, it is straightforward to see that Eqs. (4.12) and (4.15) reduce to and Thus, for such configurations, the evolution of a system of 2 dipoles decouples from that of the quadrupole even without invoking separability. By also using Eq. (4.8) for the dipole S-matrix in the Gaussian approximation at finite N c , we can easily solve Eq. (4.32) to find which simplifies furthermore to provided the latter holds in the initial condition at Y 0 (as is indeed the case within the MV model, as one can similarly check). Then it is clear that Eq. (4.31) can be solved as an inhomogeneous first order differential equation and it gives .

(4.35)
Again, when the initial condition is given by the MV model (where Eq. (4.36) below holds indeed, as one can explicitly check), the above becomes very simple:  This result is truly remarkable: within the class of configurations at hand, the quadrupole factorizes into two dipoles independently of how small or large the various distances r ij are, and for any N c , so long as the two constraints r 13 = r 14 and r 23 = r 24 are satisfied. It would be interesting to check this factorization via numerical solutions to the JIMWLK equation, as a non-trivial test of the present MFA. Now let us proceed to our second example and consider the operator The choice is of direct phenomenological interest, since the operator above is the most complicated quantity appearing in the calculation of di-hadron production in proton-nucleus collisions 8 [2][3][4][5]. Considering the same configuration as before, that is, taking r 13 = r 14 and r 23 = r 24 , we see that the evolution couples the operator in Eq. (4.37) toŜ x 1 x 2Ŝ x 1 x 2Ŝ x 2 x 1 . After a straightforward calculation, similar to the one leading at Eq. (4.36), we arrive at (4.38) (Once again, we have assumed this equation to hold already in the initial condition at Y 0 , which is in particular true within the MV model, as one can check.) It is interesting that for this configuration, and similar to Eq. (4.36), the result depends only on r 12 and r 34 , but not on r 14 and r 23 . Also, by keeping N c finite, we notice that the last term in the r.h.s. (the single dipole S-matrix) dominates deeply at saturation. This is in agreement with our earlier discussion in 3.4, since the linear combination is a special case of the operator in Eq. (3.37) for this particular configuration. In fact, the operator that appears in the di-hadron cross section is the combination 40) which is also the quantity studied numerically in [11]. Using the previous results one finds that, for our special configuration, the expectation value ofŜ 6 is given by the relatively simple expression It is amusing to note that for this particular configuration and for large N c , the 4-point function relevant for di-hadron production factorizes as Ŝ 6 Y . This is precisely the factorization formula used (for a generic configuration) in the phenomenological study in [6] -at that time, by lack of a better formula. Such a factorization however has no deep justification and is merely a property of the configuration at hand. As our next example will show, this 'factorization' can badly fail for other, equally simple, configurations.
Specifically, let us consider the expectation value of the operator (4.37) for the 'line' configuration studied in [11,12] and shown in Fig. 3.d. Note that the two quarks of the quadrupole and the antiquark of the dipole are put in a same point (x 1 = x 3 ), and similarly for the two antiquarks of the quadrupole and the quark in the dipole (x 2 = x 4 ). Thus, only one non-trivial distance r ≡ r 12 = r 23 = r 34 = r 14 characterizes the configuration. Then one can easily check that the evolution ofQ x 1 x 2 x 1 x 2Ŝ x 2 x 1 couples again toŜ x 1 x 2Ŝ x 1 x 2Ŝ x 2 x 1 , leading to a 2×2 inhomogeneous system of equations. Expressing γ Y (r) in terms of the dipole Ŝ (r) Y and using an obvious shorthand notation we have The solution to this system is straightforward to obtain; so long as QŜ Y is concerned, one finds where we assumed that the above is already valid at Y 0 , as is the case in the MV model. Using this result together with Eq. (4.40), it is straightforward to evaluate Ŝ 6 Y for this particular configuration and compare with the numerical results in Ref. [11]. We shall find it rewarding to plot Ŝ 6 Y in two different ways; first as a function of rQ s and then as a function of 1 − Ŝ Y . To be in accordance with [11], the saturation momentum is defined by the condition Ŝ Y = 1/ √ e for rQ s = √ 2. We show this comparison in Fig. 4, where for some of the curves we have used the numerical data of [11]. On the left we show the correlator of interest as a function of rQ s and for two for illustrative purposes we also show Ŝ 3 Y . JIMWLK curves are constructed from the numerical solution given in [11]. MFA curves are analytical expressions in terms of Ŝ , which is again provided by the numerical solution in [11] for the purposes of the left figure. values of the rapidity: Y = 0, which is where one starts the evolution (with initial conditions of the MV type), and Y = 5.18, which is large enough for the effects of the evolution to be fully developed. The curves denoted as 'MFA' represent our present results, cf. Eq. (4.44) and (4.40), whereas the 'JIMWLK' curves follow from the numerical solution to the JIMWLK equation. For Y = 0, the two types of curves overlap with each by construction, as they both reduce to the respective prediction of the MV model. What is remarkable though, is that a very good agreement between the numerical solution and the MFA persists for Y = 5.18. In the limiting regimes of weak and respectively strong scattering, the respectives curves are practically indistinguishable. In the transition region around rQ s ∼ 1, the agreement is not that perfect anymore, but the two curves are still very close to each other, confirming that the MFA is also an excellent global approximation.
From Fig. 4 (left) we also notice that the shape of the curve changes as we evolve from Y = 0 to higher values of rapidity (Y = 5.18 in the figure). However, this change is mostly attributed to the evolution of Ŝ Y as a function of rQ s as the rapidity grows. Indeed, in the right panel of Fig. 4 we show the correlator of interest as a function of 1 − Ŝ Y . One can see that the curves obtained from the numerical solution to JIMWLK for various values of Y form a very thin "band" whose borderline on the "lower" side is the MFA. This "band" is practically a line perfectly overlapping with the MFA when the scattering is either weak or strong, and becomes just a bit wider in the transition region. Thus, as we evolve in rapidity, the shape of the curve is barely changing. In fact, if one expands out the plot in order to better disentangle the various steps in the evolution, one can see that the high-Y curve stabilizes very close to the Y = 0 curve after just a few units in rapidity.
Still in the right plot we also show two different approximations for Ŝ 6 Y , the one is the large-N c result, while the other consists of factorizing QŜ Y into Q Y Ŝ Y , as it would be justified at large N c , but then using the finite-N c Gaussian approximation for Q Y . This latter approximation takes into account some 1/N 2 c corrections, but not in a systematic way, and was used in [11] since the (MV-like) expression (4.44) was not available at the time. Comparing with the numerical findings in [11], we already saw that the complete result Eq. (4.44) at finite-N c is the one which shows the best agreement. Even though it is not very significant, we note that the large-N c expression is the next one closer to the numerical data, perhaps because it is at least a systematic approximation. For illustrative purposes we also show Ŝ 3 Y , which is based neither on a large-N c approximation nor on a mean field one, but simply corresponds to a 'naive' counting of Wilson lines. It fails badly even in the BFKL regime and clearly it has no chance to describe properly the correlator of interest.