Higher-point correlations from the JIMWLK evolution

We develop a new approximation scheme aiming at extracting higher-point correlation functions from the JIMWLK evolution, in the limit where the number of colors is large. Namely, we show that by exploiting the structure of the 'virtual' terms in the Balitsky-JIMWLK equations, one can derive functional relations expressing arbitrary n-point functions of the Wilson lines in terms of the 2-point function (the scattering amplitude for a color dipole). These approximations are correct not only in the regime of strong scattering, where the evolution is indeed controlled by the 'virtual' terms, but also in the regime of weak scattering, where they reduce to the corresponding BFKL solutions. This last feature follows from the fact that the JIMWLK Hamiltonian is a linear combination of the pieces responsible for the 'real' and 'virtual' terms, respectively. We apply this scheme to two examples: the 'color quadrupole' (the 4-point function of the Wilson lines which enters the cross-section for the production of a pair of jets at forward rapidities) and the 'color sextupole' (the 6-point function). For particular configurations of the quadrupole, our general formula reduces to relatively simple expressions that have been previously proposed on the basis of the McLerran-Venugopalan model and which were recently shown to agree quite well with exact, numerical, solutions to the JIMWLK equation.


Introduction
Since the derivation, more than a decade ago, of the equations describing the high-energy evolution in QCD to leading logarithmic accuracy -the infinite hierarchy of Balitsky equations [1] for the n-point functions of the Wilson lines (= scattering amplitudes in the eikonal approximation) and the functional JIMWLK 1 equation [2][3][4][5][6][7][8][9] for the color glass condensate (CGC = the small-x part of the wavefunction of an energetic hadron) -, there has been little progress towards solving these equations for correlations which are more complicated than the 2-point function (the scattering amplitude of a color dipole). The progress has been mostly hindered by the extreme complexity of these equations (an issue that we shall shortly return to), but also by the fact that, for quite some time, the dipole amplitude was the only such a correlation to be directly relevant to phenomenology (e.g., for deep inelastic scattering and for single-inclusive particle production in hadron-hadron collisions; see e.g. [10][11][12] and references therein).
However, the situation has changed in the recent years, with the advent of new, less inclusive, data, which probe higher-point functions via the multiparticle correlations in the particle production at high energy. In particular, the recent data at RHIC, on the azimuthal correlations in the forward di-hadron production in deuteron-gold collisions [13], are sensitive to a 4-point function of the Wilson lines known as the 'color quadrupole' [14][15][16][17]. These data, together with their successful description by a phenomenological analysis using the general ideas of the CGC [18], have revived the interest in the higher-point correlations and triggered new efforts in that sense, mostly in relation with the quadrupole [19][20][21][22]. But these efforts have been hampered by the main difficulty alluded to above -the extreme complexity of the Balitsky-JIMWLK equations -, which so far has precluded the obtention of any analytic solution, even approximate. It is our purpose in this paper to present a major step in that sense, in the form of an approximation scheme for the higher-point functions with n ≥ 4, that we believe to be new. But before describing our approach, let us briefly recall the state of the art.
The only analytic estimate for the quadrupole which is currently available within the CGC framework is that obtained within the McLerran-Venugopalan (MV) model [23,24], which refers to a large nucleus at not too high energy. This estimate has first been computed in the limit of a large number of colors N c in Ref. [14], and more recently for generic N c in Ref. [20]. The evolution equation for the quadrupole is also known: this has first been derived in Ref. [14] at large N c , by using the dipole picture [25], and then by several authors [21,26,27] including ourselves, for generic N c , using the JIMWLK Hamiltonian. (Our version of the derivation is presented in Appendix A.) But this is only one out of an infinite hierarchy of coupled equations, which moreover exhibit a complicated non-local structure in the transverse coordinates. This non-locality becomes more and more severe with increasing the number n of external points.
The Balitsky-JIMWLK hierarchy considerably simplifies in the large-N c limit, where it reduces to a triangular hierarchy. It is then enough to solve a finite number of equations in order to determine a given correlation function. The dipole obeys a closed, non-linear, equation, known as the Balitsky-Kovchegov (BK) equation [1,28]. The quadrupole evolves according to an inhomogeneous equation, in which the dipole plays the role of a source term. The equation for the sextupole (a color trace of six Wilson lines) also involves the dipole and the quadrupole, and so on. But even in that limit, the non-locality of the equations is such that it renders prohibitive their (numerical or analytic) study, except in the simplest case: the BK equation, which has been studied at length [10][11][12]. As a matter of facts, the most promising approach towards numerical studies refers to the functional JIMWLK equation by itself, and not to the ordinary evolution equations in the hierarchy: this functional equation can be reformulated as a Langevin process [29], which is well suited for the numerics. The feasibility of such numerical solutions has been demonstrated in Refs. [30][31][32], which however focused on the 2-point function (the dipole amplitude) alone. It was only very recently -when the present study was under progress -that this numerical method has been extended to the quadrupole [22].
The analysis in Ref. [22] follows the evolution of the quadrupole with increasing energy, starting with initial conditions of the MV type and for two special configurations of the four external points in the transverse plane: the 'line' and the 'square'. (These configurations are illustrated in Figs. 1.a and b below.) Interestingly, the results show a good agreement with a heuristic extrapolation to high energy of a formula, in particular an expression for the quadrupole in terms of the dipole, that was derived within the MV model [14,20]. Such an agreement is hard to explain without a more fundamental understanding of the relation between the extrapolation allude to above and the actual JIMWLK dynamics. Moreover, the continuation of the numerical solutions in [22] to more systematic studies (say, in view of the phenomenology) can be very difficult and demanding in computer power. Thus, it is necessary to remedy such shortcomings by completing the numerical solutions with controlled analytic approximations -a task that we shall accomplish here.
Specifically, our strategy will consist in solving a simplified version of the Balitsky-JIMWLK equations at large N c , in which we keep the 'virtual' terms, but we drop the 'real' ones. The distinction between 'real' and 'virtual' is meant here in the same sense as in the dipole picture [25], that is, it refers to the evolution of the projectile (dipole, quadrupole, etc): the 'real' terms describe processes in which the projectile splits via the emission of an additional, soft, gluon (at large N c , this yields e.g. a color dipole splitting into two dipoles), whereas the 'virtual' terms refer to the complementary probability that the projectile survive without splitting. By themselves, the 'virtual' terms control the evolution of the projectile S-matrix in the vicinity of the unitarity (or 'black disk') limit, so our approximations are a priori justified in the strong scattering regime deeply at saturation. But our scope in this paper is more general than that: by exploiting the structure of the 'virtual' terms, we shall derive explicit expressions for the quadrupole and the sextupole in terms of the dipole, which represent global approximations, valid for both weak and strong scattering. We shall indeed verify that, in the regime where the scattering is weak, our general results reduce to the expected, linear, relations between the quadrupole, or the sextupole, and the dipole scattering amplitude. This is not an accident: it follows from the fact that a linear relation is always preserved by an evolution equation derived from a Hamiltonian -in particular the 'virtual' part of the JIMWLK Hamiltonian. (We recall here that the JIMWLK Hamiltonian is the direct sum of the pieces responsible for the 'real' and 'virtual' terms, respectively.) Hence, by using the BFKL approximation for the dipole scattering amplitude within our general expression, we recover the correct BFKL results for the quadrupole and the sextupole (although the BFKL evolution is sensitive to both 'real' and 'virtual' terms, of course).
Our method is general and systematic: by using similar techniques, it is possible to derive expressions for all the n-point functions of the Wilson lines in terms of the dipole S-matrix. Furthermore, our general expressions are also consistent with initial conditions of the McLerran-Venugopalan type (at large N c ). Hence, they provide a unified description of the initial conditions and of their high-energy evolution, in both the dilute (BFKL) and dense (saturation) regime. In our opinion, the ultimate reason why such a strategy can work is the fact that, within the JIMWLK evolution, the multi-gluon correlations are built exclusively via high-density effects, that is, via gluon recombination in the approach towards saturation.
The numerical analysis in Ref. [22] provides already a test of our approximation scheme: for the particular quadrupole configurations investigated there (the 'line' and the 'square'), our general result coincides with the high-energy extrapolation of the respective MV formulae, that in Ref. [22] were found to agree quite well with the numerical solutions to the JIMWLK equation. The agreement found there is even better when using the finite-N c version of the MV formulae, showing that the finite-N c corrections to the evolution of the quadrupole yield somewhat sizeable effects when N c = 3. This should be contrasted with the corresponding situation for the dipole, where one has numerically found [30,31] that the finite-N c corrections in the JIMWLK evolution are surprisingly small -at the level of the percent accuracy instead of the 10% (≈ 1/N 2 c with N c = 3) that would be normally expected. Let us finally mention some possible limitations of our present analysis, which could be improved by further studies. A priori, our approximations are better justified for relatively symmetric configurations of the n-point functions, which are such that the transverse separations between the opposite-sign charges (i.e. between quarks and antiquarks) are not very different from each other. It would be interesting to test this limitation by comparing our predictions for very asymmetric, 'small-large', configurations -which turn out to be very interesting (see the discussion in Sect. 4) -to numerical solutions to the JIMWLK equation. We also note that, by construction, our approximations are guaranteed to be correct in the weak scattering regime and in the approach towards the black disk limit, but not necessarily also in the transition between these two regimes, as occurring around the saturation scale. Still, our results provide a smooth (infinitely differentiable) interpolation between the two limiting regimes and as such we believe them to be qualitatively and even quantitatively correct including in the transition region. This is again in agreement with the numerical analysis in Ref. [22].
The paper is organized as follows. In Sect. 2 we present the evolution equations for the dipole and the quadrupole that we shall use as examples to illustrate general properties of the Balitsky-JIMWLK hierarchy, like the relation between 'real' and 'virtual' terms and the simplifications which appear in the large-N c limit. (The respective equations will be explicitly derived in Appendix A and this will give us the opportunity to describe the method that we shall later use, in Appendix B, to construct the corresponding equation for the sextupole.) In Sect. 3 we proceed with a study of the evolution deeply as saturation, as described by the virtual terms, and introduce our approximation scheme on a simple, pedagogical, example: the 'line' configuration for a quadrupole. In Sect. 4 we construct our global approximation for a generic configuration of the quadrupole and study some limiting cases corresponding to simple, but interesting, configurations. In particular, we shall identify configurations for which the quadrupole is exactly factorizing in the product of two dipoles. The corresponding analysis of the sextupole, including its evolution equation, is given in Appendix B.

