1 Introduction

While the focus of the physics programme at the large hadron collider (LHC) is the discovery and understanding of the properties of the previously missing piece in the Standard Model – the Higgs boson – and the search for its eventual failure, it has also shown very surprising and unexpected aspects of quantum chromodynamics (QCD), particularly in small collisions systems, pp and pA. One of the most exciting observations made in high multiplicity pp collisions by the CMS collaboration during the first LHC run is the discovery of the correlations between produced particles over large intervals of rapidity, peaking at zero relative azimuthal angle [1]. This phenomenon was dubbed ridge due its shape in the azimuthal angle-rapidity plot, and constitute one of the key findings at the LHC (see Fig. 1).

Later on, this structure was found by other collaborations and for smaller multiplicities [2,3,4,5] and in association with Z boson production [6]. A similar ridge structure was also observed in pPb collisions at the LHC by the four large collaborations [7,8,9,10]. A maximum in the correlations also appears at azimuthal angle \(\pi \), called the away side ridge in contrast to the near side ridge peaked at zero azimuthal angle. They have also been observed in PbPb collisions, see e.g. [11,12,13] for PbPb results and a comparison with those in pPb. Similar correlations were observed in AuAu, dAu and \(^3\)HeAu collisions at the relativistic heavy ion collider (RHIC) [14,15,16,17,18,19]. They have also been observed in photoproduction on Pb in ultraperipheral collisions (UPCs) at the LHC [20]. Their existence in smaller systems like \(e^+e^-\) collisions [21] at the large electron-positron collider and deep inelastic scattering (DIS) events in ep at the Hadron-Elektron-Ringanlage [22] has been scrutinised, but the results are not conclusive. The ridge is the most striking feature in the long list of similarities between small and large collision systems in the observed results for many observables [23,24,25,26,27].

The standard explanation for such azimuthal asymmetries in heavy ion collisions (HICs) is the existence of strong final state interactions that lead to a situation where viscous relativistic hydrodynamics can be applied, see the reviews [28, 29]. The dynamics leading to such a situation, called hydrodynamisation, is unclear [30] and both strong and weak coupling explanations have been proposed, see e.g. [31] and refs. therein. Furthermore, the hydrodynamic description seems to hold for large anisotropies, i.e. rather far from local equilibrium.

Fig. 1
figure 1

Two particle correlations in pp and pPb collisions at the LHC measured by the ATLAS Collaboration [2], for different energies and particle multiplicities in the event. Taken from [2]

The hydrodynamic description of the azimuthal asymmetries in pp and pPb collisions at the LHC is successful [29, 32]. There is a hot ongoing discussion on the explanation of such success and of the fact that hydrodynamics seems to be the effective description for long wavelength modes in any field theory, see e.g. [33] and refs. therein. But it demands a very careful choice of initial conditions, specifically that the proton is modelled as a collection of constituent quarks or hot spots. The description seems to be pushed to the limit of small collision areas and low particle densities where non-hydrodynamic modes play a very important role, as seen in both hydrodynamic studies [29] and in those that consider a weak coupling quasiparticle picture in transport frameworks [34].

Therefore, the hydrodynamic explanation for the azimuthal asymmetries in small systems looks tenuous. Besides, causality arguments show that long range correlations in rapidity must come from the very early stages of the collision [35]. And hydrodynamic calculations demand initial conditions that contains long range rapidity correlations, initial energy and particle density and flow profiles and, unless they are assumed to be completely washed out by final state interactions, correlations. So, beyond addressing the obvious fundamental question: is the strong interaction dynamics capable to lead to collectivity through final state interactions even in small systems, or is the origin of the ridge correlations different in pp and pPb collisions than in HICs?, searching for correlations coming from the initial state is linked to the understanding of the dynamics prior to the use of hydrodynamics and to the provision of well grounded initial conditions for hydrodynamic calculations.

This contribution is devoted to the description and discussion of those frameworks which lead to correlations among partons in the initial stage that, if not washed out by strong final state interactions and hadronisation – that we will assume in the following, may lead to azimuthal asymmetries as observed in data. We start by those based on the weak coupling but non-perturbative description of dense partonic systems offered by the the Color Glass Condensate (CGC) effective theory, see the reviews [36,37,38] and the book [39]. This will be the subject of Sects. 2 and 3. We will then review other explanations inspired in QCD in Sect. 4, to end with a summary and discussions in Sect. 5.

Our focus will be on recent formal developments and we will base the presentation in our own works and formalism, trying to make connection with the other formalisms which differ in notation. We will comment briefly on the status of the comparison to experimental data in the summary.

2 Two particle correlations from the CGC

The observation of the ridge correlations in small size systems has triggered a lot of efforts to understand whether the structure of the initial state itself can lead, in pp and pA collisions, to such correlations without resourcing to final state interactions. Over the last decade, several mechanisms have been suggested to explain the ridge correlations in the CGC framework. The two most successful ones are the “domain structure of the target” developed in [40,41,42] and the “glasma graph approach” introduced in [35, 43, 44].Footnote 1

The underlying mechanism for the domain structure of the target can be summarised as follows: the hadronic target is assumed to contain domains of oriented chromoelectric fields in the transverse plane. When two partons (normally assumed to be gluons when the scattering takes place at high energies and the probed values of momentum fraction of the partons, x, is small) from the projectile are close enough to scatter on the same domain, they get a common final momentum that reflects the correlated structure of the fields in the target. As gluons belong to the adjoint, thus real, representation of the \(SU(N_c)\) group, the correlation holds for both parallel and antiparallel momenta, thus justifying the near and away side structures. The size of the domain in the target is assumed to be of order \(1/Q_s\), with \(Q_s\) being the saturation momentum which is the characteristic transverse momentum for the partons in the saturated target wave function described by the CGC [36,37,38]. Projectile partons lying closer than \(1/Q_s\) contribute mainly to particle production in the region of transverse momentum \(p_T\gtrsim Q_s\). Therefore, this mechanism should mainly be applicable in that transverse momentum region.

Note that this model implies a non-trivial target structure that goes beyond the usual isotropic averages employed in CGC calculations, see below. While still lacking justification from first principles (although indeed CGC numerical calculations indicate that field correlations in the hadron wave functions are characterised by length scales \(\sim 1/Q_s\) [49,50,51]), this explanation is often used for qualitative discussion and understanding of numerical results, and may have further implications on e.g. spin or Transverse Momentum Distributions (TMDs) physics. Numerical studies based on models containing this domain structure have been performed in [52,53,54]. They show correlations that go beyond leading number of colours, see the discussion below, and lead to odd harmonics, see Sect. 3.

On the other hand, the glasma graph approach to two particle correlations is very successful to describe many features of the data as shown in [55,56,57,58,59,60], but the physics behind this approach was not clear. This issue has been studied in [61] and it has been shown that a genuine quantum effect, Bose enhancement of the gluons in the projectile wave function, leads to final state correlations in the glasma graph approach.Footnote 2

The concept of Bose enhancement for a generic quantum system can be understood by considering a state with fixed occupation number, \(\{n_i(p)\}\), of N species of bosons at different momenta which, up to some normalisation factor, can be written as

$$\begin{aligned} \big | \{ n_i(p)\}\big \rangle \propto \prod _{i,p}\left[ a_i^\dagger (p)\right] ^{n_i(p)}|0\rangle , \end{aligned}$$
(1)

with \(a^\dagger _i(p)\) the creation operator of the boson and \(i=1,2,\ldots , N\). The mean particle density \({\tilde{n}}\) is defined as the expectation value of the number operator in this state:

$$\begin{aligned} {\tilde{n}}\equiv \big \langle \{ n_i(p)\} \big |\sum _j a^\dagger _j(x)a_j(x)\big | \{ n_i(p)\}\big \rangle =\sum _{i,p} n_i(p). \end{aligned}$$
(2)

The two particle correlator in momentum space C(pk) is defined in a similar way and can be calculated in a trivial manner:

$$\begin{aligned} C(p,k)= & {} \Big [ \sum _in_i(p)\Big ]\, \Big [ \sum _jn_j(k)\Big ]\nonumber \\&+\delta (p-k)\sum _i\big [n_i(p)\big ]^2. \end{aligned}$$
(3)

The first term on the right hand side of Eq. (3) is the square of the mean particle density and the second term is the Bose enhancement term. It vanishes when the momenta of the two bosons are different and gives an enhancement when the momenta of two bosons coincide which is \({\mathcal O}(1/N)\), due to the fact that it contains a single sum over the species index. The physics behind this is the fact that only bosons of the same species are correlated with each other.

Let us now describe how Bose enhancement arises in the CGC and leads to final state correlations by considering the double inclusive gluon production within the glasma graph approach. In this approach each gluon is assumed to come from a different colour charge density in the projectile wave function that is rapidity invariant.Footnote 3 For our purposes, these colour charge densities can be conveniently represented in terms of gluon creation and annihilation operators in the incoming projectile wave function. After averaging over the target fields the glasma graphs can be written as sum of three types of diagrams (see Fig. 2).

Fig. 2
figure 2

Glasma graphs for two gluon inclusive production before averaging over the projectile colour charge density \(\rho \). Black blobs denote vertices and dashed lines the cuts separating the amplitudes (to the left of the cut) from the complex conjugate amplitudes (to the right). Taken from [61]

Type A diagrams describe the case when two gluons with transverse momenta \(\mathbf{k}_1\) and \(\mathbf{k}_2\) scatter independently on the target, acquiring transfer of transverse momentum \(\mathbf{p-k}_1\) and \(\mathbf{q-k}_2\) so that the outgoing gluons have transverse momenta \(\mathbf{p}\) and \(\mathbf{q}\). Type B and Type C diagrams include interference contributions which are also interesting to study but, for now, let us focus on the Type A contribution and show how the Bose enhancement effect can be observed by studying these diagrams alone. The Type A contribution to the double inclusive gluon production can be written as

