Particle correlations from the initial state

The observation in small size collision systems, $pp$ and $p$A, of strong correlations with long range in rapidity and a characteristic structure in azimuth, the ridge phenomenon, is one of the most interesting results obtained at the Large Hadron Collider. Earlier observations of these correlations in heavy ion collisions at the Relativistic Heavy Ion Collider are standardly attributed to collective flow due to strong final state interactions, described in the framework of viscous relativistic hydrodynamics. Even though data for small size systems are well described in this framework, the applicability of hydrodynamics is less well grounded and initial state based mechanisms have been suggested to explain the ridge. In this review, we discuss particle correlations from the initial state point of view, with focus on the most recent theoretical developments.


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 π, 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]. They have also been observed in photopro-Send offprint requests to: The ridge: 3 • Two-particle correlations in pp and pPb at the LHC show features that in AA are attributed to final state interactions describable by viscous relativistic hydrodynamics and interpreted as a signal of equilibration.
• EKT and AdS/CFT: hydro works even for large momentum anisotropies.
• What about a non-hydro initial-state explanation? (anyway long range rapidity correlations must come from the very early times…). 1609.06213 N. Armesto,18.04.2018 -Multi gluon correlations in the CGC: 1. Introduction. Fig. 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. arXiv:2004.08185v1 [hep-ph] 17 Apr 2020 duction on Pb in ultraperipheral collisions (UPCs) at the LHC [19]. Their existence in smaller systems like e + e − collisions [20] at the Large Electron-Positron collider and deep inelastic scattering (DIS) events in ep at the Hadron-Elektron-Ringanlage [21] 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 [22,23,24,25,26].
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 [27,28]. The dynamics leading to such a situation, called hydrodynamisation, is unclear [29] and both strong and weak coupling explanations have been proposed, see e.g. [30] and refs. therein. Furthermore, the hydrodynamic description seems to hold for large anisotropies, i.e. rather far from local equilibrium.
The hydrodynamic description of the azimuthal asymmetries in pp and pPb collisions at the LHC is successful [28,31]. 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. [32] 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 [28] and in those that consider a weak coupling quasiparticle picture in transport frameworks [33].
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 [34]. 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 Con-densate (CGC) effective theory, see the reviews [35,36,37] and the book [38]. This will be the subject of Sections 2 and 3. We will then review other explanations inspired in QCD in Section 4, to end with a summary and discussions in Section 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 [39,40,41] and the "glasma graph approach" introduced in [34,42,43] 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 [35,36,37]. Projectile partons lying closer than 1/Q s contribute mainly to particle production in the region of transverse momentum p T Q s . Therefore, this mechanism should mainly be applicable in that transverse momentum region. 1 Apart from these two approaches, there are also other CGC-based mechanisms to describe the two particle correlations from the initial state. In [44,45], it is argued that long range rapidity correlations can be explained by the spatial variation of the partonic density in the target. On the other hand, in [46], the origin of two particle correlations is explained through the rapidity evolution of dipole operators by breaking the mean field approximation. Correlations in the hadron wave functions as described in the CGC have been recently considered in [47] but with the aim of providing initial conditions for hydrodynamic evolution beyond simple energy, flow and particle density transverse profiles. 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 ∼ 1/Q s [48,49,50]), 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 [51,52,53]. They show correlations that go beyond leading number of colours, see the discussion below, and lead to odd harmonics, see Section 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 [54,55,56,57,58,59], but the physics behind this approach was not clear. This issue has been studied in [60] 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.
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 with a † i (p) the creation operator of the boson and i = 1, 2, . . . , N . The mean particle densityñ is defined as the expectation value of the number operator in this state: The two particle correlator in momentum space C(p, k) is defined in a similar way and can be calculated in a trivial manner: 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 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 2 . 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).
Type A diagrams describe the case when two gluons with transverse momenta k 1 and k 2 scatter independently on the target, acquiring transfer of transverse momentum p − k 1 and q − k 2 so that the outgoing gluons have transverse momenta p and 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 where |in is the wave function of the incoming projectile and N (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, ∆η 1/α s . The evaluation of the expectation value of any operator in the incoming projectile state requires a two averaging procedure in the CGC. In [60], averaging over the valence colour charge density is performed first which leads to the density matrix operatorρ 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: where S ⊥ 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, 4 Tolga Altinoluk, Néstor Armesto: Particle correlations from the initial state

Expansion of the background propagator
Perturbative expansion around free classical path: Tolga Altinoluk Next-to-eikonal corrections in the CGC Fig. 2. Glasma graphs for two gluon inclusive production before averaging over the projectile colour charge density ρ. 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).
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 ∼ Q s , i.e. |p| ∼ |q| ∼ Q s , then the production amplitude is dominated by the contributions |k 1 | ∼ |k 2 | ∼ 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 (p − k 1 )N (q − k 2 ). Since, in this kinematics, the momentum transfers from the target ( |p − k 1 | ∼ |q − k 2 | 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 [61] where the correlations between the produced quarks were studied. The results in [60] 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 [61] 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. η 1 − η 2 1. Moreover, a large contribution comes from the situation where the transverse momenta of the produced quarks p and 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. |p| ∼ |q| Q s 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 [61] for the full expressions): 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 quarkquark 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 gluons 3 [64]. 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 e background propagator sion around free classical path: Tolga Altinoluk Next-to-eikonal corrections in the CGC Fig. 3. Glasma graph diagrams (after averaging over the projectile colour charges) that lead to HBT correlations.
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: |p ± q| ∼ 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 [64] 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 (η 1 η 2 ) or when the difference between them is large (|η 1 − η 2 | 1). To sum up, the correlation function C(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: 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(p, q) 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(p, q) HBT represents the HBT correlations in the glasma graph approach which directly introduces correlations between the final state gluons. Both C(p, q) BE and C(p, q) 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 2 s S ⊥ ). 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.
In the explicit calculations in the glasma graph approach to double inclusive particle production 4 , the averaging over the target configurations that leads to the 4 Three and four gluon inclusive production are considered in [65,66] within the glasma graph approach. Fig. 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. 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 (p − k 1 ) and N (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 collisions 5 . In [67], 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 6 .
Apart from taking into account the multiple scattering effects in [67], 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: ρ a1 (x 1 )ρ a2 (x 2 )ρ b1 (y 1 )ρ b2 (y 2 ) P . Here, (x i , a i ) and (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 [70,71] 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 The function µ 2 (k, p) characterises the structure of the projectile. It can be written as where F (k + p)R 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 p and q, while the projectile gluons carry transverse momenta k 1 and k 2 .
In this case, µ 2 (p, q) gives a maximal contribution when p + q = 0 which clearly can be identified as the HBT correlations of the produced gluons. µ 2 (k 1 , k 2 ) is peaked when 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: at the amplitude level. Here, T a R is the SU (N c ) generator in the representation 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 which have to be averaged over the target field distributions. The dipole and the quadrupole operators are defined as 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 2 s R 2 ) to the integral on transverse coordinates 7 -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 Then, by using the function µ 2 (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 8 .

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 Section 2, double inclusive gluon production is computed in pA collisions in [67]. The production cross section of two gluons with rapidities η 1 and η 2 , and transverse momenta k 1 and k 2 (see Fig. 5) reads i rown-Twiss (HBT) correlations between gluons far separated in rapidity. where Here, function µ 2 defines the structure of the projectile, Eq. (10). On the other hand, function L i (k, q) is the eikonal Lipatov vertex which is defined as Note that, as discussed in Section 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): A convenient way to study two particle correlations is through a Fourier decomposition into harmonics defined in for the azimuthal angle ∆φ between the produced gluons with transverse momenta k 1 and k 2 . When Fourier expanded, the double inclusive gluon spectrum Eq. (17) can be written as 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 .
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 [40], and it is shown to be due to the symmetry (k 2 → −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 ways 9 have been proposed to solve this problem 10 . On the one hand, the projectile and target can be characterised by a more involved structure than that considered in the usual MV averages [39,40,41,51,52,53].
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 [77,78] 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 [79,80].
Finally, inclusive gluon production is usually studied within the eikonal approximation in the CGC framework. In recent studies [81,82], 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 nonvanishing odd harmonics.

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 µ a (x), and the target by the colour field A µ a (x) that is given as 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: Neglecting the (+) and transverse components of the colour field of the target.
Neglecting the x − dependence in the colour field of the target. 3. A µ (x) ∝ δ(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 [83,84], 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 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 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 µ 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 µ 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 B: Here, λ, a and k = (k + , k) are the polarization, colour and momentum of the produced gluon 11 . For a target with finite longitudinal width, the gluon production amplitude M a λ (k, 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 [85]. 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 (x, y) as 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 G µν k + (x, y). Then, the (i−) component of the retarded background gluon propagator reads The background scalar propagator G ab k + (x, y) satisfies the scalar Green's equation whose solution formally can be written as a path integral × e ik + 2 x + y + dz +ż2 (z + ) U ab x + , y + ; z(z + ) ,

with the Wilson line
U ab x + , y + ; z(z + ) = P + exp ig 11 Hereafter, we use the underline notation to indicate that for coordinates x = (x + , x) and for momentum k = (k + , k).
following the Brownian trajectory z(z + ). In the limit of vanishing longitudinal width, x + − y + → 0, the background scalar propagator G ab k + (x, 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 G ab k + (x, 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.
In order to perform an eikonal expansion of the background scalar propagator G ab k + (x, 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, where the transverse positions at step n are on the straight line between the initial and final points, and the perturbation u n satisfies the boundary conditions u 0 = 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 x − y, since x − 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 G ab k + (x, y) reads The first term on the right hand side of Eq. (36) is the standard Wilson line defined in Eq. (11). O (x + −y + )/k + terms are the first order corrections to the strict eikonal limit which we refer to as next-to-eikonal (NEik) corrections. Similarly, O (x + − y + ) 2 /(k + ) 2 terms are the second order corrections and they are referred to as nextto-next-to-eikonal (NNEik) corrections. The terms that are denoted as U [α,β] (x + , y + ; y) are the decorated Wilson lines which only appear beyond strict eikonal order. The first subscript α in the decorated Wilson lines stands for the order of expansion around the classical path while the second subscript β stands for the order of the expansion around the initial transverse position 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 × U ac (x + , z + ; y) ig T e cd ∂ y i A − e (z + , y) U db (z + , y + ; y).
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 [84]). 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 (x, y) ab and, therefore, the scalar background propagator G ab k + (x, y) are the main building blocks for computing cross sections in high energy pA collisions. In [83,84], 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 [81,82].
In [86], 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 + , y). For example, the first decorated Wilson line given in Eq. (37) reduces to 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 The first term on the right hand side of Eq. (39) corresponds to the strict eikonal limit, thus it gives the eikonal O(g 2 ).
To proceed, we analyze gluon production in p-A collisions in the quark initiated channel and compute the Lipatov rtex, which is an e↵ective vertex that takes into account all the real contributions to gluon production. For that e needs to sum the amplitudes where the gluon is emitted before, during and after the interaction with the field as own in Fig. 1. Our setup is such that the right moving quark with momentum p+k q is generated by some function J(p+k q) = p + + k + q + ) at x + 0 = 1 and (x 0 , x 0? ) = 0, and then interacts with the classical gluon field A µ (x) generated one scattering source located at x 1 , picking up a momentum q. However, since we are interested in non-eikonal rrections, we consider A µ (x) with an x + dependence which has a finite support instead of treating it as a shockwave x + = 0, but we still assume that there is no dependence on x . That is, the new form of Eq. (1) is , in momentum space, rthermore, we assume that the outgoing quark has a large momentum p + compared to all other momenta in the ocess. The general strategy in this case is to keep the leading terms in +-momenta in the numerator algebra, while king the full phase corrections coming from the integration of the denominators, see below, as done in the Furry proximation and its non-abelian generalization [75]. We start by computing diagram A where the gluon is emitted with momentum k before the quark interaction with e target field as shown in Fig. 2. Using the Feynman rules, we find that the amplitude for fixed gluon and final ark momenta is th t a the SU (N c ) generators in the fundamental representation. Since p + is the largest momentum in our problem, we approximate / p / q ⇡ / p and / p + / k / q ⇡ / p and write 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 [81] 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 where k − ≡ k 2 /(2k + ). This structure was observed in the context of jet quenching in [87,88,89] previously, however the identification of this building block for its use to include non-eikonal corrections in CGC calculations is done in [81] for the first time, further illustrating the close relation between CGC and jet quenching calculations [85].
It is now straightforward to compute the non-eikonal single inclusive gluon production in pp (i.e. dilute-dilute) collisions which formally reads 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 λ + . In this case, the two target field correlator reads where function n(x + ) defines the one-dimensional target density that we take as constant, n 0 = n(x + ) for 0 ≤ x + ≤ L + , and 0 elsewhere. Moreover, function a(q) is the potential in momentum space which can be taken to be of Yukawa type, i.e. |a(q)| 2 = m 2 /(q 2 + m 2 ) 2 with m being the Debye screening mass or inverse colour correlation length. In the eikonal limit, when λ + → 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 where we assume that the longitudinal width of the target is much larger than the colour correlation length, λ + L + . In the cross section, Eq. (43), all the non-eikonal effects are encoded in the function G NE 1 (k − ; λ + ) which is defined as with k − ≡ k 2 /2k + . In the limit of vanishing (k − λ + ) we have lim and the well known eikonal limit for the single inclusive gluon production in the dilute target limit is recovered. Therefore, function G NE 1 (k − ; λ + ) 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.
In Fig. 7, the ratio of the non-eikonal to eikonal single inclusive gluon production cross sections, i.e. function G NE 1 (k − ; λ + ), is plotted as a function of rapidity for different values of the transverse momenta of the produced gluon at correlation length λ + = 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 η = 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 of non-eikonal to eikonal single inclusive gluon production cross sections, (33), as a function of the transverse roduced gluon for di↵erent values of the correlation length + , at fixed pseudorapidity ⌘ = 2. tion length + . In the limit of vanishing transverse momenta of the produced gluon, the non-eikonal sections coincide and the ratio becomes one as expected. The ratio shows up to 20% relative weight l corrections for + = 1 fm, for smaller values of + the results show a suppression from a few to up ave plotted the ratio of the non-eikonal to eikonal single inclusive gluon production cross sections, n of pseudorapidity for di↵erent values of the transverse momenta of the produced gluon at a fixed h + = 0.5 fm. The ratio of the non-eikonal to eikonal cross sections goes to one with increasing s expected, since the relative importance of the non-eikonal corrections should vanish for large values show that up to pseudorapidity ⌘ = 2.5, depending on the value of the transverse momenta of the the relative weight of the non-eikonal corrections can vary roughly between 15% and 2%. These ur analytical predictions for the importance of the non-eikonal corrections in certain kinematical B. Double inclusive gluon production beyond the eikonal approximation ion we consider double inclusive gluon production beyond the eikonal approximation. Our strategy n is the same as the calculation performed for single inclusive gluon production in the previous Fig. 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 λ + = 0.5 fm.
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 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: where the subscripts denote the single trace terms (I (i) 1tr ) or the double trace term (I (i) 2tr ) which are analogue to double dipole and quadrupole operators in pA scattering discussed previously. The explicit expressions for these contributions read and, finally, In this setup (see Fig. 5), k 1 − q 1 and k 2 − q 2 are the transverse momenta of the two gluons in the projectile, k 1 and k 2 are the transverse momenta of the produced gluons in the final state, and q 1 and q 2 are the transverse momenta transferred from the target to the projectile during the interaction. By using the definition of function µ 2 (p, 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) 2tr corresponds to the square of the single inclusive production and does not give any contribution to the correlations. The terms in I (1) 2tr corresponds to the Bose enhancement of the target gluons since the soft form factor is peaked around q 1 = q 2 . Finally, the terms in I (1) 1tr 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: G NE 1 (k − 1 ; λ + ) presented in Eq. (44) and a new function G NE 2 (k − 1 , k − 2 ; L + ). This new function is defined as which in the eikonal limit, i.e. L + → 0, goes to unity, Different from the eikonal double inclusive gluon production cross section, in the non-eikonal expression the mirror images are given by k 2 → −k 2 where k 2 ≡ (k + 2 , k 2 ). The mirror images of the terms that are accompanied by the function G NE 2 (k − 1 , k − 2 ; L + ) are now accompanied by G NE 2 (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 G NE 2 (k − 1 , k − 2 ; L + ) differs completely from G NE 2 (k − 1 , −k − 2 ; L + ). Particularly, in a kinematic region where k − 1 ∼ k − 2 , one gets This creates an asymmetry between the terms (k 1 , k 2 ) and their partners (k 2 → −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, noneikonal 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 [82] where it is assumed that: 1. the colour sources inside the projectile have a Gaussian distribution such that µ 2 (k, q) = µ 2 (2π) 2 δ (2) (k + q), with µ being the width of the Gaussian; 2. the Yukawa type potential that defines the target field correlators is given by |a(q)| 2 = µ 2 T /(q 2 + µ 2 T ) 2 , with µ T being an infrared regulator analogous to a Debye mass; 3. the transverse area of the projectile S ⊥ is defined thro- . 7 we plot the cross section eq. (13) without the prefactors outside the curly brackets (that we call normalized city) against ⌘ and using ⌘ 1 = 0, k 1 = 1 GeV and k 2 = 1.2 GeV. We can see again that the di↵erences the forward and backward peaks are visible up to 2.5 pseudorapidity units.
clusions is manuscript we have analyzed the e↵ect on the non-eikonal corrections stemming from relaxing the shockproximation for the target which, therefore, acquires a finite length, on the two gluon inclusive cross section GC. We work in the Glasma graph approximation suitable for collisions between dilute objects (pp). While esponding expressions were derived in a previous publication [76], here we focus on the numerical implemenor which several model assumptions are made. We make no attempt to compare with experimental data but dress the existence and size of the non-eikonal e↵ects on the azimuthal structure. In this analysis, function G NE 2 (k − 1 , k − 2 ; L + ) that encodes the non-eikonal effects defined in Eq. (51) is rewritten as using k − = k 2 /2k + , k + = ke η / √ 2, k = |k| 12 . In Fig. 8 the azimuthal harmonics up to v 5 are computed by using the definition given in Eq. (25). p ref T is taken to be 1 GeV for different values of √ s N N and η 1 = η 2 = η. The plot shows that the value of the odd harmonics decreases with increasing centre-of-mass energy at fixed rapidity η. 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 √ s N N ≤ 200 GeV. On the other hand, at any fixed energy the value of odd azimuthal harmonics decreases with increasing rapidity η. 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 T L + e −η 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 13 . . The parameters used for this plot are µ P = 0.2 GeV and ⌘ 1 = ⌘ 2 = 1.5. The dashed lines are the result using µ 2 (k 1 , k 2 ) / (2⇡) 2 (2) (k 1 k 2 ) and the co We explore how the non-eikonal corrections break the accidental forward-backward symmetr CGC calculations, and thus lead to sizeable odd harmonics. We discuss the di↵erent contributions: of the projectile and target wave functions and HBT, and check the stability of the qualitative beh against variations in the functional forms and parameters in the model assumptions. We find all even and all odd harmonics with respect to the length of the target, with even harmonics b odd ones growing with increasing length. The non-eikonal corrections vanish with increasing ener In Fig. 9, the ratio v n (L + )/v n (1.5 fm) is plotted as a function L + which reveals a very interesting feature of the 12 Assuming that L is the size of the target in its rest frame, in the centre-of-mass frame L + is taken as where A is the mass number of the nucleus and γ √ sNN /(2mN ) accounts for the Lorentz contraction in the centre-of-mass frame. Moreover, for the numerics the gluonic size of the projectile is taken to be Bp = 4 GeV −2 , the transverse size of the projectile is assumed to be S ⊥ = 2πBp ≈ 9.8 mb and the size of the target in its rest frame for a Pb nucleus is taken to be L = 12 fm. Finally, the number of colours is taken to be Nc = 3 and the colour correlation is set to be λ + = 0. 13 In Fig. 8, the unrealistic peaks that account for the HBT contributions are due to the use of µ 2 (k, q) ∝ δ (2) (k + q). In a more realistic treatment, µ 2 can be chosen as Gaussian which would peak around k + q = 0 and show a bell shape behaviour. 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 [79] where it is shown that the value of v 3 increases with the increasing multiplicity.

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 [90,91].
This approach was pushed forward in string models for multiparticle production, see e.g. [92]. Later on, several models that consider string interactions were argued to lead to azimuthal asymmetries: string percolation [93,94,95] with the creation of azimuthally anisotropic strong chromoelectric fields, colour reconnection [96] that is able to produce some of the QGP-like features observed in pp, and string repulsion (shoving) [97]. 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 [98] 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 [99,100]. 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 [101] 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 [102,103] and subsequent papers of this group 14 .
14 There are also attempts to describe the near side ridge as a consequence of the momentum kick given by the leading parton 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 [106,107]. Azimuthal asymmetries then arise in final observables like dijet production in DIS [108,109,110]. 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 [111], it goes beyond the standard CGC context to link with other physics like spin.

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 (k 2 → −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.
to medium constituents [104], with a medium already present in pp collisions, or to minijets [105].
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 nonleading 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. In the glasma graph approximation, several studies were done that describe pp data in a reasonable manner [54,55,56,57]. Later on, these studies were extended to pA with diverse degree of modelling [58,59]. Then, odd azimuthal harmonics were introduced following [77,78], with the corresponding numerical studies and a comparison with data performed in [79,80]. In these latter studies a successful comparison with RHIC and LHC data is claimed, which has been subject to intense scrutiny [112]. 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 [58,59] about the second Fourier coefficient defined through four particle correlations (v 2 {4}).
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 [113,26,114,115,116] 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.