Multi-particle production in proton–nucleus collisions in the color glass condensate

Abstract

We compute multi-gluon production in the Color Glass Condensate approach in dilute-dense collisions, $\hbox {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üsthoff 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.

1 Introduction

Small collision systems, proton-proton (pp) and proton-nucleus (pA), studied at the Large Hadron Collider (LHC) show many of the characteristics [1,2,3,4,5] that in heavy ion collisions are considered as signatures of the formation of hot deconfined partonic matter, the Quark Gluon Plasma. The most prominent example is the existence of azimuthal correlations in the two-particle inclusive distributions that are extended in pseudorapidity and show maxima when the particle transverse momenta are either parallel or antiparallel. This finding, named the ridge, was first observed in high multiplicity pp collisions [6], and then for smaller multiplicities [7,8,9,10], in pPb collisions [11,12,13,14,15,16] and in association with Z boson production [17]. It was also observed in pAu, dAu and $^3$HeAu collisions at the Relativistic Heavy Ion Collider (RHIC) [18,19,20,21]. Azimuthal asymmetries in particle production have also been searched for in even smaller systems: in $\gamma$Pb through ultraperipheral collisions at the LHC [22] where they were found, and in $e^+e^-$ collisions [23] at the Large Electron-Positron collider and deep inelastic scattering in ep at the Hadron–Elektron–Ringanlage [24] with inconclusive results.

The key open question nowadays is the clarification of the origin of such collective behaviour. In heavy ion collisions where the partonic density is very large, a natural explanation is that collectivity is built through the strong final state interactions of the created system. Such explanation looks justified by the success of viscous relativistic hydrodynamics [25, 26] for describing the observed experimental features in soft particle production. The open questions at this moment are how the conditions for hydrodynamics to be applicable are reached from an initial state that is very far from equilibrium [27] – the emergence of the macroscopic description given by hydrodynamics from Quantum Chromodynamics (QCD) –, for which both strong and weak coupling explanations have been proposed (see, e.g., [28]), and why hydrodynamics seems to work even for large anisotropies, outside its presumed range of applicability. Hydrodynamics appears as the effective theory for describing the soft modes of any field theory, see e.g. [29] and references therein.

The success of the application of hydrodynamics for describing azimuthal asymmetries in small systems, pp and pPb collisions [26, 30], while requiring careful choices of the initial conditions, pushes this description to small collision areas and low particle densities where non-hydrodynamic modes play a very important role [26, 31, 32]. Therefore, it seems sensible to explore other alternatives. The color glass condensate (CGC) [33, 34], as weak coupling non-perturbative effective theory for QCD at high energies and partonic densities, offers a framework where azimuthal asymmetries can be calculated from first principles, see the review [35] and references therein. Correlations in the final state reflect those found in the wave function of the projectile and target hadrons or nuclei, assuming that final state effects, including hadronisation, do not wash them out.

The initial versus final state origin of azimuthal correlations in small systems has been subject to intense scrutiny in recent years [36,37,38]. At present, no CGC-based model is able to fully describe the existing experimental data. Still, the search for observables that may discriminate initial from final effect continues, e.g., the correlation of $v_2$ with the mean transverse momentum of the particles produced in the collisions [39,40,41,42] that has also been analysed in the CGC [43]. Also many particle cumulants are expected to be crucial. For example, four particle cumulants $c_2\{4\}$, with $v_2\{4\}=[-c_2\{4\}]^{1/4}$ (definitions of all these quantities will be provided below), change sign from positive to negative with increasing particle multiplicity in the event, with a smooth behaviour from small to large systems and from smaller to larger energies. This change of sign is associated with the onset of true collective flow of final state origin because higher order cumulants are less sensitive to non-flow contributions than those computed from two-particle correlations. In the CGC numerical implementation in [44, 45] the change of sign of $c_2\{4\}$ was interpreted as the transition from a dilute-dilute situation, described by the glasma graph approach [46, 47] where azimuthal correlations correspond to the Bose enhancement of the gluons in the wave function of the colliding hadrons and to the Hanbury-Brown-Twiss (HBT) effect for the final gluons [48,49,50,51], to a dilute-dense situation where multiple scattering dominates (for a discussion on density correlations to the dilute-dense situation, see [52] and references therein).

The goal of this work is the extension of the calculations of multiparticle production in the CGC in the dilute-dense situation (suitable for pA collisions) performed in [53] to four gluon production (see [54] for inclusive cross sections involving final state quark-antiquark pairs), and the computation of the two and four particle cumulants.¹ Note that up to four gluon production was previously computed in the glasma graph approach [55], and arguments 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 Sect. 2 we introduce the formalism to compute particle production in the CGC. Section 3 is devoted to the calculation of projectile and target averages required to obtain the final results. Then, in Sect. 4 we present our results and in Sect. 5 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.

2 Theoretical background on multi-particle correlation

2.1 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 $\rho ^a(\mathbf{x })$, with superindex a denoting color and $\mathbf{x }$ the transverse position. The target is characterised by a strong field ${\mathcal {A}}^\mu (\mathbf{x })={\mathcal {A}}^{a\mu }(\mathbf{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 ${\mathcal {A}}^+=0$ and neglecting the transverse components of the field, the amplitude for producing a gluon with transverse momentum $\mathbf{k }$, pseudorapidity $\eta$, polarization $\lambda$ 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

$$\begin{aligned} {\mathcal {M}}^a_\lambda (\eta ,\mathbf{k })=g \int \frac{d^2 \mathbf{q }}{(2 \pi )^2} \overline{{\mathcal {M}}}^{ab}_\lambda (\eta ,\mathbf{k },\mathbf{q }) \rho ^b(\mathbf{k }-\mathbf{q }), \end{aligned}$$

(1)

with $\rho ^a(\mathbf{q })$ the Fourier transform of the color charge density of the projectile which is defined (see e.g. [68]) as

$$\begin{aligned} \rho ^a(\mathbf{x })=\int \frac{dp^+}{2\pi } a^{\dagger }_i( p^+,\mathbf{x }) T^a a_i(p^+,\mathbf{x }), \end{aligned}$$

(2)

where $a^{\dagger }_i( p^+,\mathbf{x })$ and $a_i(p^+,\mathbf{x })$ are the creation and annihilation operators for gluons with longitudinal momentum $p^+$ at transverse position $\mathbf{x }$, respectively. The reduced matrix amplitude $\overline{{\mathcal {M}}}^{ab}_\lambda (\eta ,\mathbf{k },\mathbf{q })$ is derived in [66, 69, 70] and it reads

$$\begin{aligned} \overline{{\mathcal {M}}}^{ab}_\lambda (\eta ,\mathbf{k },\mathbf{q })&=\epsilon ^{i*}_\lambda (\mathbf{k }) i e^{i k^- L^+} \Bigg \{ 2\frac{\mathbf{k }^i}{\mathbf{k }^2} \int _\mathbf{y } e^{-i\mathbf{q }\mathbf{y }} {\mathcal {U}}_{\mathbf{y }}^{ab}(L^+,0)-2\frac{(\mathbf{k }-\mathbf{q })^i}{(\mathbf{k }-\mathbf{q })^2} \int _{\mathbf{y },\mathbf{x }} e^{i (\mathbf{k }-\mathbf{q }) \mathbf{y }-i \mathbf{k } \mathbf{x }} {\mathcal {G}}^{ab}_{k^+}(L^+,\mathbf{x };0,\mathbf{y }) \Bigg . \nonumber \\ \Bigg .&\quad + \int _{\mathbf{x },\mathbf{y }} e^{i (\mathbf{k }-\mathbf{q }) \mathbf{y }} \frac{1}{k^+} \int _0^{L^+} dy^+ e^{-i\mathbf{k }\mathbf{x }} [\partial _{\mathbf{y }^i} {\mathcal {G}}^{ac}_{k^+}(L^+,\mathbf{x };y^+,\mathbf{y })] {\mathcal {U}}_{\mathbf{y }}^{cb}(y^+,0) \Bigg \}. \end{aligned}$$

(3)

In this equation $\int _\mathbf{x }\equiv \int d^2\mathbf{x }$, $\epsilon ^{i*}_\lambda (\mathbf{k })$ is the polarisation vector, $k^-=\mathbf{k }^2/(2k^+)$, $k^+=e^\eta /\sqrt{2}$, $L^+$ is the longitudinal length of the target, $\mathbf{q }$ is the transverse momentum transferred from the target during the interaction and $\mathbf{k}-q$ is the transverse momenta of the projectile color charge density, and

$$\begin{aligned} {\mathcal {G}}^{ab}_{k^+}(x^+,\mathbf{x };y^+,\mathbf{y })= & {} \Theta (x^+-y^+) \int _{\mathbf{y }}^{\mathbf{x }} {\mathcal {D}} \mathbf{z } \exp \left[ \frac{i k^+}{2} \int _{y^+}^{x^+} d z^+ \dot{\mathbf{z }}^2(z^+) \right] {\mathcal {U}}^{ab}_{\mathbf{z }} \left( x^+,y^+ \right) \end{aligned}$$

(4)

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 $\mathbf{z }(x^+)=\mathbf{x }$, $\mathbf{z }(y^+)=\mathbf{y }$. We use light-cone coordinates $x^\pm =(x^0\pm x^3)/\sqrt{2}$.

$$\begin{aligned} {\mathcal {U}}_{\mathbf{x }}^{ab}( x^+,y^+)={\mathcal {P}} \exp \left\{ i g \int _{y^+}^{x^+} dz^+ {\mathcal {A}}^-(z^+,\mathbf{x }) \right\} ^{ab} \end{aligned}$$

(5)

is the Wilson line that accounts for the multiple gluon exchanges with the target.

It is necessary to mention that the reduced matrix amplitude $\overline{{\mathcal {M}}}^{ab}_\lambda (\eta ,\mathbf{k },\mathbf{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^+ \rightarrow 0$, $k^+ \rightarrow \infty$ and assuming that the target field is a local shock-wave, ${\mathcal {A}}^-(z^+,\mathbf{x }) \propto \delta (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

$$\begin{aligned} {\mathcal {G}}^{ab}_{k^+}({\underline{x}},{\underline{y}})\rightarrow {\mathcal {U}}^{ab}_{\mathbf{x }} \left( x^+,y^+ \right) \delta ^{(2)}(\mathbf{x }-\mathbf{y }). \end{aligned}$$

(6)

Consequently, the reduced matrix amplitude given in Eq. (3) simplifies as well and it reads

(7)

where we have introduced the Lipatov vertex

(8)

and changed the notation of the Wilson lines, $U^{ab}(\mathbf{y })={\mathcal {U}}_{\mathbf{y }}^{ab}(L^+,0)$.

Fig. 1

Physical interpretation of the Lipatov vertex, with the vertical dashed line denoting the interaction with the target

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 $\rho ^a(\mathbf{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 $\mathbf{k }_i$, pseudorapidity $\eta _i$, color $a_i$ and polarization $\lambda _i$ is given, in terms of the amplitude matrix that is leading for $g\rho ^a (\mathbf{q })\sim 1$, as

$$\begin{aligned} 2^n (2 \pi )^{3n} \frac{d^n N}{\prod _{i=1}^n d \eta _i d^2 \mathbf{k }_i}&= \Big \langle {\mathcal {M}}_{\lambda _1}^{a_1}(\eta _1,\mathbf{k }_1)\cdots {\mathcal {M}}_{\lambda _n}^{a_n}(\eta _n,\mathbf{k }_n)\left( {\mathcal {M}}_{\lambda _n}^{a_n}(\eta _n,\mathbf{k }_n)\right) ^{\dagger }\cdots \left( {\mathcal {M}}_{\lambda _1}^{a_1}(\eta _1,\mathbf{k }_1)\right) ^{\dagger } \Big \rangle _{p,T}, \end{aligned}$$

(9)

where 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 $\eta _i$.

Using Eq. (1) and dropping the dependence on $\eta$ due to the eikonal approximation, we can write this expression as (see Fig. 2)

$$\begin{aligned} 2^n(2 \pi )^{3n} \frac{d^n N}{\prod _{i=1}^{n}d^2 \mathbf{k }_i}=&\,g^{2n} \int \left( \prod _{i=1}^{2n} \frac{d^2 \mathbf{q }_i}{(2 \pi )^2}\right) \Big \langle \rho ^{b_1}(\mathbf{k }_1-\mathbf{q }_1)\rho ^{b_2 \dagger }(\mathbf{k }_1-\mathbf{q }_2)\cdots \rho ^{b_{2n-1}}(\mathbf{k }_n-\mathbf{q }_{2n-1})\rho ^{b_{2n} \dagger }(\mathbf{k }_n-\mathbf{q }_{2n}) \Big \rangle _{p} \nonumber \\&\times \Big \langle \overline{{\mathcal {M}}}^{a_1 b_1}_{\lambda _1}(\mathbf{k }_1,\mathbf{q }_1) \overline{{\mathcal {M}}}^{b_2 a_1 \dagger }_{\lambda _1}(\mathbf{k }_1,\mathbf{q }_2)\cdots \overline{{\mathcal {M}}}^{a_{n} b_{2n-1}}_{\lambda _{n}}(\mathbf{k }_n,\mathbf{q }_{2n-1}) \overline{{\mathcal {M}}}^{b_{2n} a_{n} \dagger }_{\lambda _{n}}(\mathbf{k }_n,\mathbf{q }_{2n}) \Big \rangle _{T}. \end{aligned}$$

(10)

Solving this equation is the main point of this work and will be the focus of the discussion in the next sections.

Fig. 2

Diagram showing Eq. (10) with its momentum assignments. Each group of three vertical straight lines represents the rescattering with the target, the black blobs Lipatov vertices and the vertical dashed line the cut

Equation (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 $|\mathbf{k }|/Q_s \gg 1$ with $Q_s$ the saturation momentum of the gluons in the target [33, 34], or equivalently $g A^- \ll 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.

2.2 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:

$$\begin{aligned}&\Big \langle \overline{{\mathcal {M}}}^{a_1 b_1}_{\lambda _1}(\mathbf{k }_1,\mathbf{q }_1) \overline{{\mathcal {M}}}^{b_2 a_1 \dagger }_{\lambda _1}(\mathbf{k }_1,\mathbf{q }_2) \cdots \overline{{\mathcal {M}}}^{a_{n} b_{2n-1}}_{\lambda _{n}}(\mathbf{k }_n,\mathbf{q }_{2n-1}) \overline{{\mathcal {M}}}^{b_{2n} a_{n} \dagger }_{\lambda _{n}}(\mathbf{k }_n,\mathbf{q }_{2n}) \Big \rangle _{T} \nonumber \\&\quad \propto \int _{\mathbf{y }_1 \cdots \mathbf{y }_{2n}} e^{-i \mathbf{q }_1 \cdot \mathbf{y }_1+i \mathbf{q }_2 \cdot \mathbf{y }_2 \cdots +i \mathbf{q }_{2n} \cdot \mathbf{y }_{2n} } \Big \langle U(\mathbf{y }_1)^{a_1 b_1}U^\dagger (\mathbf{y }_2)^{b_2 a_2} \cdots U^\dagger (\mathbf{y }_{2n})^{b_{2n} a_{2n}}\Big \rangle _{T}, \end{aligned}$$

(11)

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 operator

$$\begin{aligned} {\tilde{Q}}(\mathbf{q }_1,\mathbf{q }_2,\mathbf{q }_3,\mathbf{q }_4)&=\int _{\mathbf{y }_1 \mathbf{y }_2 \mathbf{y }_3 \mathbf{y }_4} e^{-i \mathbf{q }_1 \cdot \mathbf{y }_1+i \mathbf{q }_2 \cdot \mathbf{y }_2 -i \mathbf{q }_3 \cdot \mathbf{y }_3 +i \mathbf{q }_{4} \cdot \mathbf{y }_{4} } \frac{1}{N_c^2-1}\Big \langle Tr \left[ U(\mathbf{y }_1) U^\dagger (\mathbf{y }_2) U(\mathbf{y }_3) U^\dagger (\mathbf{y }_{4})\right] \Big \rangle _{T} \nonumber \\&\equiv \int _{\mathbf{y }_1 \mathbf{y }_2 \mathbf{y }_3 \mathbf{y }_4} e^{-i \mathbf{q }_1 \cdot \mathbf{y }_1+i \mathbf{q }_2 \cdot \mathbf{y }_2 -i \mathbf{q }_3 \cdot \mathbf{y }_3 +i \mathbf{q }_{4} \cdot \mathbf{y }_{4} } Q(\mathbf{y }_1,\mathbf{y }_2,\mathbf{y }_3,\mathbf{y }_4). \end{aligned}$$

(12)

Following the arguments in [58, 59], the configuration of the transverse coordinates $\mathbf{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_s^{-1}$, 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 $\mathbf{y }_i$ and $\mathbf{y }_j$, will have a vanishing correlation when $|\mathbf{y }_i-\mathbf{y }_j| \gg Q_s^{-1}$. 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

$$\begin{aligned} Q(\mathbf{y }_1,\mathbf{y }_2,\mathbf{y }_3,\mathbf{y }_4) \approx&\,\mathop {\mathop {\mathop {\lim }\limits _{|\mathbf{y }_1 -\mathbf{y }_2| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_3-\mathbf{y }_4| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_1-\mathbf{y }_3| \lesssim Q_s^{-1}} Q(\mathbf{y }_1,\mathbf{y }_2,\mathbf{y }_3,\mathbf{y }_4) + \mathop {\mathop {\mathop {\lim }\limits _{|\mathbf{y }_1-\mathbf{y }_2| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_3-\mathbf{y }_4| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_1-\mathbf{y }_3| \gg Q_s^{-1} } Q(\mathbf{y }_1,\mathbf{y }_2,\mathbf{y }_3,\mathbf{y }_4) \nonumber \\&+ \mathop {\mathop {\mathop {\lim }\limits _{|\mathbf{y }_1-\mathbf{y }_3| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_2-\mathbf{y }_4| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_1-\mathbf{y }_2| \gg Q_s^{-1} } Q(\mathbf{y }_1,\mathbf{y }_2,\mathbf{y }_3,\mathbf{y }_4) + \mathop {\mathop {\mathop {\lim }\limits _{|\mathbf{y }_1-\mathbf{y }_4| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_2-\mathbf{y }_3| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_1-\mathbf{y }_2| \gg Q_s^{-1}} Q(\mathbf{y }_1,\mathbf{y }_2,\mathbf{y }_3,\mathbf{y }_4). \end{aligned}$$

(13)

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

$$\begin{aligned}&\mathop {\mathop {\mathop {\lim }\limits _{|\mathbf{y }_1 -\mathbf{y }_2| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_3-\mathbf{y }_4| \lesssim Q_s^{-1}}} \limits _{|\mathbf{y }_1-\mathbf{y }_3| \gg Q_s^{-1}} \Big \langle U(\mathbf{y }_1)^{a_1 b_1}U^\dagger (\mathbf{y }_2)^{b_1 a_2} U(\mathbf{y }_3)^{a_2 b_2} U^\dagger (\mathbf{y }_{4})^{b_{2} a_{1}}\Big \rangle _{T} \approx \Big \langle U(\mathbf{y }_1)^{a_1 b_1}U^\dagger (\mathbf{y }_2)^{b_1 a_2} \Big \rangle _{T} \Big \langle U(\mathbf{y }_3)^{a_2 b_2} U^\dagger (\mathbf{y }_{4})^{b_{2} a_{1}}\Big \rangle _{T} \end{aligned}$$

(14)

and analogously to the other terms. Thus we can write the quadrupole operator as a sum of products of 2-point functions,

$$\begin{aligned}&\Big \langle U(\mathbf{y }_1)^{a_1 b_1}U^\dagger (\mathbf{y }_2)^{b_1 a_2} U(\mathbf{y }_3)^{a_2 b_2} U^\dagger (\mathbf{y }_{4})^{b_{2} a_{1}}\Big \rangle _{T} \approx \Big \langle U(\mathbf{y }_1)^{a_1 b_1}U^\dagger (\mathbf{y }_2)^{b_1 a_2} \Big \rangle _{T} \Big \langle U(\mathbf{y }_3)^{a_2 b_2} U^\dagger (\mathbf{y }_{4})^{b_{2} a_{1}}\Big \rangle _{T} \nonumber \\&\qquad +\Big \langle U(\mathbf{y }_1)^{a_1 b_1} U(\mathbf{y }_3)^{a_2 b_2} \Big \rangle _{T} \Big \langle U^\dagger (\mathbf{y }_2)^{b_1 a_2} U^\dagger (\mathbf{y }_{4})^{b_{2} a_{1}}\Big \rangle _{T}+\Big \langle U(\mathbf{y }_1)^{a_1 b_1}U^\dagger (\mathbf{y }_{4})^{b_{2} a_{1}} \Big \rangle _{T} \Big \langle U^\dagger (\mathbf{y }_2)^{b_1 a_2} U(\mathbf{y }_3)^{a_2 b_2} \Big \rangle _{T}, \end{aligned}$$

(15)

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_\perp Q_s^{-2}$ with $S_\perp$ 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:

$$\begin{aligned}&\Big \langle U(\mathbf{y }_1)^{a_1 b_1}U(\mathbf{y }_2)^{a_2 b_2} \cdots U(\mathbf{y }_{2n})^{a_{2n} b_{2n}}\Big \rangle _{T} =\sum _{\sigma \in \Pi (\chi )} \prod _{\{\alpha ,\beta \}\in \sigma } \Big \langle U(\mathbf{y }_\alpha )^{a_\alpha b_\alpha } U(\mathbf{y }_\beta )^{a_\beta b_\beta } \Big \rangle _{T}\, , \end{aligned}$$

(16)

being $\chi =\{1,2,\dots ,2 n\}$ and $\Pi (\chi )$ the set of partitions of $\chi$ 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).

2.3 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_\perp \sim Q_s$. The correlation length of the fields is roughly $Q_s^{-1}$ 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.

Fig. 3

Picture of the chromoelectric fields inside the target

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

$$\begin{aligned} c_n\{2\}&= \Big \langle e^{i n ( \phi _1 -\phi _2)} \Big \rangle , \end{aligned}$$

(17)

$$\begin{aligned} c_n\{4\}&= \Big \langle e^{i n ( \phi _1 +\phi _2-\phi _3-\phi _4)} \Big \rangle - 2 \Big \langle e^{i n ( \phi _1 -\phi _2)} \Big \rangle ^2, \end{aligned}$$

(18)

where $\langle \cdots \rangle$ denotes the average over all events and particles. For convenience we define the $n{\text {th}}$-order $\kappa$-function

$$\begin{aligned} \kappa _n\{m\}= \int \left( \prod _{i=1}^m \frac{d^2 \mathbf{k }_i}{(2\pi )^2}\right) \frac{d^m N}{\prod _{i=1}^m d^2 \mathbf{k }_i} e^{i n (\phi _1+\cdots +\phi _{m/2}-\phi _{m/2+1}-\cdots -\phi _m)}, \end{aligned}$$

(19)

in such a way that the event average can be written as

$$\begin{aligned} \Big \langle e^{i n (\phi _1+\cdots +\phi _{m/2}-\phi _{m/2+1}-\cdots -\phi _m)} \Big \rangle = \frac{\kappa _n \{m\}}{\kappa _0 \{m\}}\ . \end{aligned}$$

(20)

Given this definition of the cumulants we can write the 2- and 4-particle Fourier harmonics of order n as

$$\begin{aligned} v_n\{2\}&=(c_n\{2\})^{1/2}, \end{aligned}$$

(21)

$$\begin{aligned} v_n\{4\}&=(-c_n\{4\})^{1/4}. \end{aligned}$$

(22)

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:

$$\begin{aligned} d_n\{2\}(p_\perp )&=\frac{{\tilde{\kappa }}_n \{2\}(p_\perp )}{{\tilde{\kappa }}_0 \{2\}(p_\perp )}\ , \end{aligned}$$

(23)

$$\begin{aligned} d_n\{4\}(p_\perp )&=\frac{{\tilde{\kappa }}_n \{4\}(p_\perp )}{{\tilde{\kappa }}_0 \{4\}(p_\perp )}-2 \frac{{\tilde{\kappa }}_n \{2\}(p_\perp )}{{\tilde{\kappa }}_0 \{2\}(p_\perp )} \ \frac{\kappa _n\{2\}}{\kappa _0\{2\}}, \end{aligned}$$

(24)

where $p_\perp =|\mathbf{p }|$ and we have defined the differential $\kappa$-functions as

$$\begin{aligned} {\tilde{\kappa }}_n \{m\}(p_\perp ) \equiv \frac{d\kappa _n \{m\}}{p_\perp d p_\perp }=\int _0^{2 \pi } d\phi _1 \int \left( \prod _{i=2}^m \frac{d^2 \mathbf{k }_i}{(2\pi )^2}\right) \frac{d^m N}{\prod _{i=1}^m d^2 \mathbf{k }_i} \Bigg |_{\mathbf{k }_1=\mathbf{p }} e^{i n (\phi _1+\cdots +\phi _{m/2}-\phi _{m/2+1}-\cdots -\phi _m)}. \end{aligned}$$

(25)

With this prescription, the differential azimuthal harmonics are given by²

$$\begin{aligned} v_n\{2\}(p_\perp )&= [d_n\{2\}(p_\perp )]^{1/2}, \end{aligned}$$

(26)

$$\begin{aligned} v_n\{4\}(p_\perp )&= [-d_n\{4\}(p_\perp )]^{1/4}. \end{aligned}$$

(27)

3 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.

3.1 Setting up the notation

We write the reduced matrix amplitude as

$$\begin{aligned} \Lambda _i \equiv \overline{{\mathcal {M}}}^{a_i b_i}_{\lambda _i}(\mathbf{k }_i,\mathbf{q }_i), \end{aligned}$$

(28)

where $i={1,\dots ,2 n}$. 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:

$$\begin{aligned}&\mathbf{k }_{2m}=\mathbf{k }_{2m-1}\ , \end{aligned}$$

(29)

$$\begin{aligned}&\lambda _{2m}=\lambda _{2m-1}\ , \end{aligned}$$

(30)

$$\begin{aligned}&a_{2m}=a_{2m-1}\ , \end{aligned}$$

(31)

with $m=1,\dots ,n$. Thus this notation also changes the usual labelling of the gluon final momenta, $\mathbf{k }_i$, since now they are labeled by only odd numbers ($1,3,5,\dots$) instead of $1,2,3,\dots$. With these convention we can write the 2n-point function of the reduced matrix amplitudes in the simplified form

$$\begin{aligned}&\Big \langle \overline{{\mathcal {M}}}^{a_1 b_1}_{\lambda _1}(\mathbf{k }_1,\mathbf{q }_1) \overline{{\mathcal {M}}}^{b_2 a_1 \dagger }_{\lambda _1}(\mathbf{k }_1,\mathbf{q }_2) \cdots \overline{{\mathcal {M}}}^{a_{n} b_{2n-1}}_{\lambda _{n}}(\mathbf{k }_n,\mathbf{q }_{2n-1}) \overline{{\mathcal {M}}}^{b_{2n} a_{n} \dagger }_{\lambda _{n}}(\mathbf{k }_n,\mathbf{q }_{2n})\Big \rangle _{T} = \Big \langle \Lambda _1 \Lambda _2 \cdots \Lambda _{2n-1} \Lambda _{2n} \Big \rangle _{T}. \end{aligned}$$

(32)

As we have pointed out in Sect. 2.2, 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

$$\begin{aligned} \Big \langle \Lambda _1 \Lambda _2 \cdots \Lambda _{2n-1} \Lambda _{2n} \Big \rangle _{T}=\sum _{\sigma \in \Pi (\chi )} \prod _{\{\alpha ,\beta \}\in \sigma } \Big \langle \Lambda _\alpha \Lambda _\beta \Big \rangle _{T}, \end{aligned}$$

(33)

with $\chi =\{1,2,\dots ,2 n\}$ and $\Pi (\chi )$ the set of partitions of $\chi$ with disjoint pairs. On the other hand, in order to evaluate the 2-point function we use Eq. (7),

$$\begin{aligned} \Big \langle \Lambda _\alpha \Lambda _\beta \Big \rangle _{T}&= 4 \epsilon ^{i*}_{\lambda _\alpha }(\mathbf{k }_\alpha ) L^i(\mathbf{k }_\alpha ,\mathbf{q }_\alpha ) \epsilon ^{j*}_{\lambda _\beta } (\mathbf{k }_\beta ) L^j(\mathbf{k }_\beta ,\mathbf{q }_\beta ) \int _{\mathbf{y }_\alpha } e^{i(-1)^\alpha \mathbf{q }_\alpha \mathbf{y }_\alpha + i(-1)^\beta \mathbf{q }_\beta \mathbf{y }_\beta } \Big \langle U^{a_\alpha b_\alpha }(\mathbf{y }_\alpha ) U^{a_\beta b_\beta }(\mathbf{y }_\beta ) \Big \rangle _{T}, \end{aligned}$$

(34)

where we do not write the overall sign of the equation which should be $-(-1)^{\alpha +\beta }$ 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 $\mathbf{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

$$\begin{aligned} \epsilon ^{i*}_{\lambda }(\mathbf{k }) \epsilon ^{j}_{\lambda }(\mathbf{k }) = \delta ^{ij}, \end{aligned}$$

(35)

which implies a contraction of the two Lipatov vertices with the same $\mathbf{k }$-momentum. Thus we will write directly $L^\lambda (\mathbf{k },\mathbf{q })$ in Eq. (34) instead of $\epsilon ^{i*}_{\lambda }(\mathbf{k }) L^i(\mathbf{k },\mathbf{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

$$\begin{aligned} \Big \langle U^{a_\alpha b_\alpha }(\mathbf{y }_\alpha ) U^{a_\beta b_\beta }(\mathbf{y }_\beta ) \Big \rangle _{T}&=\frac{\delta ^{a_\alpha a_\beta }\delta ^{b_\alpha b_\beta }}{(N_c^2-1)^2} \Big \langle Tr \left[ U(\mathbf{y }_\alpha ) U(\mathbf{y }_\beta )\right] \Big \rangle _{T} \equiv \frac{\delta ^{a_\alpha a_\beta }\delta ^{b_\alpha b_\beta }}{N_c^2-1} D(\mathbf{y }_\alpha ,\mathbf{y }_\beta ), \end{aligned}$$

(36)

where we have introduced the dipole operator $D(\mathbf{x },\mathbf{y })$.

Therefore, making the change of variables $\mathbf{y }_{\alpha ,\beta }=\mathbf{b } \pm \mathbf{r }/2$ we can write Eq. (34) as

$$\begin{aligned} \Big \langle \Lambda _\alpha \Lambda _\beta \Big \rangle _{T}&= 4 L^{\lambda _\alpha }(\mathbf{k }_\alpha ,\mathbf{q }_\alpha )L^{\lambda _\beta }(\mathbf{k }_\beta ,\mathbf{q }_\beta ) \frac{\delta ^{a_\alpha a_\beta }\delta ^{b_\alpha b_\beta }}{N_c^2-1} \int _{\mathbf{r },\mathbf{b }} e^{i\mathbf{b }[(-1)^\alpha \mathbf{q }_\alpha +(-1)^\beta \mathbf{q }_\beta ]+i\mathbf{r }/2[(-1)^\alpha \mathbf{q }_\alpha -(-1)^\beta \mathbf{q }_\beta ]}D(\mathbf{r },\mathbf{b }). \end{aligned}$$

(37)

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(\mathbf{r },\mathbf{b }) = D(|\mathbf{r }|)$. Thus, defining the Fourier transform of the dipole operator

$$\begin{aligned} d(\mathbf{q })=\int _{\mathbf{r }} e^{-i \mathbf{q } \cdot \mathbf{r }} D(|\mathbf{r }|), \end{aligned}$$

(38)

we obtain our final expression for the 2-point function of the reduced matrix amplitude:

$$\begin{aligned} \Big \langle \Lambda _\alpha \Lambda _\beta \Big \rangle _{T}&=4\frac{\delta ^{a_\alpha a_\beta }\delta ^{b_\alpha b_\beta }}{N_c^2-1} (2 \pi )^2 \delta ^{(2)}[\mathbf{q }_\alpha +(-1)^{\alpha +\beta }\mathbf{q }_\beta ] L^{\lambda _\alpha }(\mathbf{k }_\alpha ,\mathbf{q }_\alpha )L^{\lambda _\beta }(\mathbf{k }_\beta ,\mathbf{q }_\beta ) d(\mathbf{q }_\alpha ). \end{aligned}$$

(39)

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

$$\begin{aligned} g \rho ^{b_i}(\mathbf{k }_i-\mathbf{q }_i) \equiv \Omega _i\ , \end{aligned}$$

(40)

we can write the 2n-point function as

$$\begin{aligned}&g^{2n} \Big \langle \rho ^{b_1}(\mathbf{k }_1-\mathbf{q }_1)\rho ^{b_2 \dagger }(\mathbf{k }_1-\mathbf{q }_2) \cdots \,\rho ^{b_{2n-1}}(\mathbf{k }_n-\mathbf{q }_{2n-1})\rho ^{b_{2n} \dagger }(\mathbf{k }_n-\mathbf{q }_{2n}) \Big \rangle _{p} = \Big \langle \Omega _1 \Omega _2 \cdots \Omega _{2n-1} \Omega _{2n} \Big \rangle _{p}. \end{aligned}$$

(41)

Here, as in Eq. (32), even indices correspond to complex conjugates.

This correlator, Eq. (41), has the following Wick expansion:

$$\begin{aligned} \Big \langle \Omega _1 \Omega _2 \cdots \Omega _{2n-1} \Omega _{2n} \Big \rangle _{p}=\sum _{\omega \in \Pi (\chi )} \prod _{\{i,j\}\in \omega } \Big \langle \Omega _i \Omega _j \Big \rangle _{p}. \end{aligned}$$

(42)

In the generalised MV model this 2-point function can be written as

$$\begin{aligned} \Big \langle \Omega _i \Omega _j \Big \rangle _{p}=\frac{\delta ^{b_i b_j}}{N_c^2-1} \mu ^2\left[ \mathbf{k }_i-\mathbf{q }_i,(-1)^{i+j}(\mathbf{k }_j-\mathbf{q }_j)\right] , \end{aligned}$$

(43)

where $\mu ^2(\mathbf{k },\mathbf{q })$ is a function peaked around $\mathbf{k }+\mathbf{q }=0$. In the strict MV model we have that $\mu ^2(\mathbf{k },\mathbf{q })\propto \delta ^{(2)}(\mathbf{k }+\mathbf{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:

$$\begin{aligned} 2^n (2 \pi )^{3n} \frac{d^n N}{\prod _{i=1}^{n}d^2 \mathbf{k }_i}=&\,\int \left( \prod _{i=1}^{2n} \frac{d^2 \mathbf{q }_i}{(2 \pi )^2}\right) \left( \sum _{\sigma \in \Pi (\chi )} \prod _{\{i,j\}\in \sigma } \Big \langle \Omega _i \Omega _j \Big \rangle _{p}\right) \left( \sum _{\omega \in \Pi (\chi )} \prod _{\{\alpha ,\beta \}\in \omega } \Big \langle \Lambda _\alpha \Lambda _\beta \Big \rangle _{T}\right) , \end{aligned}$$

(44)

that, together with Eqs. (39) and (43), provides the full expression that will be used along this work.

3.2 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:

(45)

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³:

(46)

As an illustrative example let us select one of the $5!!^2=225$ terms that appear in Eq. (44) when $n=3$:

$$\begin{aligned}&\int \left( \prod _{i=1}^{6} \frac{d^2 \mathbf{q }_i}{(2 \pi )^2}\right) \Big \langle \Omega _1 \Omega _6 \Big \rangle _{p} \Big \langle \Omega _3 \Omega _4 \Big \rangle _{p} \Big \langle \Omega _2 \Omega _5 \Big \rangle _{p} \Big \langle \Lambda _1 \Lambda _5 \Big \rangle _{T} \Big \langle \Lambda _3 \Lambda _4 \Big \rangle _{T} \Big \langle \Lambda _2 \Lambda _6 \Big \rangle _{T}. \end{aligned}$$

(47)

This term will be represented by the following diagram:

(48)

Note that the integration over the $\mathbf{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 $\mathbf{k }_i$ with odd indices, and Eqs. (29) to (31))

(49)

with $\int _\mathbf{q }\equiv \int d^2\mathbf{q }/(2\pi )^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_c^2-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, $2 m \leftrightarrow 2m-1$, of a given diagram is equivalent to make the change of variables $\mathbf{k }_{2m-1} \rightarrow - \mathbf{k }_{2m-1}$. For example,

  

  (50)

  In order to prove this property it is enough to evaluate $\Big \langle \Omega _{2m} \Omega _i \Big \rangle _{p} \Big \langle \Omega _{2m-1} \Omega _j \Big \rangle _{p} \Big \langle \Lambda _{2m} \Lambda _\gamma \Big \rangle _{T} \Big \langle \Lambda _{2m-1} \Lambda _\beta \Big \rangle _{T}$, being $i,j,\gamma$ and $\beta$ 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 $\mathbf{q }_{2m} \rightarrow - \mathbf{q }_{2m-1}$ and $\mathbf{q }_{2m-1} \rightarrow - \mathbf{q }_{2m}$ and uses the fact that $L^\lambda (\mathbf{k },-\mathbf{q })=-L^\lambda (-\mathbf{k },\mathbf{q })$ and $\mu ^2(-\mathbf{k },-\mathbf{q })=\mu ^2(\mathbf{k },\mathbf{q })$, we see that

  $$\begin{aligned}&\Big \langle \Omega _{2m} \Omega _i \Big \rangle _{p} \Big \langle \Omega _{2m-1} \Omega _j \Big \rangle _{p} \Big \langle \Lambda _{2m} \Lambda _\gamma \Big \rangle _{T} \Big \langle \Lambda _{2m-1} \Lambda _\beta \Big \rangle _{T} \nonumber \\&\quad = \Big \langle \Omega _{2m-1} \Omega _i \Big \rangle _{p} \Big \langle \Omega _{2m} \Omega _j \Big \rangle _{p} \Big \langle \Lambda _{2m-1} \Lambda _\gamma \Big \rangle _{T} \Big \langle \Lambda _{2m} \Lambda _\beta \Big \rangle _{T} (\mathbf{k }_{2m-1} \rightarrow - \mathbf{k }_{2m-1}). \end{aligned}$$

  (51)

• (ii) Interchanging two columns $(2 m,2 m-1)$ and $(2 k, 2 k-1)$ of a given diagram is equivalent to make the change of variables $\mathbf{k }_{2m-1} \leftrightarrow \mathbf{k }_{2k-1}$.

  For example,

  

  (52)

  The proof of this property is trivial since, by definition, each column of dots corresponds to a momentum $\mathbf{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_c^2-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_c^2-1)^{-1}$ of this diagram by writing the Kronecker deltas

  

  (53)

  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

  

  (54)

  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)),

  

  (55)

  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_c^2-1)$ of a given diagram for general n is

  $$\begin{aligned} \frac{(N_c^2-1)^{n_p}(N_c^2-1)^{n_T}}{(N_c^2-1)^{2n}}=(N_c^2-1)^{n_p+n_T-2n}. \end{aligned}$$

  (56)

  As we will see through this work, this property is useful for organising the terms of Eq. (44) in powers of $(N_c^2-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.

  For example,

  

  (57)

  in this diagram the 3rd produced gluon is uncorrelated, gluons 1-2 and 2-4 have an HBT correlation and gluons 1-4-5 have a Bose enhancement correlation.

  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

  $$\begin{aligned} \Big \langle \Omega _{a} \Omega _b \Big \rangle _{p} \Big \langle \Lambda _{a} \Lambda _b \Big \rangle _{T}, \end{aligned}$$

  (58)

  where $a,b=1,...,2n$ are generic dots. The objects of Eqs. (39) and (43) that contain the information of the quantum interference correlations are the Dirac deltas and the functions $\mu ^2(\mathbf{k },\mathbf{q })$ respectively (the Lipatov vertices and the dipole functions give a different kind of correlation). Thus, we can write

  $$\begin{aligned} \Big \langle \Omega _{a} \Omega _b \Big \rangle _{p} \Big \langle \Lambda _{a} \Lambda _b \Big \rangle _{T}&\propto \mu ^2\left[ \mathbf{k }_a-\mathbf{q }_a,(-1)^{a+b} (\mathbf{k }_b-\mathbf{q }_b) \right] \delta ^{(2)}\left[ \mathbf{q }_a+(-1)^{a+b}\mathbf{q }_b \right] \nonumber \\&= \mu ^2\left[ \mathbf{k }_a-\mathbf{q }_a,(-1)^{a+b} \mathbf{k }_b+\mathbf{q }_a \right] \delta ^{(2)}\left[ \mathbf{q }_a+(-1)^{a+b}\mathbf{q }_b \right] . \end{aligned}$$

  (59)

  Since $\mu ^2(\mathbf{k },\mathbf{q })$ is peaked around $\mathbf{k }=-\mathbf{q }$ this implies that we have a peak around $\mathbf{k }_a=-(-1)^{a+b} \mathbf{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 $\mu ^2$ – in fact we loose any kind of correlation since in this case the Lipatov vertices and the dipole function can be factorized.

4 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) || (43) which contain two functions that need to be modelled, $\mu ^2(\mathbf{k },\mathbf{q })$ and $d(\mathbf{q })$.

As indicated before, in the strict MV model $\mu ^2(\mathbf{k },\mathbf{q })$ is proportional to a Dirac delta. However, in order to be more realistic, we choose a smoother function that is also peaked around $\mathbf{k }+\mathbf{q }=0$, such as a Gaussian⁴:

$$\begin{aligned} \mu ^2(\mathbf{k },\mathbf{q })=e^{-\frac{(\mathbf{k }+\mathbf{q })^2}{4 B_p^{-1}}}, \end{aligned}$$

(60)

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]:

$$\begin{aligned} d(\mathbf{q })=\frac{4 \pi }{Q_s^2} e^{-\frac{\mathbf{q }^2}{Q_s^2}}. \end{aligned}$$

(61)

We should also account for the infrared divergences of the Lipatov vertices. The product of two Lipatov vertices is

$$\begin{aligned} L^i(\mathbf{k },\mathbf{q }_1)L^i(\mathbf{k },\mathbf{q }_2)=\left[ \frac{\mathbf{k }^i}{\mathbf{k }^2}-\frac{(\mathbf{k }-\mathbf{q }_1)^i}{(\mathbf{k }-\mathbf{q }_1)^2}\right] \left[ \frac{\mathbf{k }^i}{\mathbf{k }^2}-\frac{(\mathbf{k }-\mathbf{q }_2)^i}{(\mathbf{k }-\mathbf{q }_2)^2}\right] . \end{aligned}$$

(62)

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:

$$\begin{aligned} L^i(\mathbf{k },\mathbf{q }_1)L^i(\mathbf{k },\mathbf{q }_2)=\frac{(2 \pi )^2}{\xi ^2} \exp \left\{ -\frac{[\mathbf{k }-(\mathbf{q }_1+\mathbf{q }_2)/2]^2}{\xi ^2}\right\} , \end{aligned}$$

(63)

where $\xi ^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

$$\begin{aligned} W^{b_1 b_2 b_3 b_4}(\mathbf{b }_1,\mathbf{p }_1,\mathbf{b }_2,\mathbf{p }_2)&=\frac{1}{(N_c^2-1)^2}\frac{1}{\pi ^4 \xi ^4 B_p^2} e^{-(\mathbf{p }_1^2+\mathbf{p }_2^2)/\xi ^2} e^{-(\mathbf{b }_1^2+\mathbf{b }_2^2)/B_p} \Big [ \delta ^{b_1 b_2} \delta ^{b_3 b_4}+\delta ^{b_1 b_3} \delta ^{b_2 b_4} 2 \pi B_p \delta ^{(2)}(\mathbf{b }_1-\mathbf{b }_2)\nonumber \\&\quad \times e^{-(\mathbf{p }_1+\mathbf{p }_2)^2/(2 B_p^{-1})} +\delta ^{b_1 b_4} \delta ^{b_2 b_3} 2 \pi B_p \delta ^{(2)}(\mathbf{b }_1-\mathbf{b }_2) e^{-(\mathbf{p }_1-\mathbf{p }_2)^2/(2 B_p^{-1})} \Big ], \end{aligned}$$

(64)

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 $\rho$’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, $\mathbf{k }_i-\mathbf{q }_i$, and not on the final momentum, $\mathbf{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⁵. 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.

4.1 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

(65)

where we have grouped the 9 diagrams by their powers in $(N_c^2-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 $\mathbf{q }_i$ integrals are Gaussian functions and, therefore, they can be trivially solved. The result is

(66)

(67)

(68)

(69)

(70)

With these five 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 $\mathbf{k }_i$ obtaining

$$\begin{aligned} N&=\kappa _0\{2\}\nonumber \\&=1+\frac{2}{N_c^2-1}\left[ \frac{2}{1+B_p \xi ^2} \right] +\frac{2}{(N_c^2-1)^2}\left[ \frac{1}{1+B_p Q_s^2}+\frac{1}{1+B_p \xi ^2} \right] . \end{aligned}$$

(71)

We can see from this equation that, apart from the suppression in powers of $(N_c^2-1)^{-1}$, the correlated terms contain suppression factors $(1+B_p Q_s^2)^{-1}$ and $(1+B_p \xi ^2)^{-1}$. Following the domain picture that we have discussed in Sect. 2.3, $B_p Q_s^2 \equiv 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 \rightarrow \infty$ since the probability of two gluons scattering off the same domain vanishes in this limit. Therefore, to fix $\xi ^2$ it makes sense to choose a value that is proportional to $Q_s^2$ in order to preserve decorrelation in the limit $n_D \rightarrow \infty$. For this reason we will choose $\xi ^2 = \alpha Q_s^2$, being $\alpha$ a real number, in the rest of this work.⁶

The 2-particle azimuthal harmonics, Eq. (17), can be obtained by performing the integration over $\mathbf{k }_i$ with the help of Eq. (B12). The result for the second order $\kappa$-function is

$$\begin{aligned} \kappa _{2n} \{ 2 \}&= \frac{8 \Gamma \left( n+1\right) ^2}{ \Gamma (2n+1)} \Bigg \{\frac{\alpha (1+\alpha )}{N_c^2-1} \times \Bigg [ \frac{(1+\alpha +\alpha n_D)}{\alpha ^4 n_D^2} \left( \frac{\alpha ^2 n_D}{2+\alpha ^2 n_D+2\alpha (1+n_D)}\right) ^{2(n+1)} \nonumber \\&\quad \times \, _2F_1\left( n+1,n+1;2n+1;\left( \frac{\alpha ^2 n_D}{2+\alpha ^2 n_D+2\alpha (1+n_D)}\right) ^2 \right) \nonumber \\&\quad +\frac{1}{(1+\alpha ^2 n_D+\alpha (2+n_D))^2} \left( \frac{1+\alpha n_D+\alpha ^2 n_D}{1+\alpha ^2 n_D+\alpha (2+n_D)}\right) ^{2n} \nonumber \\&\quad \times \, _2F_1\left( n+1,n+1;2n+1;\left( \frac{1+\alpha n_D+\alpha ^2 n_D}{1+\alpha ^2 n_D+\alpha (2+n_D)}\right) ^2 \right) \Bigg ] \nonumber \\&\quad +\frac{1}{(N_c^2-1)^2}\Bigg [ \frac{\alpha (1+\alpha +\alpha n_D)}{(1+\alpha ^2 n_D+\alpha (2+n_D))^2} \left( \frac{1+\alpha n_D-\alpha ^2 n_D}{1+\alpha ^2 n_D+ \alpha (2+n_D)}\right) ^{2 n} \nonumber \\&\quad \times \, _2F_1\left( n+1,n+1;2n+1;\left( \frac{1+\alpha n_D-\alpha ^2 n_D}{1+\alpha ^2 n_D+ \alpha (2+n_D)} \right) ^2 \right) \nonumber \\&\quad +\frac{\alpha (1+\alpha )}{1+n_D} (1+2\alpha )^{-2(n+1)} \, _2F_1\left( n+1,n+1;2n+1;\frac{1}{(1+2\alpha )^2}\right) \Bigg ] \Bigg \}, \end{aligned}$$

(72)

where we have defined $\alpha =\xi ^2/Q_s^2$, we have taken $n>0$ and due to the symmetry $\mathbf{k }_3 \rightarrow -\mathbf{k }_3$ of Eq. (65) all odd harmonics vanish. Using Eqs. (71) and (72) we can evaluate the 2-particle azimuthal harmonics as

$$\begin{aligned} v_{2n} \{2\} = \sqrt{\frac{\kappa _{2n} \{2\}}{\kappa _{0} \{2\}}}. \end{aligned}$$

(73)

In Fig. 4 we plot the dependence of $v_{2n}\{2\}$ with respect to $n_D$ and $\alpha$ by fixing $N_c=3$. The value of the even azimuthal harmonics grows rapidly as both $n_D$ and $\alpha$ 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 $\alpha$ must be taken with care because $\alpha$ gives the ratio between the momentum transfers from projectile and target. The dilute-dense approximation that we are using makes sense only for $\alpha$ sizeably smaller than 1.

Fig. 4

Dependence of the even 2-particle azimuthal harmonics, $v_{2n}\{2\}$, on $\alpha \equiv \xi ^2/Q_s^2$ (left, for $B_p Q_s^2=12$) and $n_D \equiv B_p Q_s^2$ (right, for $\xi ^2/Q_s^2=1/4$)

As we have seen in Sect. 2.3, we can also compute the azimuthal harmonics as a function of transverse momentum by using the differential $\kappa$-function defined in Eq. (25). The $\mathbf{k }_i$ integral can be solved with the help of Eq. (B11) and the result is

(74)

and

(75)

where $n>0$. The fact that ${\tilde{\kappa }}_{2n}(p_\perp )$ is proportional to $(p_\perp ^2/Q_s^2)^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

$$\begin{aligned} v_{2n}\{2\}(p_\perp )=\sqrt{\frac{{\tilde{\kappa }}_{2n}\{2\}(p_\perp )}{{\tilde{\kappa }}_0\{2\}(p_\perp )}} \end{aligned}$$

(76)

and the result is plotted in Fig. 5 for $n=1,2,3$ and $B_p= 6$ $\hbox {GeV}^{-2}$, $\xi =Q_s/2$, $Q_s^2=2$ $\hbox {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_\perp$ sizeably larger than $Q_s$.

Fig. 5

Dependence of the differential 2-particle even azimuthal harmonics, $v_{2n}\{2\}$, on transverse momentum $p_\perp$. In this graph we have used $B_p= 6$ $\hbox {GeV}^{-2}$, $\xi =Q_s/2$, $Q_s^2=2$ $\hbox {GeV}^2$ and $N_c=3$

4.2 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_c^2-1)^{-1}$ as

$$\begin{aligned}&8(2 \pi )^{9} \frac{d^3 N}{d^2 \mathbf{k }_1 d^2 \mathbf{k }_3d^2 \mathbf{k }_5} =N_3^{(0)}+N_3^{(1)}+N_3^{(2)}+N_3^{(3)}+N_3^{(4)}. \end{aligned}$$

(77)

In order to obtain each one of these terms we have to use the property iii) of Sect. 3.2. In this case the suppression of each diagram is given by $(N_c^2-1)^{n_p+n_T-6}$ and we can have three kind of configurations on each side of the diagram

(78)

(79)

(80)

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_c^2-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) on the right side of the diagram and configurations on the left side that have one link in common with the right side and the other ones different. There are $6\times 6$ possibilities in this case:

  

  and the 5 permutations of Eq. (79) for each diagram.

• (iii) $n_p=3$ and $n_T=1$.

  This implies that we have to have the configuration of Eq. (80) on the right side of the diagram and configurations on the left side that have all the links in common with the right side. There are $1 \times 8$ possibilities in this case:

  

  and the 7 permutations of Eq. (80).

All in all, we can write all the 52 Wick diagrams that have a suppression of $(N_c^2-1)^{-2}$ as

(81)

This procedure, although tedious, is straightforward to implement on a computer. Repeating it we find that there is 1 diagram suppressed by $(N_c^2-1)^{0}$:

(82)

12 diagrams suppressed by $(N_c^2-1)^{-1}$:

(83)

96 diagrams suppressed by $(N_c^2-1)^{-3}$:

(84)

and 64 diagrams suppressed by $(N_c^2-1)^{-4}$:

(85)

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_c^2-1)^{-1}$. These equations were computed in [53] up to order $(N_c^2-1)^{-2}$. In order to write Eq. (77) just as a function of $\mathbf{k }_i$ one just has to employ Eqs. (39), (43), (60), (61) and (63) and then perform the $\mathbf{q }_i$ integrals. These integrals are straightforward if the arguments of the integrals are Gaussian functions, as we assumed before.

4.3 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

(86)

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:

(87)

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 three 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_c^2-1)^{-2}$.

The Wick diagrams of Eq. (87) can be computed in the same fashion as in Sect. 4.1. However, since we are only interested in computing the 4-particle cumulants, Eq. (18), we will exploit the $\mathbf{k }_i \leftrightarrow \mathbf{k }_j$ and $\mathbf{k }_i \rightarrow -\mathbf{k }_i$ symmetries in order to simplify the calculation. When evaluating the 4-particle $\kappa$-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 $\kappa$-function when $n>0$ and therefore we can write

(88)

where we have written schematically $\text {perm}_4$, which encodes all the factorizable connected diagrams, and $\text {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\pi$ 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 $\text {perm}_5$, that all the diagrams that are related by a change of variables $\mathbf{k }_i \rightarrow -\mathbf{k }_i$ will give the same result for the integral in Eq. (88) since this transformation is equivalent to making $\phi _i \rightarrow \phi _i + \pi$ 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 $\mathbf{k }_1 \leftrightarrow \mathbf{k }_3$, $\mathbf{k }_5 \leftrightarrow \mathbf{k }_7$ and $\mathbf{k }_3 \leftrightarrow \mathbf{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

(89)

where the last diagram can be seen as the first one with the change of variables $\mathbf{k }_1 \leftrightarrow \mathbf{k }_7$.

Furthermore, 4 out of the 12 diagrams of the fourth term of Eq. (87)

(90)

only depend on $\phi _1-\phi _3$ and $\phi _5-\phi _7$ and therefore vanish due to rotational invariance. Having this into account we can write Eq. (88) as

(91)

On the other hand, the 2-particle $\kappa$-function in the case in which the partons are initially uncorrelated in the projectile wave function is

(92)

Therefore we can write Eq. (91) as

(93)

In order to compute $\kappa _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 $\mathbf{k }_i \rightarrow -\mathbf{k }_i$ or $\mathbf{k }_i \leftrightarrow \mathbf{k }_j$ ($i \ne j$) that leave the integral invariant we can write

(94)

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)

(95)

and the four particle even azimuthal harmonics is obtained as

$$\begin{aligned} v_{2n}\{4\}=(-c_{2n}\{4\})^{1/4}. \end{aligned}$$

(96)

In Fig. 6 left we have plotted our results for Eq. (95) as a function of $Q_s^2$. The absolute values of the cumulant are very small and it even becomes positive with increasing $Q_s^2$. The reason why it is so comes from the fact that we are not having into account all the contributions that come from the correlation of the partons in the projectile ensemble. Below, see Eq. (7), these contributions are taken into account and the values are reasonable and in the ballpark of the ones in experimental data.

Fig. 6

Dependence of the 4-particle integrated cumulants on $Q_s^2$ (left) and of the differential cumulants on $p_\perp$ (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\pi$)

(97)

and, therefore,

(98)

In order to compute ${\tilde{\kappa }}_{0}\{4\}(p_\perp )$ we cannot use the same symmetries that we employed for computing $\kappa _{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 $\mathbf{k }_i \rightarrow - \mathbf{k }_i$ still leave the integral invariant but now, since we are fixing $|\mathbf{k }_1|=p_\perp$, all the diagrams that are related by a change of variable $\mathbf{k }_1 \leftrightarrow \mathbf{k }_j$ will give a different value for the integral but the ones that are related by $\mathbf{k }_i \leftrightarrow \mathbf{k }_j$, with $i,j \ne 1$, still leave the integral invariant. Performing a simple counting of the permutations of Eq. (87) we can write

(99)

Thus, using Eqs. (98) and (99) the differential 4-particle cumulant can be written in a similar form as (95):

(100)

and the differential 4-particle azimuthal harmonics is defined as

$$\begin{aligned} v_n\{4\}(p_\perp )=(-d_{2n}\{4\}(p_\perp ))^{1/4}. \end{aligned}$$

(101)

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_\perp$ 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_c^2-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 $\kappa _n\{4\}$ can be written as $2 \kappa _n\{2\}^2$. Finally, we solved numerically these integrals for given values of $Q_s^2$ 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_c^2-1)^{-1}$, by using the property iii) of Sect. 3.2. 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_s^2$, the differential cumulants as a function of $p_\perp$ 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_\perp$ 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_s^2 \rightarrow 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. 7

Dependence of the 4-particle cumulant (left) and azimuthal harmonic (right) of second and fourth order with $Q_s^2$. In these graphs we have used $N_c=3$ and the values of the remaining parameters are indicated on the plots

Fig. 8

Dependence of the differential 4-particle cumulant (left) and azimuthal harmonic (right) of second and fourth order with $p_\perp$. 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 Sect. 3.2 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.

5 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 10,000, 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_\perp$. 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_\perp$ 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,$\dots$. 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.

Appendices

Appendix A: On the validity of the area enhancement argument

In this section we will study the validity of the area enhancement argument, from now on AE model, introduced in Sect. 2.2. 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

$$\begin{aligned} \langle {\mathcal {D}}(\mathbf{x },\mathbf{y }) {\mathcal {D}}(\mathbf{u },\mathbf{v }) \rangle _T^{\text {AE}} = D(\mathbf{x },\mathbf{y }) D(\mathbf{u },\mathbf{v })+\frac{1}{N_c^2} D(\mathbf{x },\mathbf{v }) D(\mathbf{u },\mathbf{y }), \end{aligned}$$

(A1)

where we have introduced the dipole operator ${\mathcal {D}}(\mathbf{x },\mathbf{y })=\frac{1}{N_c}\text {Tr}[U(\mathbf{x })U^\dagger (\mathbf{y })]$ and $D(\mathbf{x },\mathbf{y })$ is its target average.

As discussed in Sect. 2.2, 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

$$\begin{aligned} \langle {\mathcal {D}}(\mathbf{x },\mathbf{y }) {\mathcal {D}}(\mathbf{u },\mathbf{v }) \rangle _T^{\text {MV}}&= D(\mathbf{x },\mathbf{y }) D(\mathbf{u },\mathbf{v }) \left[ \left( \frac{F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })+\sqrt{\Delta }}{2 \sqrt{\Delta }} - \frac{1}{N_c^2} \frac{F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })}{\sqrt{\Delta }} \right) e^{\frac{N_c}{4} \mu ^2 \sqrt{\Delta }}\right. \nonumber \\&\quad -\left. \left( \frac{F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })-\sqrt{\Delta }}{2 \sqrt{\Delta }} - \frac{1}{N_c^2} \frac{F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })}{\sqrt{\Delta }} \right) e^{-\frac{N_c}{4} \mu ^2 \sqrt{\Delta }} \right] e^{-\frac{N_c}{4}\mu ^2 F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })+\frac{1}{2 N_c} \mu ^2 F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })}, \end{aligned}$$

(A2)

where

$$\begin{aligned} \Delta = F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })^2+\frac{4}{N_c^2} F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v }) F(\mathbf{x },\mathbf{v };\mathbf{u },\mathbf{y }) \end{aligned}$$

(A3)

and function $F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })$ is defined in [81]. In the GBW model this function reads simply