Evolution equations
In this section, we shall review the Balitsky-JIMWLK evolution equations for the dipole and the quadrupole and explain the origin and the physical significance of the various terms which appear in these equations. This discussion will be useful in view of the construction of an approximation scheme later on.
Within the CGC effective theory, valid to leading logarithmic accuracy at high energy, the evolution of a gauge invariant operatorÔ with increasing energy is determined by where the brackets refer to the average over the color fields in the target, as computed with the CGC weight function, and H is the JIMWLK Hamiltonian. The latter, when acting on gauge invariant observables, is most conveniently given in [33] and reads (2.2) In the above equations, α a (x − , x) is the target color gauge field in the covariant gauge, where and V † and V are Wilson lines in the adjoint representation: with P denoting path-ordering in x − . Our conventions are such that the nuclear target (the CGC) is a right-mover, whereas the projectile (to which refers the operatorÔ) is a left-mover.
The functional derivative in Eq. (2.2) are understood to act at the largest value of x − , that is, at the upper end point of the path-ordered exponential in Eq. (2.4). The operators that we shall deal with are all built with products of Wilson lines (in the fundamental representation, for definiteness) and represent the S-matrices for systems of partons (quark and antiquarks) scattering off the nuclear target, in the eikonal approximation. Specifically, we shall focus on the color dipole -a quark-antiquark pair in an overall color singlet state, with S-matrix operator the color quadrupole -a system of two quarks and two antiquarks, for whicĥ and the color sextupole, witĥ Higher-point correlators can be similarly considered, including those which involve the product of several color traces (see e.g. Eq. (2.12) below). In Eqs. (2.5)-(2.7), V † and V are Wilson lines in the fundamental representation. The action of the functional derivative on these Wilson lines reads (2.8) By using these rules within Eqs. (2.1)-(2.2), it is straightforward to derive the evolution equations satisfied by the expectation values of the three operators introduced above. A streamlined derivation of the respective equations for the dipole and the quadrupole will presented in Appendix A, while the corresponding equation for the sextupole will be shown (without the details of the derivation) in Appendix B. The ensuing equation for the dipole looks relatively simple (we denoteᾱ ≡ α s N c /π): ∂ Ŝ but that for the quadrupole is considerably more involved: We shall now discuss the role of the various terms in the above equations, in order to explain, in particular, the difference between 'real' and 'virtual' terms. Consider first Eq. (2.9) for the dipole S-matrix. The 'real' term is the first term in the right hand side, which is quadratic inŜ. When interpreted in terms of projectile evolution, this term 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. The 'virtual' term, that is, the negative term which is linear in S, describes the reduction in the probability that the dipole survive in its original state -that is, the probability for the dipole not to split. The relative size of these two types of terms ('real' and 'virtual') is constrained by probability conservation, as correctly encoded in the JIMWLK Hamiltonian. Remarkably, the 'real' term has been fully generated by the last two terms in the Hamiltonian (2.2), whereas the 'virtual' term originated from the first two terms there. All the terms appearing in Eq. (2.9) are of leading order in N c . Subleading terms, of relative order 1/N 2 c , had been also generated, at intermediate steps, by each of the four terms in the Hamiltonian, but they have exactly canceled after summing up all the contributions. Note the cancelation of 'ultraviolet' (i.e. short-distance) singularities between the 'real' and the 'virtual' terms: the dipole kernel (2.3) has poles at z = u and z = v, which within the integral over z generate logarithmic divergences, separately for the 'real' and for the 'virtual' term. However, these divergences cancel in the overall equation, because of probability conservation together with the property of 'color transparency' :Ŝ x 1 x 2 → 1 when x 1 → x 2 . The latter is merely the statement that a 'color dipole' of zero transverse size is a totally neutral object, which cannot interact.
An entirely similar discussion applies to the evolution equation (2.10) for the quadrupole. 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 generated by the individual pieces of the Hamiltonian) have canceled in the final equation. Furthermore, all possible ultraviolet divergencies due to the poles of the dipole kernel cancel between the 'real' and the 'virtual' terms, due to color transparency.
All the above features are generic: they also apply to the Balitsky-JIMWLK evolution of the sextupole and, more generally, to any (gauge-invariant) operator which involves a single trace of products of pairs of Wilson lines, of the form As for the multi-trace operators, of the form the only new feature is that the respective evolution equations 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 [34] for an example).
As manifest by inspection of Eqs. (2.9) and (2.10), these equations are generally not closed: 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, so altogether one has to deal with an infinite hierarchy of equations. But important simplifications occur in the large-N c limit, as anticipated in the Introduction. Then, for a multi-trace operator like (2.12), it is not hard to show that the hierarchy admits the factorized solution provided this factorization is already satisfied by the initial conditions. Then the hierarchy becomes triangular and therefore finite: Eq. (2.9) becomes a closed equation for Ŝ Y (the BK equation), while Eq. (2.10) becomes an inhomogeneous equation for Q Y with coefficients which depend upon Ŝ Y . This is the equation originally derived in [14] and that we shall further study in the next section.