$$\begin{aligned} \mathrm{Type\, A}= & {} \int \frac{d^2\mathbf{k}_1}{(2\pi )^2}\frac{d^2\mathbf{k}_2}{(2\pi )^2} \, \langle \mathrm{in}| a^{\dagger i}_a(\mathbf{k}_1) a^{\dagger j}_b(\mathbf{k}_2) a^{k}_a(\mathbf{k}_1)a^{l}_b(\mathbf{k}_2)|\mathrm{in}\rangle \nonumber \\&\times \bigg [ \delta ^{ik}-\frac{\mathbf{k}_1^i\mathbf{k}_1^k}{\mathbf{p}^2}\bigg ]\, \bigg [ \delta ^{jl}-\frac{\mathbf{k}_2^j\mathbf{k}_2^l}{\mathbf{q}^2}\bigg ] N(\mathbf{p}-\mathbf{k}_1)\, N(\mathbf{q}-\mathbf{k}_2),\nonumber \\ \end{aligned}$$
(4)

where \(|\mathrm{in}\rangle \) is the wave function of the incoming projectile and \(N(\mathbf{p-k})\) is the dipole scattering amplitude – the scattering amplitude for a two gluon system to scatter on the target. Moreover, the rapidity dependence of the gluon creation and annihilation operators is integrated over. The explicit dependence on rapidity becomes important only when the rapidity difference between the observed particles is parametrically large, \(\varDelta \eta \gtrsim 1/\alpha _s\).

The evaluation of the expectation value of any operator in the incoming projectile state requires a two averaging procedure in the CGC. In [61], averaging over the valence colour charge density is performed first which leads to the density matrix operator \(\hat{\rho }\) on the soft gluon Hilbert space. Then, the second averaging over the soft gluons is performed using this density matrix operator. The two particle correlator that appears in the Type A contribution calculated with this procedure leads to the following result:

$$\begin{aligned} C(\mathbf{k}_1,\mathbf{k}_2)= & {} S_{\perp }^2(N_c^2-1)^2 \frac{\mathbf{k}_1^i\mathbf{k}_1^k}{\mathbf{k}_1^2} \frac{\mathbf{k}_2^j\mathbf{k}_2^l}{\mathbf{k}_2^2} \frac{g^4\mu ^2(\mathbf{k}_1)\mu ^2(\mathbf{k}_2)}{\mathbf{k}_1^2\mathbf{k}_2^2} \nonumber \\&\quad \times \bigg \{1+\frac{1}{S_\perp (N_c^2-1)}\Big [ \delta ^{(2)}(\mathbf{k_1-k_2})\nonumber \\&\quad +\delta ^{(2)}(\mathbf{k_1+k_2})\Big ]\bigg \},\nonumber \\ \end{aligned}$$
(5)

where \(S_\perp \) is the transverse area of the projectile. The first term on the right hand side of Eq. (5) is the classical term which corresponds to the square of the number of gluons, while the second term is the typical Bose enhancement term, relatively suppressed by the number of states in the adjoint colour representation.

If we consider a situation where the incoming projectile has intrinsic saturation momentum \(Q_s\) and the momenta of the produced gluons are also \(\sim Q_s\), i.e. \(|\mathbf{p}|\sim |\mathbf{q}|\sim Q_s\), then the production amplitude is dominated by the contributions \(|\mathbf{k}_1|\sim |\mathbf{k}_2|\sim Q_s\). The initial state correlations are encoded in the Bose enhancement terms in Eq. (5), which are delta functions. The interaction with the target is obtained by convoluting the two particle correlator with the dipole amplitudes \(N(\mathbf{p-k_1})N(\mathbf{q-k_2})\). Since, in this kinematics, the momentum transfers from the target ( \(|\mathbf{p-k_1}|\sim |\mathbf{q-k_2}|\ll Q_s\) ) are small and the Bose enhancement terms involve delta functions, these initial state correlations naturally transform into angular correlations between the produced gluons in the final state. In a more general case, the delta functions, which are an artefact of considering a translationally invariant projectile, are smeared when convoluted with the dipole scattering amplitudes but this should not completely destroy the final state angular correlations.

The immediate question that arises after the study of gluons is whether quarks are subject to correlations in the CGC. This question has been posed in [64] where the correlations between the produced quarks were studied. The results in [61] show that the origin of the correlations between the produced gluons is the Bose enhancement of the projectile gluons. Due to their fermionic nature, one expects quarks to experience Pauli blocking which effectively amounts to a suppression of the probability of finding two quarks with the same quantum numbers in the CGC state. Therefore, one should expect a negative correlation between the final state quarks that originate from the initial state ones. On the other hand, the correlation between the gluons is found to be long range in rapidity since the CGC wave function is dominated by the rapidity integrated soft gluon field. Thus, another important question to answer is: are the (anti)correlations between the final state quarks long or short range in rapidity? The answer to this question is not obvious a priori. In the projectile wave function, quarks are produced via splitting of the rapidity invariant gluons into quark-antiquark pairs. However, the splitting amplitude itself depends on the rapidity of the quark and antiquark. Moreover, due to this splitting in the projectile wave function the expression for the production cross section of quarks is much more complicated compared to the one for gluons. These questions are answered in [64] where it was shown that the initial state correlations between the quarks in the projectile wave function are not distorted by the small momentum transfer from the target in specific kinematics. In these kinematics, the rapidity difference between the produced quarks is relatively large, i.e. \(\eta _1-\eta _2\gg 1\). Moreover, a large contribution comes from the situation where the transverse momenta of the produced quarks \(\mathbf{p}\) and \(\mathbf{q}\) are of the same order and much larger than the saturation scale of the projectile \(Q_s\), and the saturation scale of the projectile is much larger than saturation scale of the target \(Q_T\), i.e. \(|\mathbf{p}|\sim |\mathbf{q}|\gg Q_s\gg Q_T\). Then the contribution to the production cross section that is sensitive to correlations has the following behaviour (see Eqs. (3.19) and (3.20) in [64] for the full expressions):

$$\begin{aligned} \frac{d\sigma }{d^2\mathbf{p}d\eta _1\, d^2\mathbf{q}d\eta _2}\bigg |_\mathrm{corr.}\propto -\Big [e^{-(\eta _1-\eta _2)}(\eta _1-\eta _2)^2\Big ]. \end{aligned}$$
(6)

The negative sign of this contribution shows that it suppresses the classical term as opposed to the gluon case. This is the result of the Pauli blocking effect in quark-quark production. Moreover, this effect decays exponentially with the rapidity difference between the two produced quarks, which shows that it is short range in rapidity. However, this exponential decrease is tempered by two powers of the rapidity difference.

Besides, there is another physical effect present in the glasma graph approach which is referred to as the Hanbury–Brown–Twiss (HBT) correlations between the produced gluonsFootnote 4 [67]. The diagrams in the glasma graph approach that lead to HBT correlations are those in Fig. 3 after performing pair wise contraction of the colour charges in the projectile wave function. Assuming a translationally invariant projectile wave function, the contribution from Type B and Type C diagrams to the production cross section is

$$\begin{aligned} \mathrm{Type\, B}\propto \delta ^{(2)}(\mathbf{p-q}),\quad \mathrm{Type\, C}\propto \delta ^{(2)}(\mathbf{p+q}). \end{aligned}$$
(7)
Fig. 3
figure 3

Glasma graph diagrams (after averaging over the projectile colour charges) that lead to HBT correlations. Taken from [67]

If the translational invariance condition is relaxed, then the delta functions are smeared over a scale of the inverse size \(R^{-1}\) of the projectile: \(|\mathbf{p\pm q}|\sim R^{-1}\). This size R represents the radius of the gluon cloud inside the proton and its inverse is smaller than the saturation scale, \(R^{-1}< Q_s\). Moreover, it is also shown in [67] that the HBT correlations are long range in rapidity just as the Bose enhancement effect. Thus, the strength of the HBT correlations is equal when the rapidities of the two produced gluons are similar (\(\eta _1\simeq \eta _2\)) or when the difference between them is large (\(|\eta _1-\eta _2|\gg 1\)).

To sum up, the correlation function \(C(\mathbf{p,q})\), formally defined as the ratio of double inclusive gluon production cross section to the square of the single one, in the glasma graph approach contains two physical effects and can be written as follows:

$$\begin{aligned} C(\mathbf{p,q})=1+C(\mathbf{p,q})\Big |_{\mathrm{BE}}+C(\mathbf{p,q})\Big |_\mathrm{HBT}. \end{aligned}$$
(8)

The first term on the right hand side of Eq. (8) is the classical contribution which originates from the square of single inclusive production. \(C(\mathbf{p,q})\big |_{\mathrm{BE}}\) represents the effect of Bose enhancement of the gluons in the projectile wave function. As described above, this effect leads to a correlation of the final state gluons. On the other hand, \(C(\mathbf{p,q})\big |_\mathrm{HBT}\) represents the HBT correlations in the glasma graph approach which directly introduces correlations between the final state gluons. Both \(C(\mathbf{p,q})\big |_{\mathrm{BE}}\) and \(C(\mathbf{p,q})\big |_{\mathrm{HBT}}\) are rapidity independent, therefore long range in rapidity. The Bose enhancement contribution is suppressed by the transverse area of the projectile with respect to the HBT contribution (actually by the number of “sources” \(Q_s^2 S_\perp \)). However, it leads to correlations whose width in momentum space is determined by the saturation momentum \(Q_s\). On the other hand, the HBT contribution is not suppressed but it gives a narrow peak in momentum space with width \(R^{-1}\). This comparison is shown in Fig. 4.

