{'fraction_non_alphanumeric': 0.03716233359833556, 'fraction_numerical': 0.010216681608281602, 'mean_word_length': 4.626725040449224, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New observables to test the Color Glass Condensate beyond the large-N_c limit

The JIMWLK framework offers a powerful tool to calculate the energy dependence of QCD observables at high energies. Despite a growing number of observables considered for phenomenological analysis, few features of JIMWLK evolution beyond its evolution speed are yet well constrained by experiment. We argue that meson production cross-sections have the potential to provide qualitatively new insights and allow to address issues both beyond the large-N_c limit and at higher twist. These cross-sections generically contain four point functions whose evolution is shown to follow from the JIMWLK framework. The Gaussian truncation is used to provide an efficient and practical means of calculating the evolution of four point correlators beyond the large-N_c limit.


Introduction
Hard processes in hadronic collisions, which resolve the partonic structure of hadrons, are well described by the leading-twist approximation of QCD. In this weak-coupling regime, partons in the hadronic wave function scatter independently, this is the essence of collinear factorization. However, since the parton densities grow with decreasing energy fraction x, the hadronic wave function also features a low-x part where the parton density has become large and partons scatter coherently, invalidating the leading-twist approximation. This weak-coupling regime, where non-linearities are important, is called saturation, and it can be probed at high-energies since increasing the energy of a collision allows to probe lower-energy partons.
The Color Glass Condensate (CGC) has been proposed as the effective theory to describe this small-x part of the hadronic wave function in QCD . Rather than using a standard Fockstate decomposition, it is more efficient to use collective degrees of freedom, more adapted to describe the behavior of the small-x partons, which are mainly gluons. The CGC approach uses classical color fields: the long-lived large-x partons are represented by a strong color source ρ ∼ 1/g S which is static during the lifetime of the short-lived small-x gluons, whose dynamics are described by a color field A ∼ 1/g S . The color source distribution W Y [ρ] depends on the rapidity separation Y = ln(1/x) between the source and the field, and contains information on parton densities and parton correlations. These color field configurations are naturally probed in high energy collisions with a rapidity separation of Y ∼ ln(s) between projectile and target.
Using a "dilute" projectile to probe the dense target, it is possible to expand the projectile into its leading Fock state components which then interact eikonally with the target field, due to the high Equivalently, the JIMWLK equation can be cast into a coupled hierarchy of equations for n−point Wilson-line correlators, the Balitsky hierarchies [23][24][25]. If the lowest order Fock-state of the projectile is a qq pair, like in DIS at HERA, the corresponding total cross-section is driven (at zeroth order) by a two-point function or dipole correlator Ŝ qq xy Y := tr(U x U † y )/N c Y , with the average taken with the distribution Z Y [U ].
Most of the phenomenology uses a mean-field approximation which significantly simplifies the highenergy QCD evolution: it reduces the hierarchy to a non-linear equation for the two-point function, the Balitsky-Kovchegov (BK) equation [21][22][23][24][25]. This equation reduces to the BFKL equation when the amplitude is small, and contains saturation effects as the amplitude reaches unity. Although this equation is not exact, it has the advantage of being a closed equation for the dipole scattering amplitude. In addition, closer inspection, both numerically and semi-analytically [36,37], reveals that the difference between BK and JIMWLK solutions for Ŝ qq xy Y is small -much smaller than the order 1/N 2 c difference expected from a superficial analytical argument. For the energy dependence of the total cross-section, Next-to-Leading-Order (NLO) corrections are much more important: scaleinvariance breaking running-coupling contributions [30,38,39] change exact scaling behavior with an energy dependent saturation scale Q s (Y ) at asymptotically high energies into pseudo-scaling, which is reached much earlier. Together with the remaining conformal corrections at NLO [32] evolution rates are drastically reduced. This is crucially important for successful phenomenology [40][41][42][43][44][45].
There are other observables in DIS for which it has been shown that the Y -evolution is given by the same JIMWLK equation, through more complicated correlators [46][47][48][49]. For instance, diffractive structure functions can be computed in the CGC framework with the same level of success as inclusive ones [43,44,50,51]. Diffractive gluon production has also been investigated [52][53][54] as well as semi-inclusive DIS [55]. Exclusive vector-meson production was also considered to study how a finite momentum transfer affects the way the Y -dependence predicted by evolution equations is mapped onto the actual energy dependence of cross-sections [56][57][58].
Despite this recent extension of scope beyond total cross-sections, phenomenological treatments constrain the initial conditions to JIMWLK evolution only in a relatively mild manner -evolution "speed" is most tightly constrained, precise details of correlator shapes such as UV anomalous dimensions only meet weak experimental tests. This is probably best illustrated by the fact that the HERA-fits of [43,44] use correlators in the pseudo-scaling region of evolution while the authors of [42] base their analysis on a pre-asymptotic initial condition. The calculational tools and crosssections suggested here are intended as a starting point to begin exploring such issues.
It is therefore desirable to broaden the scope of observables and search for quantities that directly probe subtler properties of correlators and at the same time remain calculable in a practical sense, preferably in a suitable truncation of full JIMWLK to sidestep the numerical obstacles posed by a full JIMWLK simulation.
The selection of suitable observables is not a trivial task. To focus on yet unconstrained features of JIMWLK evolution and its initial condition, one must by necessity resort to more and more differential observables. This must be done in a way that avoids infrared problems that would invalidate a perturbative treatment via an evolution equation. At the same time one would like to limit the amount of complexity, both to keep cross-sections large enough and to keep numerical effort down: As will become obvious from the example considered below, the price to be paid for more detailed information is the appearance of more complicated n-point functions in the analytical expressions for the cross-section, whose correspondingly involved configuration space structures result in costly, since numerically delicate, integrals.
The discussion below will mainly focus on meson production cross-sections. They serve to illustrate how to calculate the energy dependence of both momentum-transfer-(t-) dependent and t-integrated cross-sections with specific restrictions on the final state. Properly chosen, such restrictions can serve as a filter for information that is unaccessible in the total cross-section. The treatment provides a model for other observables.
A central part of this paper is played by a set of diffractive observables in DIS -based on the target dissociative part of vector-meson production or deeply virtual Compton scattering (DVCS)-that have the potential to expose unexplored features of JIMWLK evolution without excess complexity. As will be explained in the text, these observables directly probe the correlator difference [59] Ŝ qq a slight generalization of the "correlator factorization violations" extensively discussed in [37]. Such correlators strictly vanish in the BK-mean-field approximation (due to 1/N c suppression of their leading contributions), and are also twist suppressed compared to the total cross-section (they require the exchange of at least four gluons between projectile and target in the t-channel). The corresponding cross-sections are necessarily small, but should probe features of JIMWLK evolution that remain invisible in the global observables considered to date. Our main achievement is the calculation of the correlator difference (1) (in the Gaussian truncation of JIMWLK-evolution introduced next), which directly enters the formulation of the diffractive dissociation cross-section (57). Indeed, thanks to the hard scale provided by the photon virtuality, we are able to calculate in QCD (in the high-energy limit), the diffractive dissociation introduced by Good and Walker in the context of soft diffraction [60].
Since (1) vanishes identically in the BK-approximation, it becomes imperative to find an approximation capable of capturing the relevant contributions. By adopting the Gaussian truncation (GT) [27,37,61] of JIMWLK-evolution in place of the BK-approximation, it turns out to be possible to efficiently calculate the energy dependence of both two-and four-point correlators in these observables. In the large-N c limit, the Gaussian truncation reduces to the BK truncation, but in general it allows to calculate correlator-factorization violations. Instead of using Ŝ qq xy Y directly as its degree of freedom, the GT uses a two-gluon t-channel exchange correlator G Y,xy that enters the qq correlator as Ŝ qq xy Y = e −C f G Y,xy . Unlike the BK-approximation it can naturally be used to parametrize more complicated n-point functions in a consistent manner. This is used below to derive an expression for the four Wilson-line correlator Ŝ qq y x Ŝ qq xy Y (including its Y -dependence) in the Gaussian truncation.
The dynamical content driving the energy dependence in the GT approximation is in fact identical to that of the BK approximation [37]. As with the BK approximation, its evolution equation emerges directly from the JIMWLK equation and can therefore be extended to (in principle) arbitrary loop order. It is, however, possible to systematically extend the treatment by including genuine multi-t-channel gluon-correlations into the formalism, building on the non-abelian exponentiation theorem [62][63][64].
The derivation of evolution equations for the various cross-sections considered below remains at a one loop level. This is done in part to keep the arguments compact. A full treatment at NLO would first require a translation of the results of Balitsky and Chirilli into JIMWLK-language, a task left for a separate publication.
The plan of the paper is as follows. Section 2, outlines the proof at leading order of how JIMWLKevolution determines the small-x dependence of certain diffractive cross-sections in DIS, including vector-meson production. The calculation clarifies how, four-point Wilson line correlators enter the more differential cross-section in addition to the two-point functions already present in the total cross-section. The generalization to NLO is discussed qualitatively. Section 3 serves to recall the Gaussian truncation approximation of JIMWLK. It is contrasted with another extension beyond the BK-approximation to JIMWLK evolution as advocated in [65,66]. A re-derivation of the 2point function Ŝ qq xy Y in the Gaussian truncation provides the starting point to address the case of non-trivial higher-n-point correlators -an expression for Ŝ qq y x Ŝ qq xy Y is explicitly worked out in Section 4. Section 5 finally discusses how, in diffractive vector-meson production or deeply virtual Compton scattering (DVCS), the target-dissociating part of the cross-section is related to the correlator difference Ŝ qq (1). Section 6 attempts to put the results into perspective with a review of the tools developed thus far and a short discussion of systematic improvements left for future study.
2 Evolution of DIS cross-sections with restrictions on the projectile final state As already indicated in the introduction, all cross-sections in the high energy limit of the CGC framework will be expressed in terms of Wilson-line operators that reflect the Fock-state content of the projectile. The non-perturbative information about the target wave function probed at a given Y is encoded in averages . . . Y taken with that target wave function at Y (denoted graphically by in the following). In the JIMWLK formalism, this is phrased in terms of an average over the where the dots ". . ." stands for some generic U −content that depends on the observable being considered. What can be obtained by weak-coupling methods within the JIMWLK formalism is the evolution of . . . Y with the factorization rapidity Y . The derivations of JIMWLK evolution have mostly focused on the total cross-section, the literature on evolution of diffractive observables in a JIMWLK context is very terse [47,49]. For this reason a detailed explanation of how to extend JIMWLK evolution to observables beyond the total cross-section is in order. Below, this will be done starting from an amplitude or wave function picture for the total cross-section, in analogy to [61]. In this framework it is easy to illustrate the effects of restrictions on the final state (as constituted by the rapidity gap events observed at HERA or the meson production cross-sections forming the core of this exploration) on the structure of the evolution equation.
Meson production cross-sections will be discussed below mainly in the case of DIS, scattering a virtual photon γ * on a hadronic target with atomic number A. Fig. 1 provides the diagrammatic notation used below to compactly summarize the contributions.

Wilson line correlators in observables within and beyond the total cross-section
The small x approximation to the DIS total cross-section at zeroth order in ln(1/x) involves the eikonal scattering of the qq-Fock component of the virtual photon on the target. The cross-section emerges from the absolute value squared of the difference between the out-state (in which the qq pair interacts with the target and picks up non-Abelian eikonal factors at fixed transverse positions) and the in-state (where this interaction is absent). Diagrammatically, one considers The unrestricted transverse momentum integrals will identify the transverse coordinates left and right of the cut, so that the the U -content of these diagrams is partially simplified: Note in particular that the U -factors left and right of the cut in the out-out overlap cancel against each other so that both contain no interaction with the target. The interaction with the target is fully encoded in the U -content of the two remaining diagrams which takes the form of two dipole operatorŝ In the expression for the total cross-section, these operators are averaged over the target wave function, an operation that involves both perturbative and non-perturbative information, that proves the most difficult part of this calculation and induces the energy (or Y -) dependence of the cross-section. The tool to extract this energy dependence is the JIMWLK equation.
The total cross-section takes the form of a convolution between a wave function overlap -involving the qq component of the virtual photon from right and left of the cut-with the Wilson line correlators that summarize the interaction with the target. The wave function overlap is the same in all four terms of (3) and consist of both transverse and longitudinal contributions.
In the present context the target averages of (5) satisfy Ŝ qq xy (Y ) = Ŝ qq yx (Y ) * , for all applications considered below they will, in fact, be real. With this assumption, the four terms in (3b) assemble in such a way that the target interaction of the total cross-section is fully summarized by the dipole amplitude and one has re-derived the general formula for the total cross-section at small x, where r = x − y and b = (x + y)/2 denote dipole size and impact parameter. At fixed Q 2 , |Ψ(α, r 2 , Q 2 )| 2 encodes the probability to find a qq pair of size |r| and longitudinal momentum fraction α inside the virtual photon. The impact-parameter-integrated dipole amplitude has the interpretation of a qq-dipole cross-section on the target.
JIMWLK-evolution has been originally derived for the total cross-section via the optical theorem [18,23] and applies directly to the evolution of the total cross-section in Eq. (3) or (7).
The derivation of the expressions for the zeroth order cross-section used here avoids using the optical theorem. This allows for easier generalization to other, more differential observables. Such observables can also be addressed using JIMWLK evolution, their treatment, however can only be found either in very compact form [47], or in a somewhat abstract formalism [49] -the treatment here will address all these quantities in a simple unified formalism. The re-derivation of the evolution equations will be done at the one loop level only, with the main goal of identifying the underlying structures and requirements that will have to be met at any loop order.
As already mentioned in the introduction, inclusive vector-meson production is an example for such more differential observables and is of particular phenomenological interest. Addressing inclusive vector-meson production cross-sections requires only a formally trivial modification of the expression in Eq. (3): one only needs to project the qq states in the final state onto the vector-meson state under consideration, without imposing any further restrictions on the final state. This amounts to replacing (3) by Expanding the product, one finds that only the operator content of the out-out overlap is changed markedly from its counterpart in the total cross-section. Since the qq in the final state are now projected onto the vector-meson wave function (instead of being integrated over all transverse momenta) the coordinates of the Wilson lines right and left of the cut will not become equal. Also the sum over all possible colors is replaced by a projection onto a singlet combination, effectively separating the single color trace appearing in the first term of Eq. (3b) into a product of traces.
Since there is no further restriction of the final state, these two traces will remain in a common average, as induced by the single pair of target states present in Eq. (8). This leaves us with a nontrivial correlator in the out-out overlap, namely For consistency of the JIMWLK-formalism, one would expect that both the two and four point correlators entering (8) evolve under JIMWLK evolution. This is indeed generally assumed in the literature. Since the Wilson line factors originate from both sides of the cut, this is not entirely obvious. Sec. 2.2 takes a closer look at how JIMWLK-evolution emerges in the present S-matrix formalism (without recourse to the optical theorem), first for the total cross-section and then for the non-trivial evolution of the out-out-overlap in Eq. (8). The argument employs "flip-back-identities" sometimes referred to as the Mueller optical theorem to confirm that (9) indeed evolves as four Wilson-line operator under JIMWLK-evolution.
The argument also clarifies that the correlators in the remaining contributions, in-in-, in-out-, and out-in-overlaps and their energy dependence do not change compared to their total cross-section counterpart.
Note that the large-N c limit (the BK-approximation) sidesteps all these issues by factoring the correlator: In calculating the analytical form for the exclusive vector-meson from the diagrams in Eq. (8) one encounters wave function factors of the virtual photons in the initial states as well as those of the vector-meson in the final states in each of the amplitudes left and right of the cut. All four terms in (8) are convoluted with two separate wave function overlap factors, originating from the right and left side of the cut respectively. They depend on the photon polarization, quark masses, phenomenological mass and size parameters of the mesons and dynamical variables, longitudinal momentum fraction α, inter-quark-distance r, and momentum transfer Q 2 and will be denoted Ψ T (α, r, Q 2 ) and Ψ L (α, r, Q 2 ) respectively. Explicit expressions are given in App. A.
With the convention t = −l 2 , the expression for the inclusive vector-meson production cross-section takes the form [67] 4π The cross-section for exclusive vector-meson production, in which one requires the target to remain intact, on the other hand, contains two separate U -field averages on the level of the amplitudes without any need for an additional approximation (such as Eq. (10)). Eq. (8) is replaced by and differs from (8) (and (11)) only in the U -field average of the out-out-overlap which now is factorized into The arguments given in Sec. 2.2 ensure that these indeed evolve independently under JIMWLKevolution as one would expect. Due to this factorization of correlators, also the full analytic expression corresponding to (12) factorizes into independent complex conjugate amplitude-factors: Note that by taking the difference between (8) and (12) one directly measures correlator factorization violations of the form a special case of which [x = x] have been extensively discussed in [37] from a theoretical perspective in the light of full JIMWLK evolution. Measuring such differences would allow a direct experimental test of rather subtle features of JIMWLK evolution. These features are of particular interest in that they evidently lie beyond the BK-approximation, for which (15)vanishes identically.
Inclusive and exclusive vector-meson production are by no means the only observables that can be addressed in this way. For example heavy meson production cross-sections are closely related as far as the Wilson line correlators are concerned, despite the change in projectile required. For the inclusive case one would consider where the state marked "in" may be some hadron with quantum numbers suitable to produce a meson with a heavy quark Q. In the phase space region where the mass m Q can be considered large, the transverse position of the corresponding Wilson lines left and right of the cut will be approximately the same, one encounters the coincidence-limit of (9) at x = x which inherits its energy dependence from the more general case, but takes a simpler form.

Evolution and the Mueller optical theorem
The main subtlety with generalizing JIMWLK evolution for the total cross-section (as initially derived using the optical theorem) to the case of more differential observables is the fact that this derivation builds on a whole set of real-virtual cancellations automatically taken care of by the optical theorem, which for more differential observables are no longer valid. Already for a description of rapidity gap events at HERA, the restriction of the final state leads to a modification of the overall evolution equation of the diffractive cross-section that at first sight appears to be only indirectly related to JIMWLK (or BK) evolution [see for example [46], Eq. (11)]. Nevertheless it can be decomposed into individual contributions akin to the terms of (8) that still follow JIMWLK (or BK) evolution, subject to specific initial conditions [46,47].
To address JIMWLK-corrections to the observables in Sec. 2.1 one needs to include gluon emission from the zeroth order expressions shown there. Fig. 2 shows the diagrammatic shorthand notations employed to that end in an exemplary fashion.
There is a strong set of regularities in the small x corrections that is intimately connected with the properties of the JIMWLK-Hamiltonian. The full pattern is already visible if one considers the leading order (LO) corrections to the in-out overlap in Eq. (3), which consists of the diagrams shown in Fig. 3 (with the photon-and target-wave-function factored -the relationships remain valid even if the qq is in an octet state) In Fig. 3, the first two columns cancel amongst themselves. Technically, this is due to the diagrammatic identities shown in detail in the left column of Eq. (17) Eq. (17b) states that there without interaction with the target, no gluons are emitted into the final state. Eq. (17a), on the other hand, expresses the fact that there is no JIMWLK evolution in the absence of interaction with the target -for constant U -fields which are gauge equivalent to U = 1. This becomes obvious once one recognizes that the remaining third column contains the JIMWLK contributions to the evolution of S qq as originally calculated via the optical theorem in [18,23]: Alternatively, one may notice that the diagrams of each of the first three lines in Fig. 3 separately sum to zero following the prototypical pattern of final state cancellations. This recasts the contributions to JIMWLK evolution in terms of the diagrams of the last line in Fig. 3 and offers an alternative interpretation with gluons emitted into the final state during evolution.
The same argument trivially ensures that all corrections to the in-in-overlap must vanish as well.
For the total cross-section, in fact, only the mixed in-out-and out-in-overlaps receive small-x corrections, also the corrections to the out-out-overlap vanish altogether: Arranging the contributions as in Fig. 3, one first notices that there are no cancellations of columns -interaction with the target creates JIMWLK-contributions in the first and last column, and also allows a gluon to appear in the final state of the middle column. However, the first three lines cancel in the sense of final state cancellations just as for the in-out-overlap. In addition, all U fields in the last line cancel (c.f Eq. (4)), mapping the contribution onto a contribution equivalent to the JIMWLK-diagrams with constant U -fields, setting the complete set of corrections to the out-out-overlap to zero.
One recovers the established result for JIMWLK evolution of the total γ * A-cross-section: If one considers the corrections to the vector-meson production cross-sections of Eq. (8) and (12), the argument for the mixed in-out-and out-in-overlaps still follow the pattern of Fig. 3 since they are driven by cancellations based on (17) that occur within one non-interacting amplitude factor. What changes are the corrections to the out-out-overlaps to be discussed next. The diagrammatic identification is straightforward and holds for specific pairs of diagrams in the sums over insertions displayed on both the left and right hand side of the equations in (21). Eq. (21) lists only one loop correlators, but the situation is prototypical and generalizes to arbitrary orders.
One may dub relations of the type shown in Eq. (21) "flip-back-identities" as they relate correlators with (anti-) quark Wilson-lines on both sides of the final state cut (on the right), to equivalent ones that can be thought of as having those lines only on one side of the cut. Note, however, that in doing so one reinterprets a quark Wilson line U † x in the complex conjugate amplitude as an anti-quark line in an amplitude and vice versa.
Also the term Mueller optical theorem has been used to refer to this type of relationship. Its diagrammatic content is by no means trivial: It is remarkable that diagrams like those in the middle column of the right hand side, with an additional gluon in the final state, map onto amplitudes in which no gluon is emitted.
The contributions sum into which, as outlined in Eqns. (20) and (21)  into the final state is prohibited, evolution of the out-out-overlap factorizes into separate evolution for the amplitude and its complex conjugate, as anticipated already upon examining the zeroth order correlator. The corresponding evolution equation is given by Eq. (23) states that the factorized correlators entering (12) and (14) indeed evolve independently with JIMWLK. The proof given in Fig. 4 is strictly one loop and needs to be reexamined already at two loop order, where two gluons can enter the final state and potentially form a color singlet.
To perform this analysis efficiently, it becomes essential to translate the NLO results of Balitsky and Chirilli [32][33][34][35] into JIMWLK language -a task well beyond the scope of questions addressed here.
relate the contributions. Note in particular that the gluon Wilson-lines in the top middle diagram cancel between the amplitude and its complex conjugate since the gluon's final state phase space remains unrestricted.
3 The Gaussian truncation of JIMWLK evolution

Evolution of correlators
Solving the JIMWLK equations to obtain solutions for the vector-meson production observables introduced above is at present curtailed by severe practical limitations -truncations of JIMWLKevolution, however, can be successfully confronted with data. Aside from the BK-approximation which imposes a large-N c limit, we have also access to the Gaussian truncation (GT) which, although sharing the dynamical content of the BK-equation, offers a more refined prescription to map this information onto correlators. This includes correlators that are not accessible in the BKapproximation itself and leads to a closer match with JIMWLK evolution, see [37] for a detailed discussion. The exposition below first recapitulates the main steps in deriving qq and qqg-correlators already given there and in [27] and extends these results to the four Wilson line correlators needed to address inclusive vector-meson production.
The starting point is a parametrization of JIMWLK-averages in terms of the functional prescription where the dots ". . ." stands for some generic U -field correlator such as tr(U x U † y ), U ab z tr(t a U x t b U † y ) or the correlators in Eq. (38). In the light-cone gauge A − = 0, the Wilson lines depend only on the A + component of the target gluon field, and this gauge choice allows to write Eq. (24) in terms of δ/δA + explicitely. We have done so for simplicity, but it is not necessary. In App. B, we show that this parametrization of the target averages . . . Y is equivalent to parametrizing W Y [ρ] as a non-local Gaussian, as proposed in [68,69].
This recasts JIMWLK-evolution in terms of a single Y -dependent two point function G, which represents two gluon exchanges between projectile and target.
Its Y -dependence can be determined from the Balitsky evolution equation for a qq-or more generally a generic RR-correlator (with R denoting a generic representation), y ) Y are known in terms of G, see [27,37] and the discussion below. 2 The result is still a truncation of the Balitsky hierarchies (and thus JIMWLK-evolution) in that the Balitsky equation for a general three point function (for example the one already present on the r.h.s. of (25)) would impose additional conflicting constraints on G -this can be only overcome by introducing additional degrees of freedom in form of higher n-point functions G Y ;x1...xn into the functional (24) to step beyond the Gaussian truncation.
Once such higher end point functions are crucial to the physics content of an observable, such generalizations might become phenomenologically indispensable. Still, the Gaussian truncation will yield an important "baseline" contribution to which the higher n-point functions G Y ;x1...xn provide corrections and must, therefore, be understood first.
Staying within the Gaussian truncation, the U -field correlators can generally be expressed in terms of or its Y -derivative G Y,xy .
The equation for G has already been derived in [27] (and in a somewhat different form earlier in [61] 3 ), starting from the qq-dipole evolution equation (25) with R chosen to be the fundamental representation. In [37] it was shown that (24) treats all dipole equations consistently: Applying (24) to the generic dipole evolution equation (25) yields one and the same equation for G, irrespective of the representation R.
The dynamical information driving the evolution in the Gaussian truncation is in fact equivalent to that of the BK equation [37]. As a consequence the kernel in Eq. (27) is in fact the BK-kernel Like the solutions to the BK equation, the solutions to Eq. (27) generically approch the black disk limit with G Y,xy → ∞ at |x − y| → ∞.
It is possible to extend the result to NLO either by using the full results of [32][33][34][35] supplemented with a resummation of running coupling corrections, or, as briefly outlined and used phenomenologically in [43,44], by combining the resummed running coupling results with "energy conservation corrections" [70] to summarize the conformal NLO contributions in a numerically advantageous manner.
The Gaussian-truncation improves upon the BK-truncation in the way it maps this information onto correlators: it allows to calculate consistent expressions for any Wilson line correlator, even those subleading in a 1/N c -expansion. The correlators obtained from (24) automatically respect group theory constraints that are inherited from field level relationships in various coincidence limits, such as to use the correlators entering (25) as an example. [d R stands for the dimension of the representation (d f = N c for the fundamental representation, d A = N 2 c − 1 for adjoint, etc.) andtr denotes the trace in the adjoint representation.] Such coincidence limits are built into full JIMWLK-evolution, but are true in the BK-truncation only in the leading N c sense. This is the main reason for the closer match between full JIMWLKevolution and GT compared to BK observed in [37]. For instance, with R chosen to be the fundamental representation, the coincidence limit (29a) used together with the Fierz identity yields the relatioñ This relates fundamental and adjoint dipole operators at the field level, prior to target averaging. Performing this averaging in the BK approximation, one finds that Eq. (31) is only true in the leading 1/N c sense. By contrast, all group theory contstraints remain exact in the Gaussian truncation, and in this respect evolution in the GT approximation is closer to full JIMWLK evolution. In the black-disk regime, one has t r(Ũ yŨ † x ) Y = 0, and therefore Eq. (31) imposes that tr(U x U † y ) 2 Y = 1 in the fully saturated regime, meaning |Ŝ qq xy | 2 Y = 1/N 2 c . The result remains finite as required by group theory -a subtle feature that is lost in the BK-approximation.
The Gaussian truncation is not the only approximation proposed to go beyond BK. In the large-N c limit, only dipole degrees of freedom are left in the JIMWLK evolution, meaning that Wilson lines can enter only through the operatorŜ qq xy . In this context, a family a solutions was found in Ref. [65,66]. It is parametrized by a real parameter c, with values between zero and one. The treatment goes beyond the BK approximation as it does not assume correlator factorization and solves the whole hierarchy of equations for n-Ŝ qq xy operators emerging from JIMWLK evolution in this approximation. The solutions of Ref. [65,66] [65,66] for details). As a consequence it leads to a "gray disk limit" in the sense that Ŝ qq xy c Y → 1 − c at |x − y| → ∞ instead of zero, the "black disk limit." By contrast, the GT approximation does not assume the large-N c limit, but solves only the first equation of the hierarchy. The correlations obtained in the GT approximation are more complex and always satisfy group theory constraints exactly (c.f. (29)) instead of up to terms subleading in 1/N c . In addition it always leads to Ŝ qq xy Y → 0 at |x − y| → ∞ i.e. the black disk limit. For the initial conditions considered in [36,37] agreement with JIMWLK-simulations is excellent.