$$\begin{aligned} \mu ^2 F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v }) = \frac{Q_s^2}{2 C_F} (\mathbf{x }-\mathbf{y }) \cdot (\mathbf{u }-\mathbf{v }). \end{aligned}$$

(A4)

Taking the large-$N_c$ limit we can simplify Eq. (A2) drastically to read, in the GBW model,

$$\begin{aligned} \langle {\mathcal {D}}(\mathbf{x },\mathbf{y }) {\mathcal {D}}(\mathbf{u },\mathbf{v }) \rangle _T^{\text {MV}}&= D(\mathbf{x },\mathbf{y })D(\mathbf{u },\mathbf{v }) +\frac{1}{N_c^2} \frac{F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })^2}{F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })^2}\Bigg [ D(\mathbf{x },\mathbf{v })D(\mathbf{u },\mathbf{y }) +D(\mathbf{x },\mathbf{y })D(\mathbf{u },\mathbf{v }) \nonumber \\&\quad \times \left( \frac{Q_s^2}{2}(\mathbf{u }-\mathbf{x }) \cdot (\mathbf{v }-\mathbf{y })-1\right) \Bigg ] + {\mathcal {O}}\left( \frac{1}{N_c^3}\right) . \end{aligned}$$

(A5)

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 $\Psi (\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v })$ computed in the MV and AE models is