From virtual terms to global approximations
In order to better appreciate the approximation scheme that we shall eventually propose, it is useful to investigate the behavior of the solutions to the evolution equations introduced in the previous section in two limiting regimes: the dilute, or weak scattering, regime, where all the transverse separations r ij ≡ |x i − x j | between the external points are very small compared to 1/Q s (Y ), and the dense, or strong scattering, regime, where all these distances are much larger than 1/Q s (Y ). Here, Q s (Y ) ∝ e ωᾱY is the saturation momentum in the target, and it is also the transverse momentum scale which marks the transition from weak to strong scattering, for both the dipole and the quadrupole. The discussion of the dilute regime is straightforward. Introducing the dipole T -matrix operatorT x 1 x 2 ≡ 1 −Ŝ x 1 x 2 and the corresponding scattering amplitude T x 1 x 2 Y , then for r = |x 1 − x 2 | ≪ 1/Q s the scattering is weak, T x 1 x 2 Y ≪ 1, and the amplitude obeys the linearized (in T Y ) version of Eq. (2.9), that is, the BFKL equation: Consider similarly the quadrupole: when r ij ≪ 1/Q s for all the six pairs (i, j) (with i, j = 1, 2, 3, 4) of external points, the scattering is necessarily weak, because there is no net color charge in the projectile over an area ∼ 1/Q 2 s (which is the area where high gluon density effects become important in the target). Then, for consistency, it is enough to keep the lowest order terms, of order g 2 α 2 , in the perturbative expansion of 1− Q Y . These are obtained by expanding the Wilson lines in Eq. (2.6) up to quadratic order. After also comparing with the corresponding expansion ofŜ in Eq. (2.5), one finds where it is understood that, in the r.h.s., T Y is the dipole amplitude in the dilute regime and obeys the BFKL equation (3.1). Eq. (3.2) also tells us that, for generic configurations at least, the quadrupole scattering can become strong when at least one (which necessarily means at least three) of the six transverse distances r ij is of order 1/Q s , or larger. However, care must be taken when discussing very asymmetric configurations (see below).
Consider now the strong scattering regime at r ij ≫ 1/Q s (Y ), where we shall study the approach towards the 'black disk' limit, i.e. the way how S-matrices like Ŝ Y and Q Y approach to zero. The corresponding analysis for the dipole is well known [35,36], but it will be briefly reviewed here, in preparation for the quadrupole and higher n-point functions. The crucial observation is that, at least for large N c , this approach is controlled by the virtual terms in the respective evolution equations 2 . Indeed, at large N c we have, e.g., ŜŜ Y ≃ Ŝ Y Ŝ Y , and then the 'real' term in the r.h.s. of Eq. (2.9), which is quadratic in Ŝ Y , is much smaller than the linear, 'virtual', term there whenever Ŝ Y ≪ 1. More precisely, when r = |x 1 − x 2 | ≫ 1/Q s (Y ), the linear term dominates over the quadratic one unless either |x 1 − z| or |z − x 2 | is as small as 1/Q s (Y ). Then, the r.h.s. of Eq. (2.9) is dominated by the logarithmic regions of integration where one has This discussion implies that the approach of the dipole S-matrix towards the 'black disk' limit is determined by the following, simplified, equation which is correct to leading logarithmic accuracy with respect to the large logarithm ln(r 2 Q 2 s ). That is, in writing the r.h.s. of Eq. (3.4), we control the coefficient in front of the logarithm but not also the constant term under the logarithm. Using the fixed coupling relation ln[ This expression is valid for Y ≥ Y 0 , where Y 0 is roughly the rapidity at which the scattering of the dipole with size r becomes strong. We now turn to the quadrupole and similarly study the strong scattering regime for generic configurations of the four external points in the transverse plane. Note however that our ultimate objective in this study is different, and actually more ambitious, than in the corresponding study of the dipole: what we are truly interested in, is not the law for the approach towards the black disk limit (although this law will emerge too from our subsequent analysis), but rather the functional relation between the quadrupole S-matrix Q Y and the dipole one Ŝ Y , as predicted by the solutions to the respective evolution equations with the virtual terms alone. Indeed, as announced in the Introduction and will be explicitly checked on the final results, this functional form is correct also outside the strong scattering regime -namely, it has the right limit at weak scattering, as shown in Eq. (3.2).
Note first that a quadrupole configuration is bound to be in the strong scattering regime whenever all the four quark-antiquark separations (r 12 , r 14 , r 23 and r 34 ) are much larger than 1/Q s , independently of the relative separations between the two quarks (r 13 ) and the two antiquarks (r 24 ). Indeed, when r qq ≫ 1/Q s , there is uncompensated color charge over distances 1/Q s , which implies that the scattering is strong. In this regime and for large N c , the evolution towards the unitarity limit is controlled by the virtual terms, by the same arguments as for the dipole. After neglecting the 'real' terms, Eq. (2.10) simplifies considerably and, in particular, it becomes 'diagonal' with respect to the quadrupole configuration: both the l.h.s. and the r.h.s. involve the same configuration (x 1 , x 2 , x 3 , x 4 ). This makes it possible to study the approach towards the unitarity limit configuration by configuration and thus obtain explicit, analytic solutions.
The case of a general configuration (which respects the conditions stated above) will be addressed in the next section. But before that general discussion, it is instructive to consider a simple, particular case as a warm-up. In what follows, we shall often use the four special configurations shown in Fig. 1 in order to illustrate our results. These configurations are useful in that they involve only few independent transverse separations r ij , due to their high degree of symmetry. Besides, two of them -the 'line' in Fig. 1.a and the first 'square' in Fig. 1.bhave been previously chosen in the numerical studies in [22], so for them we know already the respective numerical predictions of the JIMWLK equation. This information will be useful for testing our forthcoming analytic results.
We shall pick the 'line' configuration, where we place the two quarks on top of each other and the same for the antiquarks (see Fig. 1.a), as our warm-up example. For this example, we shall perform in detail all the steps to be repeated for the general case in the next section. (The three other configurations in Fig. 1 will be studied in Sect. 4, as special limits of our general result for the quadrupole.) As obvious from Fig. 1.a, the 'line' configuration involves a single transverse scale r, which is the common size of the four qq pairs: r 12 = r 14 = r 23 = r 34 ≡ r.
When applied to this configuration, the general evolution equation (2.10) reduces to When the separation r is much larger than 1/Q s , the virtual terms (the last two terms in Eq. (3.6)) dominate and give We have written the argument of the quadrupole simply asQ(r), because r is the only independent transverse scale for this particular configuration. After also using Eq. (3.4), it turns out that this inhomogeneous equation admits a relatively simple solution, which is (for Y ≥ Y 0 ) Fig. 1.a) . (3.8) In view of its derivation above, one may expect this formula to hold only deeply at saturation, that is, for Y > Y 0 , with Y 0 the rapidity at which saturation sets in on the given transverse scale r. But as a matter of facts, Eq. (3.8) is more general than that: it has also the right limit at weak scattering, in agreement with Eq. (3.2), and hence it provides a global approximation. To see that, let us study the weak scattering limit of Eq. (3.8), by writing Ŝ (r) Y = 1 − T (r) Y , with T (r) Y ≪ 1, and then expanding the r.h.s. to linear order in T (r) Y . After doing that, one finds that the following relation is satisfied at Y , provided it was already satisfied at Y 0 : As anticipated, Eq. (3.9) is the correct weak-scattering limit for this particular configuration, as predicted by Eq. (3.2). Eq. (3.8) can be further simplified if one assumes that the initial conditions at Y 0 are provided by the MV model at large N c . Indeed, it is easy to check that Eq. (3.8) reduces to provided the above formula already holds at the initial rapidity Y 0 -a property which is indeed satisfied within the MV model at large N c . To summarize, by using Eq. (3.10) together with a reasonable approximation for the dipole S-matrix Ŝ (r) Y -say, the solution to the BK equation with initial conditions provided by the large-N c version of the MV model -one obtains a unified description for (i) the initial conditions at low energy, (ii) the BFKL regime at high energy, and (iii) the approach towards the black disk limit at all the energies.
As demonstrated by the numerical analysis in Ref. [22], Eq. (3.10) agrees well with the numerical solution to the JIMWLK equation for this particular configuration. In the forthcoming section and in Appendix B, we shall generalize Eq. (3.8) to generic configurations for the quadrupole and the sextupole, and we shall also study some other particular configurations.