Fig. 4
figure 4

Schematic separation in q (here the modulus of the difference in transverse momentum between the produced gluons) between the contributions to the HBT effect (solid line) and to the Bose enhancement effect (dashed line) in the two particle correlation function. Taken from [67]

In the explicit calculations in the glasma graph approach to double inclusive particle production,Footnote 5 the averaging over the target configurations that leads to the dipole scattering amplitude is performed expanding this amplitude to the lowest non trivial order in the target field strength, corresponding to two gluon exchange between the gluons in the projectile wave function and the target. The dipole scattering amplitudes \(N(\mathbf{p-k_1})\) and \(N(\mathbf{q-k_2})\) introduced in Eq. (4) are assumed to originate from single pairs of target fields and therefore this approach does not take into account the effects of multiple scatterings in a dense target. Therefore, this approach is only valid for pp collisionsFootnote 6. In [70], the inclusive production of two and three gluons is computed beyond the glasma graph approach by including the multiple scattering effects, which extends the validity of the glasma graph approach from pp to pA collisions.Footnote 7

Apart from taking into account the multiple scattering effects in [70], a systematic way to identify each term in the double inclusive gluon production cross section and to characterise whether it is a Bose enhancement or HBT contribution is introduced. This identification is performed by adopting the following strategy. When calculating the double inclusive gluon production, one has to average over four colour charges (two in the amplitude and two in the complex conjugate amplitude) in the projectile wave function: \(\langle \rho ^{a_1}(\mathbf{x_1})\rho ^{a_2}(\mathbf{x_2})\rho ^{b_1}(\mathbf{y_1})\rho ^{b_2}(\mathbf{y_2})\rangle _P\). Here, \((\mathbf{x_i},a_i)\) and \((\mathbf{y_i},b_i)\) stand for the transverse position and colour indices of the colour charge densities in the amplitude and in the complex conjugate amplitude, respectively. The averaging over the colour charge distributions in the projectile is commonly performed by using a generalised McLerran–Venugopalan (MV) model [73, 74] where the weight functional is Gaussian. Then, the average of any product of colour charge densities factorises into a product of all possible pair, Wick-like contractions. The correlator of two colour charge densities in momentum space can be defined as

$$\begin{aligned} \langle \rho ^a(\mathbf{k})\rho ^b(\mathbf{p})\rangle _P=\delta ^{ab}\mu ^2(\mathbf{k, p}). \end{aligned}$$
(9)

The function \(\mu ^2(\mathbf{k, p})\) characterises the structure of the projectile. It can be written as

$$\begin{aligned} \mu ^2(\mathbf{k, p})=T\bigg ( \frac{\mathbf{k-p}}{2}\bigg ) F\big [ (\mathbf{k+p})R\big ], \end{aligned}$$
(10)

where \(F\big [ (\mathbf{k+p})R\big ]\) is a soft form factor with maximal value F(0), and R is the radius of the projectile. Function T defines the transverse momentum dependent distribution of the valence charges. The soft form factor identifies whether a term is a contribution to the Bose enhancement of the projectile gluons or a contribution to the HBT correlations of the produced gluons. For example, in our set up the produced gluons have momenta \(\mathbf{p}\) and \(\mathbf{q}\), while the projectile gluons carry transverse momenta \(\mathbf{k_1}\) and \(\mathbf{k_2}\). In this case, \(\mu ^2(\mathbf{p,q})\) gives a maximal contribution when \(\mathbf{p+q}=0\) which clearly can be identified as the HBT correlations of the produced gluons. \(\mu ^2(\mathbf{k_1,k_2})\) is peaked when \(\mathbf{k_1+k_2}=0\) which is a contribution to the Bose enhancement of the gluons in the projectile wave function.

On the other hand, multiple scattering effects on the dense target are taken into account by introducing the standard Wilson lines in the CGC framework. In this framework, the interaction between the projectile and the target is assumed to be eikonal which amounts to the situation where each parton produced by the projectile colour charge scatters on the target by picking up a colour rotation described by a Wilson line which is defined as an exponential of the target field ordered in the \(x^+\) coordinate:

$$\begin{aligned} U_{{\mathcal {R}}}(\mathbf{x})={{\mathcal {P}}}_+\, e^{ig\int dx^+\, T^a_{{\mathcal {R}}}\, A^-_a(x^+,\mathbf{x})}, \end{aligned}$$
(11)

at the amplitude level. Here, \(T^a_{{\mathcal {R}}}\) is the \(SU(N_c)\) generator in the representation \({{\mathcal {R}}}\) which can be the fundamental one for a quark and the adjoint one for a gluon. This leads to the appearance in the cross section of double dipole and quadrupole amplitudes (in the adjoint representation) of the type

$$\begin{aligned} \langle s(\mathbf{x,y})s(\mathbf{z,w})\rangle _T,\quad \langle Q(\mathbf{x,y,z,w})\rangle _T, \end{aligned}$$
(12)

which have to be averaged over the target field distributions. The dipole and the quadrupole operators are defined as

$$\begin{aligned}&s(\mathbf{x},\mathbf{y})=\frac{1}{N_c^2-1} \mathrm{tr}\big [U(\mathbf{x})U^\dagger (\mathbf{y})\big ], \end{aligned}$$
(13)
$$\begin{aligned}&Q(\mathbf{x,y,z,w})=\frac{1}{N_c^2-1}\mathrm{tr}\big [ U(\mathbf{x})U^\dagger (\mathbf{y})U(\mathbf{z})U^\dagger (\mathbf{w})\big ]. \end{aligned}$$
(14)

The cross section has to be integrated over four transverse coordinates. In principle, the maximal contribution should come from the area in coordinate space, i.e. when all the four coordinates are far away from each other. However, all four points cannot be far away from each other since the target field ensemble has to be colour neutral, and colour neutralisation in the CGC happens on scales of order \(1/Q_s\). Therefore, the maximal contribution to the integral must come from the configurations where the four points are combined into pairs, such that each pair is a singlet and the distance between the pairs is large. This is the leading contribution in \(1/(Q_s^2 R^2)\) to the integral on transverse coordinatesFootnote 8 – not, by any means, a good representation of the target averages of ensembles of Wilson lines by themselves. Taking into account only such configurations is equivalent to calculating the target averages of products of any number of Wilson lines by factorising them into averages of pairs with basic Wick contraction. In this case, target averaging of the double dipole and the quadrupole amplitudes can be written as

$$\begin{aligned} \langle Q(\mathbf{x,y,z,w})\rangle _T\approx & {} d(\mathbf{x,y})d(\mathbf{z,w}) + d(\mathbf{x,w})d(\mathbf{z,y})\nonumber \\&+\frac{1}{N_c^2-1}d(\mathbf{x,z})d(\mathbf{y,w}), \end{aligned}$$
(15)
$$\begin{aligned} \langle s(\mathbf{x,y})s(\mathbf{z,w})\rangle _T\approx & {} d(\mathbf{x,y})d(\mathbf{z,w}) \nonumber \\&+\, \frac{1}{(N_c^2-1)^2}\Big [ d(\mathbf{x,w})d(\mathbf{z,y})\nonumber \\&+d(\mathbf{x,z})d(\mathbf{y,w})\Big ], \end{aligned}$$
(16)

where we have defined \(d(\mathbf{x,y})\equiv \langle s(\mathbf{x,y})\rangle _T\). Then, by using the function \(\mu ^2(\mathbf{k,p})\) given in Eq. (10) for the projectile colour charge density correlators and using the factorisation ansatz described above for the double dipole and quadrupole amplitudes, the double inclusive gluon production cross section is computed and the nature of all the terms is identified. Moreover, it is also shown that the contributions to final state correlations comes from the quadrupole terms in two gluon production.Footnote 9

Fig. 5
figure 5

Momentum assignment for the double inclusive gluon production. The grey blobs represent the colour charge densities in the projectile wave function. \(\mathbf{q}_1\) and \(\mathbf{q}_2\) are the transverse momenta transferred from the target during the interaction

3 Odd azimuthal harmonics from the CGC

To describe one key existing problem in usual CGC calculations, namely the absence of odd harmonics, let us briefly discuss double inclusive particle production in more depth. Within the approximations described in Sect. 2, double inclusive gluon production is computed in pA collisions in [70]. The production cross section of two gluons with rapidities \(\eta _1\) and \(\eta _2\), and transverse momenta \(\mathbf{k}_1\) and \(\mathbf{k}_2\) (see Fig. 5) reads

$$\begin{aligned}&\frac{d\sigma }{d^2\mathbf{k}_1d\eta _1 d^2\mathbf{k}_2d\eta _2}=\alpha _s^2(4\pi )^2(N_c^2-1)^2 \!\!\int \frac{d^2\mathbf{q}_1}{(2\pi )^2} \frac{d^2\mathbf{q}_2}{(2\pi )^2} d(\mathbf{q}_1)\nonumber \\&\quad \times \, d(\mathbf{q}_2) \bigg \{ I_0+\frac{1}{N_c^2-1}I_1+\frac{1}{(N_c^2-1)^2}I_2\bigg \} + (\mathbf{k}_2\rightarrow -\mathbf{k}_2),\nonumber \\ \end{aligned}$$
(17)

where