$$\begin{aligned}&\frac{\int _{\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }} \langle {\mathcal {D}}{\mathcal {D}} \rangle ^{\text {MV}} \Psi (\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v })}{\int _{\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }} \langle {\mathcal {D}}{\mathcal {D}} \rangle ^{\text {AE}} \Psi (\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v })}\nonumber \\&\quad =1 + \frac{1}{N_c^2} \int _{\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }} \Psi (\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }) \left[ \left( \frac{F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })^2}{F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })^2} - 1 \right) D(\mathbf{x },\mathbf{v })D(\mathbf{u },\mathbf{y })+ \frac{F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })^2}{F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })^2} \right. \nonumber \\&\qquad \times \left. \left( \frac{Q_s^2}{2}(\mathbf{u }-\mathbf{x }) \cdot (\mathbf{v }-\mathbf{y }) - 1 \right) D(\mathbf{x },\mathbf{y })D(\mathbf{u },\mathbf{v })\right] {\bigg /} \int _{\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }} D(\mathbf{x },\mathbf{y })D(\mathbf{u },\mathbf{v })\Psi (\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }) + {\mathcal {O}}\left( \frac{1}{N_c^3}\right) . \end{aligned}$$

(A6)

Using the saddle point approximation, and noting that $F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })\rightarrow 0$ when $\mathbf{x }\rightarrow \mathbf{y }$ or $\mathbf{u }\rightarrow \mathbf{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 $\propto 1/Q_s$, to behave $\delta$-like with respect to the integration area. Therefore, corrections must be order $1/Q_s^2$ 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

$$\begin{aligned} d\Omega = d^2\mathbf{x } d^2\mathbf{y } d^2\mathbf{u } d^2\mathbf{v } \Theta \left( \sqrt{2 B_p}-|\mathbf{x }|\right) \Theta \left( \sqrt{2 B_p}-|\mathbf{y }|\right) \Theta \left( \sqrt{2 B_p}-|\mathbf{u }|\right) \Theta \left( \sqrt{2 B_p}-|\mathbf{v }|\right) ; \end{aligned}$$

(A7)

that is, we integrate over a 4-sphere of radius $\sqrt{2 B_p}$ in such a way that the integral over the phase space leads to the expected result

$$\begin{aligned} \int d\Omega = (2 \pi B_p)^4 = S_\perp ^4. \end{aligned}$$

(A8)

In order to compare the MV and AE models, we perform a Fourier transform over the phase space measure defined as

$$\begin{aligned} \langle {\mathcal {D}}(\mathbf{q }_1,\mathbf{q }_2) {\mathcal {D}}(\mathbf{q }_3,\mathbf{q }_4) \rangle _T&= \int d\Omega e^{i \mathbf{q }_1 \cdot \mathbf{x }-i \mathbf{q }_2 \cdot \mathbf{y }+i \mathbf{q }_3 \cdot \mathbf{u }-i \mathbf{q }_4 \cdot \mathbf{v }} \langle {\mathcal {D}}(\mathbf{x },\mathbf{y }) {\mathcal {D}}(\mathbf{u },\mathbf{v }) \rangle _T. \end{aligned}$$

(A9)

In Fig. 9 we show the result for the ratio of the Fourier transforms of Eqs. (A2) and (A1) for different values of $B_p$, taking $Q_s^2=1$ $\hbox {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_s^2$ 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$.

Fig. 9

Ratio of the Fourier transforms of Eqs. (A2) and (A1) at different values of $B_p$. The values of the ratio were computed using four sets of random momenta with moduli between 0.5 and 1.5 GeV. We present both the norm (blue lines) and the argument (yellow lines). We have suppressed the values where the estimated error in the Monte Carlo integration becomes larger than 10%

An analogous discussion can be performed to the expected value of the quadrupole operator. In the AE model it can be written as

$$\begin{aligned} \langle {\mathcal {Q}}(\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }) \rangle _T^{\text {AE}} = D(\mathbf{x },\mathbf{y }) D(\mathbf{u },\mathbf{v })+ D(\mathbf{x },\mathbf{v }) D(\mathbf{u },\mathbf{y }), \end{aligned}$$

(A10)

where ${\mathcal {Q}}(\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v })=\frac{1}{N_c}\text {Tr}[U(\mathbf{x })U^\dagger (\mathbf{y })U(\mathbf{u })U^\dagger (\mathbf{v })]$. In the MV model it reads [96]

