Double and triple inclusive gluon production at mid rapidity: quantum interference in p-A scattering

We compute double and triple inclusive gluon production in p-A scattering beyond the so-called “glasma graph” approximation. We consider quantum interference effects and identify in this general setup the terms responsible for the gluon HBT and initial wave function Bose enhancement which lead to correlations in particle production. Both of these terms originate from the factorizable part of the quadrupole and sextupole terms in the production cross section. We also show that the target Bose enhancement in this regime is suppressed at large number of colors.


Introduction
The interest in correlated particle production in the recent years has been triggered by the observation at the Large Hadron Collider (LHC) of the so called ridge correlations in p-p and p-Pb scattering [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Many features of the ridge correlations are shared in these reactions with similar observations in Pb-Pb collisions at the LHC and earlier observations of the same structure in Au-Au collisions at the Relativistic Heavy Ion Collider, later observed in d-Au and 3 He-Au [15][16][17][18][19]. In heavy ion collisions the accepted explanation of the origin of correlations is due to a collective behavior of the final state of the dense system of gluons, which follows a hydrodynamic evolution starting a short time after scattering. A similar explanation has been put forward for the observed ridge in p-A and p-p as well [20][21][22]. Nevertheless, the question remains of whether the mechanism that leads to these correlations in small systems (p-p and p-A) is of the same origin.
It has been suggested that the structure of the wave function of the highly energetic proton can be nontrivial and cona e-mail: tolga.altinoluk@ncbj.gov.pl tain preexisting correlations which are reflected in the final state of the scattering [23,24]. Motivated by this idea, a reasonable description of large parts of the data has been provided by calculations based on the Color Glass Condensate (CGC) approach in [25][26][27][28]. This work used the so called "glasma graphs" approximation which is based essentially on the dilute-dilute limit of the CGC framework [23,24].
The physics of "glasma graph" approximation used in [25][26][27][28] was elucidated later in [29,30]. There it was shown that the physical origin of the correlations in the "glasma graph" approach where the quantum interference effects that unavoidably appear in a system of identical bosons, i.e. gluons. The two distinct contributing quantum effects identified in [29,30] were the Bose enhancement of gluons in the incoming projectile wave function and the Hanbury Brown-Twiss (HBT) interference effect in the emission of gluons in the final state. Diagrammatics of the HBT effect in QCD was discussed in [31], and in the CGC framework more recently in [32,33]. Similar effects (albeit with opposite sign) have been later identified in the quark production in the CGC based framework in [34], see also [35,36]. The glasma graph approach was also used in [37,38] to study triple and quadruple gluon correlations. However, these studies were performed in the dilute-dilute limit of glasma graphs which should be only applicable for p-p collisions at midrapidity, while we aim to go beyond this limit and consider p-A.
We mention another recent paper [39] where quantum interference effects for gluons were studied in a model calculation from the point of view of multi parton interactions. More recently a calculation of quantum interference effects in forward quark production in the framework of multi parton interactions has been performed in [40,41]. There it was shown that these effects are ubiquitous in identical particle production, and that in the framework of the eikonal multiple scattering approach are contained in the quadrupole amplitude.
It is desirable to have a similar detailed understanding of quantum interference effects in gluon production at mid rapidity in the framework of the unabridged CGC formulation, i.e., beyond the "glasma graph" approximation. 1 The purpose of the present paper is precisely to provide such a calculation for two and three particle inclusive production. We perform it within the dilute-dense CGC approach which is a fully consistent scheme for calculation of particle production in p-A collisions. In this respect we follow the same line as [54] and [32,33]. We generalize the recent work along the same lines [55] (which is based on the approach of [32,33]) by identifying more generally the origin of the various terms in the double inclusive gluon production as well as calculating quantum interference effects in three gluon production for the first time.
In order to identify various physical effects, in the present paper we adapt the approach of [40,41] to production of particles in the adjoint representation. We use the insight of [32,33,40,41] regarding the approximation of the quadrupole amplitude in terms of dipoles and explain similarly to [40,41] that this approximation correctly accounts for the quantum interference effects. This considerably simplifies the color algebra, which in general leads to several different Wilson line ensembles [32,33,56] whose target averages have to be modelled. It also allows us to identify in full generality the Bose enhancement (both in the projectile and the target wave function) as well as HBT contributions to particle production within this general framework, besides providing results beyond any approximation of large number of colors N c .
We show that just like for fundamentally charged particles, the effects due to the quantum statistics that are encoded in the quadrupole amplitude are leading correlation effects at large N c . They contribute to the correlated production at order 1/N 2 c , while similar terms that arise from the dipole squared (no color exchange) term are suppressed as 1/N 4 c . We also observe that in the dilute-dense framework the terms in the production amplitude that arise due to Bose enhancement in the target are suppressed by 1/N 2 c relative to the Bose enhancement terms in the projectile. This N c counting is quite different from that in the glasma graph approximation, where both effects are of the same order in 1/N c . 1 Note that several other attempts exist in the literature to go beyond glasma graphs in multi particle production [42][43][44]. These works however do not include the quantum interference effects, but rather deal with a more careful evaluation of non factorizable contributions to products of dipoles within the McLerran-Venugopalan model [45,46]. There is also a body of literature devoted to the evaluation of higher multipoles within the same framework [47][48][49][50][51][52][53], as well as discussing their evolution with energy, however the relevance of these objects to multi gluon production has not been elucidated.
The plan of this paper is the following. In Sect. 2 we perform the calculation of the double inclusive gluon production at central rapidities in p-A. We describe the formalism and set out the approximation that we are using for calculating the target averages for a saturated target with confinement radius given by the inverse of saturation momentum. Besides, we discuss the meaning of various contributions to correlated production, identifying the Bose enhancement and HBT ones. We also show that the interesting correlated terms in this expression arise from the contribution of the quadrupole in the cross section. Section 3 contains the calculation for the triple inclusive gluon production. In Sect. 4 we provide a short discussion of our results.

The double inclusive gluon production in dilute-dense scattering
The aim of this section is to calculate the double inclusive gluon production in p-A collisions in terms of the dipole scattering amplitude. Our starting point is the general expression [57,58] for production of two gluons with pseudorapidities η 1 , η 2 and transverse momenta k 1 , k 2 Here, ρ a (x) is the color charge density in the projectile, U (x) is the adjoint Wilson line in the color field of the target representing the scattering matrix of a gluon at transverse coordinate x, a is the color index running from 1 to N 2 c − 1, z ≡ d 2 z and the Weiszacker-Williams field A i is given by Equation (1) is graphically illustrated in Fig. 1.
The averaging over ρ, · · · P , in Eq. (1) will be performed using the McLerran-Venugopalan (MV) model [45,46]. As for averaging over U (x), · · · T , it will not be performed explicitly but some properties of the distribution of U 's will be used to simplify and interpret our general expressions. We start with averaging of the cross section with respect to the projectile color charge distribution. We use the generalized MV model where the weight functional is Gaussian. Thus the average of any product of the color charge densities factorizes into a product of all possible Wick contractions: For the average of two projectile color charges we take a general form: We do not assume translational invariance of the projectile wave function. This means that the function μ 2 (x, y) depends both on the difference x − y and the center-of-mass coordinate x+y 2 . The finite transverse size of the projectile R is reflected in vanishing of Then the average of four projectile color charges reads Implementing this projectile averaging procedure we get We define the dipole and the quadrupole amplitudes in the standard way as and the corresponding Fourier transforms as The cross section is then written most conveniently as the sum of three terms: and where, for convenience, we have defined the Lipatov vertex The different contractions of the projectile color charge density leading to the three terms are illustrated in Figs. 2, 3 and 4. 2.2 Target averaging in double inclusive gluon production Our next step is to understand the general features of the target averages. Here we follow the logic of [40,41].
The cross section involves integration over the four coordinates of the product of the eikonal matrices U (x). One therefore expects that the main contribution will come from the region of the transverse plain where as many points are far away from each other as possible. Configurations in which points come close to each other in the transverse plain will give contributions suppressed by powers of area of the projectile.
On the other hand, one cannot have all four points far away from each other. This follows since the target field ensemble has to be color invariant. As a mater of fact, it is reasonable to assume that the color neutralization in the target ensemble is achieved on transverse distance scales of order 1/Q s . In order for the S-matrix on such a target to be non vanishing, the objects that scatter must be color singlets of size of order or smaller than 1/Q s . Thus, the maximal contribution must come from the configurations where the four points are combined into pairs, such that each pair is a singlet and the distance between the pairs is large. Taking into account only such configurations is equivalent to the calculation of target averages in which one factorizes the average of a product of any number of U matrices into averages of pairs with the basic "Wick contraction" given by where Note that only one color structure appears in this expression, the one where the left and right indexes of the two U -matrices are in color singlets. The physical reason for this is the following. Recall that the left index of U specifies the color of the incoming gluon, while the right index the color of the scattered gluon. We expect that on a saturated target the S-matrix of any nonsinglet colored state vanishes (black disk limit). The structure of Eq. (15) encodes precisely this property.
With these physical assumptions on the target field ensemble we have It is now straightforward to rewrite the double gluon inclusive production cross section entirely in terms of the dipole averages assuming translational invariance, We find Finally let us organize the terms in powers of N 2 c −1. Then dσ with 2.3 Identifying terms in double inclusive gluon production Given the expressions above it is quite straightforward to identify the physical origin of the various terms above. First of all we note, that the cross section is symmetric under the transformation (k 1 , k 2 ) → (k 1 , −k 2 ). This property of the dilute-dense approximation is well known in the literature. It is also known that the symmetry is "accidental" and it disappears once one includes higher order perturbative corrections [59][60][61]. We will therefore only consider half of the terms in Eq. (24), namely those that are written out explicitly.
To understand the meaning of the various terms we have to first of all specify μ 2 (k, p) a little further. Recall that it is defined as thee correlation function of the color charge density, which operationally reads In the hypothetical translational invariant limit, where the distribution of ρ is invariant under translations in the transverse plain, we would have The translationally invariant approximation is quite reasonable if we are interested in production of high transverse momentum gluons, however in a more accurate calculation the transverse size of the projectile should be reflected in μ 2 .
A reasonable way to include it is to introduce a form factor of the type where R is the radius of the projectile. Here T is roughly speaking the transverse dependent distribution of the valence charges, and F(x) is a soft form factor which is maximal at x = 0 and rapidly decreases to zero at |x| > 1. The exact form of the function F(x) does not matter for our purposes, it can be taken as a Gaussian F G (x) = exp (−x 2 ) , or as a Lorentzian F L (x) = 1/(1 + x 2 ), or any other function with these properties. The important property is that F(x) vanishes when the sum of the transverse momenta is not soft. Let us now examine the various terms in Eq. (24): • First off, the term I 0 obviously is just an uncorrelated production cross section, which is equal to the square of single gluon emission probability. It is not interesting from the point of view of correlated production. • The expression for I 1 contains two distinct terms, and it is easy to see that they have quite different origin. The first term is proportional to Note that the momenta k 1 −q 1 and k 2 −q 2 are the momenta the two gluons have in the projectile wave function, since k 1 and k 2 are the momenta in the final state, while q 1 and q 2 , being the arguments of the dipole scattering amplitudes, are the momentum transfers imparted to the two gluons during the propagation through the target. Due to the properties of the form factor F, this term is sharply peaked when the momenta of the two gluons in the projectile wave function are very close to each other, i.e., within the inverse projectile radius. This term therefore embodies the Bose enhancement in the incoming projectile wave function. • The nature of the second term in I 1 in Eq. (24) is defined by the factor This term enhances production of pairs of gluons with equal (up to 1/R) transverse momenta in the final state. This is a typical HBT contribution. • The first term in I 2 is proportional to Here the momentum exchange in the scattering of two gluons is the same. Naturally this term is associated with Bose enhancement in the target wave function. This term is somewhat different from the others in that it looks like the target gluon distribution here is regulated by the projectile size R. This is in fact natural. The Bose enhancement in the target wave function should be certainly regulated by the target and not the projectile size. However, very long transverse wave length gluons of the target are not probed by a smaller projectile, and thus do not contribute to the cross section. This is the reason why even though the enhancement is due to the properties of the target wave function, in the expression for the cross section the regulator is the projectile size.
It is interesting to note that in the glasma graph approximation the projectile and the target Bose enhancement terms appear at the same order in 1/N c . On the other hand in the proper dilute-dense treatment the target Bose enhancement terms come with further suppression factors. It is indeed a well known fact that some aspects of 1/N c counting are different in the dilute and dense limits [62,63]. • Finally the second term in I 2 has the same structure, as the projectile Bose enhancement term and constitutes an 1/N 2 c suppressed correction to this effect.

The triple inclusive gluon production in dilute-dense scattering
In this section, we study inclusive triple gluon production. The set up of the process is illustrated in Fig. 5. The formal expression of the three gluon production cross section, for gluons with pseudorapidities η i and transverse momenta k i , i = 1, 2, 3, can be simply written as where the coordinate of each Wilson line is written as a subscript.

Projectile averaging in triple inclusive gluon production
We first perform the averaging over the projectile color charge densities. As in the previous section we adopt the generalized MV model for the average of two projectile color charge densities and write down all possible Wick contractions of the product of the color charge densities. Then, the average of six projectile color charges reads  Fig. 6 Three-dipole (ddd) contribution to the three-gluon production cross section Here, we introduce a compact notation and write the coordinate of each color charge density as a subscript (and also omitted the subscript P from the averages). These terms can be categorized in three main groups. We have labeled the first term as three-dipole (ddd) contribution since it leads to the product of three dipoles when multiplied by the target scattering matrices U . The graphical illustration of this term is presented in Fig. 6. The next three terms are named as dipole-quadrupole (dQ) contribution since they generate a dipole and a quadrupole term after multiplication by U 's. One of these terms is illustrated in Fig. 7. The last four terms in Eq. (35) are labeled as sextupole (X) contribution since these terms generate the trace of all six Wilson lines. One of these terms (the one that is proportional to ρ a 1 x 1 ρ a 2 x 2 ) is illustrated in Fig. 8.
We start with the three-dipole contribution to the threegluon production cross section. We use Eq. (4) and substitute the first term of Eq. (35) in the three-gluon production cross section. A straightforward algebra gives the three-dipole contribution as We can now Fourier transform the three-dipole contribution. Using the standard definition of the dipole amplitude, Eq. (7), we can write the three-dipole contribution to the three-gluon production as where the function L i , defined in Eq. (14), gives the transverse momentum structure. Next, we consider the dipole-quadrupole (dQ) contribution to the three-gluon production cross section. The second, Fig. 7 Graphical illustration of the first term of the dQ-contribution to the three-gluon production cross section. It corresponds to the independent emission of gluon k 1 and interference of gluon k 2 and gluon k 3 Fig. 8 Graphical illustration of the term that is proportional to ρ a1 x1 ρ a2 x2 in the sextupole contribution third and fourth terms in Eq. (35) fall into this category and, using these three terms, we can write the dQ-contribution as The graphical illustration of the first term in Eq. (38) is shown in Fig. 7. This term represents the independent emission of the gluon k 1 while gluons k 2 and k 3 interfere. Indeed, the interference of the gluons k 2 and k 3 are exactly the type A and type C contributions introduced in the double inclusive gluon production calculation. The remaining two terms in Eq. (38) correspond to independent emission of the gluon k 2 with interference of gluons k 1 and k 3 , and independent emission of the gluon k 3 with interference of gluons k 1 and k 2 , respectively.
Fourier transform of Eq. (38) can be performed in the same way as before. Then, one can use the dipole (Eq. (9)) and quadrupole (Eq. (10)) amplitudes in momentum space and write the dipole-quadrupole contribution to the three-gluon production cross section as where the function L dQ is defined as Finally, let us consider the sextupole (X) contribution to the three-gluon production cross section. This contribution stems from the last four terms in Eq. (35). These are terms that include interference of all the three gluons . We have shown the illustration of the term that is proportional to ρ a 1 x 1 ρ a 2 x 2 in Fig. 8. After contracting the color indexes of these four terms with the Wilson line structure from the target side, we can write the X-contribution to the three-gluon production cross section as The sextuple amplitude is defined in the usual way as and its momentum space expression can be written as Finally, we can write the sextupole contribution the threegluon production as where and

Target averaging in triple inclusive gluon production
In this Subsection we perform the target averaging for the three-gluon production cross section by adopting the same procedure introduced in Sect. 2.2. We use Eq. (15) for the product of two Wilson lines and calculate the pairwise factorized expressions for the three-dipole, dipole-quadrupole and sextuple contributions to the triple inclusive gluon production cross section. In order to identify the "irreducible" quantum interference effects that involve all three gluons, we need to calculate the triple inclusive gluon production cross . Therefore, we will present the results for the pairwise contraction of three-dipoles, dipolequadrupole and sextuple amplitudes to all orders in number of colors but we will only take into account the relevant terms when calculating the explicit expressions for each contribution. We will assume translational invariance for the dipole, Eq. (19), when computing the cross sections.
Let us start with the three-dipole contribution. The result for the pairwise factorization of a generic three-dipole amplitude to all orders in N c reads Substituting this factorized expression into Eq. (36) we write the three-dipole contribution to the triple inclusive gluon production cross section as where we have defined I ddd,0 as Moreover, for the O 1 (N 2 c −1) 2 terms we have introduced a compact notation with The remaining terms can be defined by using the explicit expression of I ddd,1 and the symmetry properties: The pairwise contraction of a generic dipole-quadrupole term can be written as Using Eqs. (54) in (38), one can organize this contribution to the triple inclusive gluon production cross section as where we have introduced the same notation used in the threedipole contribution and we define with The remaining terms can again be written by using the symmetry properties and they are defined as In a similar manner, one can calculate the pairwise contraction of a generic sextupole term: We can now substitute Eq. (62) into the sextupole contribution to the triple inclusive gluon production cross section given by Eq. (41). The result reads × I X,1 + I X,2 + I X,3 + I X,4 + I X,5 where withĨ X,1 andĨ X,1 defined as The terms I X,2 and I X,3 can again be defined by using the symmetry properties as The explicit expressions for the remaining two terms read Finally, the triple inclusive gluon production cross section can be organized according to the powers in the number of colors and the result reads This is our final result for the triple inclusive gluon production cross section explicitly written up to order (N 2 c −1) −3 . It is straightforward to calculate the (N 2 c −1) −3 and (N 2 c −1) −4 terms with all the ingredients introduced in this Subsection. However, as mentioned earlier, in order to observe the quantum interference effects it is enough to calculate the triple inclusive gluon production cross section to order (N 2 c −1) −2 . Thus, we have not written the higher order terms explicitly.

Identifying terms in triple inclusive gluon production
Now, we can take a closer look at each term in the triple inclusive gluon production cross section separately and identify them one by one. We will follow the same logic that was introduced in Sect. 2.3 to perform this analysis, i.e., we will use the fact that where function F is a soft form factor that is peaked around zero and R is the radius of the projectile. • The only contribution that we get at O(1) is the I ddd,0 term. It is clear that this term is the classical contribution to triple inclusive gluon production cross section and it is equal to the product of three single inclusive gluon production cross sections. It is responsible for independent emission of all three gluons, with transverse momenta k 1 , k 2 and k 3 , and does not generate any correlations.
• I dQ,1 term has two contributions: The first one is proportional to The form factor is peaked around k 2 −q 2 = k 3 −q 3 while the gluon k 1 does not interfere with the remaining two gluons. Thus, it is clear that this term is a contribution to forward peak of Bose enhancement of gluons k 2 − q 2 and k 3 − q 3 in the projectile with the third gluon emitted independently. The mirror image of this term, given by the transformation k 2 → −k 2 , contributes to the backward peak of the Bose enhancement of the gluons k 2 −q 2 and k 3 − q 3 in the projectile. Obviously, in this case the third gluon is emitted independently as well.
The second term of I dQ,1 is proportional to In this contribution the form factor is peaked around k 2 = k 3 . Therefore, it is clear that this term is a contribution to the forward peak of the HBT correlations of the gluons k 2 and k 3 while the third gluon is emitted independently. The mirror image of this term is given by the transformation k 2 → −k 2 and will be a contribution to the backward peak of the HBT correlations of the gluons k 2 and k 3 . Since I dQ,2 and I dQ,3 are related to I dQ,1 by symmetry, it is obvious that these terms exhibit the same behavior as I dQ,1 but with gluons interchanged (1 ↔ 2 and 1 ↔ 3 respectively).
(iii) O(1/N 4 c ) terms: • I ddd,1 term: This term is proportional to In this term, the form factor is peaked around q 2 = q 3 where q 2 and q 3 are the momenta of the gluons in the target wave function. Therefore, this term is clearly a contribution to the forward peak of the Bose enhancement of the gluons q 2 and q 3 in the target wave function while the third gluon is emitted independently. Its mirror image given by the transformation k 3 → −k 3 is a contribution to the backward peak of the Bose enhancement of the gluons q 2 and q 3 in the target wave function. The remaining two terms that stem from the threedipole contribution at O(1/N 4 c ) are I ddd,2 and I ddd,3 . These terms can be obtained by exchanging (1 ↔ 2) and (1 ↔ 3). Hence, they exhibit the same behavior as I ddd,1 . Namely, I ddd,2 is a contribution to the (forward/backward peaks) Bose enhancement of the gluons q 1 and q 3 in the target wave function while gluon q 2 is emitted independently, and I ddd,3 is a contribution to the (forward/backward peaks of the) Bose enhancement of the gluons q 1 and q 2 in the target wave function while the gluon q 3 is emitted independently.
• I dQ,1 term: This term is proportional to The form factor is again peaked around k 2 −q 2 = k 3 −q 3 . Therefore, this term contributes to the forward peak of the Bose enhancement of the gluons k 2 −q 2 and k 3 −q 3 in the projectile wave function with the third gluon emitted independently. Clearly, the mirror image is a contribution to the backward peak. I dQ,2 and I dQ,3 exhibit the same behavior with the exchange of (1 ↔ 2 ) and (1 ↔ 3). We would like to emphasize that these three terms are N csuppressed corrections to the (forward/backward peaks of the) Bose enhancement of the two gluons in the projectile wave function while the third gluon is emitted independently, which was the behavior that we have observed in the first part of the I dQ,1 , I dQ,2 and I dQ,3 terms. • I X,1 term: This is the first term that we are analyzing which stems from the sextupole contribution. Unlike the previous terms that we have analyzed, all the terms that originate from the sextupole contribution lead to the interference of all three gluons and there are no independent emissions. After this short comment, let us take a closer look at I X,1 term. It is composed of four different contributions. The first contribution is proportional to The first form factor is peaked around k 1 = k 2 while the other two form factors peaked around k 1 − q 1 = k 3 − q 3 . Thus, this term is a contribution to the forward peak of the HBT correlations of gluons k 1 and k 2 together with the forward peak of the Bose enhancement of gluons k 1 − q 1 and k 3 − q 3 in the projectile wave function. Its mirror image, given by the transformation k 3 → −k 3 , is a contribution to backward peak of the Bose enhancement of k 1 − q 1 and k 3 − q 3 in the projectile wave function together with a contribution to the forward peak of the HBT correlations of gluons k 1 and k 2 .
The second term in I X,1 is proportional to The first form factor is peaked around k 1 = −k 2 and the remaining two form factors are the same as the previous contribution. Therefore, one can identify this term as a contribution to the backward peak of the HBT correlations of the gluons k 1 and k 2 together with the forward peak of the Bose enhancement of the gluons k 1 − q 1 and k 3 − q 3 in the projectile wave function. Obviously, its mirror image, given by the transformation k 3 → −k 3 , is a contribution to the backward peak of HBT correlations of gluons k 1 and k 2 together with the contribution to the backward peak of the Bose enhancement of the gluons k 1 − q 1 and k 3 − q 3 in the projectile wave function. The third term in I X,1 is proportional to The first form factor in this term is peaked around k 1 = k 2 while the other two form factors are peaked around k 1 − q 1 = q 3 − k 3 . Thus, it is clear that this term is a contribution to the forward peak of HBT correlations of gluons k 1 and k 2 together with the backward peak of the Bose enhancement of gluons k 1 − q 1 and k 3 − q 3 in the projectile. Its mirror image, given by the transformation k 1 → −k 1 , is a contribution to the backward peak of HBT correlations of gluons k 1 and k 2 together with the forward peak of the Bose enhancement of gluons k 1 − q 1 and k 3 − q 3 . The last term in I X,1 is proportional to The first form factor is peaked around k 1 = −k 2 and the remaining two form factors are peaked around k 1 − q 1 = q 3 − k 3 . Therefore, it is a contribution to the backward peak of HBT correlations of gluons k 1 and k 2 together with the backward peak of the Bose enhancement of gluons k 1 − q 1 and k 3 − q 3 in the projectile wave function. Its mirror image, given by the transformation k 1 → −k 1 , is a contribution to the forward peak of HBT correlations of gluons k 1 and k 2 together with the forward peak of the Bose enhancement of gluons k 1 − q 1 and k 3 − q 3 in the projectile wave function. It can easily be shown that I X,2 and I X,3 exhibit a similar behavior as I X,1 with gluons interchanged. I X,2 contributes to (backward/forward) HBT correlations of gluons k 2 and k 3 together with (backward/forward) Bose enhancement of gluons k 1 − q 1 and k 2 − q 2 in the projectile wave function. Finally, I X,3 contributes to (backward/forward) HBT correlations of gluons k 1 and k 3 together with (backward/forward) Bose enhancement of gluons k 2 −q 2 and k 3 −q 3 in the projectile wave function. • I X,4 term: This term has four different contributions and they all contribute to the Bose enhancement of all three gluons. The first contribution in I X,4 is proportional to The first form factor is peaked around k 1 − q 1 = k 2 − q 2 , the second one is peaked around k 1 −q 1 = k 3 −q 3 and the last one is peaked around k 2 −q 2 = k 3 −q 3 . Therefore, it is clear that this term is a contribution to the forward peak of the Bose enhancement of gluons k 1 − q 1 and k 2 − q 2 , gluons k 1 −q 1 and k 3 −q 3 , and gluons k 2 −q 2 and k 3 −q 3 in the projectile wave function. Its mirror image, given by the transformation k 3 → −k 3 , is a contribution to the forward peak of the Bose enhancement of the gluons k 2 − q 2 and k 1 − q 1 together with the backward peak of the Bose enhancement of the gluons k 1 − q 1 and k 3 − q 3 , and gluons k 2 − q 2 and k 3 − q 3 in the projectile wave function. The second contribution in I X,4 is proportional to Clearly, this is a contribution to the forward peak of Bose enhancement of gluons k 2 − q 2 and k 3 − q 3 together with the backward peak of the Bose enhancement of gluons k 1 − q 1 and k 2 − q 2 , and gluons k 1 − q 1 and k 3 − q 3 in the projectile wave function. Its mirror image is given by the transformation k 1 → −k 1 . Therefore, it contributes to the forward peak of Bose enhancement of gluons k 2 − q 2 and k 3 −q 3 , gluons k 1 −q 1 and k 2 −q 2 , and gluons k 1 −q 1 and k 3 − q 3 in the projectile wave function. The third contribution in I X,4 is proportional to By looking at the peaks of the form factors, it is straightforward to see that this is a contribution to the forward peak of the Bose enhancement of the gluons k 1 − q 1 and k 3 − q 3 together with a backward peak of the Bose enhancement of the gluons k 1 − q 1 and k 2 − q 2 , and gluons k 3 −q 3 and k 2 −q 2 in the projectile wave function. Its mirror image is given by the transformation k 3 → −k 3 . Therefore, it contributes to the forward peak of the Bose enhancement of the gluons k 2 − q 2 and k 3 − q 3 together with a backward peak of the Bose enhancement of the gluons k 1 − q 1 and k 3 − q 3 , and gluons k 1 − q 1 and k 2 − q 2 in the projectile wave function. The last contribution in I X,4 is proportional to Clearly, this term is a contribution to the forward peak of the Bose enhancement of gluons k 1 − q 1 and k 2 − q 2 together with a backward peak of the Bose enhancement of the gluons k 1 − q 1 and k 3 − q 3 , and gluons k 2 − q 2 and k 3 − q 3 in the projectile wave function. Its mirror image is given by the transformation k 1 → −k 1 . Thus, it contributes to the forward peak of the Bose enhancement of gluons k 1 − q 1 and k 3 − q 3 together with a backward peak of the Bose enhancement of the gluons k 1 − q 1 and k 2 − q 2 , and gluons k 2 − q 2 and k 3 − q 3 in the projectile wave function • I X,5 term: The last term that stems from the sextupole contribution at O(1/N 4 c ) is the I X,5 term. It has four different contributions and all of them contribute to the HBT correlations of all three gluons. The first contribution in I X,5 is proportional to In this term, the form factors are peaked around k 2 = k 1 , k 1 = k 3 and k 3 = k 2 . Thus, it is a contribution to the forward peak of the HBT correlations of gluons k 1 and k 2 , gluons k 1 and k 3 , and gluons k 2 and k 3 . The mirror image of this term is given by the transformation k 3 → −k 3 and it contributes to the forward peak of the HBT correlations of k 1 and k 2 together with backward peak of the HBT correlations of the gluons k 1 and k 3 , and gluons k 2 and k 3 . The second contribution in I X,5 is proportional to The form factors in this term are peaked around k 2 = k 3 , k 1 = −k 2 and k 1 = −k 3 . Hence, this term contributes to the forward peak of the HBT correlations of gluons k 2 and k 3 together with the backward peak of the HBT correlations of gluons k 1 and k 2 , and gluons k 1 and k 3 . The mirror image of this term is given by the transformation k 1 → −k 1 . Therefore, it contributes to the forward peak of the HBT correlations of gluons k 1 and k 3 , gluons k 2 and k 3 , and gluons k 1 and k 2 . The third contribution in I X,5 is proportional to In this term the form factors are peaked around k 2 = −k 1 , k 3 = k 1 and k 2 = −k 3 . It contributes to the forward peak of HBT correlations of gluons k 3 and k 1 together with a backward peak of the HBT correlations of gluons k 1 and k 2 , and gluons k 2 and k 3 . On the other hand, its mirror image (given y the transformation k 3 → −k 3 ) contributes to the forward peak of the HBT correlations of gluons k 3 and k 2 together with the backward peak of the HBT correlations of gluons k 1 and k 2 , and gluons k 1 and k 3 . Finally, the last contribution in I X,5 is proportional to Clearly, this term contributes to the forward peak of the HBT correlations of the gluons k 1 and k 2 together with the backward peak of the HBT correlations of the gluons k 2 and k 3 , and gluons k 3 and k 1 . Its mirror image is given by the transformation k 1 → −k 1 . Therefore, it contributes to the forward peak of the HBT correlations of the gluons k 1 and k 3 together with the backward peak of the HBT correlations of the gluons k 2 and k 3 , and gluons k 1 and k 2 .

Conclusions
To conclude, we have calculated the inclusive two and three gluon production in p-A collisions at mid rapidity in the CGC formalism. We use the generalized McLerran-Venugopalan model to perform the projectile averaging. This model allows for accounting for the finite transverse size of the projectile. We observe that in the full dilute-dense limit that goes beyond the glasma graphs, the expression for the cross section double inclusive gluon production contains two types of terms: a product of two dipole amplitudes and a quadrupole. The origin of the quadrupole term is the contribution to scattering where the two incoming gluons exchange color while propagating through the target.
We further used simple physical assumptions about the target structure to express the quadrupole average in terms of products of averages of two dipoles. We stress that this factorization is not a result of the large N c limit and is, in fact, completely unrelated with the large N c expansion. It is rather the result of the physical expectation that the color neutralization of the target ensemble happens on transverse scales of order 1/Q s . As discussed in [40,41], this approximation does not take into account the "classical" correlated term arising from the contributions to transverse integrals where the produced particles have to be close to each other in the incoming projectile wave function.
The resulting expression for double inclusive production is quite simple, see Eq. (23). It exhibits very similar terms to those obtained for double inclusive quark production in [40,41]. Just like there, we find that all the quantum interference effects that constitute the genuine multiparticle correlations (as can be extracted e.g. using the cumulant method [64][65][66][67]) at order 1/N 2 c , given by the I 1 term in (23), originate from the quadrupole. Our results in this part of the paper are consistent with [55]. We observe two types of quantum interference effects -the Bose enhancement of gluons in the projectile wave function and the Hanbury-Brown-Twiss interference effect. In the case of gluons (differently from the production of identical quarks) the two effects enhance the production of gluons which are both collinear and anticollinear in the transverse plane. These same effects have been observed in the glasma graph calculation earlier [25][26][27][28][29][30]32,33].
Additionally to [55] we identify the Bose enhancement terms in the target wave function. There is an interesting difference in this aspect between the full dilute-dense calculation presented here, and the glasma graphs of [25][26][27][28][29][30]. In the glasma graph calculation, which is based on the dilutedilute limit, and is therefore symmetric between the projectile and the target, the Bose enhancement in the target wave function is of the leading order in 1/N c , just like the Bose enhancement in the projectile. In the complete dilute-dense framework utilized in the present paper this effect, although present, is suppressed as 1/N 2 c relative to the projectile Bose enhancement effect.
We have also computed the triple inclusive gluon production in the dilute-dense limit for the first time. Our result is given in Eq. (71). Although the expressions are lengthy, the basic physics is very simple. One observes terms in the production which involve the interference between two of the gluons, and independent emission of the third one. Such terms would be subtracted if one would calculate the three particle cumulant rather than write up the three particle inclusive cross section. Additionally, there are genuine three particle correlation terms which appear at order 1/N 4 c -these are the I X,i terms in (71). These terms originate solely from the sextupole in the cross section. They include interference due to Bose enhancement and HBT contributions of all three particles, as well as "mixed" terms where two of the particles are coupled via Bose enhancement, and other two due to the HBT effect.
These features are very similar to those observed in three quark production in [40,41]. There is one interesting albeit expected difference. In production of three quarks the correction due to interference of all three quarks has an opposite sign to the term where one quark is emitted independently of the other two (which interfere). It therefore has a flavor of "unitarization correction" as discussed in [40,41]. For gluon production this is not the case. All terms in the cross section are positive, and thus the genuine three gluon interference term enhances the correlations rather than weakening them.
Correlations among more than two particles have been studied at RHIC and the LHC, e.g., in the context of the study of properties of the produced medium [68][69][70], for the study of HBT correlations [71] or for extraction of the azimuthal asymmetries [72,73]. Without some simplifica-tions and implementation of models for the proton and nuclei, it is very difficult to provide definite predictions beyond the long range pseudorapidity nature of our correlations. Such study and the computation of several gluon inclusive production that is paved by the present work, with the obvious extension to the four gluon inclusive case, are left for future work.