$$\begin{aligned} I_0= & {} \frac{1}{2}\mu ^2(\mathbf{k}_1-\mathbf{q}_1, \mathbf{q}_1-\mathbf{k}_1) \, \mu ^2(\mathbf{k}_2-\mathbf{q}_2, \mathbf{q}_2-\mathbf{k}_2) \nonumber \\&\times \, L^i(\mathbf{k}_1, \mathbf{q}_1)L^i(\mathbf{k}_1, \mathbf{q}_1) \, L^j(\mathbf{k}_2, \mathbf{q}_2)L^j(\mathbf{k}_2, \mathbf{q}_2), \end{aligned}$$
(18)
$$\begin{aligned} I_1= & {} \mu ^2\ (\mathbf{k}_1-\mathbf{q}_1, \mathbf{q}_2-\mathbf{k}_2) \, \mu ^2(\mathbf{k}_2-\mathbf{q}_2, \mathbf{q}_1-\mathbf{k}_1)\nonumber \\&\times \, L^i(\mathbf{k}_1, \mathbf{q}_1)L^i(\mathbf{k}_1, \mathbf{q}_1) \, L^j(\mathbf{k}_2, \mathbf{q}_2)L^j(\mathbf{k}_2, \mathbf{q}_2) \nonumber \\&+\, \mu ^2\ (\mathbf{k}_1-\mathbf{q}_1, \mathbf{q}_1-\mathbf{k}_2) \, \mu ^2(\mathbf{k}_2-\mathbf{q}_2, \mathbf{q}_2-\mathbf{k}_1)\nonumber \\&\times L^i(\mathbf{k}_1, \mathbf{q}_1)L^i(\mathbf{k}_1, \mathbf{q}_2) \, L^j(\mathbf{k}_2, \mathbf{q}_1)L^j(\mathbf{k}_2, \mathbf{q}_2), \end{aligned}$$
(19)
$$\begin{aligned} I_2= & {} \mu ^2(\mathbf{k}_1-\mathbf{q}_1, \mathbf{q}_2-\mathbf{k}_1) \, \mu ^2(\mathbf{k}_2-\mathbf{q}_2, \mathbf{q}_1-\mathbf{k}_2) \nonumber \\&\times L^i(\mathbf{k}_1, \mathbf{q}_1)L^i(\mathbf{k}_1, \mathbf{q}_2) \, L^j(\mathbf{k}_2, -\mathbf{q}_1)L^j(\mathbf{k}_2, -\mathbf{q}_2) \nonumber \\&+ \mu ^2(\mathbf{k}_1-\mathbf{q}_1, \mathbf{q}_2-\mathbf{k}_2) \, \mu ^2(-\mathbf{k}_1-\mathbf{q}_2, \mathbf{q}_1+\mathbf{k}_2) \nonumber \\&\times L^i(\mathbf{k}_1, \mathbf{q}_1)L^i(\mathbf{k}_1, -\mathbf{q}_2) \, L^j(\mathbf{k}_2, -\mathbf{q}_1)L^j(\mathbf{k}_2, \mathbf{q}_2).\nonumber \\ \end{aligned}$$
(20)

Here, function \(\mu ^2\) defines the structure of the projectile, Eq. (10). On the other hand, function \(L^i(\mathbf{k},\mathbf{q})\) is the eikonal Lipatov vertex which is defined as

$$\begin{aligned} L^i(\mathbf{k}, \mathbf{q})\equiv \frac{(\mathbf{k}-\mathbf {q})^i}{(\mathbf{k}-\mathbf{q})^2}-\frac{\mathbf{k}^i}{\mathbf{k}^2}. \end{aligned}$$
(21)

Note that, as discussed in Sect. 2, all correlations are subleading in \(1/N_c\) which is due to the use of Gaussian averages – the MV model. Moreover, the double gluon inclusive production cross section is written in terms of the dipole averages assuming translational invariance of the target (a standard approximation which is reasonable in pA for large nuclei):

$$\begin{aligned} d(\mathbf{x}_1,\mathbf{x}_2)= & {} \int \frac{d^2\mathbf{q}_1}{(2\pi )^2}\frac{d^2\mathbf{q}_2}{(2\pi )^2} e^{-i\mathbf{q}_1\cdot \mathbf{x}_1+i\mathbf{q}_2\cdot \mathbf{x}_2}\nonumber \\&\times d\bigg (\frac{\mathbf{q}_1+\mathbf{q}_2}{2}\bigg )\delta ^{(2)}(\mathbf{q}_1- \mathbf{q}_2). \end{aligned}$$
(22)

A convenient way to study two particle correlations is through a Fourier decomposition into harmonics defined in for the azimuthal angle \(\varDelta \phi \) between the produced gluons with transverse momenta \(\mathbf{k}_1\) and \(\mathbf{k}_2\). When Fourier expanded, the double inclusive gluon spectrum Eq. (17) can be written as

$$\begin{aligned}&N(k_1,k_2, \varDelta \phi )\nonumber \\&\quad =a_0(k_1,k_2)\left[ 1+\sum _{n=0}^{\infty }2V_{n\varDelta }(k_1,k_2)\cos (n\varDelta \phi )\right] , \end{aligned}$$
(23)

where

$$\begin{aligned} V_{n\varDelta }(k_1,k_2)=\frac{\int _0^{\pi } N(k_1,k_2,\varDelta \phi )\cos (n\varDelta \phi )d\varDelta \phi }{\int _0^{\pi }N(k_1,k_2,\varDelta \phi )d\varDelta \phi }. \end{aligned}$$
(24)

One way of defining the \(p_T\) dependence of the Fourier coefficients is by fixing one of the momenta (\(k_1=p^{ref}_T\)), and treating the other one as a free variable (\(k_2=p_T\)). With this choice, the azimuthal harmonics are defined as

$$\begin{aligned} v_n(p_T)=\frac{V_{n\varDelta }(p_t,p^{ref}_T)}{\sqrt{V_{n\varDelta }(p^{ref}_t,p^{ref}_T)}}. \end{aligned}$$
(25)

The key theoretical problem for the description of the two particle correlations within the CGC is the absence of the odd harmonics. This problem is analyzed in [41], and it is shown to be due to the symmetry \((\mathbf{k}_2\rightarrow -\mathbf{k}_2)\) of the double inclusive gluon production cross section Eq. (17) which is also referred to as the “accidental symmetry of the CGC”.

Three waysFootnote 10 have been proposed to solve this problem.Footnote 11 On the one hand, the projectile and target can be characterised by a more involved structure than that considered in the usual MV averages [40,41,42, 52,53,54].

On the other hand and as discussed previously, within the CGC framework each produced gluon originates from a separate colour charge density in the projectile wave function as shown in Fig. 5. The contributions to the projectile wave function that emerge from merging of the gluons before the interaction with the target or splitting of a gluon into two gluons emitted from the same colour charge density, are not taken into account in the standard CGC calculations. Recently, in [80, 81] it is shown that the accidental symmetry of the CGC can be broken by including such corrections to the projectile wave function. The even and odd parts of the double inclusive gluon production cross section under the accidental symmetry are computed separately, and finally, the azimuthal harmonics are calculated. The corresponding numerical studies and a comparison with data are performed in [82, 83].

Finally, inclusive gluon production is usually studied within the eikonal approximation in the CGC framework. In recent studies [84, 85], it is shown that the accidental symmetry can be broken by going beyond this eikonal approximation. In the next subsection, we introduce a systematic way to include subeikonal corrections in CGC calculations and show how these corrections give rise to non-vanishing odd harmonics.

3.1 Subeikonal corrections in the CGC

In inclusive gluon production at central rapidity in pA collisions, both the projectile and the target are highly energetic since they are boosted from their initial rapidity to the central rapidity where the collision occurs. Therefore, in this case both colliding objects can be treated in the CGC framework. This corresponds to defining the projectile by the colour charge \(J^\mu _a(x)\),

$$\begin{aligned} J^\mu _a(x)=\delta ^{\mu +}\delta (x^-)\rho _a(\mathbf{x}), \end{aligned}$$
(26)

and the target by the colour field \(A^\mu _a(x)\) that is given as

$$\begin{aligned} A^{\mu }_a(x)=\delta ^{\mu -}\delta (x^+)A^-_a(\mathbf{x}). \end{aligned}$$
(27)

Let us recall that these expressions of the colour charge of the projectile and the colour field of the target are defined within the eikonal approximation, which is justified by the large energy of both colliding objects. If for the dilute projectile the eikonal approximation can be trusted at a given energy, the same approximation for a large target can be true only for larger energies. The eikonal approximation for the target amounts to the following three conditions:

  1. 1.

    \(A^{\mu }_a(x)\simeq \delta ^{\mu -}A^-_a(x)\): Neglecting the (+) and transverse components of the colour field of the target.

  2. 2.

    \(A^{\mu }_a(x)\simeq A^{\mu }_a(x^+,\mathbf{x})\): Neglecting the \(x^-\) dependence in the colour field of the target.

  3. 3.

    \(A^\mu (x)\propto \delta (x^+)\): Assuming that the target field is peaked around \(x^+=0\) due to Lorentz contraction, which is also known as the shockwave approximation.

In realistic kinematical conditions under which the experiments are performed, the energies are not asymptotic and the eikonal approximation is not always justified. While for a dilute projectile it is usually valid even for high energy collisions, this is not necessarily true for a large nucleus. Relaxing any of the above approximations accounts for corrections to the eikonal limit. In [86, 87], a systematic method to compute the corrections to the eikonal limit by relaxing the third approximation is developed. This corresponds to treating the colour field of the target with a finite longitudinal support \(L^+\) along the \(x^+\) direction, thus replacing Eq. (27) by

$$\begin{aligned} A^{\mu }_a(x)\simeq \delta ^{\mu -}A^-_a(x^+,\mathbf{x}). \end{aligned}$$
(28)

Such subeikonal corrections are thus subleading with respect to the infinite Lorentz contraction of the target.