$$\begin{aligned} \langle {\mathcal {Q}}(\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }) \rangle _T^{\text {MV}}&= D(\mathbf{x },\mathbf{y }) D(\mathbf{u },\mathbf{v }) \left[ \left( \frac{F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })+\sqrt{\Delta }}{2 \sqrt{\Delta }} - \frac{F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })}{\sqrt{\Delta }} \right) e^{\frac{N_c}{4} \mu ^2 \sqrt{\Delta }}\right. \nonumber \\&\quad \left. -\left( \frac{F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })-\sqrt{\Delta }}{2 \sqrt{\Delta }} - \frac{F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })}{\sqrt{\Delta }} \right) e^{-\frac{N_c}{4} \mu ^2 \sqrt{\Delta }} \right] e^{-\frac{N_c}{4}\mu ^2 F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v }) +\frac{1}{2 N_c} \mu ^2 F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })}. \end{aligned}$$

(A11)

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_c^2$.

In the large-$N_c$ limit, Eq. (A11) can be simplified to read

$$\begin{aligned} \langle {\mathcal {Q}}(\mathbf{x },\mathbf{y },\mathbf{u },\mathbf{v }) \rangle _T^{\text {MV}} = D(\mathbf{x },\mathbf{y })D(\mathbf{u },\mathbf{v }) - \frac{F(\mathbf{x },\mathbf{y };\mathbf{u },\mathbf{v })}{F(\mathbf{x },\mathbf{u };\mathbf{y },\mathbf{v })} \Big [D(\mathbf{x },\mathbf{y })D(\mathbf{u },\mathbf{v }) - D(\mathbf{x },\mathbf{v })D(\mathbf{u },\mathbf{y })\Big ] + {\mathcal {O}}\left( \frac{1}{N_c^2}\right) . \end{aligned}$$