Efficient construction of correlators
To address the four Wilson line correlator encountered in inclusive vector-meson production, one has to step beyond the correlators already obtained in [37] -to prepare for that it helps to review how this was achieved with a number of simple examples.
While direct calculation of specific correlators using the averaging procedure (24) is in many cases straightforward, the calculation can be often simplified significantly by using differential equations. This relies on the observation that which can be used to advantage if the right-hand side turns out to be in some form proportional to . . . (Y ) itself. Note that the derivatives δ/δA a+ x,Y in (32) -contrary to those in (24)-only touch the largest value of x − when acting on the U 's, which are x − -ordered exponentials of t a A a+ . Our convention is that larger values of x − appear further to the left in a U , and therefore further to the right in a U † .
Two point projectile R-R correlators: Using the notation (26) and a prime to denote a Y derivative, straightforward algebra leads to which is readily solved to obtain R tr( The freedom in the initial condition was used to accommodate the normalization factor d R .
Three point projectile adjoint-R-R correlators: These involve several distinct color structures.
where we have used The nontrivial point here is, that this holds for any representation R. Integrating (35) one finds again with the free initial condition used to set the normalization properly.

Four Wilson lines in the Gaussian truncation: correlators for inclusive vector-meson production
To address inclusive vector-meson production in the Gaussian truncation, one needs to step beyond two-and three-point functions and to derive an expression for the four point correlator Ŝ qq in Eq. (9).
As it turns out, this does not take the form of a simple exponentiation as the examples encountered earlier. The reason is the richer color structure associated with the four Wilson lines: with two quarks and two anti-quarks available, one can form not only one, but two non-equivalent singlets, which will mix under evolution in Y . Where (9) forms a singlet directly from qq pairs, one may equally well have each of them in an octet state which then together form a singlet -transition elements between the two singlets in general are also non-vanishing. One is led to consider a 2 × 2matrix of correlators, of which only the (1, 1)-component enters the meson-production cross-sections of Sec. 2.1 directly: where d A = N 2 c −1 is the dimension of the adjoint representation. The diagrammatic representation in the second line employs standard birdtrack notation to clarify the color structures (see [71] for a textbook introduction). The decomposition of a qq state into singlet and octet is written as and used to form a color basis to span the space of two non-equivalent singlets according to 1 and 1 with normalization factors given by with the 1 4 related to the the normalization of generators tr(t a t b ) = 1 2 δ ab .
This notation facilitates a more immediate structural identification of the meaning of the analytical expressions given in the first line of Eq. (38).
The reason one can not simply focus on A 11 (Y ) to the exclusion of the remaining three components is that, in general, they mix under JIMWLK-evolution even if the off diagonal elements in (38) vanish at some Y 0 . To see this, one notes that already the linear contributions to JIMWLK-evolution (the terms in Eq. (20) with non-interacting gluons), such as couple the four components in A, since M is a 2×2-matrix with non-vanishing off diagonal elements. Its entries depend on all four coordinates x, y, x , and y present in A via the BK-kernel. 4 These linear terms are the qq-analogues to the second term on the right hand side of (25). The non-linear term (the analogue to the first term on the right hand side of (25)) involves the entries of a matrix of entirely new correlators in a space spanned by all singlets in a (qq) 2 g tensor product: there are altogether 6 independent such singlets, and therefore 6 × 6 (qq) 2 g-correlators contributing to the nonlinear term of full JIMWLK evolution of (38).
Even if one restricts oneself to the Gaussian truncation -and therefore ignores the full complications of the nonlinear structure of the JIMWLK-equation-one can not escape the channel mixing already seen on the linear level in Eq. (41).
The correlator matrix A is hermitian even in the most general case of four independent coordinates x, y, x y , and thus can in principle be diagonalized -the diagonalizing transformation, however, will in general be Y -dependent. This means that the procedure outlined in Sec. 3.2 to derive an expression for the four Wilson line correlator in the Gaussian truncation leads to a solution in terms of a Y -ordered exponential. Applying (32) to (38) (which chooses (40) as a basis) leads to a matrix equation of the form where matrix entries of M can be obtained by expanding A(Y ) to first order in G . One finds with (suppressing the Y -dependence on the G for compactness) [C d = (N 2 c − 4)/N c is defined via the totally symmetric tensors through d abc d abc = C d δ cc -it vanishes at N c = 2, where d abc → 0 itself.] The equation for A(Y ) can be integrated to yield (with P Y denoting path ordering in rapidity) which relates A(Y 0 ) to A(Y ) solely in terms of G , the degrees of freedom in the Gaussian truncation.
As it stands, the initial condition A(Y 0 ) can not be determined from that of a dipole correlatorit contains new information that can not be extracted from dipoles alone. The freedom to choose A(Y 0 ) and its coordinate dependence is not unlimited, however, it is rather strongly constrained by coincidence limits.
Direct inspection of Eq. (38) reveals that the off-diagonal elements vanish when either y → x or y → x (simply since the generators are traceless) in a Y -independent manner. In these limits, the diagonal elements reduce to the two and three point correlators encountered earlier, whose initial conditions contain no freedom whatsoever, they are fully determined in terms of G(Y 0 ) alone. This then also determines A(Y 0 ) in this limit without further freedom. The diagonal entries of A(Y 0 ) follow from Eqns. (34) and (37), together with the Fierz identity (30).
From the above, it is clear that if one swaps coordinates x and y before forming singlets and octets, one obtains a matrix of correlators, related to the above by using The only set of pairwise coincidence limits not yet discussed is the situation in which either both U -factors or both U † -factors are taken at the same point. Y -independent diagonalizability of A(Y ) for this pair of coincidence limits is exposed in a basis where one first decomposes both quarkand anti-quark lines into symmetric and antisymmetric contributions that then map into singlets. Details, including the limiting correlators on the diagonal are given in App. C.
Evidently, this full set of coincidence limit constraints on A(Y 0 ) does not leave much quantitative freedom. The nontrivial freedom on the initial condition is restricted to the configuration space regions away from the coincidence limits, it resides in the difference between path-ordered-exponentiation as in Eq. (45) and "rigid" exponentiation as defined by 5 and A rigid (Y ), contrary to A(Y ), can be expressed exclusively in terms of G(Y ) instead of involving Y -ordered expression in terms of G (Y ) and differs from A(Y ) only away from the coincidence limits. 6 The large-N c limit is of no help in understanding the freedom left. It simplifies the situation so drastically that all expressions for A(Y ) in both the bases (40) and (46) become completely diagonal and ordering plays no role at all (c.f. both (38) and (44)). 7 The simplifications obtained, however, not only affect configuration space away from the pairwise coincidence limits, also the result for the limits looses subleading 1/N c contributions correctly encoded into the rigid exponentiation expression (47) at finite N c : The diagonal entries in basis (40) appear as the product of two large-N c dipoles in keeping with Eq. (10) with the two diagonal entries swapped in basis (46).
Perhaps more useful is the observation that A(Y ) reduces to A rigid (Y ) identically if one approximates the Y -dependence of G by an ansatz that factorizes coordinate-dependence from Ydependence according to with arbitrary functions f and g. The widely used Golec-Biernat-Wüsthoff parametrization as well as the McLerran-Venugopalan model are of that type. They replace G xy (Y ) by respectively. The factorization property (51) implies that the Y -dependence factors out of a Y , b Y , and c Y ; eigenvectors of (38) become Y -independent and Y -ordering is no longer an issue. (47) is indeed identical to the result obtained in [59] for the Ŝ qq y x Ŝ qq xy Y correlator, where the McLerran-Venugopalan model was employed. The 4-point correlator calculated in [72] can also be recovered using a different color basis.
While this type of factorized ansatz is not compatible with the GT-evolution equation (27), it has met with much phenomenological success as an ansatz for the initial condition to evolution for dipoles, sometimes supplemented with a short evolution interval to modify the initial condition away from the simple models before one starts comparison with data.
Exploiting evolution to at least partially erase features of the initial condition is useful in situations in which one expects most gluons in the system to be perturbatively produced. In the process details of the initial condition are erased and supplemented by universal properties imposed by the nonlinearities of the evolution equation as the solutions of that equation approach the asymptotic scaling regime. For this situation, an analogous strategy can also applied: Aiming at a phenomenological comparison with data at Y > Y 0 one uses to parametrize the evolution of the four point function.
For sufficiently large ∆Y , this procedure should erase much of the arbitrariness introduced into A(Y 0 ) (and A(Y ) at Y > Y 0 ) by using any reference to A rigid at all. In this sense, Eq. (54) should provide a reasonable ansatz for the 1/N c -corrections characteristic of the Gaussian truncation also away from the coincidence limits (where they are correctly implemented by construction). The quality of the result should improve with Y .
The strategy to calculate a correlator in the Gaussian truncation, as laid out in the above, is in fact completely general. For any combination of quarks, anti-quarks and gluons, one first needs to determine a full set of non-equivalent singlet projections (if non-equivalent, their basis elements will automatically be orthogonal to each other). From this one forms the analogue of A(Y ) in Eq. (38) as a correlator matrix. The next step is to find M(Y ) by expanding the correlators to lowest order in G . With these ingredients one can then construct (numerically) A(Y ) via (54), using the solutions to (27) as input.
The dimensionality of the problem generically grows with the number of quarks and gluons: As already mentioned above, the (qq) 2 g-correlator gives rise to 6-nonequivalent singlets and their transition elements. Adding more Wilson lines will generically increase the number of singlets even more. The main exception to this rule are the baryon-or anti-baryon-correlators which allow only one singlet to be formed via the totally antisymmetric tensor in N c -dimensions [52,73]. As a consequence the baryon correlators can be evaluated by direct exponentiation to give 5 Measuring correlator-factorization violations with target diffractive dissociation As already indicated in Sec. 2.1, taking the difference of inclusive and exclusive vector-meson production cross-sections, Eqns. (8) and (12), one directly probes correlator factorization violations (15), since all but the out-out-overlap contributions cancel and one is left with Using t = −l 2 , one finds the following expressions for target-dissociating parts of the transverse and longitudinal cross-sections: Let us first discuss a subtlety of this formula. The fully differential vector-meson production crosssections also depend on the fraction z of the meson longitudinal momentum with respect to the photon. This dependence has been integrated out in (57), in the regime where the eikonal approximation we are using to describe the γ * A scattering is valid, meaning with 1 − < z < 1 (the result does not depend on but should be small for the result to hold). Over this kinematical z range, even inclusive events feature a rapidity gap in the final state, between the vector-meson and the system X, of invariant mass M X , coming from the dissociation of the target A: indeed from kinematics a small implies M X M A . The difference with exclusive events is, that in these the target A really escapes the collision intact. Obtaining the z dependence of the fully differential vector-meson cross-sections for arbitrary z, as well as a non-trivial (or maximum-M X ) dependence of the z integrated cross-section (57), implies to go beyond the eikonal approximation.
Using the change of variables r = x − y, b = αx + (1 − α)y and r = x − y , b = α x + (1 − α )y allows to write the t−integrated cross-sections as where in the target averages . . .
These quantities directly probe corrections beyond the leading 1/N c approximation of the BK truncation of JIMWLK evolution. The Gaussian truncation, however, offers consistent finite expressions for such correlator differences, that can be constructed from the solution to the evolution equation for the Gaussian truncation (27). The unfactorized term follows numerically from the (1, 1)-component of (54), the factors in the factorized term can be read off from Eq. (34) to yield exp −C f G y x + G xy .
After t integration, also heavy meson production cross-sections are of particular interest: In the phase space region where the four Wilson line correlator is dominated by the conicidence limit of the coordinates of the heavy quark in amplitude and complex conjugate amplitude (x = x), the Gaussian truncation provides an analytical expression in terms of G: This set of observables provides a very attractive tool to test nontrivial features of JIMWLKevolution in direct comparison with data. In the case of vector-meson production, the cross-section difference we are interested in could be measured at a future electron-ion collider [74]. In the case of heavy meson production, one may not be able to measure the cross-sections directly. Perhaps the inclusive cross-section can be extracted in heavy-ion collisions, from studying the propagation of mesons through cold-nuclear matter. In this case the Gaussian truncation provides 1/N 2 c corrections to the BK result.
To compare with data, however, one must find a way to treat the b-dependence in Eqns. (57) and (58) or any other meson production cross-section. JIMWLK evolution already fails to correctly describe the b-dependence of the total cross-section formula (7), since gluon emission in JIMWLK is perturbative and has 1 |z| 2 power law tails that lead to exponential growth of an initially finite size target. Comparison with data already in this case is done by modeling the b-dependence using a factorized b-profile according to N r 2 ;Y is then taken from JIMWLK (or one of its truncations) while the profile choice leaves only a normalization constant d 2 b T (b 2 ) ≈ 2πR 2 target . The situation for the factorization violation measurements as well as the separate inclusive and exclusive meson production cross-sections is more sensitive to profiles: They explicitly enter the integrands of expressions such as (14) and will affect and distort how Y -dependence as calculated from evolution equations is mapped onto the actual behavior of cross-sections. Still, there are models on the market that have been successfully applied to data (see for instance [75]) which one may use as starting point for measuring correlator factorization violations, as in [67]. A thorough evaluation of sensitivity on model details and model independent statements will be left to a future phenomenological study.

Conclusions
We have shown that the high energy limit of inclusive vector-meson production cross-sections in DIS is determined by two-and four-point Wilson-line correlators [Eqns. (8), and (11)], whose energy dependence in that limit is entirely determined by JIMWLK evolution (22)). The situation for exclusive meson production is different only in the contribution with four Wilson-lines, which are factorized into a product of two point functions (Eqns. (12), and (13)) with accordingly factorized (independent) JIMWLK evolution (23).
Unlike particle multiplicities or forward particle production cross-sections, only a minimal set of phenomenological assumptions is needed to compare with data -modelling is needed only with respect to impact parameter profiles and meson wave functions. No additional assumptions about the applicability of k t -factorization or the twist expansion are needed as long as the energies are high enough so that the no-recoil approximation remains valid to justify a description solely in terms of Wilson line degrees of freedom.
Already the inclusive and exclusive vector-meson production cross-sections probe new information about JIMWLK evolution not visible in the total cross-section. This is induced by their sensitivity to the energy dependence of four point correlators (which cancels in simpler observables) both at finite momentum transfer and for the t-integrated cross-section (the latter being both easier to calculate and to measure).
The difference of the two cross-sections -the target dissociative contributions to vector-meson production-is directly driven by the correlator difference Ŝ qq y x Ŝ qq xy Y − Ŝ qq y x Y Ŝ qq xy Y . This quantity, although expected to be small, is sensitive to entirely unexplored contributions which are not accessible in the large-N c limit and even at low densities start only at higher twist, with four gluons in the t-channel necessary to yield a non-vanishing contribution.
We have explained how to efficiently implement JIMWLK-evolution of all the correlators involved beyond the large-N c limit using the Gaussian truncation (24). Such a step beyond the BK approximation to JIMWLK-evolution is mandatory, since the BK approximation is built on the large-N c limit and assumes Ŝ qq xy Y from the outset. As a consequence, differences between inclusive-and exclusive-vector-meson production become inaccessible -they are beyond the scope of the approximation. The GT approximation, on the other hand, allows to calculate this difference and to obtain a finite result for the target-dissociating part of the cross-section (57).
The GT approximation can be applied to even more general correlators -the method to obtain the four point correlator from a matrix evolution equation outlined in Sec. 4 readily generalizes to arbitrary n-point functions -at the price of increasing numerical cost with each added point. Constraints on the coordinate hyper-volumes help mitigate that cost -t-integrated cross-sections, for example, require less numerical effort already at the four-point level.
While we have been mainly concerned with vector-meson production in DIS, there exist other observables in hadron-hadron collisions known to involve complicated correlators, to which our method can be applied. Heavy quark-antiquark pair production [68] and inclusive dijet production [69,76] which have already been studied in the BK approximation can be reexamined in the Gaussian truncation.
The rich structure of the full JIMWLK equation and the Balitsky hierarchies beyond the simplest truncations contain an immense amount of information that can only be studied by widening the pool of observables considered.