Before discussing the results, let us give a brief sketch of the method employed to derive the non-eikonal corrections. Let us consider the production of a single gluon with transverse momenta \(\mathbf{k}\) and longitudinal momenta \(k^+\) in pA collisions at central rapidity. The dilute projectile is still treated in the eikonal approximation and defined with the charge density \(J^\mu _a(x)\) given in Eq. (26). On the other hand, the eikonal approximation is relaxed for the dense target that is defined by the colour field \(A^\mu _a(x)\) given in Eq. (28) with a finite support from 0 to \(L^+\) in the longitudinal direction. In this case, the production cross section can be written as the square of the gluon production amplitude averaged over the projectile and target distributions and integrated over impact parameter \(\mathbf{B}\):

$$\begin{aligned} 2k^+\frac{d\sigma }{dk^+d^2\mathbf{k}}= \int d^2\mathbf{B}\sum _{\lambda } \bigg \langle \Big \langle |{{\mathcal {M}}}^a_\lambda (\underline{k}, \mathbf{B})|^2\Big \rangle _P\bigg \rangle _T. \end{aligned}$$
(29)

Here, \(\lambda \), a and \({\underline{k}}=(k^+, \mathbf{k})\) are the polarization, colour and momentum of the produced gluonFootnote 12. For a target with finite longitudinal width, the gluon production amplitude \(\mathcal{M}^a_\lambda (\underline{k}, \mathbf{B})\) is composed of three different contributions: gluon production before, while and after the projectile propagates through the target. At leading order, it is possible to relate the total gluon production amplitude and the background retarded gluon propagator by using the LSZ reduction formula and the perturbative expansion of the colour field of the target [88]. In the light cone gauge \(A^+=0\), the total gluon production amplitude can be written in terms of the \((i-)\) component of the background retarded gluon propagator \(G_R^{\mu \nu }(x,y)\) as

$$\begin{aligned} {{\mathcal {M}}}^a_{\lambda }(\underline{k},\mathbf{B})= & {} \epsilon ^{i*}_\lambda (2k^+)\lim _{x^+\rightarrow 0}\int d^2\mathbf{x}\int dx^- e^{ik\cdot x} \nonumber \\&\times \, \int d^4y\, G^{i-}_R(x,y)_{ab} \, J^+_b(y). \end{aligned}$$
(30)

Since the colour field of the target is independent of \(x^-\), one can introduce the one-dimensional Fourier transform of the background retarded gluon propagator and write it in terms of of the background scalar propagator \({{\mathcal {G}}}_{k^+}^{\mu \nu }({\underline{x}}, {\underline{y}})\). Then, the \((i-)\) component of the retarded background gluon propagator reads

$$\begin{aligned} G_R^{i-}(x,y)_{ab}= & {} \int \frac{dk^+}{2\pi }e^{-ik^+(x^- -y^-)}\nonumber \\&\times \frac{i}{2(k^++i\epsilon )^2}\partial _{\mathbf{y}^i}\mathcal{G}_{k^+}^{ab}(\underline{x}, \underline{y}). \end{aligned}$$
(31)

The background scalar propagator \({{\mathcal {G}}}_{k^+}^{ab}(\underline{x}, \underline{y})\) satisfies the scalar Green’s equation whose solution formally can be written as a path integral

$$\begin{aligned} {{\mathcal {G}}}_{k^+}^{ab}(\underline{x}, \underline{y})= & {} \theta (x^+-y^+)\int _{\mathbf{z}(y^+)=\mathbf{y}}^{\mathbf{z}(x^+)=\mathbf{x}}\big [ {{\mathcal {D}}}{} \mathbf{z}(z^+)\big ] \, \nonumber \\&\times \, e^{\frac{ik^+}{2}\int _{y^+}^{x^+} dz^+ \mathbf{\dot{z}}^2(z^+)} U^{ab}\Big ( x^+,y^+;\big [\mathbf{z}(z^+)\big ]\Big ),\nonumber \\ \end{aligned}$$
(32)

with the Wilson line

$$\begin{aligned}&U^{ab}\Big ( x^+,y^+;\big [\mathbf{z}(z^+)\big ]\Big )\nonumber \\&\quad ={{\mathcal {P}}}_+\, \exp \bigg \{{ig\int _{y^+}^{x^+}\, d\tilde{ z}^+\, T^c\, A^-_c\Big ( \tilde{z}^+,\mathbf{z}(z^+)\Big )}\bigg \}^{ab} \end{aligned}$$
(33)

following the Brownian trajectory \(\mathbf{z}(z^+)\). In the limit of vanishing longitudinal width, \(x^+-y^+\rightarrow 0\), the background scalar propagator \({{\mathcal {G}}}_{k^+}^{ab}(\underline{x}, \underline{y})\) reduces to the standard Wilson line introduced in Eq. (11) and one recovers the eikonal limit. Therefore, it can be safely concluded that all the non-eikonal effects that are due to the finite longitudinal width of the target are encoded in the background scalar propagator. This also means that an expansion of \({{\mathcal {G}}}_{k^+}^{ab}(\underline{x}, \underline{y})\) can be performed in terms of an eikonal parameter, with the first term in this expansion corresponding to the eikonal limit and higher order terms to the corrections to this limit.

Fig. 6
figure 6

Diagrams that contribute to the computation of the Lipatov vertex. The black blob represents the Lipatov vertex which is the sum of all real diagrams for gluon production shown on the right hand side of the equation. Taken from [84]

In order to perform an eikonal expansion of the background scalar propagator \({{\mathcal {G}}}_{k^+}^{ab}(\underline{x}, \underline{y})\), one should first discretise the scalar background propagator. In the eikonal limit, \(k^+/(x^+-y^+)\) is much larger than any transverse scale in the problem. In the large \(k^+\) limit, it is natural to consider a generic path as a perturbation around the classical free path,

$$\begin{aligned} \mathbf{z}_n=\mathbf{z}_n^\mathrm{cl}+\mathbf{u}_n, \end{aligned}$$
(34)

where the transverse positions at step n are on the straight line

$$\begin{aligned} \mathbf{z}_n^\mathrm{cl}=\mathbf{y}+\frac{n}{N}(\mathbf{x}-\mathbf{y}) \end{aligned}$$
(35)

between the initial and final points, and the perturbation \(\mathbf{u}_n\) satisfies the boundary conditions \(\mathbf{u}_0=\mathbf{u}_N=0\) with N being the number of discretised steps. Once the expansion around the free classical path is performed for fixed initial and final positions, one should perform another expansion for small \(\mathbf{x}-\mathbf{y}\), since \(\mathbf{x}-\mathbf{y}\) is parametrically small in the large \(k^+\) limit. After performing these two expansions up to second order in \((x^+-y^+)\) – the finite longitudinal width of the target, the scalar background propagator \(\mathcal{G}_{k^+}^{ab}(\underline{x}, \underline{y})\) reads

$$\begin{aligned}&\int d^2\mathbf{x} \, e^{-i\mathbf{k}\cdot \mathbf{x}} \, \mathcal{G}_{k^+}^{ab}(\underline{x}, \underline{y}) \nonumber \\&\quad =\theta (x^+-y^+)e^{-i\mathbf{k}\cdot \mathbf{y}}e^{-k^-(x^+-y^+)} \Big \{U(x^+,y^+; \mathbf{y}) \nonumber \\&\qquad +\, \frac{(x^+-y^+)}{k^+}\Big [\mathbf{k}^i U^i_{[0,1]}(x^+,y^+; \mathbf{y})\nonumber \\&\qquad + \frac{i}{2}U_{[1,0]}(x^+,y^+; \mathbf{y})\Big ] \nonumber \\&\qquad + \frac{(x^+-y^+)^2}{(k^+)^2} \Big [ \mathbf{k}^i\mathbf{k}^j U^{ij}_{[0,2]}(x^+,y^+; \mathbf{y}) \nonumber \\&\qquad +\frac{i}{2}{} \mathbf{k}^iU^i_{[1,1]}(x^+,y^+; \mathbf{y}) -\frac{1}{4}U_{[2,0]}(x^+,y^+; \mathbf{y})\Big ]\Big \}^{ab}.\nonumber \\ \end{aligned}$$
(36)

The first term on the right hand side of Eq. (36) is the standard Wilson line defined in Eq. (11). \({{\mathcal {O}}}\big [(x^+-y^+)/k^+\big ]\) terms are the first order corrections to the strict eikonal limit which we refer to as next-to-eikonal (NEik) corrections. Similarly, \({{\mathcal {O}}}\big [ (x^+-y^+)^2/(k^+)^2\big ]\) terms are the second order corrections and they are referred to as next-to-next-to-eikonal (NNEik) corrections. The terms that are denoted as \(U_{[\alpha ,\beta ]}(x^+,y^+; \mathbf{y})\) are the decorated Wilson lines which only appear beyond strict eikonal order. The first subscript \(\alpha \) in the decorated Wilson lines stands for the order of expansion around the classical path while the second subscript \(\beta \) stands for the order of the expansion around the initial transverse position \(\mathbf{y}\). The reason why these objects are referred to as decorated Wilson lines is related with their structure. These objects involve a background field insertion into the standard Wilson lines along the \(+\)-direction in a given \(+\)-coordinate. For example, the first decorated Wilson line is defined as

$$\begin{aligned}&\Big [U_{[0,1]}^{i}(x^+,y^+; \mathbf{y})\Big ]^{ab}=\int _{y^+}^{x^+} dz^+ \frac{z^+-y^+}{x^+-y^+} \nonumber \\&\quad \times \, U^{ac}(x^+,z^+; \mathbf{y}) \big [ ig\, T^e_{cd}\, \partial _{\mathbf{y}^i}A^-_e(z^+,\mathbf{y})\big ]\, U^{db}(z^+,y^+; \mathbf{y}).\qquad \end{aligned}$$
(37)