(A12)

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_c^2$ 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. (A11) and (A10) for three sets of random momenta with moduli between 0.5 and 1.5 GeV for different values of $B_p$, taking again $Q_s^2=1$ $\hbox {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_c^2$ 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_s^2$. Therefore, the tendency is not as clear as in Fig. 9.

Fig. 10

Ratio of the Fourier transform of Eqs. (A11) and (A10) at different values of $B_p$. The values of the ratio were computed using three sets of random momenta with moduli between 0.5 and 1.5 GeV. We present both the norm (blue lines) and the argument (yellow lines). We have suppressed the values where the estimated error in the Monte Carlo integration becomes larger than 10%

Appendix B: Table of integrals

In this section we present all the closed-form solutions of the integrals that have been used in this work. First, let us introduce the well known definitions

$$\begin{aligned}&I_n(x)=\sum _{k=0}^{\infty } \left( \frac{x}{2}\right) ^{2k+n}\frac{1}{k!\Gamma (k+n+1)}, \end{aligned}$$

(B1)

$$\begin{aligned}&{}_{p}F_{q}(a_{1},\ldots ,a_{p};b_{1},\ldots ,b_{q};z) =\sum _{n=0}^{\infty }\frac{(a_{1})_{n}\cdots (a_{p})_{n}}{(b_{1})_{n}\ldots (b_{q})_{n}}\,\frac{z^{n}}{n!}, \end{aligned}$$

(B2)

$$\begin{aligned}&(a)_{n}={\frac{\Gamma (a+n)}{\Gamma (a)}}, \end{aligned}$$

(B3)

$$\begin{aligned}&\Gamma (k)=\int _0^{\infty } dx e^{-x}x^{k-1}, \end{aligned}$$

(B4)

where $I_n(x)$ is the modified Bessel function of the first kind, $\,{}_{p}F_{q}(a_{1},\ldots ,a_{p};b_{1},\ldots ,b_{q};z)$ is the generalised hypergeometric function, $(a)_{n}$ is the rising factorial (or Pochhammer symbol) and $\Gamma (k)$ is the gamma function.

We will also introduce for convenience the Jacobi-Anger expansion

$$\begin{aligned} e^{x \cos \phi }=\sum _{n=-\infty }^{\infty }I_n(x)e^{-i n \phi }. \end{aligned}$$

(B5)

The first integral and the one that will be the most used in this work is the well known Gaussian integral

$$\begin{aligned} \int d^2 \mathbf{k }\ e^{-A \mathbf{k }^2+\mathbf{B } \cdot \mathbf{k }}=\frac{\pi }{A}e^{\frac{\mathbf{B }^2}{4A}}. \end{aligned}$$

(B6)

When computing the 2-particle differential cumulant we will have to deal with the following integral

$$\begin{aligned}&\int _0^{2\pi } d \phi _1 \int d^2 \mathbf{k }_2 \ e^{i2n(\phi _1-\phi _2)}e^{-A_1 \mathbf{k }_1^2-A_2\mathbf{k }_2^2 +A_{12}\mathbf{k }_1 \cdot \mathbf{k }_2} \Big \vert _{|\mathbf{k }_1|=p_\perp } \nonumber \\&\quad = \int _0^{2\pi } d \phi _1 d\phi _2\ e^{i2n(\phi _1-\phi _2)} \int _0^{\infty } dk_2 k_2 \ e^{-A_1 p_\perp ^2-A_2 k_2^2+A_{12} p_\perp k_2 \cos (\phi _1-\phi _2)} \nonumber \\&\quad =2 \pi e^{B} \int _0^{\infty } dk_2 k_2 \ e^{-A_2 k_2^2} \int _0^{2 \pi } d\phi \ e^{i2n\phi +C k_2 \cos \phi },\nonumber \\ \end{aligned}$$

(B7)

where we have defined $B=-A_1 p_\perp ^2$ and $C=A_{12} p_\perp$. Taking into account the Jacobi-Anger expansion Eq. (B5) and the fact that n is an integer we can write

$$\begin{aligned}&2 \pi e^{B} \int _0^{\infty } dk_2 k_2 \ e^{-A_2 k_2^2} \int _0^{2 \pi } d\phi \ e^{i2n\phi +C k_2 \cos \phi }\nonumber \\&\quad = 2 \pi e^{B} \sum _{m=-\infty }^{\infty } \int _0^{\infty } dk_2 k_2 \ e^{-A_2 k_2^2} I_m(C k_2) \int _0^{2 \pi } d\phi \ e^{i(2n-m)\phi }= (2 \pi )^2 e^{B} \int _0^{\infty } dk_2 k_2 \ e^{-A_2 k_2^2} I_{2n}(C k_2). \end{aligned}$$

(B8)

Using the definition of Eq. (B1) and then making the change of variable $x=A_2 k_2^2$ we have that

$$\begin{aligned}&(2 \pi )^2 e^{B} \int _0^{\infty } dk_2 k_2 \ e^{-A_2 k_2^2} I_{2n}(C k_2) =(2 \pi )^2 e^{B} \sum _{k=0}^{\infty } \int _0^{\infty } dk_2 k_2 \ e^{-A_2 k_2^2} \left( \frac{C k_2}{2}\right) ^{2k+2n}\frac{1}{k!\Gamma (k+2n+1)} \nonumber \\&\quad =(2 \pi )^2 e^{B} \sum _{k=0}^{\infty } \frac{1}{k!\Gamma (k+2n+1)} \int _0^{\infty } \frac{dx}{2 A_2} e^{-x} \left( \frac{C^2 x}{4 A_2}\right) ^{k+n} \nonumber \\&\quad = \frac{(2 \pi )^2}{2 A_2} e^{B} \left( \frac{C^2}{4 A_2}\right) ^{n} \sum _{k=0}^{\infty } \left( \frac{C^2}{4 A_2}\right) ^{k} \frac{1}{k!\Gamma (k+2n+1)} \int _0^{\infty } e^{-x} x^{k+n} \nonumber \\&\quad = \frac{(2 \pi )^2}{2 A_2} e^{B} \left( \frac{C^2}{4 A_2}\right) ^{n} \sum _{k=0}^{\infty } \left( \frac{C^2}{4 A_2}\right) ^{k} \frac{1}{k!\Gamma (k+2n+1)} \Gamma (k+n+1), \end{aligned}$$

(B9)

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

$$\begin{aligned}&\sum _{k=0}^{\infty } \left( \frac{C^2}{4 A_2}\right) ^{k} \frac{1}{k!\Gamma (k+2n+1)} \Gamma (k+n+1) = \frac{\Gamma (n+1)}{\Gamma (2n+1)} {}_1F_1\left( n+1;2n+1;\frac{C^2}{4 A_2}\right) \end{aligned}$$

(B10)

and, therefore, the result for the integral is

$$\begin{aligned}&\int _0^{2\pi } d \phi _1 \int d^2 \mathbf{k }_2 \ e^{i2n(\phi _1-\phi _2)} e^{-A_1 \mathbf{k }_1^2 -A_2\mathbf{k }_2^2+A_{12}\mathbf{k }_1 \cdot \mathbf{k }_2} \Big \vert _{|\mathbf{k }_1|=p_\perp } \nonumber \\&\quad = \frac{(2 \pi )^2}{2 A_2} e^{-A_1 p_\perp ^2} \left( \frac{A_{12}^2 p_\perp ^2}{4 A_2}\right) ^{n} \frac{\Gamma (n+1)}{\Gamma (2n+1)} {}_1F_1\left( n+1;2n+1;\frac{A_{12}^2 p_\perp ^2}{4 A_2}\right) . \end{aligned}$$

(B11)

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

$$\begin{aligned}&\int d^2 \mathbf{k }_1 d^2 \mathbf{k }_2 \ e^{i2n(\phi _1-\phi _2)} e^{-A_1 \mathbf{k }_1^2 -A_2\mathbf{k }_2^2+A_{12}\mathbf{k }_1 \cdot \mathbf{k }_2} = \frac{(2\pi )^2}{4A_1 A_2} \left( \frac{A_{12}^2}{4 A_1 A_2}\right) ^{n} \frac{\Gamma (n+1)^2}{\Gamma (2n+1)} \,_2F_1\left( n+1,n+1;2n+1;\frac{A_{12}^2}{4 A_1 A_2}\right) .\nonumber \\ \end{aligned}$$

(B12)

Appendix C: The Wigner function approach

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 $\mathbf{p }$ scattering on a dense target and leaving with a momentum $\mathbf{k }$ is given, at leading order, by

(C1)

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 $\mathbf{k }$ and color a coming from the projectile after the interaction with the target can be written as

(C2)

where ${\hat{\rho }}$ is the single gluon density matrix.

Using the completeness relation for the initial state we can write this equation as

(C3)

Doing the change of variables $\mathbf{q }_{1,2}=\mathbf{p } \pm \mathbf{q }/2$ and $\mathbf{x },\bar{\mathbf{x }}=\mathbf{b } \pm \mathbf{r }/2$ we get

(C4)

Defining

$$\begin{aligned} \Phi ^{b_2 b_1}(\mathbf{b },\mathbf{k }-\mathbf{p })&=\int _{\mathbf{r }} e^{-i(\mathbf{k }-\mathbf{p })\cdot \mathbf{r }} U^{ab_1} \left( \mathbf{b }+\frac{\mathbf{r }}{2}\right) U^{b_2 a}\left( \mathbf{b }-\frac{\mathbf{r }}{2}\right) ^\dagger \end{aligned}$$

(C5)

and realizing that

(C6)

is the Weyl transform of the density matrix, that is, the Wigner function, we can write the single inclusive gluon spectrum as

$$\begin{aligned}&\frac{dN}{d^2 \mathbf{k }}=\int _{\mathbf{b }}\int _{\mathbf{p }} W^{b_1 b_2}(\mathbf{b },\mathbf{p }) \Phi ^{b_2 b_1}(\mathbf{b },\mathbf{k }-\mathbf{p }). \end{aligned}$$

(C7)

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):