A global approximation for the quadrupole
We shall now extend the analysis in the previous section to a generic configuration for the quadrupole -that is, we shall study the approximation obtained by keeping only the virtual terms in Eq. (2.10). In line with our general strategy, we shall rely on this approximation, which is strictly speaking valid only in the vicinity of the black disk limit, in order to derive a functional relation between the quadrupole and the dipole, that will be then promoted to a global approximation.
After neglecting the 'real' terms and restricting the integrations over z to the regions yielding large logarithmic contributions, like in Eq. (3.3), we get where in writing the r.h.s. we have also factorized the 2-dipole S-matrix, as appropriate at large N c . Now we can use Eq. (3.4) in order to express the logarithms in terms of the dipole S-matrix. The merit of that is three-fold. First, we can express the quadrupole in terms of the dipole without worrying about the explicit solution for the latter. Second, in order to get the dominant logarithmic contribution we do the same kind of approximation in both the dipole and the quadrupole equations since their structures are the same. Therefore we probably have better control than originally announced, that is, we have better accuracy than the one given by the leading logarithmic terms. Third, we do not really need to specify the energy dependence of Q s and thus the approach is fully valid even if we consider a running coupling scenario (at least in the case whereᾱ is evaluated at Q s ). Thus, writing all logarithms in Eq. (4.1) terms of S we have This is an ordinary first order inhomogeneous differential equation whose general solution is straightforward to find; it reads As explained in the Introduction, the above expression is guaranteed to be correct also for weak scattering, even though it has been obtained here via manipulations which are well justified only deeply at saturation. This is so since, at weak scattering, the relation between the quadrupole and the dipole becomes linear, cf. Eq. (3.2). Such a linear relation is then preserved by any evolution equation generated by a Hamiltonian -in particular, by the equations including only the virtual terms, which are generated, as we have seen, by a truncated version of the JIMWLK Hamiltonian (the first two terms in Eq. (2.2)). This argument about linearity is simple, but at the same time important for the present construction, so let us make a small digression in order to explain it better. Consider a set of three operators,Â,B, andĈ, which obey a linear relation:Â = c 1B +c 2Ĉ . Also assume that the evolution of these operators with 'time' Y is governed by some Hamiltonian H; e.g. ∂ YÂ = HÂ. Then one can successively write This equation implies that the linear relationÂ = c 1B + c 2Ĉ is preserved at any Y provided it was satisfied at the initial 'time' Y 0 . Note that specific form of H was irrelevant for the above argument; all that matters is the fact that this is a linear operator.
Returning to the physics problem at hand, the relation (4.3) between the quadrupole and the dipole has been obtained by solving a Hamiltonian evolution equation. The corresponding Hamiltonian is incorrect in the dilute regime, where it would predict a wrong evolution for the quadrupole, or the dipole, taken separately. However, as a linear operator, it must preserve the linear relation (3.2) between the two operatorsQ andŜ which is valid in that regime. And indeed, it can be easily checked -by rewriting Ŝ x i x j y = 1 − T x i x j y and then expanding the r.h.s. of Eq. (4.3) to linear order in the various T x i x j y 's -that Eq. (4.3) predicts the correct relation, Eq. (3.2), between the quadrupole and the dipole at weak scattering and at rapidity Y provided this relation was already satisfied at the initial rapidity Y 0 .
In view of the above, we conclude that Eq. (4.3), when used with a numerical solution of the BK equation (or some good approximation to it), provides a reliable approximation to the solution to the quadrupole equation which is valid in all regimes (at least for generic configurations, which are not very asymmetric). This solution provides a smooth (infinitely differentiable) interpolation between the BFKL solution in the weak scattering regime at r ij ≪ 1/Q s and the correct approach towards the black disk limit at r ij ≫ 1/Q s .
The expression (4.3) has some other good properties. It is symmetric under the interchange of the two quarks (or the two antiquarks) at any Y , provided that this is true at the initial rapidity Y 0 . This is a true property of the JIMWLK evolution at large-N c , as it can be directly checked at the level of the evolution equations in Sect. 2. Furthermore, there is an interesting relation between Eq. (4.3) and the corresponding expression in the MV model [14,20]. To see this connection, it is useful to write the dipole S-matrix in the form Now, (4.3) becomes formally identical with the quadrupole formula in the MV model (for large N c ) provided one assumes Γ Y (x i , x j ) to be a separable function of Y and the transverse coordinates, plus an arbitrary function of Y . This property is certainly not satisfied in general, but it is locally satisfied under some approximations -e.g., in the window for extended geometric scaling, where Γ Y (r) ≃ (r 2 Q 2 s (Y )) γs with γ s ≈ 0.63 (modulo more slowly varying logarithmic terms) [41][42][43], and deep at saturation where Γ Y (r) ≃ (1/2ω) ln 2 [r 2 Q 2 s (Y )] (cf. Eq. (3.5)). This is also fulfilled in some dipole models, like the GBW model [44,45], where Γ Y (r) = r 2 Q 2 s (Y )/4. Assuming this to be the case, then it is possible to explicitly perform the integral over y in Eq. (4.3) and thus obtain, after some algebraic manipulations, which is indeed the same as the corresponding expression in the MV model at large N c [14,20]. Notice that the quadrupole at Y 0 does not appear anymore since in writing Eq. (4.6) we have assumed that the latter is already valid at Y 0 , as correct for initial conditions of the MV type.
Returning to the general formula (4.3), let us notice that there are special quadrupole configurations, like those shown in Fig. 1, which have enough symmetry for the integral over y to be exactly doable, without additional assumptions. For instance, consider the case where all transverse separations between a quark and an antiquark are the same, that is r 12 = r 34 = r 14 = r 32 ≡ r. (This includes the 'line', Fig. 1.a, and also the first 'square', Fig. 1.b.) Then the integrand in Eq. (4.3) is a total derivative, so one easily finds The arguments of Q Y , which here is a function of r, r 13 , and r 24 , are kept implicit. The 'line' configuration in Fig. 1.a corresponds to r 13 = r 24 = 0, and then Eq. (4.7) reduces to Eq. (3.10), as expected. As for the 'square' configuration in Fig. 1.b, one has r 13 = r 24 = √ 2r, so Eq. (4.7) gives Fig. 1.b) . (4.8) Notice that the above reduce to provided Eq. (4.9) is also valid at Y 0 , as it is indeed the case in the MV model at large-N c . More generally, Eq. (4.9) is valid at Y even if it was not valid at Y 0 provided the scattering of a dipole with the given size r was weak at Y 0 . (But of course, the scattering can be strong at Y > Y 0 for the same size r.) Indeed, assuming that T (r) Y 0 ≪ 1, and therefore also T ( √ 2r) Y 0 ≪ 1, it is straightforward to show that the pieces linear in T Y 0 coming from the first and the last term of Eq. (4.8) mutually cancel and then Eq. (4.9) is again obtained.
In Ref. [22], Eq. (4.9) has been compared to the numerical solution to JIMWLK equation and the agreement appears to be remarkably good -actually, even better than for the 'line' configuration. Clearly, expressions like (4.9) or (3.10), which involve a single transverse scale r, exhibit geometric scaling in the same regime where the dipole does so. This feature too is in agreement with the numerical results in Ref. [22] A second class of configurations for which the integrand in Eq. (4.3) is a total derivative are those which are constrained by r 13 = r 14 and r 23 = r 24 , whereas the two other separations r 12 and r 34 are left arbitrary. The two last configurations in Fig. 1 -the second 'square', Fig. 1.c, and the 'double triangle', Fig. 1.d -belong to this category. For any configuration in this class, Eq. (4.3) yields a very simple, factorized, expression which involves the unconstrained separations r 12 and r 34 , but it is independent of the other ones (those which are constrained). As before, Eq. (4.10) reduces to an even simpler expression 11) provided this last formula holds already at Y 0 (a property which, once again, is satisfied within the MV model at large N c ). In particular, for the second 'square' configuration in Fig. 1.c, one has r 13 = r 14 = r 23 = r 24 ≡ r and r 12 = r 34 = √ 2 r, and Eq. (4.11) gives Fig. 1.c) . But the most interesting configuration in our opinion is the last one in Fig. 1: the 'double triangle'. What is remarkable about this configuration is the fact that the separation between the qq pair (x 1 , x 2 ) and the second qq pair (x 3 , x 4 ) can be made arbitrary, yet the corresponding S-matrix in Eq. (4.10) is independent of this separation. In particular, so long as r 12 and r 34 are much smaller than 1/Q s , the scattering of the quadrupole as a whole remains weak independently of the distance ∼ r 14 separating the two pairs -including in the case where r 14 ≫ 1/Q s (Y ). This may look counterintuitive since normally one expects the scattering to be strong whenever the separation between color charges is of order 1/Q s or larger. However, we believe that this result can be physically understood as follows: when the two quark-antiquark pairs (x 1 , x 2 ) and respectively (x 3 , x 4 ) are individually small, whereas the distance ∼ r 14 between them is of order 1/Q s or larger, the color exchanges between the two pairs are strongly suppressed by saturation; then, the color charges compensates locally, within each of the small qq pairs, which therefore behave like two individual dipoles. In view of this argument, the factorized structure of Eq. (4.10) looks quite natural. It would be very interesting to verify this formula via numerical simulations of the JIMWLK evolution, along the lines in Ref. [22].
For pedagogy, we shall first derive the evolution equation for a color dipole, that is, the 2-point function of the Wilson lines shown in Eq. (2.5). The action of the functional derivative on these Wilson lines is computed according to Eq. (2.8), thus yielding where in writing the r.h.s. we have dropped some terms proportional to δ uv since the kernel M uvz vanishes when u = v. Consider the first term in the parenthesis of Eq. (2.2), that is, the identity matrix. For it, both terms arising from Eq. (A.1) contribute equally and together yield the following contribution to HŜ x 1 x 2 : where we have also used t a t a = (N 2 c − 1)/2N c . Next consider the contribution coming from the second term in the parenthesis of Eq. (2.2). After expressing adjoint Wilson lines in terms of fundamental ones according to it becomes straightforward to show that this second term yields the same contribution as the first one, that is, Eq. (A.2). The third term in the Hamiltonian leads to contributions like By using the Fierz identity tr t a A t a B = 1 2 trA trB − 1 2N c tr(AB), (A. 5) we find that the corresponding contribution to HŜ x 1 x 2 , which also is equal to the one coming from the fourth term in Eq. (2.2), is which after simple manipulations using the Fierz identity in Eq. (A.5) leads tō Acting with the second term of the Hamiltonian we obtain and by also making use of Eq. (A.3) we are lead tō The action of the third term of the Hamiltonian gives Finally the action of the fourth term in the Hamiltonian gives yielding the same result as shown in Eq. (A.9). Putting together all the contributions, we see that the terms explicitly suppressed by 1/N 2 c cancel each other, as anticipated. Thus, the action of the Hamiltonian on the pair made with the two quarks at x 1 and x 3 gives Similarly, when acting on the quark and the antiquark at x 1 and x 2 respectively, we obtain