The other decorated Wilson lines have similar structure with one or more background field insertions. We do not present the structure of all the decorated Wilson lines due their complexity and lengthy expressions (see [87]). One can easily get the expression for the gluon production amplitude at NNEik accuracy given in Eq. (30) by using the expression of the retarded background gluon propagator Eq. (31) and the expression derived for the background scalar propagator Eq. (36).

As discussed above, the retarded background gluon propagator \(G_R^{\mu \nu }(x,y)_{ab}\) and, therefore, the scalar background propagator \({{\mathcal {G}}}_{k^+}^{ab}(\underline{x}, \underline{y})\) are the main building blocks for computing cross sections in high energy pA collisions. In [86, 87], these propagators were used to calculate the single inclusive gluon production cross section in pA collisions at NNEik accuracy. The same formalism can be adopted to compute double inclusive gluon production and hence the azimuthal harmonics in pA collisions [84, 85].

In [89], the results of the single inclusive gluon production cross section at NNEik accuracy in pA collisions are used to study the weak field limit (i.e. glasma graph approximation) of this result which corresponds to single inclusive production in pp collisions. In this limit, the decorated Wilson lines are expanded to first order in the background field of the target \(A_a^-(z^+,\mathbf{y})\). For example, the first decorated Wilson line given in Eq. (37) reduces to

$$\begin{aligned} \Big [U_{[0,1]}^{i}(x^+,y^+; \mathbf{y})\Big ]^{ab}\rightarrow & {} \int _{y^+}^{x^+} dz^+\, \frac{z^+-y^+}{x^+-y^+}\nonumber \\&\times \, \left[ ig\, T^c_{ab}\, \partial _{\mathbf{y}^i}A^-_c(z^+,\mathbf{y})\right] . \end{aligned}$$
(38)

This simplification allows us to calculate the Lipatov vertex at NNEik accuracy. After expanding the eikonal and non-eikonal terms to first order in powers of the background field, which corresponds to the glasma graph approach in usual CGC calculations, the Lipatov vertex at NNEik accuracy can be written as

$$\begin{aligned} L^i_{\mathrm{NNEik}}({\underline{k}},\mathbf{q};x^+)= & {} \bigg [\frac{(\mathbf{k}-\mathbf {q})^i}{(\mathbf{k}-\mathbf{q})^2}-\frac{\mathbf{k}^i}{\mathbf{k}^2}\bigg ] \nonumber \\&\times \, \bigg \{1+i\frac{\mathbf{k}^2}{2k^+}x^+-\frac{1}{2}\bigg (\frac{\mathbf{k}^2}{2k^+}x^+\bigg )^2\bigg \}.\nonumber \\ \end{aligned}$$
(39)

The first term on the right hand side of Eq. (39) corresponds to the strict eikonal limit, thus it gives the eikonal Lipatov vertex defined in Eq. (21). The second and the third terms are the NEik and NNEik corrections respectively. The structure of the vertex suggests that the corrections to the amplitude due to finite width of the target may exponentiate.

This observation is further studied recently in [84] where it was shown that indeed the non-eikonal corrections due to finite longitudinal width of the target in the weak field limit exponentiate and can be written as modified Lipatov vertex. By computing the corresponding three diagrams (see Fig. 6) at the amplitude level and keeping the phase \(e^{ik^-x^+}\) which is taken to be one in the eikonal limit, the non-eikonal Lipatov vertex that accounts for all order corrections to the eikonal limit due to finite longitudinal width of the target in the weak field limit reads

$$\begin{aligned} L^i_{\mathrm{NonEik}}({\underline{k}}, \mathbf{q};x^+)=\bigg [\frac{(\mathbf{k}-\mathbf {q})^i}{(\mathbf{k}-\mathbf{q})^2}-\frac{\mathbf{k}^i}{\mathbf{k}^2}\bigg ] e^{ik^-x^+}\ , \end{aligned}$$
(40)

where \(k^-\equiv \mathbf{k}^2/(2k^+)\). This structure was observed in the context of jet quenching in [90,91,92] previously, however the identification of this building block for its use to include non-eikonal corrections in CGC calculations is done in [84] for the first time, further illustrating the close relation between CGC and jet quenching calculations [88].

It is now straightforward to compute the non-eikonal single inclusive gluon production in pp (i.e. dilute-dilute) collisions which formally reads

$$\begin{aligned}&\frac{d\sigma }{d^2\mathbf{k}d\eta }\bigg |_{\mathrm{dilute}}^\mathrm{NonEik}= 4\pi \alpha _s\, C_A \, g^2\int dx_1^+dx_2^+ \int \frac{d^2\mathbf{q}_1}{(2\pi )^2} \frac{d^2\mathbf{q}_2}{(2\pi )^2} \nonumber \\&\quad \times \, \delta ^{c{\bar{c}}}\big \langle A^-_c(x_1^+, \mathbf{q}_1)A^-_{{\bar{c}}}(x_2^+, \mathbf{q}_2) \big \rangle _T \; \mu ^2\big [\mathbf{k}-\mathbf{q}_1, \mathbf{q}_2-\mathbf{k}\big ] \nonumber \\&\quad \times \, L^i_{\mathrm{NonEik}}({\underline{k}}, \mathbf{q}_1;x_1^+)\, L^i_{\mathrm{NonEik}}({\underline{k}}, \mathbf{q}_2;x_2^+). \end{aligned}$$
(41)

An additional modification that is needed to account for the finite longitudinal width of the target is to adopt a modified expression for the correlator of two target fields. Motivated by the non zero longitudinal extent of the target, the fields can be located at different positions that are separated by the colour correlation length in the target \(\lambda ^+\). In this case, the two target field correlator reads

$$\begin{aligned}&\langle A^-_c(x_1^+, \mathbf{q}_1)A^-_{{\bar{c}}}(x_2^+, \mathbf{q}_2) \big \rangle _T=\delta ^{c\bar{c}}\, n(x_1^+)\, \frac{1}{2\lambda ^+}\nonumber \\&\quad \times \, \theta \big (\lambda ^+-|x_1^+-x_2^+|\big ) \, (2\pi )^2 \delta ^{(2)}(\mathbf{q}_1-\mathbf{q}_2)\, |a(\mathbf{q}_1)|^2,\nonumber \\ \end{aligned}$$
(42)

where function \(n(x^+)\) defines the one-dimensional target density that we take as constant, \(n_0=n(x^+)\) for \(0\le x^+\le L^+\), and 0 elsewhere. Moreover, function \(a(\mathbf{q})\) is the potential in momentum space which can be taken to be of Yukawa type, i.e. \(|a(\mathbf{q})|^2=m^2/(\mathbf{q}^2+m^2)^2\) with m being the Debye screening mass or inverse colour correlation length. In the eikonal limit, when \(\lambda ^+\rightarrow 0\) for a constant potential and constant one-dimensional target density, one recovers the standard MV expression for the two target field correlator. Using Eq. (42) one can integrate over the longitudinal coordinates that appear in the phases in non-eikonal Lipatov vertices. The final result of the non-eikonal single inclusive gluon production cross section in pp collisions then reads

$$\begin{aligned}&\frac{d\sigma }{d^2\mathbf{k}d\eta }\bigg |_{\mathrm{dilute}}^\mathrm{NonEik}= 4\pi \alpha _s\, C_A (N_c^2-1)\, g^2 {{\mathcal {G}}}_1^\mathrm{NE}(k^-;\lambda ^+)\nonumber \\&\quad \times \int \frac{d^2\mathbf{q}}{(2\pi )^2} \, \mu ^2\big [\mathbf{k}-\mathbf{q}, \mathbf{q}-\mathbf{k}\big ] \, L^i(\mathbf{k},\mathbf{q}) L^i(\mathbf{k},\mathbf{q})\, |a(\mathbf{q})|^2\ ,\nonumber \\ \end{aligned}$$
(43)

where we assume that the longitudinal width of the target is much larger than the colour correlation length, \(\lambda ^+\ll L^+\). In the cross section, Eq. (43), all the non-eikonal effects are encoded in the function \({{\mathcal {G}}}_1^\mathrm{NE}(k^-;\lambda ^+)\) which is defined as

$$\begin{aligned} {{\mathcal {G}}}_1^\mathrm{NE}(k^-;\lambda ^+)=\frac{1}{k^-\lambda ^+}\sin (k^-\lambda ^+), \end{aligned}$$
(44)

with \(k^-\equiv \mathbf{k}^2/2k^+\). In the limit of vanishing \((k^-\lambda ^+)\) we have

$$\begin{aligned} \lim _{k^-\lambda ^+\rightarrow 0}{{\mathcal {G}}}_1^\mathrm{NE}(k^-;\lambda ^+)=1, \end{aligned}$$
(45)

and the well known eikonal limit for the single inclusive gluon production in the dilute target limit is recovered. Therefore, function \({{\mathcal {G}}}_1^\mathrm{NE}(k^-;\lambda ^+)\) can be interpreted as the function that accounts for the relative importance of the non-eikonal effects with respect to the eikonal limit of the single inclusive gluon production in the dilute target limit.

Fig. 7
figure 7

The ratio of non-eikonal to eikonal single inclusive gluon production as a function of the rapidity of the produced gluon for different values of its transverse momenta at fixed correlation length \(\lambda ^+=0.5\) fm. Taken from [84]

In Fig. 7, the ratio of the non-eikonal to eikonal single inclusive gluon production cross sections, i.e. function \({{\mathcal {G}}}_1^\mathrm{NE}(k^-;\lambda ^+)\), is plotted as a function of rapidity for different values of the transverse momenta of the produced gluon at correlation length \(\lambda ^+=0.5\) fm. The results show that with increasing rapidity of the produced gluon, the effects of the non-eikonal corrections vanish as expected from analytical predictions. Up to rapidity \(\eta =2.5\), the relative importance of the corrections varies between 15 and 2% depending on the value of the transverse momenta.