$$\begin{aligned} 2 (2\pi )^3 \frac{dN}{d^2 \mathbf{k }}&=\int _{\mathbf{q }_1\mathbf{q }_2} \Big \langle \rho ^{b_1}(\mathbf{q }_1)\rho ^{b_2}(\mathbf{q }_2)^* \Big \rangle _p 4 L^i(\mathbf{k },\mathbf{k }-\mathbf{q }_1) L^i(\mathbf{k },\mathbf{k }-\mathbf{q }_2) \int _{\mathbf{x } \bar{\mathbf{x }}} e^{-i(\mathbf{k }-\mathbf{q }_1) \cdot \mathbf{x } +i(\mathbf{k }-\mathbf{q }_2) \cdot \bar{\mathbf{x }}} U^{a b_1}(\mathbf{x }) U^{b_2 a}(\bar{\mathbf{x }})^\dagger . \end{aligned}$$

(C8)

Doing the same change of variables that we did before can write this expression as

$$\begin{aligned} \frac{dN}{d^2 \mathbf{k }}&=\frac{1}{2 (2\pi )^3} \int _{\mathbf{b }}\int _{\mathbf{p }} \int _{\mathbf{q }} e^{i \mathbf{q } \cdot \mathbf{b }} \Big \langle \rho ^{b_1}\left( \mathbf{p }+\frac{\mathbf{q }}{2}\right) \rho ^{b_2}\left( \mathbf{p }-\frac{\mathbf{q }}{2}\right) ^* \Big \rangle _p 4\nonumber \\&\quad \times L^i\left( \mathbf{k },\mathbf{k }-\mathbf{p } -\frac{\mathbf{q }}{2}\right) L^i\left( \mathbf{k },\mathbf{k }-\mathbf{p } +\frac{\mathbf{q }}{2}\right) \Phi ^{b_2 b_1}(\mathbf{b },\mathbf{k }-\mathbf{p }). \end{aligned}$$