B The sextupole
In order to demonstrate the generality of our procedure, we shall succinctly present here the corresponding analysis for the next non-trivial high-point correlator, which is the sextupole defined in Eq. (2.7). The respective evolution equation has not been yet presented in the literature, thus we have derived it on this occasion. We shall not give here the details of the derivation, since this is merely a lengthy but straightforward repetition of the manipulations shown in the previous Appendix. This is the final result: We notice that all the subleading terms of relative order 1/N 2 c have canceled in the final equation, as expected for a single trace operator. As a cross check, one can see that when choosing, for example, x 5 = x 6 , the above equation reduces to the quadrupole equation (2.10).
The subsequent manipulations follow the general procedure outlined in Sect. 4. First, one separates the 'virtual' terms in the r.h.s. of Eq. (B.1), that is, the terms in which the integration variable z appears only in the dipole kernel, but not in the correlation functions; these are the terms proportional to Ŝ (6) Y and to ŜQ Y ≃ Ŝ Y Q Y . Then in evaluating the integrals over z, one keeps only the dominant logarithmic contributions and use Eq. (3.4) to express the logarithms in term of the dipole derivative. We thus arrive at an ordinary first order inhomogeneous differential equation whose solution is Again, one can check that choosing, for example, x 5 = x 6 , the above solution reduces to the quadrupole one given in Eq. (4.3). Expanding Eq. (B.2) in the weak scattering region, we find where we have also used Eq. (3.2) for Q Y and we have assumed that Eq. (B.3) is already valid at Y 0 . It is straightforward to confirm that the same result is obtained when we expand Eq. (2.7) to order g 2 α 2 .
As a particular example, we shall consider the 'line' configuration, already studied in the case of the quadrupole. This is obtained by putting all the quarks at the same point, x 1 = x 3 = x 5 , and similarly for the anti-quarks, x 2 = x 4 = x 6 , with the two points separated by a distance r. (One can trivially visualize this configuration by looking at Fig. 1.a and adding there x 5 and x 6 at the left and right ends, respectively.) Then by also using Eq. (3.8) we see that the y-integration in Eq. (B.2) can be exactly performed and gives where clearly all quantities are evaluated at r. The above simplifies to if we assume that the above and Eq. (3.10) are already valid at Y 0 or if we assume that the scattering at Y 0 for the given r is weak. The sextupole in Eq. (B.5) exhibits geometric scaling in the same regime where the dipole does so.