JIMWLK evolution for multi-particle production in Langevin form

Within the effective theory for the Color Glass Condensate, we study multi-particle production with rapidity correlations in proton-nucleus collisions at high energy. The high-energy evolution responsible for such correlations is governed by a generalization of the JIMWLK equation which describes the simultaneous evolution of the (strong) nuclear color fields in the direct amplitude and the complex conjugate amplitude. This functional equation can be used to derive ordinary evolution equations for the cross-sections for particle production (a generalization of the Balitsky hierarchy). However, the ensuing equations appear to be too complicated to be useful in practice, including in the limit where the number of colors is large. To circumvent this problem, we propose an alternative formulation of the high-energy evolution as a Langevin process, which is better suited for numerical implementations. This process is directly oriented towards the calculation of the cross-sections, so its detailed structure depends upon the nature of the final state. We present the stochastic equations appropriate for two gluon production, and also for three gluon production, with generic rapidity differences.


Introduction
The study of multi-particle production in hadron-hadron collisions at very high energy, as currently pursued at RHIC and the LHC, provides strong evidence in favour of collective phenomena associated with high parton densities. In the case of proton-proton or proton-nucleus collisions, such phenomena are generally attributed to 'initial-state interactions', that is, to the existence of high gluon densities in the wavefunctions of the incoming hadrons, as generated via QCD evolution (i.e. parton branching) with increasing energy, and also -in the case of a large nucleus with atomic number A 1 -via radiation from a large number of valence quarks. This paradigm has recently been challenged by new data for p+Pb collisions at the LHC, which refer to long-range rapidity correlations (the 'ridge effect') and exhibit strong flow components [1][2][3][4], suggestive of 'final-state interactions' -that is, interactions between the partons liberated by the collision and which might form a dense fireball at the intermediate stages of a collision. For nucleus-nucleus collisions, the importance of final-state interactions is by now well established, via phenomena like jet quenching or hydrodynamical flow, as observed at both RHIC and the LHC. But even in that case, it is quite clear that the final multi-particle distribution as measured in the detectors is the result of a complex interplay between initial-state and final-state interactions. In order to distinguish between the various possible mechanisms and arrive at a global and unambiguous understanding of the data, it is essential and urgent to have reliable calculations of multi-particle production, from first principles.
The present study provides a further step in that sense, by proposing a new and hopefully efficient method for computing multi-particle production with rapidity correlations in a 'densedilute' set-up, such as proton-nucleus collisions. Some physics problems that we have in mind are two-hadron production at central-forward rapidities 1 and the 'ridge effect' alluded to above. Our method includes the effects of 'initial-state interactions' in the framework of the Color Glass Condensate (CGC) effective theory [5][6][7][8][9][10][11]. In this framework, correlations in rapidity and transverse momentum (or azimuthal angle) are built in the nuclear wavefunction, via the JIMWLK evolution [12][13][14][15][16][17][18] of the respective gluon distribution, and get transmitted to the produced partons via multiple scattering during the collision. Our approach generalizes previous CGC calculations of particle production which referred to simpler observables -either single inclusive hadron production [19][20][21][22][23], or a pair of particles with similar rapidities (e.g. di-hadron production at 'forward rapidities', i.e. in the fragmentation region of the projectile) [24][25][26][27][28][29][30][31]. It builds up on previous extensions of the CGC formalism towards multi-particle production [25,[32][33][34], which brought important conceptual clarifications but failed to provide a tractable calculational scheme. In fact, the method that we shall develop here is the extension to particle production of the only method that has proven so far to be useful for actually solving the JIMWLK equation: the Langevin approach pioneered in Ref. [35] (see also Refs. [36][37][38] for numerical implementations of this method).
The main difficulty with this kind of problems is the treatment of the high-energy evolution at intermediate rapidities, between those of the produced particles. This cannot be factorized as 'initial-state' (i.e. prior to the collision) evolution of the (proton and nucleus) wavefunctions and greatly complicates the calculation of the cross-sections. To better explain this difficulty, we start with a problem for which the CGC factorization is by now well established: the production of one or several partons with similar rapidities in proton-nucleus [19][20][21][22][23][24][25][26][27][28][29][30][31], or nucleus-nucleus collisions [39][40][41][42]. The CGC factorization scheme, which holds within the leading logarithmic approximation in perturbative QCD at high energy, separates the high energy evolution of the incoming wavefunctions from the cross-section for the partonic subprocess, while taking into account the high-density effects (gluon saturation and multiple scattering) on the nuclear side. For definiteness, consider two gluon production in p+A collisions. So long as the rapidity difference ∆η between the produced gluons is not too large, such that α s ∆η 1, one can ignore the intermediate evolution. It is then convenient to divide the total high-energy evolution between the BFKL evolution [43][44][45] of the dilute projectile (the proton) and the JIMWLK evolution of the dense target (the nucleus), each of them extending over the respective rapidity difference w.r.t. the produced gluons. Moreover, when computing these evolutions, one finds that only the initial-state emissions contribute to the final result: the effects of the final-state evolution -emissions of unresolved gluons which occur after the collision -cancel between 'real' and 'virtual' emissions [46]. This is a consequence of causality: within the approximations at hand, the evolution gluons are significantly faster than the measured ones, so the respective formation times are also much larger, by Lorentz time dilation. Hence, the final-state emissions associated with the high-energy evolution occur long after the particle production has been completed and cannot affect the cross-section for the latter.
We now move to the more interesting case where the rapidity separation between the produced gluons is relatively large, α s ∆η 1, and the intermediate evolution cannot be neglected anymore. The (unresolved) gluons associated with this evolution can in particular be emitted by the fastest among the two measured gluons, that is, the one which is closer in rapidity to the projectile. Such emissions can occur both before and after the scattering between the emitter and the target. Whenever such an emission occurs, it modifies the state (color and kinematics) of the fast measured gluon, hence it has an observable effect on the cross-section -including for the final-state emissions. Accordingly, the effects of the final-state evolution do not cancel anymore. Moreover these effects depend upon the color field of the target: they describe the multiple scattering of the evolution gluons off the strong target field. We see that the complications associated with the final-state emissions are generally twofold: not only they cannot be simply factorized like the initial-state evolution, but they can neither be treated as the standard BFKL evolution of the fast measured gluon from its own rapidity down to that of the softer measured one 2 . Rather, one has to deal with BFKL evolution in the presence of a strong background field. A first method in that sense, which applies in the limit where the number of colors N c is large and exploits the 'color dipole' formulation [47] of the BFKL evolution, has been proposed in [25]. A more general method, which works for any N c and involves a suitable generalization of the CGC formalism, has been introduced in Ref. [32], in a somewhat different context -the study of diffraction. Subsequently, this method has been generalized to particle production in [33,34]. This second method will be here given an equivalent Langevin formulation, which is better suited for numerical calculations.
The basic idea is that one needs to follow the background-field evolution of the wavefunction squared of the fastest measured gluon while keeping trace of the difference between the background field in the direct amplitude (DA) and that in the complex conjugate amplitude (CCA). Indeed, at the intermediate stages of the calculations, this 'background' field does not reduce itself to the classical field of the target (which is of course the same in the DA and the CCA), but also includes the quantum fields representing the softer gluons to be emitted -the 'evolution' gluons at intermediate rapidities and the soft measured gluons. (This doubling of the quantum degrees of freedom is reminiscent of the Keldysh-Schwinger formulation of quantum field theory in real time, which involves a similar doubling of the time axis, with the upper side corresponding to the DA, and the lower side to the CCA.) Hence, the wavefunction squared plays the role of a generating functional for both evolution and particle production. Its evolution can be viewed from either the projectile side, or from the target side, and it is the last perspective which will be more useful for us here: as in the standard CGC formalism, it is the target perspective which allows for a Langevin reformulation of the JIMWLK evolution [35].
To describe our method, let us first remind that, in a Lorentz frame where the target carries most of the high-energy evolution -a 'target infinite momentum frame' -, the partons from the projectile couple to the (strong) color field of the target via Wilson lines, i.e. time-ordered exponentials of the gauge potential, where the 'time' is the light-cone time of the projectile (i.e. x + if the projectile is a right-mover). The time ordering of the fields in the Wilson lines reflects their ordering in rapidity: the inner core near the light-cone (x + = 0) corresponds to the target valence degrees of freedom and is surrounded by layers of fields at larger and larger values of |x + |, associated with quantum gluons which are softer and softer relative to the target (i.e. closer and closer to the rapidity of the projectile). This 'multi-layer' structure is manifest in the Langevin formulation of the JIMWLK evolution, where the Wilson lines are built by successively adding new layers at larger values of |x + | with increasing the rapidity difference w.r.t. to the (valence d.o.f. of the) target. In the original formulation of the CGC effective theory, which is oriented towards the calculation of the elastic S-matrix for dilute-dense scattering, all the modes of the target field are treated as classical fields, irrespective of their rapidity. This formulation can also be used for computing multi-particle production, but only so long as the produced partons 2 The second complication disappears in the limit where the target is itself dilute and one can neglect multiple scattering. In that case, one recovers the standard kT factorization of multi-particle production at high energy, based on the BFKL evolution (see also the discussion in Appendix A).
have similar rapidities. In that case, one can associate by hand two Wilson lines -one in the DA, the other one in the CCA -to each of the produced partons. These two Wilson lines live at different transverse coordinates, but they are both built with the same gauge fieldthe classical field of the target. However this construction becomes insufficient for the general problem of multi-particle production, where the measured partons have arbitrary rapidities. In that case, the measured gluons have to be put on-shell, hence they must be treated as quantum fields, which are different in the DA and respectively the CCA. In the presence of non-linear phenomena like gluon saturation and multiple scattering, this distinction must be carried over to all the 'evolution' (unresolved) gluons at intermediate rapidities, in between the produced ones. This in turn requires a generalized CGC formalism, in which Wilson lines in the DA and the CCA are built independently from each other.
The simultaneous evolution of the target fields in the DA and the CCA is governed by a generalization of the JIMWLK Hamiltonian [32] -essentially, its extension to the Keldysh-Schwinger closed time contour -that we shall refer to as the 'evolution Hamiltonian' H evol . In Refs. [33,34], this has been mostly employed to study the evolution of the (dilute) projectile at large N c , and thus deduce equations which extend the Balitsky-Kovchegov equation [48,49] to the cross-section for two gluon production 3 . This strategy will be illustrated in Appendix A, where we shall use H evol to construct coupled evolution equations (valid for any N c ), which generalize the Balitsky hierarchy [48] to multi-particle production in dense-dilute scattering. The ensuing equations are however extremely complicated, both because of the non-linear effects associated with saturation in the target, and because of the transverse non-locality inherent in the calculation of particle production. Some simplifications occur at large N c , where the hierarchy closes at the level of the 'dipole' and the 'quadrupole' -color traces of two and, respectively, four Wilson lines in the fundamental representation, which in the present context represent generating functionals for projectiles made with a single quark and, respectively, a quark-antiquark dipole. But even those equations appear to be too complicated to be useful in practice. (The simplest closed equation, which describes gluon radiation by a color dipole, has the same degree of complexity as the Balitsky-JIMWLK equation for the evolution of the quadrupole S-matrix [25,33]; see Eq. (A.9) in Appendix A.) This complexity appears to be inextricable and motivates our search for an alternative strategy. Our proposal is to use H evol in the same way as the original JIMWLK Hamiltonian is generally used in the context of the CGC: to describe the evolution of the target. This allows us to set-up a Langevin formalism in which the Wilson lines in the DA and respectively the CCA are built independently from each other, over the whole rapidity interval separating the produced particles. One additional difficulty as compared to the standard JIMWLK evolution is the fact that, as alluded to above, one needs to construct a generating functional -a functional of the Wilson lines at lower values of the rapidity, which can be used to generate an arbitrary number of gluons via functional differentiations. In terms of the Langevin process, this implies that one has to solve a stochastic equation with functional initial conditions (generic Wilson lines, rather than numerically-valued color matrices). This is conceptually well-defined, but not well suited for numerics. The solution that we propose to circumvent this difficulty, is to adapt the Langevin process to the specific cross-section of interest: for a given final state, one can build ordinary (i.e. non-functional) Langevin equations for both the Wilson lines and their functional derivatives which enter the calculation of the cross-section. For instance, in order to compute two gluon production with generic rapidity separation, one needs additional Langevin equations for the first-order functional derivatives of the Wilson lines, for three gluon production, one also needs the equations obeyed by the respective second-order functional derivatives, etc. All such equations are straightforward to write down and a priori accessible to numerical simulations. But the complexity of the numerical implementation should rapidly increase with the number of produced particles: indeed, each additional functional derivative increases the non-locality of the Langevin process in the transverse space (see Sect. 5 for details). Hopefully, this whole scheme will turn out to be tractable, at least, for the calculation of two-gluon production, with consequences for the phenomenology of di-hadron correlations in proton-nucleus collisions at RHIC and the LHC.
As mentioned earlier, the physical original of the large ridge effect seen in p+Pb collisions at the LHC is still unclear. Previous CGC-inspired calculations appear to describe these data quite well [50][51][52] (see also Refs. [53][54][55] for earlier related work), but this agreement cannot be viewed as fully conclusive in view of the many underlying approximations, whose effect is difficult to control. We hope that the method proposed in this paper will open the way towards controlled calculations from first principles, within the accuracy limits of the present formalism.
The plan of the paper is as follows. In Sect. 2 we give a brief review of the CGC formalism and the JIMWLK equation, with emphasis on the Langevin formulation of the latter. In Sect. 3 we consider multi-particle production at similar rapidities, that is, in the absence of any high energy evolution between the parent partons and the produced ones. In this context, we introduce the notion of generating functional for particle production (the 'wavefunction squared of the projectile'), together with a special operator -the 'production Hamiltonian' H prodwhich generates on-shell gluon emissions when acting on the generating functional [33,34]. This Hamiltonian looks very similar to the JIMWLK Hamiltonian, for reasons which should become clear in the next section. Specifically, in Sect. 4 we study the high energy evolution with increasing rapidity difference between the parent partons and the produced gluons, or between the produced gluons themselves. As already mentioned, this evolution is governed by H evol -a generalization of the JIMWLK Hamiltonian which generates gluon emissions in both the DA and the CCA [32]. (The 'production Hamiltonian' H prod is the 'cut' version of H evol , which produces on-shell gluons alone.) The emphasis in this section will be on the evolution of the target (the dense nucleus), up to the rapidity of the 'most forward' produced parton. (The complimentary viewpoint of the projectile evolution will be developed in Appendix A.) In this context, we shall describe the generalization of the CGC formalism to multi-particle production. This involves an 'off-diagonal' version of the CGC weight function, which allows for different field configurations in the DA and, respectively, the CCA, and encodes the high-energy evolution of the generating functional. Sect. 5 presents our main new results: a Langevin formulation for the JIMWLK evolution of the generating functional. In particular, we shall present a version of the formalism which is free of functional aspects and thus suitable for numerical implementations. In Sect. 5, we show the stochastic equations appropriate for two gluon production, while in Appendix B we describe the additional equations which appear when computing three gluon production with generic rapidity differences.

The Color Glass Condensate and the JIMWLK evolution
In the effective theory of the Color Glass Condensate (CGC) valid in the leading logarithmic approximation at high energy, the expectation value of an observableÔ which is local in rapidity and which is associated with the scattering between a dilute projectile ('proton') and a dense target ('nucleus'), is schematically given by the following functional integral Here, U ≡ U (x) is the Wilson line describing the scattering between a parton with transverse coordinate x from the projectile and the strong color field of the target, in the eikonal approximation (see Eq. (2.3) below for an explicit formula).
[DU ] is the functional group-invariant (or de Haar) measure on SU(N c ). W Y [U ] is a functional probability density describing the distribution of the Wilson lines in the target, known as the 'CGC weight-function'. Y is the rapidity difference between the projectile and the target, and it is assumed to be large, that is, α s Y > 1 with α s the QCD coupling. Eq. (2.1) is written in a frame in which most of this rapidity is carried by the target, so one can ignore the high-energy evolution of the projectile. For instance 4 , represents the S-matrix for the scattering of a right-moving gluonic dipole off the strong color field of the left-moving nucleus. In the above N g = N 2 c − 1 with N c the number of colors, and where the T a 's are the color group generators in the adjoint representation, A µ a = δ µ− α a is the color field generated by the nucleus and P stands for path ordering: with increasing x + , matrices are ordered from right to left. The integral over x + formally extends along the real axis, but in practice it is limited to the support of the target field, which is localized near x + = 0 by Lorentz contraction (a 'shockwave'). This support depends upon the rapidity difference Y : with increasing Y , the scattering probes gluon modes in the nuclear wavefunction which carry smaller and smaller fractions x = e −Y of the target longitudinal momentum (k − ), and hence are more and more delocalized in x + . The information about the support of the target field and about its correlations is encoded in the CGC weight-function W Y [U ]. When increasing Y , this evolves according to a functional renormalization group equation known as the JIMWLK equation 5 [12][13][14][15][16][17][18] : H is the JIMWLK Hamiltonian, which for our purposes is most conveniently written as 6

5)
4 From now on, we shall denote the dependence of the Wilson lines upon the transverse coordinates with an index instead of an argument. 5 The acronym stands for Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov and Kovner. 6 We shall use in general the shorthand notation uv..
uz the Weizsäcker-Williams emission kernel: In Eq. (2.5), L and R are 'left' and 'right' Lie derivatives, i.e. the generators of local color rotations of the Wilson lines on the left and, respectively, on the right. They can be formally, but unambiguously, defined via their action on the Wilson lines, which reads For the purpose of the physical interpretation, it is however useful to notice that the Lie derivatives can be explicitly realized as functional derivatives with respect to the field α a x (x + ) at the end-points of the Wilson lines: the L derivative acts at the largest value of x + (which is positive), meaning after the scattering with the shockwave, whereas the R derivative acts at the smallest value of x + (which is negative), meaning before the scattering. Hence, each step in the evolution described by Eqs. (2.4) and (2.5) adds two infinitesimal layers, one 'to the left' and one 'to the right', to the support of the target field: with increasing Y the shockwave expands in x + , as anticipated. This expansion is however not symmetric under reflexion in x + , as shown by the structure of the JIMWLK Hamiltonian, where the Wilson lines multiply the right derivatives alone. The physical meaning of this dissymmetry will be explained later.
The Lie derivatives are non-commuting objects which obey the color group algebra (with f abc the structure constants for SU(N c )) yet one does not need to worry about the particular ordering of the derivatives in Eq. (2.5) since Note also that the 'left' and 'right' Lie derivatives are not independent operators, rather they are related by a color rotation with the Wilson line: This relation is a consequence of the unitarity of the Wilson lines and can be checked by using Eq. (2.7) together with the color group identity We are interested in the evolution of the observables with increasing Y . After taking a derivative with respect to Y in Eq. (2.1) and using Eq. (2.5), there are two ways to proceed. The first one is to integrate by parts the Lie derivatives and make them act on the operatorÔ which defines the observable. This amounts to transferring the evolution step from the wavefunction of the target to that of the projectile. From this perspective, the role of the Lie derivatives is to generate soft gluon emissions from the color sources represented by the Wilson lines. This procedure yields an evolution equation for Ô Y , which however is not a closed equation, but a member in an infinite hierarchy of coupled equations for the correlations of the Wilson lines -the Balitsky hierarchy [48]. This hierarchy becomes tractable in the limit of large N c , where the first respective equation reduces to a closed, non-linear, equation -the Balitsky-Kovchegov equation [48,49]. The solution to this equation can be combined with appropriate truncation schemes, like a Gaussian approximation [23,29,[56][57][58][59][60][61], and used to construct solutions for the higher equations in the hierarchy. The second strategy (the only one to be useful at finite N c and whenever the Gaussian approximation fails to apply), is based on the fact that Eqs. (2.4) and (2.5) correspond to a functional Fokker-Planck equation, that can be given an equivalent Langevin formulation [35], which is better suited for numerical studies [36][37][38].
In the Langevin formulation, the JIMWLK evolution of the color field in the target is depicted as a random walk in the functional space of the Wilson lines, with Y playing the role of 'time'. The CGC average of any observable, in the sense of Eq. (2.1), can be then computed as an average over the noise term in the Langevin equation which governs this random walk. For instance, for the gluonic dipole in Eq. (2.2), one writes where we have discretized the rapidity interval according to Y − Y in = N , with N → ∞ and → 0. (The 'initial' rapidity Y in is where we introduce the initial conditions for the evolution; see below.) As already mentioned, each evolution step adds two new layers in the support of the Wilson lines in x + , leading to infinitesimal, left and right, color rotations. These rotations are random, reflecting the quantum nature of the fluctuations that have been 'integrated' out. This leads to the following Langevin equation [60,62] where the left and right color matrix fields read In Eqs. (2.12) and (2.13) ν n is a Gaussian white noise local in rapidity, that is ν ia m,x ν jb n,y = 1 δ ij δ ab δ mn δ xy , (2.14) and this is the meaning of the average over ν in Eq. (2.10). Physically, the noise term accounts for the color charge density and the polarization of the target gluons radiated in this particular step of the evolution, which act as sources for the new layers α R n and α L n in the target field. These sources too are left movers, but they are not as fast as the sources produced in the previous steps. Accordingly, the gluons radiated at negative x + , meaning ahead of the shockwave, can be caught by the latter, and then they are color-rotated. This is the origin of the Wilson line visible in the r.h.s. of Eq. (2.13), which in turn is responsible for generating the BFKL cascade via iterations. Alternatively, and perhaps simpler, one can interpret Eq. (2.11) as one step in the projectile evolution; then, ν represents an on-shell gluon radiated by the projectile, which is a right mover and hence it can scatter off the target field, provided it was emitted at negative x + (i.e. prior to the scattering). Fig. 1 illustrates two steps in the evolution generated by Eq. (2.11).
For the purpose of computing correlation functions, or deriving the associated evolution equations, one needs to keep terms up to order in the r.h.s. of the Langevin equation (2.11). Note however that the noise term scales like 1/ √ , as manifest on Eq. (2.14), so in order to achieve the desired accuracy, one has to expand the left and right rotations in Eq. (2.11) up to quadratic terms. Moreover, the quadratic terms can already be replaced by their average, since they cannot be multiplied by noise terms coming from other sources to the order of interest. Figure 1. Two successive steps in the evolution described by Eq. (2.11). Gluon '2' has a larger k + than gluon ' . From the viewpoint of target evolution, cf. Eq. (2.11), gluon '2' is emitted after gluon '1', so in particular it can scatter off the color field created by the latter. From the dual viewpoint of projectile evolution, gluon '2' is emitted first and can act as a source for gluon '1'. The vertical line represents the shockwave of the nucleus, located near x + = 0.
When doing that, the quadratic terms coming from the right rotation become independent of the Wilson line in the previous step (this is trivially true for left rotations) and correspond to virtual terms. Indeed one finds In both approaches to the JIMWLK evolution (the Balitsky hierarchy and the Langevin formulation), one needs an initial condition, that is, the CGC weight-function at a given rapidity Y in . The initial condition is typically given by the MV (McLerran-Venugopalan) model [63,64] in which the color charge ρ of the nucleus is distributed according to a Gaussian distribution. Then, in principle, one can construct the initial values for all the correlators that appear in the Balitsky hierarchy. In the Langevin approach, a Gaussian initial condition means that the color charges are randomly distributed on the two-dimensional transverse plane and, in turn, this determines a probability distribution for the initial values of U in and U † in needed to start developing the Wilson lines according to Eq. (2.11).
3 Generating functional and same rapidity gluon production So far, the gluonic dipole introduced in Eq. (2.2) has been interpreted as a scattering amplitude, namely the amplitude for the elastic scattering between a pair of gluons in a color singlet state and the strong color field of the target. But the same quantity also enters the cross-section for the inclusive production of a gluon with transverse momentum p and pseudo-rapidity η p in proton-nucleus collisions: where xg(x) is the usual ('integrated') gluon distribution in the proton and x is the longitudinal momentum fraction of the projectile gluon participating in the scattering. (Note that we assume collinear factorization at the level of the proton.) In this context, the Fourier transform in Eq. (3.1) describes the transverse momentum broadening of a gluon preexisting in the projectile. The two Wilson lines withinŜ xx now describe the same gluon, which has transverse coordinate x in the direct amplitude (DA) and respectivelyx in the complex conjugate amplitude (CCA). The color trace and the normalization factor 1/N g have been generated by the sum (average) over the gluon color indices in the final (initial) state.
In what follows, we will be mostly interested in the radiation of new gluons triggered by the interactions between the projectile -that will be taken to be simply a gluon -and the target. To compute the corresponding cross-section, it is convenient to introduce a generalization of the 'gluonic dipole' in Eq. (2.2), which distinguishes between the Wilson lines in the DA (for which we shall use the same notations as before, i.e. U and U † ) and those in the CCA (to be denoted with a bar:Ū andŪ † ). Specifically, we introduce the following functional of U andŪ To start with, consider the emission of a second gluon which is softer than the first one, but not much softer: the emission vertex can be described in the eikonal approximation, but the rapidity separation between the two gluons remain small enough to allow one to neglect the high-energy evolution in between. In that sense both gluons have the same rapidity difference Y w.r.t. the target. Let (η p , p) and (η k , k) be the pseudo-rapidity and transverse momentum of the parent and the emitted gluon respectively, withᾱ(η p − η k ) 1. The cross-section for inclusive two-gluon production can be computed as with the production Hamiltonian [32][33][34] H prod (k) = 1 4π 3 where in the r.h.s. there is no 'bar' anymore, neither on the Wilson lines nor in the functional derivatives. It is a straightforward exercise to perform the various derivatives, then integrate over u and v, and thus obtain (up to a factor δ(1 − x p − x k ) expressing the conservation of the plus component of the momentum; here, x p and x k are the longitudinal momentum fractions taken by the two final gluons and are such that x k 1 and x p 1) in agreement with the calculation in [31]. The above manipulations can be easily generalized to the production of several gluons having similar rapidities. For instance, in order to produce two gluons with transverse momenta k A and k B and such that gluon A is softer than gluon B (η B > η A , but such thatᾱ(η p − η A ) 1 andᾱ(η p − η B ) 1), one needs to act twice with the production Hamiltonian on the generating functional and only then setŪ = U , as follows: That is, the softer gluon (here, gluon A) is produced after the harder one (gluon B), and in particular it can be emitted by the latter: the functional derivatives inside H prod (k A ) can also act on the Wilson lines which enter the structure of H prod (k B ).
Returning to the cross-section for two gluon production, Eq. (3.6), notice that if we integrate over p and over η p , that is, we do not measure the parent gluon after the scattering, the first phase factor in the integrand of Eq. (3.6) leads to (2π) 2 δ (2) (x − x) in Eq. (3.6) and one recovers the well-known result for single gluon production out of gluon source [19] If on the other hand, one integrates Eq. (3.6) over k and η k , that is, one measures just the parent gluon, then one obtains a part of the radiative corrections to the cross-section in Eq. (3.1) for single inclusive gluon production -namely, the 'real' part associated with the emission of a gluon which is not measured. The complete radiative corrections (to leading logarithmic accuracy at high energy) also include 'virtual' processes, i.e. diagrams where the evolution gluon is emitted and reabsorbed on the same side of the cut. For the single gluon production in Eq. (3.1) and also for the production of two gluons with roughly the same rapidity, the quantum evolution with increasing energy is controlled by the usual Balitsky-JIMWLK evolution of the Wilson line correlators which enter Eqs. (3.1) or (3.6). But in the case where the produced gluons are widely separated in rapidity, we need a generalization of the Balitsky-JIMWLK equations, to be described in the next section. between the two measured ones, as shown in Fig. 3. The first gluon is assumed to be close to the projectile rapidity (that is, its rapidity difference w.r.t. the target is equal to Y ) and thus can be viewed as the "parent" gluon. The second, "emitted", gluon has a relative rapidity Y A w.r.t. the target, such thatᾱ(Y − Y A ) > 1. In order to emit this gluon we need to act with the corresponding production Hamiltonian H A prod on the WFS of the projectile evolved from Y A up to Y . Here, H A prod is given by an expression similar to Eq. (3.4), but where the Wilson lines and the functional derivatives refer to the rapidity Y A ; we shall denote the corresponding Wilson lines as U A (in the DA) andŪ A (in the CCA). Accordingly, we need to evolve the generating functional of the projectile gluon from rapidity Y A to rapidity Y , with an initial condition at Y = Y A given by Eq. (3.2) with U → U A andŪ →Ū A . The answer reads with the superscript A standing for the functional dependence on U A andŪ A . The QCD evolution in the above is encoded in the conditional CGC weight-function, which satisfies [32] where H 11 is the JIMWLK Hamiltonian given in Eq. (2.5), H 22 is obtained by putting a bar in all Wilson lines and derivatives in Eq. (2.5), and H 12 by putting a bar only in the quantities in the second square bracket factor. This conditional weight-function obeys the initial condition Moreover, at any Y it satisfies the following, important, property states that the DA and the CCA evolve in the same way, so the respective field configurations coincide with each other at any Y provided they do so in the initial condition at Y = 0. This is of course the same property which allows one to relate cross-sections for particle production like Eqs. (3.6) and (3.7) to forward scattering amplitudes (hence, to total cross-sections).  It is easier to interpret these diagrams in terms of projectile evolution, i.e. by assuming that the evolution step is achieved by boosting the projectile gluon, and not the nuclear target. This will be shortly discussed in more detail.
The corresponding cross-section can be given a formal expression in terms of a path integral: first we act with H A prod ≡ H prod (k A )[U A ,Ū A ] on the generating functional in Eq. (4.1), then we letŪ A = U A and finally we average over where the proportionality constant (including the Fourier transform) can be read from Eq. (3.3). This equation can be suggestively and more succinctly rewritten as where the notation emphasizes the fact that the present calculation involves two types of target averaging: one involving the conditional CGC weight-function, which encodes the evolution from Y A up to Y , and one using the standard weight-function, for the evolution from Y in up to Y A . So far, we have privileged the viewpoint of target evolution, as manifest on equations like (4.1) and (4.5). This turns out to be more fruitful for our present purposes, and we shall indeed stick to it in most of our subsequent developments. But for the sake of the physical interpretation, it is useful to notice that the evolution from Y A up to Y can alternatively be viewed as the BFKL evolution of the projectile in the presence of the strong target field (evolved up to rapidity Y A ). This appears to be more complicated than the usual BFKL evolution of the gluon distribution in the projectile, in that it involves both 'initial-state' and 'final-state' emissions: the evolution gluons can be emitted and/or reabsorbed both before and after the scattering with the shockwave, as illustrated in Fig. 4. This introduces a dependence upon the background field (the Wilson line U A ), via processes in which the evolution gluon crosses the shockwave only once; such processes are generated by terms like LU R or LUR in H evol . For the processes where the gluon crosses the shockwave twice, as generated by RUŪ †R , the Wilson lines Figure 5. Cancellations in the final-state evolution of the "parent" gluon when it is not measured. The two diagrams, generated by RU L and RUL, add to zero since they can be obtained from one another by simply moving a gluon vertex attachment from the DA to the CCA. cancel between the DA and the CCA, by unitarity. (More precisely, this cancellation occurs after 'producing' the gluon at rapidity Y A , that is, after acting with H A prod and then lettingŪ A = U A .) But the usual BFKL evolution is recovered in the case where the parent gluon is not measured, as recovered by chosingx = x in the previous formulae. In that case, the final-state evolution cancels out between the DA and the CCA, as illustrated in Fig. 5. This cancellation occurs already at the level of the generating functional, i.e. before acting with H A prod : forx = x, the evolution of the generating functional in Eq. (4.1) is controlled solely by the terms involving 'right' derivatives in H evol -that is, the terms proportional to RR,RR, and RR -, which describe the initial-state evolution of the projectile. This is required by causality: as their name suggest, the 'final-state emissions' occur after the gluon at rapidity Y A has been produced and hence they cannot influence its emission.
The simplest way to render manifest the above physical picture is by constructing evolution equations for quantities like Ŝ 12 -the projectile generating functional and respectively the cross-section for producing a gluon at Y A . Such equations can be derived in the same way as the Balitsky equations: starting with the CGC representation (4.1) for the generating functional, one takes a derivative w.r.t. Y , uses the generalized JIMWLK equation (4.2), and then integrates the action of H evol by parts, to make the functional derivatives act onŜ 12 (xx). The evolution equations thus obtained generalize the Balitsky-JIMWLK hierarchy to multi-particle production. The first few equations in this new hierarchy will be presented, together with their physical interpretation, in Appendix A. (At large N c , similar equations have been previously constructed in Refs. [25,34].) But these equations appear to be even more involved than the original Balitsky equations. Although some simplifications occur at large N c (where only dipoles and quadrupoles remain as the relevant degrees of freedom; see Appendix A), it appears very difficult to make progress with them in practice. In what follows, we shall return to the viewpoint of target evolution and propose a Langevin formulation for expressions like Eq. (4.5), that we expect to be better suited for numerical solutions.

The Langevin description of multi-particle production
To arrive at a Langevin formulation for the two-gluon cross-section in Eq. (4.5), we start with Eq. (4.1) for the evolution of the projectile generating functional over the intermediate rapidity range Y − Y A . In an analogous way to Eq. (2.10), one can write Here, the rapidity interval is discretized as Y − Y A = N A and the Wilson lines are built starting at n = 0, corresponding to Y 0 ≡ Y A , with U 0 ≡ U A andŪ 0 ≡Ū A . Specifically, the Wilson lines in the DA are built according to Eqs. (2.11) -(2.14), while for the CCA we similarly writē where the left and right matrix fields now read Importantly, the noise term ν n in these stochastic equations is exactly the same as in the corresponding equations, (2.12) -(2.13), for the DA -in particular, it satisfies Eq. (2.14). This is in agreement with the fact that, as already discussed, the evolution of the DA and that of the CCA are strictly correlated with each other: a gluon which is emitted say in the DA can be then absorbed either in the DA, or in the CCA. It is furthermore consistent with the structure of the evolution Hamiltonian in Eq. , this means that the 'left' matrix field in the CCA is the same as that in the DA (ᾱ L n,x = α L n,x ), while the 'right' one in the CCA (ᾱ R n,x ) differs from the corresponding one in the DA (α R n,x ) only through its dependence on the Wilson line of the previous step. Therefore any difference between the 'barred' and 'unbarred' Wilson lines at the final rapidity Y can be traced back to a difference between the respective initial conditions (U A andŪ A ) at rapidity Y A .
Although mathematically well defined, the stochastic process that we have just described is not well suited for numerical simulations, because of the need to keep trace of the functional dependence upon the initial Wilson lines U A andŪ A . Yet, as we now explain, one can set up an alternative Langevin process, which is oriented towards the physical problem at hand -that is, it depends upon the specific structure of the observable (here, the cross-section for two gluon production) -and which circumvents this problem: in this new process, the initial Wilson lines are (random) numbers (more properly, numerically-valued color matrices), and not functions.
To introduce this alternative procedure, let us recall that the evaluation of the cross-section of interest requires the action of H A prod on the evolved generating functional in Eq. (5.1). This in turn involves the linear combination of 4 terms, such as , there is no need to distinguish between the DA and the CCA anymore. Yet, one still needs to cope with the functional dependence upon the initial conditions, which is the problem alluded to above. Our proposal to circumvent this problem is to enlarge the Langevin process, in such a way to include the building blocks of expressions like Eq. (5.5) -that is, R A U † and L A U † .
To this aim, one must supplement the Langevin equation (2.11) for the Wilson lines, which is local in the transverse plane, with the following, bilocal, recurrence formula for R A U † , with ν i n,z ≡ ν ib n,z T b the noise color matrix in the adjoint representation. This equation follows from Eq. (2.11) after using the Leibniz rule for differentiation together with a few manipulations that we shall now explain. The emergence of the first term in the r.h.s. being pretty obvious, let us concentrate on the second term, which expresses the action of R a A on the 'right' infinitesimal color rotation. To obtain this term, we have kept only the linear term in the expansion of exp[−i gα R n,x ] : indeed, the zeroth order term (the unity) gives trivially zero under the action of R a , and the same is also true for the quadratic term, which is effectively independent of U † n−1 to the accuracy of interest (since it can be averaged over the noise; recall the discussion after Eq. (2.11)). After this expansion, one is led to evaluate (cf. Eq. (2.13)) where we have also used the identities U †cb T b = U T c U † and R a U = −U (R a U † )U (with the latter arising from R a (U U † ) = 0). In writing Eq. (5.6), we have used the same discretization conventions as in relation with Eq. (5.1). That is, we have written Y − Y A = N A and used n with n = 1, 2, . . . N A to denote a generic intermediate step. The initial condition, represented by n = 0, corresponds to Y 0 ≡ Y A and U 0 ≡ U A , and reads The quantity L a A,u U † n,x obeys an equation identical to Eq. (5.6), but with a different initial condition, as also indicated above. However, there is no need to separately consider the associated equation, because of the relation between 'left' and 'right' Lie derivatives mentioned after Eq. (2.9): once that the quantity R a A,u U † N A ,x is obtained by iterating Eq. (5.6), one can imme- x . Also, one does not need to separately calculate R a A,u U n,x : it is an easy exercise to prove by induction that it can be obtained by taking the hermitian conjugate of R a A,u U † n,x . To conclude, we just need to follow two independent Langevin processes: Eq. (2.11) for the Wilson line U † n,x and Eq. (5.6) for its 'right' Lie derivative R a A,u U † n,x . Eq. (5.6) introduces a numerical complication as compared to Eq. (2.11): unlike Eq. (2.11), which is local in the transverse coordinates, Eq. (5.6) is bi-local. Actually, the corresponding initial condition is still local, because of the presence of the factor δ ux in Eq. (5.8), but already its first iteration acquires a bi-local structure, as one can easily check. In spite of this potential complication, which is purely numerical, the procedure that we have just outlined has a decisive advantage over Eq. (5.1), namely it allows for a fully numerical implementation, at least in principle. Indeed, the initial condition U A for Eq. (5.6) is not a generic color matrix anymore, but the numerical matrix that has been generated in the previous steps of the Langevin process, from Y = Y in (where U = U in is randomly selected according to the Gaussian MV distribution), up to Y = Y A (where U = U A ). Of course, during the first part of this evolution, from Y in up to Y A , there is no analog of Eq. (5.6) : this part involves just the usual Langevin process, Eq. (2.11), for the Wilson lines.
Note also an alternative form for the Langevin equation (5.6), which is perhaps more convenient for numerics, in that it involves the infinitesimal right rotation exp[i gα R n,x ] alone 7 : Let us observe here that Eq. (5.9) should be consistent with the fact that U n,x R a A,u U † n,x must be a member of the Lie algebra. This is manifest for the first term on the r.h.s. of Eq. (5.9), but it is not so obvious for the second term due to the presence of a single infinitesimal right rotation. However, one should keep in mind that the equation is valid only up to order and thus, recalling that the noise scales like 1/ √ , one can expand exp[i gα R n,x ] to first order. Then the second term of Eq. (5.9) becomes where we have been allowed to take the average for the contribution quadratic in the noise. It is now clear that both terms in Eq. (5.10) are members of the Lie algebra. To summarize, by iterating the recurrence formula Eq. (2.11) over the whole rapidity interval between Y in (where we insert the initial condition of the MV type) up to Y , and the new formula (5.6) (or (5.9)) over the upper rapidity range between Y A and Y , one can construct the r.h.s. of Eq. (5.5) via a fully numerical procedure. The only potential difficulty that we can foresee is a numerical complication associated with the bi-local structure of Eq. (5.6) in transverse coordinates. At this point it is should be clear that the method can be extended to the production of more than two gluons, but it becomes more and more tedious (since increasingly non-local) with increasing the number of measured gluons. As an illustration, we shall outline the corresponding procedure for three gluon production in Appendix B.

A Evolution equations for the projectile generating functionals
As discussed at the end of Sect. 4, one method to study the evolution of the generating functionals with Y is to construct the respective evolution equations, which generalize the Balitsky hierarchy to multi-particle production. These equations are obtained by acting on the generating functionals with the evolution Hamiltonian H evol in Eq. (4.2). Here, we are interested only in illustrating the procedure, and for this purpose we shall consider the simplest generating functional, that of a physical quark. This is represented by a mathematical dipole, similar to that in Eq. (3.2), but where the Wilson lines are now in the fundamental representation and will be denoted by V . Differentiating the quark analog of Eq. (4.1) w.r.t. Y , using Eq. (4.2), and integrating by parts, we find with F standing for the fundamental representation. A straightforward exercise, similar to the derivation of the Balitsky hierarchy from the JIMWLK equation, leads to the following equation (after also using the Fierz identities to rearrange the terms) InŜ F 11 both Wilson lines live in the DA, while inŜ F 22 they both live in the CCA, that is, The mathematical quadrupoleQ F 1221 stands for the generating functional of a physical quarkantiquark dipole, and readŝ where x and z are the coordinates of the quark and the antiquark in the DA, whilex andz similarly refer to the CCA. Below, we shall be interested in emitting a gluon at rapidity Y A out of a quark with rapidity Y > Y A . For that purpose, Eq. (A.2) has to be integrated from Y A up to Y , with a functional initial condition at Y A : is not a closed equation, rather it involves new generating functionals in the r.h.s., for which we need to construct the respective equations. This complication should not be a surprise, given our experience with the Balitsky hierarchy. But Eq. (A.2) looks considerably more involved than the corresponding Balitsky equation (that for the dipole S-matrix), both because of the complex structure of its r.h.s. and because of the functional initial conditions. Like in the Balitsky hierarchy, important simplifications occur at large N c . But before we discuss them, let us notice some special cases of Eq. (A.2) : (a) If one setsV A = V A in the initial conditions, that is, if one identifies the Wilson lines in the DA and the CCA at the initial rapidity, then this property remains true at any Y > Y A , cf. Eq. (4.4), and then there is no difference anymore betweenŜ F 11 ,Ŝ F 22 , andŜ F 12 . In that case, one can easily see that the last term in Eq. (A.2) vanishes, while the first two combine to give the 'BK equation' (more properly, the first equation in the Balitsky hierarchy). In the present context, this equation describes the evolution with Y of the cross-section for single inclusive quark production (cf. Eq. (3.1)). But of course, after identifyingV = V , one loses 8 the very meaning of a 'generating functional' : one cannot useŜ F 12 to emit gluons anymore. (b) The case where the parent quark is not resolved is obtained by settingx = x, and then the first two terms in Eq. (A.2) individually vanish. This is the expected cancellation of the final state evolution between the DA and the CCA (cf. Fig. 5). Correspondingly, the surviving term in the third line is associated with initial state evolution alone. Indeed, one can verify that this term could directly be obtained by keeping only the 'right' derivatives in H evol , that is, the three terms which involve RR,RR, and RR.
Since the 'quadrupole' generating functional in Eq. (A.4) (the WFS of a physical dipole) enters the evolution equation (A.2) for the 'dipolar' one (the WFS of a physical quark), it is useful to also have a look at the corresponding evolution equation. After some algebra, this is found as where M is the 'dipole kernel', which vanishes when u = v. The initial condition at Y = Y A is given by Eq. (A.4) with V → V A andV →V A . As before one can recognize a couple of special cases in Eq. (A.5) : (a) In the limit of where one identifiesV = V ('no quantum fluctuations'), one recovers the standard Balitsky equation for the evolution of the S-matrix of a physical quadrupole [25,33]. In the present context, this equation could refer e.g. to the cross-section for producing a quarkantiquark pair in deep inelastic scattering off the nucleus [29].
(b) If the parent quark and antiquark are not resolved, one setsx = x andz = z, and then Eq. (A.5) reduces to its very last term, which describes the evolution of the WFS of a qq color dipole via initial state emissions alone. In particular at large N c one can factorize 8 Alternatively, one can say that one loses the 'quantum fluctuations', by which we here mean the gluons which can be emitted or absorbed by the projectile and which are encoded in the (functional) difference betweenV and V . In the Keldysh-Schwinger formalism, this would refer to the difference between the quantum fluctuations on the two branches of the closed time contour. Q FQF Q F Q F (see below), and then one recovers, as expected, the evolution equation for the generating functional of an 'onium' (the system of dipoles generated via the BFKL evolution of an original dipole, in the limit where N c 1) [47]. We are now in a position to describe the simplifications which occur at large N c . First, expectation values of products of color traces can be factorized into products of averages of the individual traces. This property is well known to hold for the solutions to the Balitsky-JIMWLK equations provided it is already satisfied by the initial conditions at Y = Y in (as indeed happens within the MV model at large N c ). The respective argument can be taken over to the evolution Hamiltonian H evol in Eq. (4.2) and it implies a similar factorization for the generating functionals; e.g.
Moreover, still at large N c , the quantities Ŝ F 11 Y and Ŝ F 22 Y lose their meaning as generating functionals, because the two functional derivatives inside H prod must act on Wilson lines from a same color trace, to generate the maximal power of N c . Hence, inside these quantities one can setV = V already before acting with H prod , and then they reduce to the standard dipole S-matrix (the solution to BK equation). Similar simplifications apply to Eq. (A.5) for the quadrupole.
We conclude that, in the large N c limit, the correspondingly simplified versions of Eqs. (A.2) and (A.5) (withz = z) form a 2 × 2 system for Ŝ 12 Y and Q 1221 Y . In that sense, multiparticle production at large N c involves only dipoles and quadrupoles, as already noticed in the literature [25,34,66]. In fact, the large N c versions of Eqs. (A.2) and (A.5) are equivalent to the respective equations constructed in Refs. [25,34] directly at large N c . The coefficients in these equations also involve the dipole S-matrix Ŝ F Y , which here describes multiple scattering for the (unresolved) gluons associated with final-state evolution.
But even at large N c , the equations above have the drawback that they must be solved with functional initial conditions, which is not very convenient in practice. It is therefore interesting to notice that ordinary evolution equations -i.e. equations for the production cross-sections which describe the evolution from Y A up to Y (at large N c ) and admit numerical initial conditions at Y A -can be deduced by acting with H A prod on the above equations. We recall that the derivatives inside H A prod act on the functional dependence of the WFS upon the Wilson lines V A andV A , as introduced via the initial conditions at Y A . Also, after acting with H A prod , one has to set V A =V A (which ensures thatV = V at any Y > Y A ) and then average over To be more precise, let us consider the evolution of the cross-section for quark-gluon production. Starting with the large N c version of Eq. (A.1) and performing the manipulations explained above, one finds where Ŝ F Y is the solution to the BK equation and For more clarity, we have indicated through our notations that the generating functionals have to be evolved from Y A up to Y . Also, we have used the 'double-averaging' notation introduced in Eq. (4.6). Up to an overall factor and a Fourier transform, N k (xx) Y is the cross-section for quark-gluon production, with the quark at rapidity Y and the gluon at rapidity Y A < Y (cf. Eq. (3.3)). Similarly, N k (xxzz) Y is proportional to the cross-section for producing a gluon with momentum k at rapidity Y A out of a qq dipole at rapidity Y and such that the quark and the antiquark can be measured as well. To get closed equations, one also needs the evolution equation obeyed by the last quantity, at least for the case where the antiquark is not measured (z = z). This follows from Eq. (A.5) and reads ∂ N When one also setsx = x, i.e. when neither the quark nor the antiquark are being measured, the first two terms in the r.h.s of this equation individually vanish, while the surviving terms reproduce the BFKL equation for N (2) k (xxzz), as expected: the cross-section for producing a gluon out of the onium evolves with the energy in the same way as the dipole (or gluon) distribution of the onium, since the gluon can be emitted by any of the dipoles composing the onium. Note that the 'color transparency' property Ŝ F (xx) Y = 1 has been important too for this argument (it has been used to simplify the last term in Eq. (A.9)). In the present context, this property expresses the cancellation of the 'final-state interactions' of the unmeasured quark -i.e. its interactions with the nuclear target which occur after emitting the gluon -between the DA and the CCA.
Eqs. (A.7) and (A.9) must be integrated from Y A up to Y , with initial conditions given by the respective cross-sections (for quark-gluon production and respectively dipole-gluon production). These equations are linear w.r.t. the quantities that one needs to solve for, i.e. the cross-sections N k (xx) Y and N (2) k (xxzz) Y , but in general they include multiple scattering effects via the dipole S-matrix Ŝ F Y . Such effects are associated with the final-state evolution, or the finalstate interactions, of the parent quark and disappear whenx = x, i.e. when this quark is not measured 9 . In that limit, Eqs. (A.7) and (A.9) describe the BFKL evolution of the unintegrated gluon distribution in the quark and respectively dipole projectile, at large N c . 9 There is of course another limit in which the effects of multiple scattering can be neglected, including for the case where the parent partons are measured, i.e. for multi-particle production: this is the limit where the target itself is dilute, as e.g. in proton-proton collisions at not too high energies. The corresponding limit of Eqs. (A.7) and (A.9) is obtained by replacing Ŝ F Y → 1 for all the dipole S-matrices which appear in these equations. This amounts to neglecting the target rescattering for both the projectile partons and the evolution gluons. The only interactions with the target which survive in this limit are those of the produced gluon at the lowest rapidity YA. They enter the solutions to the simplified equations via the respective initial conditions at Y = YA, themselves computed in the single scattering approximation. After this replacement, the equations become fully linear and describe the evolution of the BFKL Green's function (for quark and dipole impact factors respectively).
To summarize, the general strategy in this case would be as follows: (i) First construct approximate solutions to the usual Balitsky-JIMWLK equations for the dipole and the quadrupole S-matrices, Ŝ F (xx) Y and Q F (xxzz) Y . For instance, one can solve the BK equation for Ŝ F (xx) Y and then use the Gaussian approximation to the JIMWLK evolution to compute the quadrupole in terms of the dipole. (ii) Then, use the aforementioned (approximate) solutions in order to compute the cross-sections N k (xx) Y A and N (2) k (xxzz) Y A at rapidity Y A . The respective expressions are well known: the cross-section for quark-gluon production looks very similar to Eq. (3.5) and can be found in [27], whereas that for dipole-gluon production is presented in [25]. (iii) The cross-sections obtained in the previous step are now used as initial conditions for Eqs. (A.7) and (A.9), to be integrated from Y A up to Y . Namely, one needs to first solve Eq. (A.9) for N (2) k Y , then plug the result into Eq. (A.7), and solve the latter for N k Y .
(iv) Finally one needs to take the Fourier transform as in Eq. (3.3), to ascribe a transverse momentum p to the measured quark.
This being said, this procedure appears as prohibitively cumbersome, including for numerical simulations, and we believe that the Langevin procedure described in the main text should be more tractable in practice.

B Three gluon production
Here we shall discuss the extension of the Langevin procedure to the production of three gluons at (Y A , k A ), (Y B , k B ) and (Y, p) with Y > Y B > Y A . The analog of Eq. (4.5) is As in the corresponding discussion in Sect. 5, this equation is easier to read by following the evolution backwards in rapidity, from Y down to Y in . For generic Wilson lines U B andŪ B at rapidity Y B , one needs to first compute the action of H B prod on the generating functional of the projectile gluon evolved from Y B up to Y : As usual, we consider a discretization of the relevant rapidity interval according to Y −Y B = N B . The result of this step is itself a functional of U B andŪ B . Next, the Wilson lines U B andŪ B must in turn be built from generic Wilson lines U A andŪ A , via iterations which cover the rapidity interval Y B − Y A = N A ; as a result, Eq. (B.2) becomes a functional of U A andŪ A . Then we can act with H A prod and subsequently setŪ A = U A , with U A constructed by evolving the initial condition (say, as given by the MV model) from Y in to Y A , via the standard JIMWLK evolution. Note that by imposingŪ A = U A , we automatically enforceŪ B = U B , and eventuallȳ U N B = U N B , because of the property (4.4) of the conditional CGC weight-function.
But although conceptually clear, the above procedure -which would involve only the standard local Langevin Eq. (2.11), but with functional initial conditions at Y B and Y A -is not convenient in practice. To achieve a fully numerical implementation for the three gluon production, one should rather proceed upwards in rapidity and construct in each step all the building blocks (Wilson lines and their Lie derivatives) which enter the cross-section in Eq. (B.1). In view of the procedure we followed in Sect. 5 we can rewrite Eq. (B.1) as At this point let us first recall that 'left' Lie derivatives can be written in terms of 'right' ones (cf. the discussions after Eqs. (2.9) and (5.8)). This means that H prod , so long as the square brackets in Eq. (3.4) are considered, can be rewritten as As in the case of the two-gluon production, it is enough to study the (double) action of H prod on the Wilson line in the DA. Denoting by y A and y B the transverse coordinates of the produced gluons in the DA at Y A and Y B respectively, one needs to calculate quantities of the form where we recall that R A and R B are the Lie derivatives w.r.t. the Wilson lines U A and U B . All the structures appearing in Eq. (B.5) can be constructed via a Langevin procedure. Specifically, at Y = Y in , one starts with U † in,x , which is randomly selected according to the MV model, and uses the Langevin equation (2.11) to build U † A at Y = Y A . Then, using the procedure described in Sect. 5, one evolves from Y A to Y B , to construct U † B and R A U † B . Finally one evolves from Y B up to Y to build R B U † and R A R B U † (during which process, one also needs to follow U † n and R A U † n ). All these quantities can thus be numerically computed, but the respective Langevin process is more involved than in the case of two-gluon production: during the uppermost rapidity interval from Y B to Y , the structure R b A,u R d B,w U † n,x is tri-local in the transverse plane. This method considerably simplifies in the case where two of the three produced gluons have similar rapidities, say Y Y B in Eq. (B.1). Indeed, this avoids the most complicated step in the above procedure, that of the evolution between Y B and Y . The corresponding calculation is similar to that of two-gluon production discussed in Sect. 5, except for the fact that the projectile emitting the soft gluon at rapidity Y A is not a bare gluon anymore, but rather a system of two gluons, one of which (the gluon B) has been emitted by the other one (the 'parent' gluon). Hence, the generating functional describing the WFS of the emitter is not justŜ 12 (xx) anymore, but rather H prod (k B )Ŝ 12 (xx), which is explicitly shown 10 in Eq. (3.6). The cross-section for producing three gluons with rapidities Y Y B Y A can therefore be computed by replacingŜ 12 (xx) → H prod (k B )Ŝ 12 (xx) in the procedure outlined in Sect. 5, that is, in equations like (5.1) and (5.5).