Non-eikonal double inclusive gluon production cross section in pp scattering can be computed in a similar manner. The main difference between the single and double inclusive production is that one needs to compute the target average of the four field correlator. This can be accomplished by factorising the average of the colour fields of the target into all possible Wick contractions which can be written

$$\begin{aligned}&\big \langle A^-_a(x_1^+,\mathbf{q}_1)A^-_b(x_2^+,\mathbf{q}_2)A^-_c(x_3^+,\mathbf{q}_3)A^-_d(x_4^+,\mathbf{q}_4) \big \rangle _T \nonumber \\&\quad = \big \langle A^-_a(x_1^+,\mathbf{q}_1)A^-_b(x_2^+,\mathbf{q}_2)\big \rangle _T \big \langle A^-_c(x_3^+,\mathbf{q}_3)A^-_d(x_4^+,\mathbf{q}_4) \big \rangle _T \nonumber \\&\qquad + \big \langle A^-_a(x_1^+,\mathbf{q}_1)A^-_d(x_4^+,\mathbf{q}_4) \big \rangle _T \big \langle A^-_c(x_3^+,\mathbf{q}_3) A^-_b(x_2^+,\mathbf{q}_2) \big \rangle _T \nonumber \\&\qquad +\big \langle A^-_a(x_1^+,\mathbf{q}_1)A^-_c(x_3^+,\mathbf{q}_3)\big \rangle _T \big \langle A^-_b(x_2^+,\mathbf{q}_2)A^-_d(x_4^+,\mathbf{q}_4) \big \rangle _T,\nonumber \\ \end{aligned}$$
(46)

where each two target field correlator is defined in Eq. (42). The integrals over the longitudinal coordinates can be performed using this definition and the final result for the non-eikonal double inclusive gluon production cross section in pp scattering can be organised in the following way:

$$\begin{aligned}&\frac{d\sigma }{d^2\mathbf{k}_1d\eta _1d^2\mathbf{k}_2d\eta _2}\bigg |_{\mathrm{dilute}}^\mathrm{NonEik}=\alpha _s^2 \, (4\pi )^2 \, g^4\, C_A^2\, (N_c^2-1) \nonumber \\&\quad \times \int \frac{d^2\mathbf{q}_1}{(2\pi )^2} \frac{d^2\mathbf{q}_2}{(2\pi )^2} |a(\mathbf{q}_1)|^2 |a(\mathbf{q}_2)|^2 {{\mathcal {G}}}_1^\mathrm{NE}(k_1^-;\lambda ^+) {{\mathcal {G}}}_1^\mathrm{NE}(k_2^-;\lambda ^+) \nonumber \\&\quad \times \bigg \{ I^{(0)}_{\mathrm{2tr}}+ \frac{1}{N_c^2-1}\Big [I^{(1)}_{\mathrm{2tr}}+I^{(1)}_\mathrm{1tr}\Big ]\bigg \}, \end{aligned}$$
(47)

where the subscripts denote the single trace terms (\(I^{(i)}_\mathrm{1tr}\)) or the double trace term (\(I^{(i)}_{\mathrm{2tr}}\)) which are analogue to double dipole and quadrupole operators in pA scattering discussed previously. The explicit expressions for these contributions read