A Photon-vector-meson wave function overlaps
Here we list the overlaps Ψ T,L (α, r, Q 2 ) of virtual photon and vector-mesons introduced in Sec. 2 and used throughout.
The general structure and the complete photon-side contributions are calculated from field theory, the vector-meson side relies on phenomenological models for the functions φ L,T entering below.
Among the models for φ L,T used in the literature, two examples are of particular interest: the boosted Gaussian (BG) wave functions [77,78] φ BG L,T = N L,T exp − used with δ = 1, and the light-cone Gauss (LCG) wave functions [79,80] φ LCG L =N L exp −r 2 /(2R 2 L ) , used with δ = 0. The parameters R and N L,T are constrained by the normalizations of the wave functions, as well as by electronic decay widths. They are given for instance in the Appendix of Ref. [58], along with the effective chargesê f , quark masses m f , and meson masses M V .
B Expressing the Gaussian truncation in terms of a nonlocal Gaussian distribution of ρ's In this Appendix, we show that the GT approximation is equivalent to using the following ansatz for the CGC wave function in terms of target color sources ρ: It is a Gaussian distribution for the color sources ρ c , whose variance µ 2 represents the transverse color charge density squared on the path of a projectile moving along the x − direction. The longitudinal extend of the target, over which µ 2 doesn't vanish, is proportional to e −Y . Rather than labeling µ 2 with this explicit Y -dependence, it is recast in the boundaries of the rapidity integration (equivalent to an x − integration).
In Eq. (2), after expanding the Wilson-line content of ". . ." in powers of the color field A + , (we work in the light-cone gauge A − = 0 where A + is a linear function of ρ), any average with the Gaussian weight W Y [ρ] can be computed using Wick's theorem with µ 2 (Y , x, y) parametrizing the ρ-correlators ρ c (x, Y )ρ d (y,Ȳ ) Y = δ cd δ(Y −Ȳ )µ 2 (Y , x, y) .
To establish the connection with (24), one first notes that it simply restates Wick's theorem as and assumes that the AA correlator is in fact local in Y to pick the leading contribution in the non-abelian exponentiation theorem.
The only step left is the translation of the AA-correlator into ρ-correlators via the Yang-Mills equation −∇ 2 ⊥ A c+ = ρ c . Inverting this relation gives and therefore G Y ,xy = d 2 zd 2 z µ 2 (Y , z, z )G 0 (x − z)G 0 (z − y) .
with normalization factors arising from The off-diagonal elements of Eq. (38) in this last basis (72) vanish when y → x or x → y, since anti-symmetrization of a symmetric object yields zero The surviving correlators that appear on the diagonals in this basis in these limits constitute the last set of coincidence limit conditions on A(Y 0 ). Using the basis (72), one obtains explicit expressions for A in the limit x = y e −(C f + Nc −1 2Nc )(Gx,y+G y,y −G x,y )−C f G x,y 0 0 e −(C f − Nc +1 2Nc )(Gx,y+G y,y )− Nc +1 2Nc G x,y , while the limit y = x leads to These structures are not completely meaningless: For N c = 3 the anti-symmetrized correlator (lower right corner) must agree with the color correlator of a proton (anti-proton) Eq. (55). This is indeed the case, since the prefactors of all three terms in the exponent become equal for N c = 3: The symmetrized channel shows a structure similar to the qqg correlator with (C f + Nc−1 2Nc ) replacing the Nc 2 present there. The two agree at N c = 2.