(C9)

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:

$$\begin{aligned} W^{b_1 b_2}(\mathbf{b },\mathbf{p })&= \frac{1}{2 (2\pi )^3} \int _{\mathbf{q }} e^{i \mathbf{q } \cdot \mathbf{b }} \left\langle \rho ^{b_1}\left( \mathbf{p } +\frac{\mathbf{q }}{2}\right) \rho ^{b_2}\left( \mathbf{p } -\frac{\mathbf{q }}{2}\right) ^* \right\rangle _p 4 L^i\left( \mathbf{k },\mathbf{k }-\mathbf{p }-\frac{\mathbf{q }}{2} \right) L^i\left( \mathbf{k },\mathbf{k }-\mathbf{p }+\frac{\mathbf{q }}{2}\right) . \end{aligned}$$

(C10)

Using the models employed through this work, Eqs. (63) and (60), we can write the single particle Wigner function as

$$\begin{aligned} W^{b_1 b_2}(\mathbf{b },\mathbf{p })&=\frac{\delta ^{b_1 b_2}}{N_c^2-1} \frac{1}{2 (2\pi )^3} \frac{(4 \pi )^2}{\xi ^2} e^{-\mathbf{p }^2/\xi ^2} \times \int _{\mathbf{q }} e^{i \mathbf{q } \cdot \mathbf{b }} e^{-\mathbf{q }^2/(4 B_p^{-1})} =\frac{\delta ^{b_1 b_2}}{N_c^2-1} \frac{1}{\pi ^2} \frac{1}{\xi ^2 B_p} e^{-\mathbf{b }^2/B_p-\mathbf{p }^2/\xi ^2}, \end{aligned}$$