$$\begin{aligned} I^{(0)}_{\mathrm{2tr}}= & {} \mu ^2\big [ \mathbf{k}_1-\mathbf{q}_1,\mathbf{q}_1-\mathbf{k}_1\big ] \mu ^2\big [ \mathbf{k}_2-\mathbf{q}_2,\mathbf{q}_2-\mathbf{k}_2\big ] \nonumber \\&\quad \times L^i(\mathbf{k}_1,\mathbf{q}_1)L^i(\mathbf{k}_1,\mathbf{q}_1)L^j(\mathbf{k}_2,\mathbf{q}_2)L^j(\mathbf{k}_2,\mathbf{q}_2), \nonumber \\\end{aligned}$$
(48)
$$\begin{aligned} I^{(1)}_{\mathrm{2tr}}= & {} \Big \{ {{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k_2^-;L^+) \mu ^2\big [\mathbf{k}_1-\mathbf{q}_1, \mathbf{q}_2-\mathbf{k}_1\big ]\nonumber \\&\times \mu ^2\big [\mathbf{k}_1-\mathbf{q}_1, \mathbf{q}_2-\mathbf{k}_1\big ] L^i(\mathbf{k}_1,\mathbf{q}_1)L^i(\mathbf{k}_1,\mathbf{q}_2)\nonumber \\&\times L^j(\mathbf{k}_2,\mathbf{q}_2) L^j(\mathbf{k}_2,\mathbf{q}_1) \Big \} +({\underline{k}}_2\rightarrow -{\underline{k}}_2) \end{aligned}$$
(49)

and, finally,

(50)

In this setup (see Fig. 5), \(\mathbf{k}_1-\mathbf{q}_1\) and \(\mathbf{k}_2-\mathbf{q}_2\) are the transverse momenta of the two gluons in the projectile, \(\mathbf{k}_1\) and \(\mathbf{k}_2\) are the transverse momenta of the produced gluons in the final state, and \(\mathbf{q}_1\) and \(\mathbf{q}_2\) are the transverse momenta transferred from the target to the projectile during the interaction. By using the definition of function \(\mu ^2(\mathbf{p},\mathbf{k})\) given in Eq. (10) and the behaviour of the soft form factor, one can easily identify each term in the non-eikonal double inclusive gluon production cross section given in Eqs. (48), (49) and (50). Clearly, the term in \(I^{(0)}_{\mathrm{2 tr}}\) corresponds to the square of the single inclusive production and does not give any contribution to the correlations. The terms in \(I^{(1)}_{\mathrm{2 tr}}\) corresponds to the Bose enhancement of the target gluons since the soft form factor is peaked around \(\mathbf{q}_1=\mathbf{q}_2\). Finally, the terms in \(I^{(1)}_{\mathrm{1 tr}}\) contribute to the HBT correlations of the produced gluons and to Bose enhancement of the projectile gluons.

In the non-eikonal double inclusive gluon production cross section, two functions appear that account for the non-eikonal effects: \({{\mathcal {G}}}_1^\mathrm{NE}(k_1^-;\lambda ^+)\) presented in Eq. (44) and a new function \({{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k_2^-;L^+)\). This new function is defined as

$$\begin{aligned}&{{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k_2^-;L^+) \nonumber \\&\quad = \bigg \{ \frac{2}{(k_1^--k_2^-)L^+}\sin \bigg [ \frac{(k_1^--k_2^-)}{2}L^+\bigg ]\bigg \}^2 \end{aligned}$$
(51)

which in the eikonal limit, i.e. \(L^+\rightarrow 0\), goes to unity,

$$\begin{aligned} \lim _{L^+\rightarrow 0}{{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k_2^-;L^+)=1. \end{aligned}$$
(52)

Different from the eikonal double inclusive gluon production cross section, in the non-eikonal expression the mirror images are given by \({\underline{k}}_2\rightarrow -{\underline{k}}_2\) where \({\underline{k}}_2\equiv (k_2^+,\mathbf{k}_2)\). The mirror images of the terms that are accompanied by the function \({{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k_2^-;L^+)\) are now accompanied by \({{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,-k_2^-;L^+)\). However, as obvious from the definition given in Eq. (51), this function is not symmetric under this transformation. Moreover, in certain kinematic regimes the behaviour of \({{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k_2^-;L^+)\) differs completely from \({{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,-k_2^-;L^+)\). Particularly, in a kinematic region where \(k_1^-\sim k_2^-\), one gets

$$\begin{aligned} {{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k_2^-;L^+)\gg {{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,-k_2^-;L^+). \end{aligned}$$
(53)

This creates an asymmetry between the terms (\({\underline{k}}_1,{\underline{k}}_2\)) and their partners (\({\underline{k}}_2\rightarrow -{\underline{k}}_2\)). This asymmetry which comes from the non-eikonal effects reminds the asymmetry between the forward and backward peaks of the ridge structure observed in two particle production. Therefore, non-eikonal corrections break the accidental symmetry present in usual CGC calculations and can give rise to odd harmonics. In the remaining of the section we briefly examine the numerical relevance of this effect.

A detailed numerical analysis of the azimuthal structures in two particle correlations based on the non-eikonal double inclusive gluon production cross section given in Eq. (47) with Eqs. (48), (49) and (50), is performed in [85] where it is assumed that:

  1. 1.

    the colour sources inside the projectile have a Gaussian distribution such that \(\mu ^2(\mathbf{k}, \mathbf{q})=\mu ^2(2\pi )^2\delta ^{(2)}(\mathbf{k}+\mathbf{q})\), with \(\mu \) being the width of the Gaussian;

  2. 2.

    the Yukawa type potential that defines the target field correlators is given by \(|a(\mathbf{q})|^2=\mu _T^2/(\mathbf{q}^2+\mu ^2_T)^2\), with \(\mu _T\) being an infrared regulator analogous to a Debye mass;

  3. 3.

    the transverse area of the projectile \(S_\perp \) is defined through \((2\pi )^2\delta ^{(2)}(\mathbf{q}-\mathbf{q})\rightarrow S_{\perp }\).

Fig. 8
figure 8

Two particle azimuthal harmonics computed in the non-eikonal Glasma graph approach using the definition Eq. (25). The values are calculated taking \(\mu _T=0.4 \; \mathrm{GeV}\), \(\mu _P=0.2 \; \mathrm{GeV}\) and \(p_T^{ref}=1\; \mathrm{GeV}\) for different values of centre-of-mass energy and different gluon rapidities \(\eta _1=\eta _2=\eta \). Taken from [85]

In this analysis, function \({{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k^-_2;L^+)\) that encodes the non-eikonal effects defined in Eq. (51) is rewritten as

$$\begin{aligned}&{{\mathcal {G}}}_2^\mathrm{NE}(k_1^-,k_2^-;L^+) \nonumber \\&\quad = \left\{ \frac{\sqrt{2}}{(k_1e^{-\eta _1}-k_2e^{-\eta _2})L^+} \sin \left[ \frac{(k_1e^{-\eta _1}-k_2e^{-\eta _2})}{\sqrt{2}}L^+\right] \right\} ^2\nonumber \\ \end{aligned}$$
(54)

using \(k^-=k^2/2k^+\), \(k^+=ke^\eta /{\sqrt{2}}\), \(k=|\mathbf{k}|\).Footnote 13

In Fig. 8 the azimuthal harmonics up to \(v_5\) are computed by using the definition given in Eq. (25). \(p_T^{ref}\) is taken to be \(1\, \mathrm{GeV}\) for different values of \(\sqrt{s_{NN}}\) and \(\eta _1=\eta _2=\eta \). The plot shows that the value of the odd harmonics decreases with increasing centre-of-mass energy at fixed rapidity \(\eta \). This behaviour is the natural outcome of the fact that non-eikonal corrections become smaller with the increasing Lorentz gamma factor. Therefore, one can conclude that the non-eikonal corrections can be negligible for collisions at high centre-of-mass energy such as the ones at the LHC but they can be important for collisions at RHIC with \(\sqrt{s_{NN}}\le 200\, \mathrm{GeV}\). On the other hand, at any fixed energy the value of odd azimuthal harmonics decreases with increasing rapidity \(\eta \). This behaviour is also expected, since the value of the odd harmonics is directly related to the non-eikonal corrections. The eikonal expansion parameter can be written as \(p_TL^+e^{-\eta }\) in terms of rapidity and, therefore, non-eikonal corrections (and thus the value of odd harmonics) decrease with increasing rapidity and vanish completely in the strict eikonal limit.Footnote 14

Fig. 9
figure 9

Scaling of \(v_n(L^+)/v_n(1.5\, \mathrm{fm})\) with \(L^+\). Taken from [85]

In Fig. 9, the ratio \(v_n(L^+)/v_n(1.5\, \mathrm{fm})\) is plotted as a function \(L^+\) which reveals a very interesting feature of the effects of non-eikonal corrections on azimuthal harmonics: odd harmonics depend strongly on the size of the target while even ones are almost independent of it. Even though the explicit relation between multiplicity and \(L^+\) requires a more dedicated study, the scaling of the odd harmonics with \(L^+\) in Fig. 9 qualitatively resembles the results of the analysis performed in [82] where it is shown that the value of \(v_3\) increases with the increasing multiplicity.

4 Non-CGC explanations

Besides explanations to the ridge phenomenon based on the CGC, there are others that address its origin in the initial state of the collision or, at least, do not demand hydrodynamics or transport at work. It must be noted that the existence of long range rapidity correlations was discussed long ago as a consequence of multiple scattering, see [93, 94].

This approach was pushed forward in string models for multiparticle production, see e.g. [95]. Later on, several models that consider string interactions were argued to lead to azimuthal asymmetries: string percolation [96,97,98] with the creation of azimuthally anisotropic strong chromoelectric fields, colour reconnection [99] that is able to produce some of the QGP-like features observed in pp, and string repulsion [100, 101]. It is not clear whether the dynamics contained in these approaches can be considered as pure initial state but they offer a mechanism to produce the ridge in collisions between small systems that does not require any explicit final state rescattering, see [102] for a model that explicitly requires parton and/or hadron rescattering to build azimuthal asymmetries even with just two strings. In all these approaches, particle production from a single string is still isotropic and the anisotropy is built after string breaking.

A string-based model is also proposed in [103, 104]. There, valence diquark-quark flux tubes or strings in the incoming protons overlap and produce more particles in the transverse than in the longitudinal direction of the flux tube. Such anisotropic particle production leads to azimuthal asymmetries and the prediction has been made that it should also be visible in photoproduction, with particle production becoming maximal in the plane of the deflected electron (in ep collisions) or proton (in UPCs).

It should be noted that all these approaches are inspired in the string behaviour of the QCD interaction in the non-perturbative domain, in contrast to the CGC that relies on perturbation theory for a small coupling constant. Indeed already in the framework of Reggeon field theory, some ideas have been pushed [105] on the spatial variation of the transverse density in the hadron that resemble those in the CGC. Or CGC arguments have been extended to the soft physics domain and applied to describe azimuthal correlations, see [106, 107] and subsequent papers of this group.Footnote 15

Finally, let us indicate that azimuthal asymmetries arise in several processes when the nucleon is studied and characterised beyond collinear parton densities, as in the framework of Wigner distributions and TMD parton densities, see [110, 111]. Azimuthal asymmetries then arise in final observables like dijet production in DIS [112,113,114]. But, although these calculations are often performed in a framework close to that of the CGC which is related with the TMD framework at small x [115], it goes beyond the standard CGC context to link with other physics like spin.

5 Summary and discussions

In this manuscript we have discussed the explanations that are currently proposed to describe particle correlations, the ridge, observed in experimental data in small collision systems, pp and pA, from the initial state point of view. Our main focus has been those weak coupling explanations based on the CGC. We have assumed that correlations among partons in the initial stage leave an imprint on those among particles in the final state, i.e. they are not washed out by final strong final state interactions and hadronisation.

First, we have focused on the standard eikonal treatment within the glasma graph approximation which is valid for collisions between dilute objects – pp. We have reviewed the studies which have shown that this approximation encodes two different type of contributions, namely the Bose enhancement of both projectile and target gluons and also HBT correlations of the produced final gluons.

We have summarised the procedure that should be adopted to extend the validity of glasma graph approximation from dilute-dilute to dilute-dense (i.e. from pp to pA) collisions by taking into account the multiple scattering effects in the dense target. We have shown that the structure of the double inclusive gluon production cross section is symmetric under (\(\mathbf{k}_2\rightarrow -\mathbf{k}_2\)), which is known as the accidental symmetry of the CGC. Since this symmetry is the reason for vanishing odd harmonics in the CGC framework, we have discussed the suggested mechanisms to break this accidental symmetry.

In particular, we have focused on a specific mechanism to break this symmetry which is based on going beyond the standard eikonal approximation and including the subeikonal corrections that are due to the finite longitudinal width of the target. We have argued that such non-eikonal corrections, when included in the glasma graph approach to two particle correlations, successfully generate non-zero odd harmonics in specific kinematics. We would like to emphasise here that we make no attempt to compare the results with experimental data but only address the existence and size of the non-eikonal effects on the azimuthal structure. As expected from non-eikonal corrections, their value and thus that of the odd harmonics decrease rapidly with increasing centre-of-mass energy. This decrease is strong since the analysis is performed for a dilute target – a slower decrease of the size of the odd harmonics with increasing energy could be expected in a dilute-dense collision. Besides, the treatment of such non-leading eikonal corrections shows explicitly the link of the formalisms used in CGC and jet quenching calculations.

At this point, we should comment briefly on the comparison with experimental data. The main characteristics concerning azimuthal asymmetries in small systems observed in experiment [1,2,3,4,5, 7,8,9,10, 13, 16,17,19] are:

  • The even and odd harmonics extracted using correlations between two and more particles, are of similar size to those found in larger systems.

  • They show the same dependence on the mass of the measured hadron as found in larger systems (see [116] for an approach in the CGC).

  • Even harmonics show a much weaker dependence on the multiplicity in the event than odd harmonics.

  • \(v_2\) and \(v_3\) found in pAu, dAu and \(^3\)HeAu collisions at RHIC show, for central (head-on) collisions, the ordering \(v_2^{p\mathrm {Au}}<v_2^{d\mathrm {Au}}\approx v_2^{^3\mathrm {HeAu}}\), \(v_3^{p\mathrm {Au}}\approx v_3^{d\mathrm {Au}}< v_3^{^3\mathrm {HeAu}}\).

  • Measurements of many particle cumulants show evidence of collectivity. For example, four particle cumulants \(c_2\{4\}=\langle e^{i2(\phi _1+\phi _2-\phi _3-\phi _4)}\rangle -2\langle e^{i2(\phi _1-\phi _2)} \rangle ^2\) (\(v_2\{4\}=[-c_2\{4\}]^{1/4}\)) change sign from positive to negative with increasing associated multiplicity, with a smooth behaviour from small to large systems and from smaller to larger energies.

In the glasma graph approximation, several studies were done that describe pp data in a reasonable manner [55,56,57,58]. Later on, these studies were extended to pA with diverse degree of modelling [59, 60]. Then, odd azimuthal harmonics were introduced following [80, 81], with the corresponding numerical studies and a comparison with data performed in [82, 83]. In these latter studies a successful comparison with RHIC and LHC data was initially claimed, which was later corrected after the criticism in [117]. Nevertheless, it must be stated that none of the numerical calculations can be considered as a full implementation of the theoretical framework and that some results are still to be clarified from an analytical point of view, e.g. those in [59, 60] about the second Fourier coefficient defined through four particle correlations (\(v_2\{4\}\)) where the mentioned change of sign in \(c_2\{4\}\) is attributed to multiple scattering beyond glasma graphs.

Finally, we have also shortly commented on approaches that are not based on, or go beyond, CGC ideas, to study the two particle correlations from the initial state.

To conclude, let us indicate that future experimental programmes and facilities [27, 118,119,120,121] will address the physics of small systems and the transition from small to large, particularly the onset and understanding of collectivity which is a central question in QCD at high energies.