Multi-particle production in proton-nucleus collisions in the Color Glass Condensate

We compute multi-gluon production in the Color Glass Condensate approach in dilute-dense collisions, p$A$, extending previous calculations up to four gluons. We include the contributions that are leading in the overlap area of the collision but keep all orders in the expansion in the number of colors. We develop a diagrammatic technique to write the numerous color contractions and exploit the symmetries to group the diagrams and simplify the expressions. To proceed further, we use the McLerran-Venugopalan and Golec-Biernat-W\"usthoff models for the projectile and target averages, respectively. We use a form of the Lipatov vertices that leads to the Wigner function approach for the projectile previously employed, that we generalise to take into account quantum correlations in the projectile wave function. We provide analytic expressions for integrated and differential two gluon cumulants and show a smooth dependence on the parameters defining the projectile and target Wigner function and dipole, respectively. For four gluon correlations we find that the second order four particle cumulant is negative, so a sensible second Fourier azimuthal coefficient can be defined. The effect of correlations in the projectile on this result results qualitatively and quantitatively large.

in [56] suggested that in such approximation c 2 {4} > 0 -a result also found in [45] where only quark scattering is considered and partons in the projectile wave function are uncorrelated. In this work we use the argument in [53,[57][58][59] that captures those contributions of the ensembles of Wilson lines to multiparticle production that are leading in the overlap area of the collision (i.e., in the number of color domains or correlated particle sources), while keeping contributions to all orders in the number of colors. We use the Golec-Biernat-Wüsthoff (GBW) model [60,61] for the target, and the generalised McLerran-Venugopalan (MV) model [62,63] for the projectile. In order to push the analytical calculations as far as possible, we employ the Wigner function ansatz used in [45,64,65] but extended to include quantum correlations in the projectile wave function. Due to the Gaussian forms that we employ for both the Wigner function and dipole, our results cannot be considered reliable for transverse momenta sizeably larger than the saturation scale.
This manuscript is organized as follows. In Section II we introduce the formalism to compute particle production in the CGC. Section III is devoted to the calculation of projectile and target averages required to obtain the final results. Then, in Section IV we present our results and in Section V we give a summary and our conclusions. Appendices A, B, C and D contain a discussion on the validity of the area enhancement argument that we use for computing target ensembles of Wilson lines, useful integrals, a discussion of the Wigner function approach and details on the calculation of four gluon production, respectively.

A. Gluon production in dilute-dense collisions
In this section we present a quick overview on multi-particle production in proton-nucleus collisions in the CGC framework. We follow [66,67] and references therein. The projectile is considered a highly boosted dilute system that is composed, mostly, by large-x partons that act each as a color source with color charge density ρ a (x), with superindex a denoting color and x the transverse position. The target is characterised by a strong field A µ (x) = A aµ (x)T a , with T a the generators of the SU(N c ) group in the adjoint representation. The nucleus ensemble is supposed to be much larger than the projectile in the transverse plane. In this picture, working in the light-cone gauge A + = 0 and neglecting the transverse components of the field, the amplitude for producing a gluon with transverse momentum k, pseudorapidity η, polarization λ and color a in the projectile-target collision is obtained by using the LSZ reduction formula (at leading order in the QCD coupling constant g) leading to M a λ (η, k) = g with ρ a (q) the Fourier transform of the color charge density of the projectile which is defined (see e.g. [68]) as ρ a (x) = dp + 2π a † i (p + , x)T a a i (p + , x), where a † i (p + , x) and a i (p + , x) are the creation and annihilation operators for gluons with longitudinal momentum p + at transverse position x, respectively. The reduced matrix amplitude M ab λ (η, k, q) is derived in [66,69,70] and it reads y,x e i(k−q)y−ikx G ab k + (L + , x; 0, y) + x,y e i(k−q)y 1 k + L + 0 dy + e −ikx [∂ y i G ac k + (L + , x; y + , y)]U cb y (y + , 0) .
In this equation x ≡ d 2 x, i * λ (k) is the polarisation vector, k − = k 2 /(2k + ), k + = e η / √ 2, L + is the longitudinal length of the target, q is the transverse momentum transferred from the target during the interaction and k-q is the transverse momenta of the projectile color charge density, and G ab k + (x + , x; y + , y) = Θ(x + − y + ) x y Dz exp ik + 2 x + y + dz +ż2 (z + ) U ab z x + , y + is the scalar gluon propagator, with the path integral taking into account the Brownian motion of the gluon in transverse plane, for fixed ends of the trajectory z(x + ) = x, z(y + ) = y. We use light-cone coordinates x ± = (x 0 ± x 3 )/ √ 2.
U ab x (x + , y + ) = P exp ig is the Wilson line that accounts for the multiple gluon exchanges with the target. It is necessary to mention that the reduced matrix amplitude M ab λ (η, k, q) given in Eq. (3) is written for a target with a longitudinal width L + and therefore goes beyond the standard eikonal approximation commonly adopted in CGC calculations. In the eikonal approximation (see [35] and references there) the target and projectile are taken as very highly boosted systems without longitudinal extent because of Lorentz contraction. This is equivalent to taking the limit L + → 0, k + → ∞ and assuming that the target field is a local shock-wave, A − (z + , x) ∝ δ(z + ). In the present work, we restrict ourselves to the standard CGC framework and adopt the eikonal approximation. Within this approximation, the scalar gluon propagator simplifies and can be written as Consequently, the reduced matrix amplitude given in Eq. (3) simplifies as well and it reads M ab λ (k, q) = 2i i * λ (k)L i (k,q) y e −iqy U ab (y), where we have introduced the Lipatov vertex and changed the notation of the Wilson lines, U ab (y) = U ab y (L + , 0). The physical interpretation of the Lipatov vertex, see Fig. 1, is such that the first element in the sum in Eq. (8) accounts for interaction of the color source ρ a (x) with the target before emitting the gluon and the second element accounts for a gluon being emitted from the source and then interacting with the target.
In [71,72] it was shown that the corrections with respect to the eikonal approximation stemming from the target having a finite length can be important for collision energies below a few hundred GeV. Corrections coming from the inclusion of transverse components of the background field have also been considered in [73][74][75][76][77], but until now no estimation is available of their quantitative impact on particle production. In the remainder of this work we restrict to the eikonal approximation.
The multiplicity for producing n-gluons with transverse momentum k i , pseudorapidity η i , color a i and polarization λ i is given, in terms of the amplitude matrix that is leading for gρ a (q) ∼ 1, as 2 n (2π) 3n d n N n i=1 dη i d 2 k i = M a1 λ1 (η 1 , k 1 ) · · · M an λn (η n , k n ) M an λn (η n , k n ) † · · · M a1 λ1 (η 1 , where · · · p,T denotes the average over the color charge density configurations of the projectile and target. The factor of 2 n on the right hand side of Eq. (9) originates from the Lorentz invariant phase space written in terms of rapidity η i . Using Eq. (1) and dropping the dependence on η due to the eikonal approximation, we can write this expression as (see Fig. 2) 2 ρ b1 (k 1 − q 1 )ρ b2 † (k 1 − q 2 ) · · · ρ b2n−1 (k n − q 2n−1 )ρ b2n † (k n − q 2n ) p × M a1b1 λ1 (k 1 , q 1 )M b2a1 † λ1 (k 1 , q 2 ) · · · M anb2n−1 λn (k n , q 2n−1 )M b2nan † λn (k n , q 2n ) Solving this equation is the main point of this work and will be the focus of the discussion in the next sections. Eq. (10) involves the 2n-point correlation functions of the color charge densities of the projectile and the target. Solving these objects is a highly non trivial task. It is known that using the MV Gaussian weight [62,63] it is possible to find a closed-form solution for these correlators [45,[78][79][80][81][82][83][84]. However, the solution can be extremely complicated for n > 2 if the large-N c limit is not taken. Corrections to the Gaussian weight for the 2-point correlator were considered in [85], and solutions based on high-energy evolution at large N c in [86], but in this work we restrict ourselves to the MV model.
When |k|/Q s 1 with Q s the saturation momentum of the gluons in the target [33,34], or equivalently gA − 1, we can expand the product of 2n Wilson lines up to order g 2n . Thus, if a Gaussian weight is chosen for the target color charge density we can apply Wick's theorem and write the 2n-point functions as the sum of (2n − 1)!! products of n 2-point functions, thus being able to solve the correlator exactly. This is the glasma graph approach previously mentioned and has been used in several works [46,55,[87][88][89][90] to produce phenomenological results. The main disadvantage of this approximation is that it works in a small kinematic range, being only suitable for dilute-dilute collisions.
Another approach for evaluating the n-point functions that keeps the simplicity of the glasma graph approximation but without having to restrict ourselves to the dilute limit is the so called area enhancement argument [53,[57][58][59]. We will explain the this argument in the next section.

B. The area enhancement argument
One of the key points to evaluate the multiplicity for multi-particle production is the calculation of the average over charge color densities of 2n matrix amplitudes: where we have used the fact that the only part of the reduced amplitude that depends on the target charge density are the Wilson lines. For the sake of simplicity, let us just consider the case where we just have 4 Wilson lines in such a way that the color indices are contracted forming a single trace. In this case the object that we have to evaluate is the quadrupole operatorQ (q 1 , q 2 , q 3 , q 4 ) = y 1 y 2 y 3 y 4 e −iq 1 ·y 1 +iq 2 ·y 2 −iq 3 ·y 3 +iq 4 ·y 4 1 e −iq 1 ·y 1 +iq 2 ·y 2 −iq 3 ·y 3 +iq 4 ·y 4 Q(y 1 , y 2 , y 3 , y 4 ).
Following the arguments in [58,59], the configuration of the transverse coordinates y i that maximises the integral is such that these legs are as far away as possible between them. On the other hand, in the CGC picture the target ensemble is composed by domains of chromoelectric field with a typical correlation length, Q −1 s , that is fixed by the saturation scale, where color neutralises. Therefore two objects that only depend on the target color charge density, sitting at two different points y i and y j , will have a vanishing correlation when |y i − y j | Q −1 s . This implies that the only way of obtaining a non vanishing correlator is by grouping the legs in, at least, pairs where the distance between the transverse points is smaller than the correlation length. Thus, for the case of the quadrupole, this is equivalent to write The first term of this equation, although it gives a non vanishing contribution to the 4-point correlator, is constrained to a smaller region of phase space than the other 3 terms. This will imply that, after performing the integration in Eq. (12), it will be suppressed by the area of the target with respect to the other ones. On the other hand, the other three terms can be just written as a product of dipoles, that is and analogously to the other terms. Thus we can write the quadrupole operator as a sum of products of 2-point functions, keeping in mind that this approximation is only good after performing the phase space integral since, otherwise the first term in Eq. (13) is non-negligible. In Appendix A we discuss the validity of this argument, that we call area enhancement argument. This result can be generalised to the case of any number or configuration of the Wilson lines by noting that the contribution of the multipole that is enhanced by the area of the target, i.e., that is leading in S ⊥ Q −2 s with S ⊥ the area of the projectile (or the overlap area in a dilute-dense collision), is always a sum over all possible combinations of 2-point functions. This is analogous to assume that the target averages of Wilson lines follow a Gaussian statistics and thus we are able to apply Wick's theorem to them: being χ = {1, 2, . . . , 2n} and Π(χ) the set of partitions of χ with disjoint pairs. Eq. (16) simplifies enormously the evaluation of multipoles and shares its simplicity with the glasma graph approach through Wick's theorem. The main difference between them is that the first does not rely on the dilute limit and thus is applicable to dilute-dense scattering. This approach has been used recently [43,53] in order to evaluate the phase space integral of 4-point and 6-point functions. We will use it in order to evaluate Eq. (10).

C. Particle correlations
In this section we summarise the main ideas behind particle correlations within the CGC effective theory and provide the general formulae that we will employ to study azimuthal correlations. For a complete review of the former aspect, we refer to [35] and references therein.
Following the argument in [91,92], we make the picture of angular correlations in dilute-dense scatterings through the interactions of the projectile partons with the strong chromoelectric fields generated by the target. The strength of the chromoelectric field in the target wave function is characterised by the saturation momentum, Q s , which is also the typical momentum of partons in the wave function, k ⊥ ∼ Q s . The correlation length of the fields is roughly Q −1 s and the target ensemble can be modelled as a compound of domains with different chromoelectric fields that change from event to event as illustrated in Fig. 3. When a parton coming from the projectile wave function hits the target it will scatter in one of these domains and will pick a momentum that is proportional to the chromoelectric field inside this domain. Thus angular correlation appears when two partons scatter in the same chromoelectric domain. As gluons belong to a real representation of SU(3), scattering with parallel and antiparallel momenta is identical. Thus, this picture is also able to explain the absence of odd azimuthal correlations in gluon production. In order to study the azimuthal harmonics appearing in multi-particle correlation it is convenient to use the cumulant method [93]. This method aims to reduce the contribution of the so-called "non-flow" correlation, i.e., contributions to the correlation function that come from other processes other than true collective flow, such as resonance decays or jet correlations, to the definition of the azimuthal harmonics. In this method we define the 2-and 4-particle cumulants of order n as where · · · denotes the average over all events and particles. For convenience we define the n th -order κ-function in such a way that the event average can be written as Given this definition of the cumulants we can write the 2-and 4-particle Fourier harmonics of order n as Similarly, the harmonics v n can also be defined as a function of the transverse momentum. In order to do that we also define the so-called "differential" cumulants: where p ⊥ = |p| and we have defined the differential κ-functions as With this prescription, the differential azimuthal harmonics are given by 2

III. EVALUATING THE TARGET AND PROJECTILE CORRELATION FUNCTIONS
In this section we will introduce the notation and the arguments followed in order to evaluate the 2n-point correlation functions for both projectile and target ensembles.
A. Setting up the notation We write the reduced matrix amplitude as where i = 1, . . . , 2n. We should note, however, that in this notation when i is even the reduced matrix element is conjugate (i.e., to the right of the cut) and when it is odd it is not conjugate (i.e., to the left of the cut). Furthermore, since the produced gluon has the same momentum, polarisation and color both in the real and conjugate spaces we have to apply the following constraints: with m = 1, . . . , n. Thus this notation also changes the usual labelling of the gluon final momenta, k i , since now they are labeled by only odd numbers (1, 3, 5, . . . ) instead of 1, 2, 3, . . . . With these convention we can write the 2n-point function of the reduced matrix amplitudes in the simplified form As we have pointed out in Section II B, in this work we will use the area enhancement argument in order to evaluate the multipole correlators. Thus, using the same arguments that we have used for obtaining Eq. (16), we apply Wick's theorem and Eq. (32) reads with χ = {1, 2, . . . , 2n} and Π(χ) the set of partitions of χ with disjoint pairs. On the other hand, in order to evaluate the 2-point function we use Eq. (7), where we do not write the overall sign of the equation which should be −(−1) α+β since the number of real and complex conjugate matrix elements are the same and therefore the net sign of Eq. (33) will always be positive. Another simplification that we can make is by noting that the Lipatov vertices will always be contracted with the one that is evaluated at the same momentum k i . This follows from the fact that two gluons with the same transverse momentum will also have the same polarisation and thus the polarisation vectors fulfill which implies a contraction of the two Lipatov vertices with the same k-momentum. Thus we will write directly L λ (k, q) in Eq. (34) instead of i * λ (k)L i (k, q) because both expressions lead to the same result. On the other hand, we should also evaluate the average of two Wilson lines. In order to do so we follow [59] and use the fact that the target ensemble is globally color invariant, which implies that the average of any tensor in this ensemble has to be proportional to a linear combination of invariant tensors. Thus where we have introduced the dipole operator D(x, y). Therefore, making the change of variables y α,β = b ± r/2 we can write Eq. (34) as This equation can be simplified further if we exploit the fact that the target ensemble has a much larger extension in the transverse plane than the projectile one and then we assume translational invariance of the dipole operator, that is, D(r, b) = D(|r|). Thus, defining the Fourier transform of the dipole operator we obtain our final expression for the 2-point function of the reduced matrix amplitude: In order to obtain a final expression for Eq. (10) we should also evaluate the 2n-point function of the projectile color charge densities. In this case we will use the generalised MV model and also use the Wick's theorem. Introducing again the simplified notation we can write the 2n-point function as Here, as in Eq. (32), even indices correspond to complex conjugates. This correlator, Eq. (41), has the following Wick expansion: In the generalised MV model this 2-point function can be written as where µ 2 (k, q) is a function peaked around k + q = 0. In the strict MV model we have that µ 2 (k, q) ∝ δ (2) (k + q). All in all, using the area enhancement argument for computing the target correlator and the MV model for computing the projectile one, we arrive at the following general result for the multiplicity of n-gluon production: that, together with Eqs. (39) and (43), provides the full expression that will be used along this work.

B. Wick diagrams
Since the expression of Eq. (44) involves the product of two Wick expansions it includes the sum of (2n − 1)!! 2 products of 2n 2-point functions. Thus, when n > 2 we will have to deal with a large number of terms and, for this reason, it is convenient to introduce a shorthand notation for each of these objects involved in the sum. Therefore we introduce in this work a diagrammatic notation for each term inside the sum of Eq. (44) analogous to the diagrams introduced in [90] within the glasma graph approach. In our case, the diagrams consist of two parts that are separated by a vertical dashed line. In both parts we draw 2 rows and n columns of dots where the dots of the upper row are labelled by odd numbers and the ones of the lower row are labelled by even numbers, and the labels are the same in both sides: Each column of both parts of the diagram corresponds to a produced gluon. The columns defined by (1,2) corresponds to gluon 1, the ones defined by (3,4) to gluon 2 and so on. The upper row (odd indices) will represent the real space and the lower row (even indices) will represent the conjugate space. As we have said, each term of the sum of Eq. (44) will have a product of n 2-point functions coming from the projectile average and n 2-point functions coming from the target average that are labelled by 2 indices that goes from 1 to 2n. We will draw these 2-point functions as lines that connect the dots in the diagram. We choose the left part of the diagram to represent the 2-point correlators of the projectile and the right part to represent the 2-point correlators of the target and schematically what we will draw is the following 3 : As an illustrative example let us select one of the 5!! 2 = 225 terms that appear in Eq. (44) when n = 3: This term will be represented by the following diagram: Note that the integration over the q s is implicit.
If we want to write the diagram in an equation form we just have to use Eqs. (39) and (43). For example, the diagram in Eq. (48) reads (remember that we are labeling k i with odd indices, and Eqs. (29) to (31) with q ≡ d 2 q/(2π) 2 .
Besides making the notation more compact we can also exploit the structure of the diagrams in order to find symmetries between them, the associated power in (N 2 c − 1) for each diagram and which kind of quantum correlations (Bose enhancement or HBT) it includes, by making use of the following properties: i) Interchanging two dots within a column, 2m ↔ 2m − 1, of a given diagram is equivalent to make the change of variables k 2m−1 → −k 2m−1 .
In order to prove this property it is enough to evaluate Ω 2m Ω i i, j, γ and β arbitrary indices, since it is the only piece of Eq. (44) that depends on the dots 2m and 2m − 1. This expression can be computed using Eqs. (39) and (43). Then, if one makes the change of variables with unit Jacobian q 2m → −q 2m−1 and q 2m−1 → −q 2m and uses the fact that L λ (k, −q) = −L λ (−k, q) and ii) Interchanging two columns (2m, 2m − 1) and (2k, 2k − 1) of a given diagram is equivalent to make the change of variables k 2m−1 ↔ k 2k−1 .
The proof of this property is trivial since, by definition, each column of dots corresponds to a momentum k i and, therefore, interchanging two columns in both sides is equivalent to interchange the label of two momenta.
iii) We can extract the powers of (N 2 c − 1) by looking at the structure of each side of the diagram. Before making a statement of the property we will start by using Eq. (48) as an illustrative example. Using Eqs. (39) and (43) we can extract the counting in powers of (N 2 c − 1) −1 of this diagram by writing the Kronecker deltas where the second group of deltas of the first line is introduced to preserve Eq. (31), that is, that the color of the produced gluons is the same in the real and the conjugate spaces. The first group of deltas of both lines accounts to the target configuration (right side of the diagram) and the last group of deltas accounts to the projectile configuration. If we organize this equation in such a way that all the indices in the deltas are closed we have that where we have written the deltas that come from the target side of the diagram in a different color by convenience. We can do the same procedure that we did in the last equation in a diagrammatic and faster way by just drawing the target (right) side of the diagram on top of the left side and counting the number of closed lines that we obtain (which is equivalent to the second line of Eq. (54)) and drawing vertical lines in the right side of the diagram and counting the number of closed lines that we obtain (which is equivalent to the first line of Eq. (54)), where we can identify the red lines in the second diagram as the red Kronecker deltas of Eq. (54). In general, if we call n p the number of closed lines that we obtain by projecting the right side of the diagram on top of the left side and n T the number of closed lines that we obtain by projecting vertical lines on top of the right side of the diagrams the counting in powers of (N 2 c − 1) of a given diagram for general n is As we will see through this work, this property is useful for organising the terms of Eq. (44) in powers of (N 2 c − 1) −1 in a systematic way, especially when n is large.
iv) The types of quantum correlation that we have in a given diagram can be obtained as follows. If the same two dots are linked in both sides of the diagram we have two possibilities: if the dots belong to the same column labelled by (2k, 2k − 1) it means that the gluon k is uncorrelated (disconnected piece) and if the dots belong to different columns it means that the gluons that define these columns have an HBT correlation. By exclusion, all the gluons involved in other kind of links have a Bose enhancement correlation either in the projectile or in the target wave function.
In order to check this property it is enough to evaluate the terms in Eq. (44) that contain the same links in both the projectile and target sides. That is, we are interested in terms that contain where a, b = 1, ..., 2n are generic dots. The objects of Eq. (39) and Eq. (43) that contain the information of the quantum interference correlations are the Dirac deltas and the functions µ 2 (k, q) respectively (the Lipatov vertices and the dipole functions give a different kind of correlation). Thus, we can write Since µ 2 (k, q) is peaked around k = −q this implies that we have a peak around k a = −(−1) a+b k b which is an HBT correlation. In the case in which a and b belong to the same column, that is, a = 2k − 1 and b = 2k (or vice-versa), it is clear that we loose the correlation in function µ 2 -in fact we loose any kind of correlation since in this case the Lipatov vertices and the dipole function can be factorized.

IV. RESULTS
In this section we present the calculation of Eq. (44) for n = 2, 3 and 4. Larger values of n can be also considered in the same fashion, contingent upon sufficient computation power. In order to compute Eq. (44) we need Eqs. (39) and (43) which contain two functions that need to be modelled, µ 2 (k, q) and d(q).
As indicated before, in the strict MV model µ 2 (k, q) is proportional to a Dirac delta. However, in order to be more realistic, we choose a smoother function that is also peaked around k + q = 0, such as a Gaussian 4 : where B p is the gluonic transverse area of the projectile. For the dipole we use the Fourier transform of the GBW saturation model [60,61]: We should also account for the infrared divergences of the Lipatov vertices. The product of two Lipatov vertices is Usually these divergences are regulated by introducing an infrared cutoff both in all the integration over the momenta. However, in this work we use the following expression for the product of two Lipatov vertices: where ξ 2 is a parameter with dimensions of momentum squared. This choice, although it does not maintain some important properties of the Lipatov vertices, it is much simpler to deal with and, as we show in Appendix C, it is equivalent to using the Wigner function approach [44,45,64,65,94] but including quantum correlations in the projectile wave function. Thus, for two partons in the projectile the joint Wigner function that we use reads where one sees the uncorrelated term as the first one in the sum on the right hand side and the four color indices correspond to the four ρ's in the projectile average for the double inclusive gluon cross section. We note that the main problem of Eq. (63) is that it only depends on the momentum of the parent parton, k i − q i , and not on the final momentum, k i . Therefore, Eq. (63) only includes the contribution in which the gluon is emitted from the source and then interacts with the target, thus missing part of the physics. The final momentum is acquired by the interaction with the target which is suitable for the projectile collinear limit. In principle, in this limit the so-called "hybrid factorization" is employed and it corresponds to forward production of partons near the proton fragmentation region [95]. The approach that we adopt in this manuscript is suitable for central production even though the approximation used for the Lipatov vertices in Eq. (63) is more appropriate for considering the forward limit. Therefore, admittedly the validity of our approach is reduced to the forward region but not yet near the proton fragmentation one. In this region, the projectile partons are defined in terms of Wigner functions (see [44,45,64,65,94]). However, we would like to emphasize that the Wigner functions adopted in these references are factorized for two partons and do not include quantum correlations in the projectile. The two parton joint Wigner function (given in Eq. (64)) that we use in our approach indeed encodes the correlations in the projectile which is one of the novelties of the present manuscript 5 . Moreover, adopting Eq. (63) for the Lipatov vertices and Eq. (64) for the joint Wigner function to describe the projectile partons, allows us to perform the computation analytically until the very end, even though they restrict the validity region of our results. In our approach, one can generalise the computation to the production of any number of particles and can perform the study analytically within its limits of the validity. Other approaches that are strictly valid for central production, such as the study performed in [43] or the one in [45], rely on final numerical integrations which would be extremely difficult in the case of four particle correlations, or the computation is performed numerically from the very beginning making it difficult to control, respectively. Finally, due to the assumed Gaussian forms, our final expressions cannot be considered reliable for transverse momenta sizeably larger than the saturation scale.

A. Double inclusive gluon production
The case n = 2, that is, the spectrum for double inclusive gluon production is the most studied case. It was well described using the exact solution for the dipole correlators in the MV model [64,65], in the glasma graph approximation [46] and using the area enhancement argument [53,57]. The result that we present in this section is the same obtained in [53] but now, with the help of Eqs. (60), (61) and (63), we are able to obtain a closed-form solution for both the multiplicity and the azimuthal harmonics.
In this case, the expansion of Eq. (44) in terms of the Wick diagrams is where we have grouped the 9 diagrams by their powers in (N 2 c − 1) −1 . The Wick diagrams of Eq. (65) can be computed using Eqs. (39), (43), (60), (61) and (63) in a straightforward way since all the arguments of the q i integrals are Gaussian functions and, therefore, they can be trivially solved. The result is With these 5 equations we have fully determined the differential multiplicity in Eq. (65). In order to obtain the value of the integrated spectrum we just have to perform again Gaussian integrations over k i obtaining We can see from this equation that, apart from the suppression in powers of (N 2 c − 1) −1 , the correlated terms contain suppression factors (1 + B p Q 2 s ) −1 and (1 + B p ξ 2 ) −1 . Following the domain picture that we have discussed in Section II C, B p Q 2 s ≡ n D is the number of color domains in the overlap area of the projectile with the target in the transverse plane. We should expect decorrelation of the produced gluons in the limit of n D → ∞ since the probability of two gluons scattering off the same domain vanishes in this limit. Therefore, to fix ξ 2 it makes sense to choose a value that is proportional to Q 2 s in order to preserve decorrelation in the limit n D → ∞. For this reason we will choose ξ 2 = αQ 2 s , being α a real number, in the rest of this work 6 . The 2-particle azimuthal harmonics, Eq. (17), can be obtained by performing the integration over k i with the help of Eq. (B12). The result for the second order κ-function is where we have defined α = ξ 2 /Q 2 s , we have taken n > 0 and due to the symmetry k 3 → −k 3 of Eq. (65) all odd harmonics vanish. Using Eqs. (71) and (72) we can evaluate the 2-particle azimuthal harmonics as In Fig. 4 we plot the dependence of v 2n {2} with respect to n D and α by fixing N c = 3. The value of the even azimuthal harmonics grows rapidly as both n D and α approach zero and it decreases slowly when these parameters are large. This decrease with n D is what we should expect in the color domain picture of particle correlation since as n D gets larger the probability of two gluons scattering in the same domain is smaller and thus the overall correlation. On the other hand, the decrease with α must be taken with care because α gives the ratio between the momentum transfers from projectile and target. The dilute-dense approximation that we are using makes sense only for α sizeably smaller than 1.
As we have seen in Section II C, we can also compute the azimuthal harmonics as a function of transverse momentum by using the differential κ-function defined in Eq. (25). The k i integral can be solved with the help of Eq. (B11) and the result is and (1 + αn D + α 2 n D ) 2 α(1 + α)(1 + α 2 n D + α(2 + n D )) n 6 Since ξ 2 is a momentum scale of the projectile wave function it should be related with B −1 p and not with Q 2 s which is a momentum scale of the target wave function. However, the choice ξ 2 = Q 2 s is the one that has given more consistent phenomenological results and for this reason we use it through all this work. In [64,94] and [65] the choices ξ 2 = B −1 p and ξ = Qs/4 have been made, respectively, and the sensitivity of the results to variations of these choices has been examined.
where n > 0. The fact thatκ 2n (p ⊥ ) is proportional to (p 2 ⊥ /Q 2 s ) n was also obtained in [65] although there a different model for the target average was employed. The differential 2-particle even azimuthal harmonics can be obtained by evaluating and the result is plotted in Fig. 5 for n = 1, 2, 3 and B p = 6 GeV −2 , ξ = Q s /2, Q 2 s = 2 GeV 2 and N c = 3. Although we do not aim for a comparison with experimental data, the obtained values are in the ballpark of them. Note that due to the Gaussian forms that we employ, our results cannot be considered reliable for p ⊥ sizeably larger than Q s .

B. Triple inclusive gluon production
In this section we show the result for Eq. (44) when n = 3, that is, the triple inclusive gluon spectrum. Since in this work we are mainly interested in computing azimuthal harmonics we will just show the expansion of the spectrum in terms of the Wick diagrams. However, it has be shown recently [43] that this result is useful for studying the correlation between the 2-particle azimuthal harmonics and multiplicity and average transverse momentum. As we did for n = 2, we can group the Wick diagrams in the expression for the n = 3 gluon spectrum in powers of (N 2 c − 1) −1 as 3 .
In order to obtain each one of these terms we have to use the property iii) of Section III B. In this case the suppression of each diagram is given by (N 2 c − 1) np+n T −6 and we can have three kind of configurations on each side of the diagram , It is easy to see that the only way of obtaining n T = 3 ,2 and 1 is having the first, second and third configuration on the right side of the diagram, Eqs. (78) to (80) respectively. On the other hand, the only way of obtaining n p = 3, 2 and 1 is having the same configuration on the left side of the diagram as the one on the right side, a configuration on the left side that has one link equal to the configuration on the right side and the other 2 links different, and at configuration on the left side that has all the links different than the configuration on the right side, respectively. One can also check that for a given configuration on the right side of the diagram the number of possibilities for n p = 3 is 1, for n p = 2 is 6 and for n p = 1 is 8. With this taken into account, let us show as an example how to find all the Wick diagrams suppressed by (N 2 c −1) −2 . In this case n p + n T = 4. There are three possibilities: i) n p = 1 and n T = 3.
This implies that we have to have the configuration of Eq. (78) on the right side and configurations on the left side that has all the links different than the one on the right side. As we have said, there are 8 possibilities for this case: ii) n p = 2 and n T = 2.
This implies that we have to have the configuration of Eq. (79)  All in all, we can write all the 52 Wick diagrams that have a suppression of (N 2 This procedure, although tedious, is straightforward to implement on a computer. Repeating it we find that there is 1 diagram suppressed by (N 2 c − 1) 0 : 12 diagrams suppressed by (N 2 c − 1) −1 : 96 diagrams suppressed by (N 2 c − 1) −3 : All in all, we get the 225 = (5!!) 2 . Eqs. (81) to (85) give the full Wick expansion of the triple inclusive gluon spectrum to all orders in (N 2 c − 1) −1 . These equations were computed in [53] up to order (N 2 c − 1) −2 . In order to write Eq. (77) just as a function of k i one just has to employ Eqs. (39), (43), (60), (61) and (63) and then perform the q i integrals. These integrals are straightforward if the arguments of the integrals are Gaussian functions, as we assumed before.

C. Four gluon inclusive production
In this section we evaluate Eq. (44) for n = 4, that is, the four gluon inclusive spectrum. Since in this case the number of diagrams involved is very large ([7!!] 2 = 11025) and thus their writing is not viable, we will start by discussing the simpler case in which the partons in the wave function of the projectile are initially not correlated. Then we will consider to the more general case that is explained in detail in Appendix D.
The case in which the partons in the projectile wave function are initially uncorrelated was discussed for scattering quarks [45], and for gluons [55] within the glasma graph approach. In our case, this implies writing Eq. (44) as and, instead of having [(2n − 1)!!] 2 terms in the sum, we just have (2n − 1)!!. Using Eq. (86) we can write the 4-gluon inclusive production as the sum of 105 diagrams (also known as rainbow diagrams [90]) in the following form: In this expression, the term in the first line corresponds to the case in which all the generated gluons are uncorrelated . The 12 terms in the second line correspond to the case in which 2 gluons are uncorrelated and 2 gluons are correlated. The 32 terms in the third line correspond to the case in which 1 gluon is uncorrelated and the remaining 3 ones are correlated. The 12 terms of the fourth line correspond to the case in which two pair of gluons are correlated independently, i.e., factorisable connected diagrams. Finally, the 48 terms of the last lines (the factor 1/2 avoids double counting of the diagrams) correspond to the case in which all the gluons are correlated between them, i.e., fully connected diagrams. Note that the first, second, third and fourth, and fifth terms in the sum on the right hand side correspond to terms with increasing powers in (N 2 c − 1) −2 . The Wick diagrams of Eq. (87) can be computed in the same fashion as in Section IV A. However, since we are only interested in computing the 4-particle cumulants, Eq. (18), we will exploit the k i ↔ k j and k i → −k i symmetries in order to simplify the calculation. When evaluating the 4-particle κ-function in Eq. (19) all the terms that contain at least one disconnected piece, i.e., two vertical lines in both sides of the diagram, will vanish trivially due to rotational invariance. For this reason the diagrams of the first three lines of Eq. (87) will not contribute to the 4-particle κ-function when n > 0 and therefore we can write where we have written schematically perm 4 , which encodes all the factorizable connected diagrams, and perm 5 , which includes all the fully connected diagrams, as all the permutations of the fourth and fifth lines of Eq. (87) respectively. Note that we have also dropped the factors 2 and 2π as they cancel in Eqs. (17), (18), (23) and (24). On the other hand, we can read from the permutations for the fully connected diagrams perm 5 , that all the diagrams that are related by a change of variables k i → −k i will give the same result for the integral in Eq. (88) since this transformation is equivalent to making φ i → φ i + π in the argument of the exponential and, thus, leaves the integral invariant. Furthermore, it is easy to check that all the diagrams of the last three lines of Eq. (87) that are related by the change of variables k 1 ↔ k 3 , k 5 ↔ k 7 and k 3 ↔ k 7 also give the same value for the integral. By inspection of the symmetries, which is detailed in Appendix D, we can write the 48 integrals defined by the permutations of the last three lines of Eq. (87) as where the last diagram can be seen as the first one with the change of variables k 1 ↔ k 7 . Furthermore, 4 out of the 12 diagrams of the fourth term of Eq. (87) only depend on φ 1 − φ 3 and φ 5 − φ 7 and therefore vanish due to rotational invariance. Having this into account we can write Eq. (88) as On the other hand, the 2-particle κ-function in the case in which the partons are initially uncorrelated in the projectile wave function is Therefore we can write Eq. (91) as In order to compute κ 0 {4} we have to have into account all the diagrams of Eq. (87). However, since all the permutations are related by the change of variables k i → −k i or k i ↔ k j (i = j) that leave the integral invariant we can write This integral can be easily performed since all the terms are just Gaussian functions. All in all, the 4-particle cumulant can be computed by using Eq. (18) and Eqs. (93) and (94) and the four particle even azimuthal harmonics is obtained as In  6: Dependence of the 4-particle integrated cumulants on Q 2 s (left) and of the differential cumulants on p ⊥ (right) in the case in which the partons in the projectile wave function are uncorrelated. In these graphs we have used N c = 3 and the values of the remaining parameters are indicated on the plots. The differential 4-particle cumulant in Eq. (24) can be computed in the same fashion but since we are fixing one of the momentums we have to be more careful with the symmetries discussed in the last paragraphs. In Appendix D we show that we can write (again dropping factors 2 and 2π) and, therefore, In order to computeκ 0 {4}(p ⊥ ) we cannot use the same symmetries that we employed for computing κ 0 {4} because now one of the momenta is fixed. All the Wick diagrams of Eq. (87) that are related by a change of variables k i → −k i still leave the integral invariant but now, since we are fixing |k 1 | = p ⊥ , all the diagrams that are related by a change of variable k 1 ↔ k j will give a different value for the integral but the ones that are related by k i ↔ k j , with i, j = 1, still leave the integral invariant. Performing a simple counting of the permutations of Eq. (87) we can writẽ Thus, using Eqs. (98) and (99) the differential 4-particle cumulant can be written in a similar form as Eq. (95): and the differential 4-particle azimuthal harmonics is defined as In Fig. 6 right we have plotted our result for Eq. (100). Again, the values are very small and they even become positive with increasing p ⊥ because we are not including the diagrams that take into account the correlation of the partons inside the projectile.
With the results of Fig. 6 we have finished our discussion of 4-gluon production in the case in which the partons are not correlated in the projectile wave function. So far, let us recapitulate what we did in this section. First, we wrote the 4-gluon spectrum in terms of the Wick diagrams by classifying them in different topologies and, thus, with a different suppression in powers of (N 2 c − 1) −1 . Then we wrote the diagrams with the same topology as just one plus a bunch of permutations, as in Eq. (87). Then we exploited the symmetries of these permutations in order to reduce the number of integrals to be performed in the 4-particle cumulant functions Eqs. (95) and (100). We also noticed that the contribution of the non vanishing factorizable connected diagrams to κ n {4} can be written as 2κ n {2} 2 . Finally, we solved numerically these integrals for given values of Q 2 s and B p . Now let us jump to the case in which we take into account all the terms of the Wick expansion of the projectile correlator. In this case we have to deal with (7!!) 2 = 11025 terms instead of 7!! = 105. While the calculation becomes more cumbersome, the approach is exactly the same. First, we group all the Wick diagrams in the 4-gluon spectrum by their topology that defines the power in (N 2 c − 1) −1 , by using the property iii) of Section III B. Then, we relate the diagrams with the same topology by permutations. Next, in order to compute the 4-particle cumulant we exploit the symmetries of these permutations and reduce as much as possible the number of integrals to be performed. Finally, we solve numerically each one of these integrals and obtain a result for the azimuthal harmonics. The detailed discussion of this procedure can be found in Appendix D.
In Figs. 7 and 8 we show our results for the four gluon cumulants as a function of Q 2 s , the differential cumulants as a function of p ⊥ and the corresponding azimuthal harmonics for n = 2 and 4. We use the same parameters that we employed in the two gluon case. Now, in contrast with the case seen above, the values obtained are negative (for the cumulants, thus real for the Fourier coefficients), larger in absolute value and in the ballpark of experimental data. Monte Carlo integration is used, yielding negligible errors except for the smallest p ⊥ for d 4 {4}. On the other hand, it is known that when the multiplicity gets low the 4-particle cumulant turns positive [10,15]. The naive assumption that the multiplicity is proportional to the saturation momentum suggests a change of sign in the cumulant as Q 2 s → 0. Indeed, in the glasma graph approach, suitable for dilute-dilute collisions and therefore for lower multiplicities, arguments [56] suggested that c 2 {4} > 0 -a result also found in [45] where a transition from positive to negative is found when multiple scattering (that goes beyond glasma graphs) is introduced. This is not seen in Fig. 7. A more detailed calculation should be done in this regime of low multiplicities where the transition for the glasma graph approach is expected. FIG. 8: Dependence of the differential 4-particle cumulant (left) and azimuthal harmonic (right) of second and fourth order with p ⊥ . For the latter we also show the results obtained from 2-particle correlations. In these graphs we have used N c = 3 and the values of the remaining parameters are indicated on the plots.
With the results for the azimuthal harmonics in the case of 4-gluon inclusive production we finish our discussion on multi-gluon production. We point out that the procedure that we have developed can be generalised for larger values of n. It implies dealing with a large number of diagrams ([(2n − 1)!!] 2 ). There is not conceptual problem for doing it since, as we have shown, we can always use the property iii) of Section III B to group all the diagrams in a systematic way and then exploit the symmetries to reduce the number of integrals to be performed. The remaining issue is dealing with a large number of 2n-dimensional integrals that must be solved numerically.
We should also note that, although the results shown in this section are consistent with experimental data, no attempt is done to compare with them. We have used the area enhancement argument that should only be valid in the case when the overlap area is large. Furthermore, we have only taken into account the scattering of gluons. For more realistic results, we should at least compute the differential multiplicities for scattering quarks, consider more involved projectile and target averages (e.g. fluctuations) and convolute the results with fragmentation functions.