(C11)

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 $\mathbf{b }$ and $\mathbf{p }$,

$$\begin{aligned} \int _\mathbf{b }\int _\mathbf{p } W^{a a}(\mathbf{b },\mathbf{p }) = 1. \end{aligned}$$

(C12)

Doing an analogous discussion we can write the 2-particle Wigner function as

$$\begin{aligned}&W^{b_1 b_2 b_3 b_4}(\mathbf{b }_1,\mathbf{p }_1,\mathbf{b }_2,\mathbf{p }_2)\nonumber \\&\quad = \frac{4}{(2\pi )^6} \int _{\mathbf{q }_1 \mathbf{q }_2} e^{i \mathbf{q }_1 \cdot \mathbf{b }_1+i \mathbf{q }_2 \cdot \mathbf{b }_2} \Bigg \langle \rho ^{b_1}\left( \mathbf{p }_1+\frac{\mathbf{q }_1}{2}\right) \rho ^{b_2}\left( \mathbf{p }_1-\frac{\mathbf{q }_1}{2}\right) ^* \rho ^{b_3}\left( \mathbf{p }_2+\frac{\mathbf{q }_2}{2}\right) \rho ^{b_4}\left( \mathbf{p }_2-\frac{\mathbf{q }_2}{2}\right) ^* \Bigg \rangle _p L^i\left( \mathbf{k }_1,\mathbf{k }_1-\mathbf{p }_1-\frac{\mathbf{q }_1}{2}\right) \nonumber \\&\qquad \times L^i\left( \mathbf{k }_1,\mathbf{k }_1-\mathbf{p }_1+\frac{\mathbf{q }_1}{2}\right) L^i\left( \mathbf{k }_2,\mathbf{k }_2-\mathbf{p }_2-\frac{\mathbf{q }_2}{2}\right) L^i\left( \mathbf{k }_2,\mathbf{k }_2-\mathbf{p }_2+\frac{\mathbf{q }_2}{2}\right) . \end{aligned}$$

(C13)

Performing the Wick expansion of the projectile correlator and using again Eqs. (63) and (60) we obtain Eq. (64). We can check that the quantity defined in that equation is not well normalised:

$$\begin{aligned}&\int _{\mathbf{b }_1\mathbf{b }_2} \int _{\mathbf{p }_1\mathbf{p }_2} W^{a a b b}(\mathbf{b }_1,\mathbf{p }_1,\mathbf{b }_2,\mathbf{p }_2) = 1+2\frac{1}{(N_c^2-1)}\frac{1}{1+B_p \xi ^2}. \end{aligned}$$

(C14)

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:

$$\begin{aligned}&W^{b_1 b_2 b_3 b_4}(\mathbf{b }_1,\mathbf{p }_1,\mathbf{b }_2,\mathbf{p }_2) =W^{b_1 b_2}(\mathbf{b }_1,\mathbf{p }_1) W^{b_3 b_4}(\mathbf{b }_2,\mathbf{p }_2). \end{aligned}$$

(C15)

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.

Appendix D: Calculation of four gluon inclusive production

In this section we analyse the four gluon inclusive spectrum by taking into account all the terms in Eq. (10). In order to do so we will follow the same arguments that we have used for writing down the triple gluon spectrum in Sect. 4.2. First we note that after performing the Wick expansion of either the target or projectile correlators we have 105 contributions on each side that can be written schematically as

(D1)

(D2)

(D3)

(D4)

(D5)

where the permutations $\mathbf{k }_i \rightarrow -\mathbf{k }_i$ and $\mathbf{k }_i \leftrightarrow \mathbf{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

(D6)

(D7)

(D8)

(D9)

(D10)

Now we generate the Wick diagrams in such a way that they are grouped by their powers of $(N_c^2-1)^{-1}$. In order to do so we exploit property iii) of Sect. 3.2. In this case the suppression of a given diagram is given by $(N_c^2-1)^{n_p+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_c^2-1)^{-1}$. As an example, let us see which are the diagrams with power suppression $(N_c^2-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. Thus we will have $1\times 48$ possibilities that are

  

  (D11)

  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 $\hbox {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\times 32$ possibilities

  

  (D12)

  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 \times 12$ possibilities:

  

  (D13)

• (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 \times 12$ possibilities for the first case,

  

  (D14)

  and $12 \times 12$ possibilities for the second case,

  

  (D15)

• (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 \times 1$ possibilities:

  

  (D16)

With (D11) to (D16) we have found all the 1152 diagrams with a suppression of $(N_c^2-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 $\mathbf{k }_i \rightarrow -\mathbf{k }_i$ or $\mathbf{k }_i \leftrightarrow \mathbf{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_c^2-1)^{0}$, 24 diagrams with a suppression of $(N_c^2-1)^{-1}$, 232 diagrams with a suppression of $(N_c^2-1)^{-2}$, 3088 diagrams with a suppression of $(N_c^2-1)^{-4}$, 4224 diagrams with a suppression of $(N_c^2-1)^{-5}$ and 2304 diagrams with a suppression of $(N_c^2-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_c^2-1)^{-3}$. Let us introduce the shorthand notation $\hat{{\mathcal {D}}}_{n_p}$ 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_c^2-1)^{-3}$ to the 4-gluon spectrum as

(D17)

with $\hat{{\mathcal {D}}}_2^{(1)}$ and $\hat{{\mathcal {D}}}_2^{(2)}$ referring to the first and second contributions to Item ii) discussed above, respectively.

In order to evaluate $\kappa _0\{4\}$ we can use the fact that all the permutations, $\hbox {perm}_i$, of (D1) to (D2) will give the same result since we are integrating over all the momentum $\mathbf{k }_i$. Thus we can write

(D18)

When we evaluate $\kappa _n\{4\}$ with $n \ne 0$ we have to integrate the spectrum times $e^{in(\phi _1+\phi _3-\phi _5-\phi _7)}$ which will break some of the symmetries encoded in the permutations, $\hbox {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 as

(D19)

By using the properties i) and ii) of Sect. 3.2 we can check that this sum has the following symmetry

$$\begin{aligned} f_2(\mathbf{k }_1,\mathbf{k }_3,\mathbf{k }_5,\mathbf{k }_7)=f_2(\mathbf{k }_1,\mathbf{k }_3,\mathbf{k }_7,\mathbf{k }_5). \end{aligned}$$

(D20)

Thus, the contribution of the $n_T=3$ diagrams to the $\kappa$-function can be written as

(D21)

where the factor 2 comes from exploiting the symmetries $\mathbf{k }_i \rightarrow -\mathbf{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

(D22)

For the $n_T=2$ diagrams that are defined by Eq. (D3) we can follow the same arguments by defining

(D23)

We can check that this function has the following symmetries

$$\begin{aligned} f_3(\mathbf{k }_1,\mathbf{k }_3,\mathbf{k }_5,\mathbf{k }_7)=f_3(-\mathbf{k }_1,-\mathbf{k }_7,-\mathbf{k }_5,-\mathbf{k }_3) =f_3(-\mathbf{k }_1,-\mathbf{k }_3,-\mathbf{k }_7,-\mathbf{k }_5)=f_3(-\mathbf{k }_1,-\mathbf{k }_5,-\mathbf{k }_3,-\mathbf{k }_7). \end{aligned}$$

(D24)

Thus, exploiting this symmetries we can write the contribution of the diagrams that have the configuration of Eq. (D3) on the right side with generic $n_p$ to the cumulant as

(D25)

For the $n_T=2$ diagrams that are defined by Eq. (D4) we define

(D26)

which has the following symmetries

$$\begin{aligned} f_4(\mathbf{k }_1,\mathbf{k }_3,\mathbf{k }_5,\mathbf{k }_7)=f_4(\mathbf{k }_3,\mathbf{k }_1,\mathbf{k }_5,\mathbf{k }_7)=f_4(\mathbf{k }_1,\mathbf{k }_3,\mathbf{k }_7,\mathbf{k }_5). \end{aligned}$$

(D27)

Thus, the contribution of the diagrams that have the configuration of Eq. (D4) on the right side with generic $n_p$ to the cumulant is

(D28)

For the $n_T=1$ diagrams that are defined by Eq. (D5) we define

(D29)

This function has the following symmetries

$$\begin{aligned} f_5(\mathbf{k }_1,\mathbf{k }_3,\mathbf{k }_5,\mathbf{k }_7)=f_5(-\mathbf{k }_5,-\mathbf{k }_3,-\mathbf{k }_1,-\mathbf{k }_7) =f_5(-\mathbf{k }_1,-\mathbf{k }_7,-\mathbf{k }_5,-\mathbf{k }_3)=f_5(-\mathbf{k }_3,-\mathbf{k }_1,-\mathbf{k }_7,-\mathbf{k }_5). \end{aligned}$$

(D30)

Therefore, the contribution of the diagrams that have the configuration of Eq. (D5) on the right side with generic $n_p$ to the cumulant is

(D31)

All in all, using Eq. (D17) and the simplifications in Eqs. (D22), (D25), (D28) and (D31), we can reduce the number of integral when computing $\kappa _n\{4\}$ at order $(N_c^2-1)^{-3}$ down to

(D32)

We compute the contribution at different orders of $(N_c^2-1)^{-1}$ in the same way. Finally, we just have to solve numerically Eqs. (D32) and (D18) and the equivalent ones at different order. By doing that we are able to obtain Fig. (7).