V. SUMMARY
In this work we have computed multi-gluon production in the CGC in dilute-dense (pA) collisions, extending the work in [53] to four gluons. Our calculation includes the contributions that are leading in the overlap area of the collision [53,[57][58][59], while keeping all orders in the expansion in the number of colors. We develop a diagrammatic technique to write the numerous color contractions and exploit the symmetries to group the diagrams and simplify the expressions. This technique reduces dramatically the number of integrals needed to compute the multiplicity distributions and integrated and differential cumulants, which results essential for the large number of diagrams, more than 10000, that appears for four gluon production. We use the GBW model [60,61] for the dipoles that result from the target averages, and the generalised MV model [62,63] for projectile averages. In order to proceed analytically as far as possible and simplify the final calculations, we use the Wigner function approach [45,64,65] that we extend to include quantum correlations in the projectile wave function. The Wigner function approach supposes that the final momenta of gluons is mainly acquired through interaction with the dense target and is thus suitable for a collinear projectile approximation.
Apart from the techniques developed and the discussions on the validity of the area enhancement argument and the Wigner function approach, our main results can be summarised in Figs. 4, 5, 7 and 8. For two gluon correlations, we provide analytic expressions for integrated and differential cumulants which show smooth dependences on the parameters defining the projectile and target Wigner function and dipole, respectively. For four gluon correlations we find that the second order four particle cumulant c 2 {4} < 0 -thus providing a sensible second order Fourier coefficient v 2 {4}, a result found in [45] (where only quark scattering is considered and partons in the projectile wave function are uncorrelated) and attributed to multiple scattering. We note that the approximation in which gluons in the projectile are uncorrelated gives results for the cumulants that are much smaller in absolute value than when correlations are included, and become positive for some values of Q s and p ⊥ . This emphasises the importance of including the full correlations in the projectile.
Our numerical results, due to the Gaussian forms that we employ for the Wigner function and dipole, cannot be considered reliable for p ⊥ sizeably larger than Q s . They lie in the ballpark of experimental data, for values of parameters that look reasonable. But we are aware that further analytic understanding is still required, and several pieces are still missing in our formalism: the contribution from quarks, more involved projectile and target averages, fragmentation functions,. . . . All these aspects should be explored before we can establish a model ready for phenomenology.
An immediate outlook of this work that we plan to address in the near future, is exploring the transition to low multiplicities, where the target should behave as a dilute object and the glasma graph approach should be valid. It has been argued that in such approximation c 2 {4} > 0 [45,56]. It would be most interesting to clarify the origin of such change of behaviour observed in data [10,15] and implement a framework that consistently goes from the dilute-dilute to the dilute-dense situation to examine the behaviour of the many particle cumulants. In this section we will study the validity of the area enhancement argument, from now on AE model, introduced in Sect. II B. For the sake of simplicity we will work in the fundamental representation of the Wilson lines instead of in the adjoint representation. Furthermore, we will only consider the expectation value of 4 Wilson lines, i.e. double quark interaction with a target. We expect the discussion presented here to be also valid for any number of Wilson lines or for a different color representation. In our case, the expected value of the double dipole operator is where we have introduced the dipole operator D(x, y) = 1 Nc Tr[U (x)U † (y)] and D(x, y) is its target average. As discussed in Sect. II B, this approximation is only valid after integration over the phase space and at leading order in the transverse size of the interaction region, B p . In order to check the validity of this model we will compare it with the result of [81,96] that was obtained by assuming multiple coherent scatterings of the quarks within the target. This result was obtained using the MV model and the result is where and function F (x, y; u, v) is defined in [81]. In the GBW model this function reads simply Taking the large-N c limit we can simplify Eq. (A2) drastically to read, in the GBW model, Thus, in the large-N c limit, the ratio between the integral of the double dipole weighted by an arbitrary smooth function of the coordinates Ψ(x, y, u, v) computed in the MV and AE models is Using the saddle point approximation, and noting that F (x, y; u, v) → 0 when x → y or u → v and the fact that the dipole functions are Gaussian functions, it is straightforward to see that, in this approximation, the MV and AE model lead to the same result. For such approximation to hold we must consider the Gaussian functions, with width ∝ 1/Q s , to behave δ-like with respect to the integration area. Therefore, corrections must be order 1/Q 2 s that, by dimensional reasons, has to be multiplied by an inverse area, with the overlap area, i.e., the size of the proton B p , being the only parameter with such dimensions.
So far the discussion in this section only relies on the dynamics of the target and for this reason B p does not appear in the expressions. We will introduce it by defining the phase space measure as that is, we integrate over a 4-sphere of radius 2B p in such a way that the integral over the phase space leads to the expected result In order to compare the MV and AE models, we perform a Fourier transform over the phase space measure defined as In Fig. 9 we show the result for the ratio of the Fourier transforms of Eq. (A2) and Eq. (A1) for different values of B p , taking Q 2 s = 1 GeV 2 . The result was generated by using four sets of random momenta with moduli between 0.5 and 1.5 GeV. We see that, as expected, as we increase the value of B p Q 2 s the results in the AE model tends to those in the MV model, being the difference between both approaches of order a few % at relatively high B p . An analogous discussion can be performed to the expected value of the quadrupole operator. In the AE model it can be written as where . In the MV model it reads [96] Q(x, y, u, v) MV It turns out that the only difference between the expectation value of the quadrupole and double dipole operators in both models is a factor 1/N 2 c . In the large-N c limit, Eq. (A11) can be simplified to read We can show again that by using the saddle point approximation and the same arguments given below Eq. (A6) that the expected value of the quadrupole is the same in both models. However, in the case of the quadrupole the difference between Eqs. (A10) and (A12) is not suppressed by any power of 1/N 2 c and, therefore, we expect a larger discrepancy between both models.
In Fig. 10 we plot the ratio between the Fourier transform, defined analogous to Eq. (A9), of Eqs. (A10) and (A11) for three sets of random momenta with moduli between 0.5 and 1.5 GeV for different values of B p , taking again Q 2 s = 1 GeV 2 . In this case the difference between both models is of order 30% at relatively high B p , being larger than in Fig. 9 due to the 1/N 2 c suppression present for the double dipole and absent for the quadrupole. It also looks that at high B p the AE model tends to the MV model but integrals become very time consuming which prevents reaching larger values of B p Q 2 s . Therefore, the tendency is not as clear as in Fig. 9. The first integral and the one that will be the most used in this work is the well known Gaussian integral When computing the 2-particle differential cumulant we will have to deal with the following integral where we have defined B = −A 1 p 2 ⊥ and C = A 12 p ⊥ . Taking into account the Jacobi-Anger expansion Eq. (B5) and the fact that n is an integer we can write Using the definition of Eq. (B1) and then making the change of variable x = A 2 k 2 2 we have that where in the last line we have used the definition of the Gamma function Eq. (B4). Using the definition of the hypergeometric function Eq. (B2) we can write the last sum as ∞ k=0 C 2 4A 2 k 1 k!Γ(k + 2n + 1) Γ(k + n + 1) = Γ(n + 1) Γ(2n + 1) 1 F 1 n + 1; 2n + 1; and, therefore, the result for the integral is n Γ(n + 1) Γ(2n + 1) 1 F 1 n + 1; 2n + 1; The solution of the integrals that we will find when we evaluate the 2-particle cumulant can be obtained in the same fashion and the result is 4A 1 A 2 n Γ(n + 1) 2 Γ(2n + 1) 2 F 1 n + 1, n + 1; 2n + 1; The Wigner function approach was used in several works [44,45,64,65,94] in order to compute multi-particle production. Here we will follow the arguments in [64]. The forward amplitude for a gluon with momentum p scattering on a dense target and leaving with a momentum k is given, at leading order, by where for simplicity we are not taking into account the longitudinal polarization of the gluons. On the other hand the distribution of gluons with momentum k and color a coming from the projectile after the interaction with the target can be written as whereρ is the single gluon density matrix.
Using the completeness relation for the initial state we can write this equation as Doing the change of variables q 1,2 = p ± q/2 and x,x = b ± r/2 we get and realizing that is the Weyl transform of the density matrix, that is, the Wigner function, we can write the single inclusive gluon spectrum as Since this expression is also dependent of the color charge density of the target we still have to perform the target average. On the other hand, in the approach used in this work, we can evaluate the single gluon spectrum, before target averaging, using Eq. (10): Doing the same change of variables that we did before can write this expression as Thus, comparing Eqs. (C7) and (C9) we see that the single particle Wigner function can be written in terms of the Lipatov vertices and the 2-point correlator of the projectile charge density: Using the models employed through this work, Eqs. (60) and (63), we can write the single particle Wigner function as which is the same function found in the literature [64]. We can also check that this function is well normalised by performing the trace and integrating over b and p, Doing an analogous discussion we can write the 2-particle Wigner function as Performing the Wick expansion of the projectile correlator and using again Eqs. (60) and (63) we obtain Eq. (64). We can check that the quantity defined in that equation is not well normalised: Therefore, in order to have a proper definition of the Wigner function we should normalise Eq. (64) by this factor. However, since in correlation studies the overall constants do not contribute to the cumulants, this normalisation factor is not important for us. We should also note that Eq. (64) breaks the factorisation assumption that is used in the literature in which the 2-particle Wigner function factorizes into a product of two single particle Wigner function: The reason for the breaking of this factorisation is that we are including in our approach quantum correlations in the projectile wave function. Thus we can interpret the terms in Eq. (64) that break factorisation as Bose enhancement contributions in the projectile wave function.
where the permutations k i → −k i and k i ↔ k j are an abuse of notation since we have not contracted the diagrams and thus we cannot apply properties i) and ii) yet. In order to make the notation lighter we write these permutations as , Now we generate the Wick diagrams in such a way that they are grouped by their powers of (N 2 c − 1) −1 . In order to do so we exploit property iii) of Section III B. In this case the suppression of a given diagram is given by (N 2 c − 1) np+n T −8 , with the values of n p and n T fixing the topology of the diagram. It is straightforward to realise that all the diagrams with n T = 4 will have the configuration of Eq. (D1) on the right side. All the diagrams with n T = 3 will have one of the 12 configurations of Eq. (D2) on the right side. All the diagrams with n T = 2 will have one of the 32 configurations of Eq. (D3) or one of the 12 configurations of Eq. (D4) on the right side and all the diagrams with n T = 1 will have one of the 48 configurations of Eq. (D5) on the right side. Therefore the value of n T is fixed by the configuration that we have on the right side of the diagram.
The value of n p , on the other hand, will depend on the configuration that we have on both sides. It is determined by the number of disconnected pieces that we obtain after drawing the right configuration of the diagram on top of the left one. Thus, the only way of obtaining n p = 4 is by having a configuration on the left that has the same links as the one on the right. The only way of obtaining n p = 3 is by having a configuration on the left that has just two links that are equal to the ones on the right. The way of obtaining n p = 2 is by having a configuration on the left that has only one link that is equal to the right configuration or by having all the links different but in such a way that, after the projection, we obtain two disconnected pieces. Finally, the only way of obtaining n p = 1 is by having a configuration on the left side that has all the links different to the right configuration in such a way that, after the projection, we have a fully connected piece. The number of possibilities for n p = 4 is 1, for n p = 3 is 12, for n p = 2 is 32 and 12, respectively, and for n p = 1 is 48.
Having this into account we can find all the diagrams with a given suppression in powers of (N 2 c − 1) −1 . As an example, let us see which are the diagrams with power suppression (N 2 c − 1) −3 . In this case we have n p + n T = 5 and we will have 4 different topologies that are fixed by this constraint: n T = 4 and n p = 1 ; n T = 3 and n p = 2 ; n T = 2 and n p = 3 ; n T = 1 and n p = 4. Let us study this situation case by case: i) n T = 4 and n p = 1. In this case we will have the configuration of Eq. (D1) on the right side of the diagram and on the left side we will have all the diagrams that have zero links in common with the one on the right in such a way that after the projection we just have one connected piece.
where in the last line we have used the fact that since on the right side of the diagram we have a fully disconnected piece we can write the sum of these 48 diagrams as just one plus perm 5 .
ii) n T = 3 and n p = 2. In this case we will have the configurations of Eq. (D2) on the right side of the diagram and on the left side we will have all the diagrams that have just one link in common with the right one. This gives a total of 12 × 32 possibilities or we can have on the left side the configurations that have no links in common with the one on the right in such a way that, after the projection, we have two connected pieces. This gives a total of 12 × 12 possibilities: iii) n T = 2 and n p = 3. This implies that we will have the configurations of Eq. (D3) or Eq. (D4) on the right side of the diagram and on the left side we will have the configurations that have two links equal to the one on the right. This gives a total of 32 × 12 possibilities for the first case, iv) n T = 1 and n p = 4. This implies that we have the configurations of Eq. (D5) on the right side of the diagram and on the left side we have the configuration that have the same links with respect to the right one. This gives a total of 48 × 1 possibilities: With Eqs. (D11) to (D16) we have found all the 1152 diagrams with a suppression of (N 2 c − 1) −3 and written them as a bunch of diagrams plus permutations -which was our goal. We should also note that some of the diagrams that are drawn in these equations can also be related by symmetries k i → −k i or k i ↔ k j , which could lead to a better optimisation of the calculation but we have not found any systematic way of finding these symmetries. Therefore, we have decided to not include them in the calculation since we see not advantage in doing this by hand.
We can find the other diagrams with a different suppression in the same fashion obtaining 1 diagram with a suppression of (N 2 c −1) 0 , 24 diagrams with a suppression of (N 2 c −1) −1 , 232 diagrams with a suppression of (N 2 c −1) −2 , 3088 diagrams with a suppression of (N 2 c − 1) −4 , 4224 diagrams with a suppression of (N 2 c − 1) −5 and 2304 diagrams with a suppression of (N 2 c − 1) −6 . The next step is to exploit the symmetries encoded within the permutations in order to evaluate the cumulants through Eq. (19). We will do as an example the calculation only for the terms that contribute, again, with a power (N 2 c − 1) −3 . Let us introduce the shorthand notationD np as the sum of all the diagrams that satisfy the topology given by n p with a given configuration on the right side. Then we can write the contribution of order (N 2 c − 1) −3 to the 4-gluon spectrum as 2 referring to the first and second contributions to Item ii) discussed above, respectively. In order to evaluate κ 0 {4} we can use the fact that all the permutations, perm i , of Eqs. (D1) to (D2) will give the same result since we are integrating over all the momentum k i . Thus we can write When we evaluate κ n {4} with n = 0 we have to integrate the spectrum times e in(φ1+φ3−φ5−φ7) which will break some of the symmetries encoded in the permutations, perm i . In order to check how we can simplify the integration let us start with the permutations of the n T = 3 case with a generic n p . In this case we can define the sum of the diagrams without the permutations asD np ≡ f 2 (k 1 , k 3 , k 5 , k 7 ). (D19) By using the properties i) and ii) of Section III B we can check that this sum has the following symmetry f 2 (k 1 , k 3 , k 5 , k 7 ) = f 2 (k 1 , k 3 , k 7 , k 5 ). (D20) Thus, the contribution of the n T = 3 diagrams to the κ-function can be written as k1k3k5k7 e in(φ1+φ3−φ5−φ7) D np + perm 2 = 2 k1k3k5k7 e in(φ1+φ3−φ5−φ7) f 2 (k 1 , k 3 , k 5 , k 7 ) + f 2 (k 5 , k 3 , k 1 , k 7 ) + f 2 (k 7 , k 3 , k 5 , k 1 ) +f 2 (k 1 , k 5 , k 3 , k 7 ) + f 2 (k 1 , k 7 , k 5 , k 3 ) + f 2 (k 7 , k 5 , k 3 , k 1 ) = 2 k1k3k5k7 e in(φ1+φ3−φ5−φ7) 2f 2 (k 1 , k 3 , k 5 , k 7 ) + 4f 2 (k 1 , k 5 , k 3 , k 7 ) , where the factor 2 comes from exploiting the symmetries k i → −k i that are in Eq. (D2) and in the last equality we have used Eq. (D20) and relabelled the variables. Therefore, we can write k1k3k5k7 e in(φ1+φ3−φ5−φ7) D np + perm 2 = k1k3k5k7 e in(φ1+φ3−φ5−φ7) 4D np + 8D np .
For the n T = 2 diagrams that are defined by Eq. (D3) we can follow the same arguments by defininĝ D np ≡ f 3 (k 1 , k 3 , k 5 , k 7 ).
We can check that this function has the following symmetries