Collectivity in small collision systems : an initial state perspective

Measurements of multi-particle correlations in the collisions of small systems such as $p+p$, $p/d/^3He+A$ show striking similarity to the observations in heavy ion collisions. A number of observables measured in the high multiplicity events of these systems resemble features that are attributed to collectivity driven by hydrodynamics. However alternative explanations based on initial state dynamics are able to describe many characteristic features of these measurements. In this brief review we highlight some of the recent developments and outstanding issues in this direction.

Measurements of multi-particle correlations in the collisions of small systems such as p + p, p/d/ 3 He + A show striking similarity to the observations in heavy ion collisions. A number of observables measured in the high multiplicity events of these systems resemble features that are attributed to collectivity driven by hydrodynamics. However alternative explanations based on initial state dynamics are able to describe many characteristic features of these measurements. In this brief review we highlight some of the recent developments and outstanding issues in this direction.

I. INTRODUCTION
Relativistic collisions of hadrons and nuclei at the modern colliders provide unique testing ground for QCD at high energies. Over a decade of experimental measurements for a wide range of energies and collision systems have been dedicated to study the properties of the matter formed in such collisions. A number of measurements have provided strong indications that the QCD matter formed in the collisions of two nuclei (A+A) behave like a strongly interacting fluid that exhibits collectivity [1][2][3][4]. Consistent measurements of strong radial and anisotropic flow, jet quenching etc. have convincingly established such properties of the medium. Small system collisions have initially been thought of providing benchmarking measurements for the observations in heavy ion collisions.
Very recently several striking observations have been made in p+p, p+A, d+A and 3 He + A collisions [5][6][7][8][9][10][11][12][13]. Observations in high multiplicity events for such small collision systems seem to resemble features that are common to A+A collisions and very often attributed to decisive signatures of collectivity due to hydrodynamic evolution. Such observations include appearance of azimuthal correlations that extend in long-range rapidity known as ridge, which is also quantified in terms of Fourier harmonic coefficient of v n , strong multiplicity dependence of mean transverse momentum, HBT radii and many others.
In this review we provide a brief overview on the experimental findings and outline a general theory perspective on the collective phenomena observed in small systems. Subsequently, we focus specifically on theoretical explanations based on the initial state dynamics. We discuss different theoretical approaches to the calculations and contrast the results with experimental findings. We conclude with a brief summary of the present status and perspectives for future studies.

A. Experimental Overview
The current debate of collectivity in small systems was triggered by the discovery of ridge like correlations that extend over a long range in rapidity in high multiplicity p+p collisions by the CMS collaboration [5]. Such long range structure of azimuthal correlations was previously seen in heavy ion collisions at RHIC [44][45][46] and LHC [47] and generally attributed to the nearly boost invariance structure of the azimuthal correlations driven by hydrodynamic flow. At the same time causality arguments suggest [48] that such correlations must develop at the very early stages of the collisions, indicating the strong influence of initial state dynamics on such observable. While the observation of similar ridge like structures made in high multiplicity p+Pb collisions at the LHC was expected [6][7][8], a surprise was that the signal strength at the same multiplicity was stronger in p+Pb compared to p+p. The latest inclusion to these measurements are the highest energy p+p collisions so far at 13 TeV [12,13]. Including the 13 TeV and 7 TeV data both ATLAS and CMS collaborations have demonstrated that the strength of ridge correlations (characterized by the near side yield) at a given multiplicity is independent of collision energy. However, a more pronounced energy dependence appears to be present in higher-order cumulants such as c 2 {4} = −v 4 2 {4} which unlike at lower energies appears to exhibit a clear sign change as a function of multiplicity in 13 TeV p + p collisions [49]. Meanwhile RHIC, being a versatile machine, has collided 3 He + Au [11] along with d + Au [9,10] and p + Au [50] to perform similar measurements of azimuthal correlations. Such measurements indicate the presence of significant v 2 and v 3 in p+Au, d+Au, 3 He+Au, the systematics of which is very similar to what is commonly seen in A+A collisions. Analyses common to A+A have been repeated in both p+Pb and p+p at LHC by triggering on high multiplicity events. Measurements of higher order harmonics of azimuthal anisotropy v n , mass dependence of v n , higher arXiv:1611.00329v2 [hep-ph] 17  tracking efficiencies using MC simulations and can be found in [23]. The combined geometriefficiency for track reconstruction exceeds 50% nd |η| < 2.4. The efficiency is greater than 90% for p T > 0.6 GeV/c. For the multiplicity range r no dependence of the tracking efficiency on and the rate of misreconstructed tracks reevel.
, pPb and peripheral PbPb collisions using the ydjet event generators, respectively, yield efactors that vary due to the different kinematic ons for the particles produced in these gene resulting correction factors from one of the ated data from one of the others gives assotions that agree within 5%. Systematic uncerquality cuts and potential contributions from (including those from weak decays) are examr tightening the track selections on dz/σ (dz) 2 to 5. The associated yields are found to be track selections within 2%.
2-D two-particle correlation functions for and high (b) multiplicity, for pairs of charged p T < 3 GeV/c. For the low-multiplicity seleche dominant features are the correlation peak , 0) for pairs of particles originating from the ngated structure at φ ≈ π for pairs of partick jets. To better illustrate the full correlation k has been truncated. High-multiplicity events show the same-side jet peak and back-touctures. However, in addition, a pronounced e emerges at φ ≈ 0 extending to | η| of at observed structure is similar to that seen in collision data at √ s = 7 TeV [17] and in AA e range of energies [3][4][5][6][7][8][9][10]. , correlation functions were also generated for ECAL photons, which originate primarily from for pairs of ECAL photons. These distributions ures as those seen in Fig. 1 [17], are also plotted (open circles).
A clear evolution of the φ correlation function as a function of both p T and N offline trk is observed. For the lowest multiplicity selection in pp and pPb the correlation functions have a minimum at φ = 0 and a maximum at φ = π , reflecting the correlations from momentum conservation and the increasing contribution from back-to-back jet-like correlations at higher p T . Results from the hijing [24] model (version 1.383), shown as dashed lines, qualitatively reproduce the shape of the correlation function for low N offline trk . For multiplicities N offline trk 35, a second local maximum near | φ| ≈ 0 emerges in the pPb data, corresponding to the near-side, long-range ridge-like structure. In pp data, this second maximum is clearly visible only for N offline trk > 90. For both pp and pPb collisions, this near-side correlated yield is largest in the 1 < p T < 2 GeV/c range and increases with increasing multiplicity. While the evolution of the correlation function is qualitatively similar in pp and pPb data, the absolute near-side correlated yield is significantly larger in the pPb case.
In contrast to the data, the hijing calculations show a correlated yield of zero at φ = 0 for all multiplicity and p T selections. The η ∆  Figure 1: Two-dimensional (2D) per-trigger-particle associated yield of cha function of | | and | | for 3 < p trig T < 3.5 GeV/c and 1 < p assoc T < 1.5 centrality ranges of PbPb collisions at s NN = 2.76 TeV. The near-side peak two most peripheral distributions to better display the surrounding structu As was done in Ref. [1], to quantitatively examine the features of short-ran azimuthal correlations, one-dimensional (1D) correlation functions are c  [6,13,42] indicates the same underlying dynamics possibly drives these two phenomena. The first step is to understand the origin of high multiplicity events that populate the long tail of experimental multiplicity distributions. The probability distribution of multiplicity P (n) contains the information of n-particle correlations, a first principle estimation of which is a challenging problem. This requires full treatment of di↵erent sources of initial state fluctuations of the wave functions of the colliding systems and a framework of multi-particle production. In the CGC EFT, the saturation scale (Q 2 S ) is the natural scale controlling sub-nucleonic scale fluctuations which combined with the knowledge of nuclear geometry accounts for major sources of such initial-state fluctuations. Significant progress has been made in recent years to model these fluctuations using the framework of IP-Glasma model. Recently the importance of additional sources of initial state fluctuations have been realized, the origin of which is intrinsically non-perturbative and not captured in the conventional framework of CGC. Such fluctuations lead to a distribution of the intrinsic saturation scale of the hadron or nucleus as shown in Fig.2(left) and has to be introduced in the IP-Glasma model [44]. With the initial state including di↵erent sources of fluctuations constrained by the HERA data, the CGC EFT can be generalized to estimate n-particle production. For a given configuration of initial color charge, the n-particle distribution is a negative binomial distribution (NBD) with mean and width related to the saturation scale. Due to fluctuation of impact parameter, a convolution of many such NBDs gives rise to the final probability distribution of multiplicity. However it has been demonstrated that such distribution is narrower compared to the data. Only after including the intrinsic fluctuations of the proton saturation scale one can describe the tails of the experimental multiplicity distributions ( Fig.2(right)). Since multiplicity is dominated by low momentum (p T < 1GeV) gluons, it is very challenging to implement a scheme of fragmentation for the production of soft hadrons. Therefore in this approach the number of produced charged Two-particle correlation function in relative pseudo-rapidity and azimuthal angle showing long range ridge-like structure in high multiplicity p+p, p+Pb as compared to peripheral Pb+Pb collisions. Figures are taken from [6,13,47] order cumulants of v 2 {n}, n ≥ 4 and many others are now available (for a comprehensive review we refer the reader to Ref [51,52]).
By now a number of intriguing results observed in p+p and p/d/ 3 He + A have accumulated, including a strong multiplicity dependence of average transverse momentum p T that is attributed to radial flow as well as the observation of sizable Fourier harmonic coefficients v n (p T ) up to n = 4 and its higher order moments of the azimuthal correlation generally attributed to anisotropic flow. Most importantly, several characteristics, such as the mass dependence of both p T and v n (p T ) have been found to be similar to what is seen in A+A collisions.
However it is worth to mention that some striking contrasts also exist. Unlike in A+A collisions, where the observation of jet-quenching has been one of the pillars of the discovery of a strongly interacting Quark Gluon Plasma (QGP) [53,54], so far no evidence of (mini) jet-quenching has been found in small systems [55][56][57][58]. Even though the standard jet-quenching analysis in small-systems is complicated due to trigger bias effects, the absence of such phenomena may provide important insights with regard to the theoretical interpretation of the observed phenomena.

B. General theorectical perspectives
It is useful to first address the question about the origin of long-range azimuthal correlations (shown in Fig.1) from a more general point of view and formulate our theoretical expectations based on previous observations in small and large systems. While causality arguments imply that any long range rapidity correlation must originate from the very early stages of the collision [48], this leaves open the question how the observed momentum space correlations are created dynamically during the space-time evolution. Specifically one can, at least from a theoretical point of view, distinguish two different mechanisms whereby momentum space correlations of hadrons produced in the final state reflect i) intrinsic momentum space correlations of the partons produced in initial (semi-) hard scatterings and/or ii) position space correlations between initial state partons, e.g. the initial state geometry, which are transformed into momentum space correlations due to final state interactions.
While in any realistic scenario, both kinds of correlations i) and ii) contribute to the long-range azimuthal correlations, their relative strength depends on the magnitude of final state effects. In low-multiplicity p + p collisions for example, the dominant source of long-range azimuthal correlations is due to the production of back to back (mini-) jets. Since in this case the density of produced partons is low, the typical (semi-) hard partons produced in the initial scattering escape the interaction region without final state effects significantly affecting their back-to-back correlation. Considering on the other hand soft particle production amidst large parton densities in nucleus-nucleus collisions, it is well established that the azimuthal anisotropy of say p T 1 GeV particles is dominated by the final state response to the initial state geometry. In this case the mean-free path of a typical (semi-) hard parton is small compared to the system size, such that the initial state momentum correlations of ∼ GeV partons are destroyed during the equilibration process. Therefore, the subsequent dynamics of the equilibrated QGP can be accurately described by relativistic hydrodynamics.
Even though it is sometimes possible to choose the kinematics such that one mechanism dominates over the other, there are various examples in-between where both initial state and final state effects are important. One prominent example includes the behavior of jets in heavyion collisions. While highly energetic jets can escape the interaction region without equilibrating, they can loose a significant part of their energy through interactions with the softer medium. Even though the dominant correlation of the leading high-p T particles is still due to the initial back-to-back correlation, the path length dependence of the energy loss in the medium also leads to an additional correlation with the initial state geometry. Such correlations are reflected e.g. by the high-momentum v n (p T ) measuring correlations between soft and hard particles.
Clearly the aforementioned examples illustrate that it is important to consider both initial state momentum space correlations and the response to the initial state geometry due to final state effects in order to describe azimuthal correlations in small systems over a wide kinematic range. Our qualitative expectation is illustrated in Fig. 2, where the azimuthal correlation strength due to initial state and final state effects is shown versus the event multiplicity e.g. in p + p collisions for a fixed transverse momentum range e.g. 1 − 3 GeV. Based on our discussion we expect that in low multiplicity or min-bias events the azimuthal correlations between 1−3 GeV particles are pre-dominantly due to back-to-back mini-jets (peaked at ∆φ = π). With increasing event-multiplicity the contribution from multi-parton processes, such as the "Glasma graphs" (Sec. II C 3), becomes increasingly important resulting in azimuthal correlations that have a symmetric structure in relative azimuthal angle ∆φ around π/2. When increasing the multiplicities even further, final state interactions in this transverse momentum region can no longer be neglected at some point and lead to a depletion of initial state correlations. Even though mini-jets do not fully equilibrate yet, the system starts to show a response to the initial state geometry, which in this low opacity region is presumably dominated by the path length dependence of the parton energy lossalso referred to as parton escape mechanism [59]. Ultimately, in the limit of very high multiplicities, mini-jets are fully quenched, resulting in the formation of a thermalized medium and the complete loss of initial state momentum space correlations. In this high opacity regime, azimuthal correlations are dominated by the response to initial geometry described by a hydrodynamic expansion of a thermalized Quark-Gluon plasma.
One can attempt to further estimate the multiplicities corresponding to the transitions from the initial state to the final state dominated regime, exploiting recent theoretical progress in the understanding of the equilibration process [61]. Since the equilibration time at weak coupling corresponds to the time scale when a semi-hard parton ∼ Q s looses all its energy to form a soft thermal bath, one naturally expects the cross-over from the initial state to final state dominated regime to occur when the associated equilibration time τ eq becomes comparable to the system size R. Conversely, as long as τ eq R typical semi-hard partons escape without encountering significant final state interactions, whereas for τ eq R semihard partons are fully quenched, equilibrium is reached early on and the dynamics is dominated by the subsequent hydrodynamic expansion. Based on the estimate of Illustration of long-range azimuthal correlations in small systems, a slightly modified version of the figure from [60].

4/3
Teq (g 2 N c ) 1/3 10 for (η/s) Teq 5/4π at realistic coupling g 2 N c 10 [62,63] and the multiplicity dN/dy ξQ 2 s πR 2 with ξ 1/4 [64] we obtain that corresponding to a cross-over at around dN/dy ∼ 100, which in fact is much larger than the min-bias multiplicities reached in p + p or p + P b collisions [65]. We caution however that the estimate in Eq. (1) is inferred from leading order weak-coupling calculations and should only serve as a ballpark figure.
Beyond simple analytic estimates probably a promising alternative approach is to directly attempt an extraction of the boundaries between the different regimes through detailed comparisons of theory and experiment. While a first principle theoretical description is complicated throughout most of the multiplicity regimes shown in Fig. 2, significant theoretical progress has been made in understanding the features of initial state correlations in the regime where final state effects can be neglected. In the following we will review the theoretical computation of initial state correlations in the Color-Glass Condensate (CGC) effective field theory of high-energy QCD and critically access to what extent these calculations are compatible with the experimental observations.

A. High multiplicity events
Experimental observations suggest that long-range ridge like correlations in small colliding systems appear in high multiplicity events. Before we turn to a more detailed discussion of possible mechanisms to produce such correlations, a first necessary step is to understand the origin of high multiplicity events that populate the long tail of experimental multiplicity distributions. Considering the most elementary case of p+p collisions, high multiplicity events are a consequence of three major sources of fluctuations 1) geometry of collisions 2) intrinsic saturation scale of the proton 3) distribution of color charge density inside the protons.
While one naturally expects a strong impact parameter dependence of the multiplicity in p + p collisions [68], the importance of additional sources of initial state fluctuations have only been realized recently. Significant progress in including all of the above into a consistent phenomenological description has been made within the IP-Glasma model which is based on the framework of CGC [69,70]. Non-perturbative large x effects, which are not captured in the conventional CGC framework, are expected to give rise to fluctuations of the intrinsic saturation scale of the proton [71][72][73][74][75][76] as illustrated in Fig.3 (middle). In [67] intrinsic fluctuations of the proton saturation scale were introduced in the IP-Glasma model [67] according to a distribution with the variance σ 2 adjusted previously to the inclusive charged particle multiplicity and rapidity distributions in p+p collisions over a wide range of energies 0.2-7 TeV [67]. An interesting consequence of such fluctuations is that even in symmetric collision systems such as p + p, event-by-event fluctuations of the saturation scale of each proton lead to an asymmetry of the rapidity distribution of the produced particles on an event-by-event basis [77]. Consequently, high-multiplicity p+p collisions in particular are always asymmetric, and in a sense expected to look more like p+A collisions. The saturation scale and the geometric profile of the proton, which are the two most important ingredients in the IP-Glasma model are obtained from the IP-Sat parameterization of the HERA DIS data [78,79]. While in the original IP-Sat model the proton shape was assumed to be round, recent modification to IP-Sat has been done by assuming the gluon density to be distributed around three valence quarks inside the proton (see Fig.3 (left) ) [66,80] and several other models have started to include similar kinds of fluctuations [81,82]. Interestingly such sub-nucleon scale fluctuations can be constrained by incoherent diffractive vector meson production and improved agreement with existing HERA data can be achieved [66]. It has also been pointed out that geometric fluctuations of the proton can provide a dynamical explanation for the "hollowness effect" observed in elastic p+p scattering at LHC energies [83] and it would further be interesting to explore to what extent the modeling of the proton geometry can be aided by first-principle lattice QCD calculations.
With the initial state including different sources of fluctuations constrained by the HERA data, the n-particle production probability P (n) can be computed within the color glass condensate framework. For a given configuration of initial color charge, the n-particle distribution is a negative binomial distribution (NBD) with mean and width related to the saturation scale [85]. Due to fluctuation of impact parameter, a convolution of many such 4. A cartoon showing the contributions of di-jet and glasma graphs in two particle correlation function Y (∆φ) integrated over a broad range of |∆η|. This is a slightly modified version of the figure from [84] NBDs gives rise to the final probability distribution of multiplicity. However it has been demonstrated that such distribution is narrower compared to the data. Only after including the intrinsic fluctuations of the proton saturation scale one can describe the tails of the experimental multiplicity distributions (Fig.3(right)). Since multiplicity is dominated by low momentum (p T < 1GeV) gluons, it is very challenging to implement a scheme of fragmentation for the production of soft hadrons. Therefore the number of produced charged particles are generally taken to be proportional to the number of gluons. The results shown in Fig.3 indicate that the high multiplicity events that populate the tail of P (n) distributions are generated due to rare high color charge density configurations of the wave functions of the colliding systems. In the next section we argue that the same underlying dynamics that leads to the origin of high multiplicity events also drives the systematics of multi-particle correlations in small collision systems.
B. Qualitative discussion of Initial state correlations Multi-particle production in Quantum Chromo Dynamics (QCD) naturally leads to correlations between particles produced in high-energy collisions. A complete theoretical understanding of these effects though is extremely challenging. Nevertheless significant progress has been achieved in recent years based on the CGC effective field theory (EFT) of high-energy QCD, which provides the basis for phenomenological applications at RHIC and LHC energies.
Let us focus our discussion on the origin and systematics of the two particle correlations seen at LHC. By far the most well established source of long-range twoparticle azimuthal correlations is due to the production of back-to-back di-jets. Such processes (also referred to as "Mueller-Navelet" jets [86]) are depicted in the right panel of Fig. 4 and can be computed within standard perturbative QCD. Di-jet production is kinematically constrained to produce only away side (peaked at ∆φ = π) collimations and dominates in low-multiplicity or minbias events. However, in high-multiplicity events one is probing rare configurations of the proton where in addition to the production of di-jets from a single hard scattering, multi-parton processes become increasingly important. A first calculation of these effects in the CGC framework was based on evaluating the associated Feynman diagrams referred to as "Glasma graphs", depicted in the left panel of Fig. 4. Such graphs give rise to nonfactorizable two particle correlations that have a symmetric structure in relative azimuthal angle ∆φ around π/2 (see Fig.4). When decomposed in terms of the Fourier coefficients of the particle distributions, they give rise to non-zero even harmonics v n . Beyond the lowest order processes depicted in the left panel of Fig. 4, further contributions to the azimuthal collimations come from the multiple scattering of partons leading to both even and odd v n . Such processes can be included in a classical Yang-Mills description and will be discussed in more detail in a following section.
Since interference effects between Glasma graphs and Jet graphs vanish to lowest order in the kinematic regime Λ QCD Q s p T , q T the resulting two-particle correlations function as a direct sum of both contributions (3) The relative strength of the di-jet production represented by the "Jet-graph" and the "Glasma-graphs" determines the features of the observed di-hadron correlations as shown in Fig.4. In high multiplicity events the Glasma graphs are enhanced by a relative factor of α −4 S compared to the "Jet-graphs", one therefore naturally expects to see a pronounced near side collimation at ∆φ ∼ 0 that extends over a wide range of rapidity referred to as the "near side ridge". The fact that near side collimation extends far in rapidity is a consequence of the nearly boost invariant nature of the glasma gluon fields. While qualitatively these features are indeed present in the experimental data, of course it requires detailed theory calculations to establish the quality of agreement. In the remainder of this section, we will outline the essential steps in the computation of initial state correlations in the CGC framework. A summary of comparisons with experimental results is presented in section III.
In the CGC framework colliding protons and nuclei are effectively described as static sources of color charge on the light-cone that generate color currents The color charge densities ρ A(B) (x ± , x ⊥ ) in each colliding hadron or nucleus fluctuate from event to event and their statistical properties are constrained by independent measurements. Computation of multi-particle production in the CGC framework is based on the calculation of the classical Yang-Mills fields created from such color currents by solving the Yang-Mills equations Different theoretical descriptions within the CGC framework, employ different levels of approximation which are discussed in more detail in the following.

C. Perturbative computation
In the perturbative framework one tries to obtain an analytical solution of the gauge fields by performing order-by-order expansion of Eq.5 in powers of the color sources ρ A(B) . In the dilute-dense framework which assumes lowest order in ρ A and all orders in ρ B or in the dilute-dilute framework one can derive analytical expressions for n-gluon production in the k ⊥ -factorized form [87]. The essential ingredients to such factorization relations are the correlator of the dilute-sources in ρ A or the Wilson lines corresponding to the dense source ρ B in momentum space. Such correlators are represented as unintegrated gluon distributions (UGDs) Φ(x, k ⊥ ) and can be expressed in terms of T (x, r ⊥ )-the forward scattering amplitude of a quark-antiquark dipole of transverse size r ⊥ on a proton/nuclear target-through the expression The x-dependence of Φ is determined by the rapidity (Y = ln(1/x)) evolution of the dipole scattering amplitude implemented in the Balitsky-Kovchegov (BK) renormalization equation [88,89] which is a simplified form of the JIMWLK renormalization equations [88,[90][91][92][93]. The leading order expression of the BK equation is given by where r 2 ≡ r − r 1 and K is the BFKL kernel. The implementation of the kernel K often used for phenomenology includes a running coupling next-to-leading-log (NLL) correction to BK and is referred to as the rcBK equation [94] given by Eq.8 requires an initial condition for T A,B (x = x 0 , r ⊥ ), one choice of which is the McLerran-Venugopalan (MV) model [95,96] with a finite anomalous dimension γ given by s0 is a non-perturbative scale. This MV-like parameterization is constrained by global fits to the DIS data [97] although it must be noted that MV along with BK evolution does not include spatial geometric structure of the proton. The scattering amplitude is known to have a strong dependence on the impact parameter [78] which is incorporated in other parameterization such as IP-Sat [79,98] or b-CGC [78] models of DIS. More recently, significant progress in consistently including the full NLL corrections to BK evolution into DIS fits [99][100][101] has been made. However, so far this has not been included in phenomenological studies of p+p and p+A collisions.
With the un-integrated gluon distribution obtained from Eq.6, one can estimate the production of n-gluons using the k ⊥ factorization approach. In the following section we describe the approach for single-inclusive, doubleinclsuive and jet production essential for the phenomenology in p+p and p+A collisions.

Single inclusive gluon distribution
The single inclusive gluon production corresponds to the simplest process that describes the emission of a single gluon of momentum p T which in the perturbative framework can be written in the k ⊥ − factorized form as where C F = (N 2 c − 1)/2N c and S ⊥ is the transverse overlap area. It must be noted that this k ⊥ -factorized form of single gluon production has a logarithmic infrared divergence which is generally regulated by putting a lower p T,min cut. It must be noted that Eq.10 assumes that the dependence on transverse geometry or the impact parameter of collisions has already been integrated out and absorbed into S ⊥ , a more general expression in such a context can be found in Ref [87].

Di-Jets
In the perturbative framework, estimation of the twoparticle correlations in p + p and p + A collisions are based on a direct computation of di-jet and the glasma graphs in the k ⊥ -factorization approximation [28-30, 32, 41]. Such approximations are valid for momenta above the saturation scale Q s and do not include multiplescattering effects.
The di-jet contribution in this framework is estimated [86,102,103] to be where G is the BFKL Green's function that generates gluon emissions between the gluons that fragment into triggered hadrons. The form of G is given by [32]. Here ω(ν, n) = −2α s Re Ψ |n|+1 2 + iν − Ψ(1) is the BFKL eigenvalue with Ψ(z) = d ln Γ(z)/dz being the logarithmic derivative of the Gamma function with the effective coupling factor α s ≡ N c α S √ q a⊥ q b⊥ /π and φ ≡ arccos q a⊥ ·q b⊥ |q a⊥ | |q b⊥ | . For a pair of hadrons with a rapidity separation of ∆y 1/α s , G does the necessary resummation of the rapidity ordered multi-gluon emissions. The effect of such gluon emission leads to angular de-correlation that affects the observed di-hadron correlation [86]. The diagrams corresponding to such BFKL emission and its impact on broadening the awayside structure of the azimuthal di-hadron correlation is shown in Fig.5. In the α S ∆y → 0 limit one can obtain the well known form of the di-jet cross-section expression in the Multi-Regge kinematics (MRK) [103,104] : As we discuss in the following section, the quantitative estimation of the di-jet cross-section is essential for the description of the di-hadron correlation observed in high multiplicity events of p+p and p+Pb collisions.

Glasma graph
In the context of two particle correlations, there are a total of eight topologies of the glasma graph (full expression can be found in Ref [32]), the contribution for one such diagram to the two-particle correlations can be written as Beyond the straightforward, perturbative result the nonperturbative factor ζ has been introduced to account for the contributions to multi-particle production below the scale Q 2 S . Constrains from the independent analysis of nparticle multiplicity distributions [105][106][107][108] suggest that ζ be in the range 0.1 − 1 and phenomenological studies typically employ the value of ζ = 1/6 obtained from [106,107].
The collimation in the glasma graph framework comes from the fact that each of the four UGDs in Eq.14 has a bell-shaped structure with a maximum corresponding to the saturation momentum. Such nature of the UGDs kinematically constrains two gluons to be produced in similar (or back-to-back) directions giving rise to the double ridge structure in two particle correlations [28,109].

D. Dilute-dense models
Besides the perturbative diagrammatic approach taken in the Glasma graph calculation, several studies have been launched to investigate the origin of structure of the long-range azimuthal correlations in the limit where a dilute projectile of individual partons scatter off the color-fields of a dense projectile. Neglecting correlations of the incoming partons in the projectile, the double inclusive distribution of the scattered partons takes the form [38,40]  where W q/g,dy1 (b 1 , k 1 ) denotes the Wigner function of quarks/gluons inside the projectile and D(x, y) = 1 Nc Tr V (x)V † (y) is the usual dipole-correlator of lightlike Wilson lines in the fundamental/adjoint representation. While the transverse momenta k 1 and k 2 of the two incoming quarks are uncorrelated, it is evident from Eq. 15 that the correlations in the momentum transfers p 1 − k 1 and p 2 − k 2 give rise to azimuthal correlations of the scattered partons p 1 and p 2 . We note that in contrast to the perturbative calculation, present calculations in the dilute-dense description usually do not take into account rapidity evolution between the produced particles and are therefore limited to the kinematic range where y 1 − y 2 1/α s where evolution effects can be neglected. While recent progress has been made in formulating evolution equations for multi-particle production in dilutedense systems [110], these have not been employed for phenomenological studies so far.
By evaluating the dipole operator in short-distance expansion (Q s |x − y| << 1) where E i (x) = iV (x)∂ i V † (x) denotes the light-cone electric field, one can establish an intuitive picture of the origin of the long-range near side correlations [25]. When a projectile parton scatters off the color field of the nucleus it receives a transverse momentum kick in the direction of the color-electric field of the target. Color fields fluctuate from event to event and are locally organized in domains of size ∼ 1/Q s as illustrated in Fig. 6. When two (or more) quarks scatter off the same domain, they will receive a similar kick whenever they are in the same color state. Naturally, this leads to a correlation which is suppressed by 1/N 2 c (in the limit of large N c ) and by the number of domains Q 2 s S ⊥ , where S ⊥ denotes the transverse area probed by the projectile.
Several studies have computed the azimuthal harmonics of the two-particle correlation function in Eq. (15), either based directly on numerical evaluations including JIMWLK evolution [38,40] or within (semi-)analytic models [34,111,112]. Generally the double-inclusive spectrum in Eq. 15 features sizable azimuthal correlations v n which are on the order of 1/ N 2 c − 1 and most pronounced when both parton momenta p 1 are p 2 are on the order of the saturation scale. While the doubleinclusive spectrum for the two incoming quarks features both even (n = 2, 4, ...) and odd (n = 1, 3, 5, ...) harmonics, it turns out that the odd harmonics for gluons vanish identically due an exact symmetry of the the correlation function under p 1 → −p 1 . Even though the dilute-dense framework described above certainly does not provide the most realistic initial state model for mid-rapidity particle production in high-multiplicity events, it turns out that the absence of odd harmonics for gluons poses a more general problem within the CGC framework. However, as discussed in Sec. III recent simulations including the early-time classical Yang-Mills dynamics have been able to resolve this puzzle to some extent.

E. Classical Yang-Mills
Beyond the approaches outlined above there have also been new theoretical developments in the study of initial state correlations in event-by-event simulations in classical Yang-Mills theory [39,69,108]. In this approach Eq.5 is solved numerically for the individual colliding hadrons or nuclei in Lorentz gauge ∂ µ A µ = 0, where and then transformed into the light-cone gauge A + (A − ) = 0, where one finds [113][114][115] The infrared regulator m in Eq. (17) is of order Λ QCD and crudely incorporates color confinement at the nucleon level by damping the Coulomb tail. The gauge field in the forward light-cone after the collision at time τ = 0 is given by the solution of the CYM equations in Fock-Schwinger gauge A τ = (x + A − + x − A + )/τ = 0 in terms of the gauge fields of the colliding nuclei [116,117]: Such gluon fields produced after the collision are evolved in time according to Classical Yang-Mills equations up to time τ ∼ 1/Q s to estimate the gluon spectrum dN g /dyd 2 k t by imposing Coulomb gauge ∂ i A i τ = 0 and extracting the equal time correlation function [118,119] dN where g µν = (1, −1, −1, −τ −2 ) denotes the Bjorken metric and λ = 1, 2 labels the two transverse polarizations. In Coulomb gauge the mode functions take the form where H (2) α denote the Hankel functions of the second type and order α (see [118] for details).
The calculations in the framework of Classical Yang-Mills naturally include the Glasma graphs and extend the reach of perturbative calculations towards lower p T by consistently including multiple-scattering effects as well as coherent re-scattering in the final state. Since eventby-event simulations in the classical Yang-Mills theory also allow for an improved treatment of the impact parameter dependence, they can be used in the future to systematically study initial state effects across different collision geometries e.g. in p + A, d + A and 3 He + A collisions at RHIC [10,11,120].

F. Hadronization
So far we have outlined the computation of initial state correlations at the parton level. The mechanism of hadronization converts the partonic correlations driven by initial state dynamics into correlated production of final-state particles. Implementation of a realistic scheme of hadronization is therefore essential for the phenomenology of small systems collisions. A first principle QCD based approach to such a problem is very challenging. One therefore resorts to several available approximation schemes for fragmentations of partons. The most commonly used approach is the standard parton-hadron independent hadronization scheme [121] in which e.g. the single inclusive hadron distributions are obtained by convoluting the gluon distributions with fragmentation functions as where D g→h (z, µ 2 ) denotes the probability that a gluon fragments into a hadron carrying z fraction of its momentum at a scale µ 2 . The lower limit of the integral is determined from the kinematic requirement that the momentum fraction of the gluons x ≤ 1. Commonly used form for D g→h (z, µ 2 ) such as KKP or DSS fragmentation functions are obtained from fits to the inclusive hadron production data in e + +e − and p+p collisions [122,123]. In case of double inclusive production relevant for the study of two-particle correlations, one assumes an ansatz of the form [28] A limitations of is standard approach of hadronization is that the applicability of the fragmentations functions are questionable at low virtuality µ 2 . Therefore such scheme of hadronization can not be used for bulk particle production which is dominated by soft processes with typical virtuality scale µ 2 1 GeV. In this regime an alternative approach such as the Local Parton Hadron Duality (LHPD) [124] can be used to describe bulk particle production driven by initial sate dynamics [125].
The state-of-the-art scheme of hadronization to describe bulk particle production used in several event generators like PYTHIA [126,127] is based on the Monte-Carlo implementation of the Lund string fragmentation function where m T and z denote the transverse mass and the light cone momentum fraction of the fragmenting hadron. The default parameters a and b are constrained by global analysis of p+p data. The hadronization scheme in PYTHIA also includes decays of hadron resonances and final-state hadronic interactions. A quantitative study of the effect of hadronization on initial-state parton level correlations using different schemes of fragmentation has been done in Ref [128]. The study indicates that the structure of the final azimuthal correlations are very much sensitive to the choice of fragmentation. For example the resonance decay and hadronic interactions can distort the initial-state correlations and introduce artificial final-state correlations in the produced hadrons. Such mechanism will complicate the interpretation of experimental data. In such a context it has been demonstrated that hardronization effects in PYTHIA combined with the scheme of color reconnection produce effects that can mimic collective flow like pattern leading to a strong mass ordering of p T [129] in p+p collisions. In addition, the hadronic transport models that implements the Lund string fragmentation of PYTHIA have also shown to qualitatively reproduce the systematics of v 2 measured in p+Pb collisions [130]. Clearly more phenomenology in this direction will further improve our understanding. Meanwhile it is essential to implement a state-of-the art framework of hadronization to test different features of multi-particle correlations in the CGC  In the framework of gluon saturation models, a long-range correlation structure is predicted to arise from initial collimated gluon emissions [40][41][42]. The energy dependence of associated yields observed in the data is qualitatively in agreement with this model at p s = 13 TeV [39], as shown in Fig. 3 (b). However, although the model calculation quantitatively describes the associated yields over the multiplicity range covered by the previous 7 TeV data, significant deviations are observed at the higher multiplicities probed by the present 13 TeV data. The associated yields predicted by this model exhibit a much faster increase with N offline trk than that seen in the data, suggesting that other other mechanisms may be active in this region. Hydrodynamic models also predict no energy dependence: they reproduce the collective flow effect in heavy-ion collisions, which is nearly unchanged from the RHIC to the LHC center-of-mass energies, although they differ by more than an order of magnitude [43][44][45]. However, it remains to be seen whether hydrodynamic models can quantitatively describe the behavior of the observables presented here.
Long-range near-side yields have also been measured for pPb and PbPb collisions by CMS [14]. Figure 4 compares the associated yields in pp, pPb, and PbPb collisions for 1 < p T < 2 GeV/c as a function of the track multiplicity. The various data sets were collected at different centerof-mass energies, but this should have negligible effect on the results, as discussed above. In all three systems, the ridge-like correlations become significant at a multiplicity value of about 40, and exhibit a nearly linear increase for higher values. For a given track multiplicity, the

FIG. 7.
Figure taken from Ref. [13] showing the measurements by the CMS collaboration on the transverse momentum dependence of the near side ridge yield (left) and the multiplicity dependence of the same quantity shown for two energies in p+p collisions (right). Model calculations (shown by dashed lines) are from Ref. [41] using the Glasma graph and BFKL approach discussed in section II C.
at the level of hadrons; work in this direction is in progress [131].

A. p+p collisions
We begin with a comparison of model calculations with the latest data in p + p collisions at LHC energies [12], where the quantity of experimental interest is the rapidity and momentum integrated correlation function defined as which so far has only been computed in the perturbative CGC approach. The quantity of experimental interest is the near-side yield Y int defined as the ZYAM subtracted integrated associated yield per trigger given by Here ∆φ = ∆φ min. corresponds to the minimum of the di-hadron correlation function (ZYAM) and N trig is the number of trigger particles. A quantitative comparison of the near-side yield in p + p collisions between initial state calculations [41] using the Glasma graph + BFKL approach and the experimental data at 7 and 13 TeV from the CMS collaboration [13] is shown in Fig. 7. One can see a good agreement between the perturbative calculations and the data at low multiplicities. At higher multiplicities one sees systematic deviations where the calculations over predict the data. An improved calculations in the framework of Classical Yang-Mills which consistently include multiple-scattering effects [39] might lead to a better description of the data. The striking feature of the data model comparison shown in Fig.7 is that the near-side yield is approximately energy independent. Such energy independent scaling is naturally explained within the framework of CGC [41]. Origin of such scaling can be understood as follows. In the CGC framework a single scale Q 2 S determines both the single and double inclusive production in p+p collisions. The dependence on the center of mass energy in both multiplicity N rec ch and the near side yield Y int enters only through Q 2 S , i.e. N rec Fixing the multiplicity N rec ch (Q 2 S ) therefore naturally fixes the value of Y int (Q 2 S ) giving energy independent scaling of the ridge-like correlations. Such scaling provides strong indication of multiparticle production driven by a single semi-hard scale as captured in the CGC.

B. p+A collisions
A detailed quantitative estimation of the di-hadron correlations in p+Pb collisions at the LHC using the perturbative di-jet and glasma-graph framework has been performed in Ref [32]. One of the interesting results from Ref [32] shown in Fig.8 [39] compared to data from ATLAS [132] and CMS collaboration [133] for inclusive hadrons.
yield in p+Pb collisions as compared to p+p collisions for events with thesame multiplicity measured by the CMS collaboration. Results shown in Fig.8 (right) indicate that such framework also very well reproduces the peripheral subtracted correlation function (1/N trig ) dN/d∆φ measured by the ALICE collaboration in p+Pb collisions. Aforementioned due to the limitations of the perturbative framework such estimations can be performed only for trigger and associated momentum p T , q T ≥ 1 GeV.
Recent simulations in the framework of classical Yang-Mills that go beyond such limitations of the perturbative regime have already opened the path for improved phenomenology. Even though at present these calculations do not yet include di-jet graphs [134] and hadronization effects -and thus do not allow for a direct comparison with experimental data -simulations in this regard have lead to new insight into the correlations at the parton level concerning in particular the dynamics of the correlations during the very early stages. The results from the recent classical Yang-Mills simulations performed in p+Pb collisions [39] are shown in Fig.9. While gluons are produced with a significant momentum space anisotropy at τ = 0 + , initially the two-particle correlation function is symmetric around ∆φ = π/2 and features only even harmonics in accordance with the perturbative result [39]. However, including the effects of the classical Yang-Mills evolution up to τ = 0.4fm/c leads to the build up of a sizable v 3 on the parton level, while the initial state v 2 remains intact [39]. These results indicate that at least at the gluon level the non-perturbative dynamics of the glasma fields gives rise to large values of both v 2 and v 3 at time scales τ ∼ 1/Q s after the collision that are comparable to what is seen in the data. Interestingly sizable v 2 and v 3 in this framework also extends to large transverse momentum probably beyond the regime of the applicability of hydro in p+Pb collisions [135].
While a systematic comparison of the initial state calculations to experimental data in p + P b collisions yields a similar level of quantitative agreement of twoparticle correlations for momenta p T > 1 GeV [30,32], there has been enormous progress on the experimental side to identify additional signatures of collective motion which could be indicative of the onset of significant final state effects. Even though some observables, such as e.g. mass ordering properties observed in correlations between identified hadrons [136] are not necessarily sensitive to the origin of the correlations, more promising directions including e.g. correlations between more than two particles have also been explored [137,138]. One of the most striking observations in this regard is the sign change of the four-particle cumulant c 2 {4} observed around N ch (|η lab | < 1) ∼ 60 by ALICE [138]. However, most of the recent measurements have focused on low p T observables where the perturbative calculations [28-30, 32, 41] outlined in Sec. II C do not necessarily apply, hence complicating the comparison between theory and experiment.
Nevertheless, there have been first attempts to extend the theoretical framework to understand whether certain features of the low p T data such as the sign change of c 2 {4} can also be explained from initial state effects. Based on the dilute-dense approximation in Eq. (16) first attempts have been made to study initial-state correlations of more than two particles. Generalizing the dilutedense formalism to n-particle scattering, it was pointed out in [36,40,139] that the four particle cumulant c 2 {4} is sensitive to non-Gaussian correlations of the colorelectric fields inside the target. While the lowest order perturbative contribution to c 2 {4} was found to be positive, it was argued [36] that non-Gaussian correlations of the domains of color-electric fields may explain the experimentally observed sign change of c 2 {4} as a function of multiplicity. Even though a dynamical explanation for the existence of sizable non-Gaussian correlations is still lacking, this topic is an active subject of further investigations. However, one should also caution that in contrast to the two particle cumulant c 2 {2, |∆η > 2}, where short-range correlations are suppressed by introducing a large rapidity gap, the four-particle cumulant c 2 {4} also receives contributions from short-range correlations and interference diagrams. While a complete theoretical calculation of higher order n-particle correlations is desirable, further progress is needed to tackle this problem.
Similarly, preliminary results from event-by-event simulations in classical Yang-Mills theory a la [39] also suggest a negative sign of c 2 {4} when all particles have low momenta, while at high p T the four particle cumulant is always positive [140]. While the first results are interesting, further theoretical progress is needed to decide unambiguously whether the sign change in c 2 {4} can be explained in terms of initial state and early-time effects. It would also be interesting to extend the experimental measurements of higher cumulants towards higher momenta to achieve convergence on this issue.

IV. SUMMARY
Experimental observations of collective correlations in small systems challenge our current understanding of the space-time evolution of high-energy collisions. While at low multiplicities one expects dominance of initial state effects, it is also expected that at sufficiently high multiplicities initial state correlations are destroyed due to final state interactions and the hydrodynamic response to the initial state geometry provides the dominant source of correlations. However, it is theoretically challenging to predict the transition from initial state to final state dominated dynamics pointing to the importance of developing a unified theoretical framework where both effects are consistently taken into account.
In this brief review we specifically outlined recent developments in the approach based on initial state dynamics. While calculations based purely on initial state correlations are able to quantitatively describe various features of the experimental data in p + p and p+Pb collisions up to the highest multiplicity windows, a simultaneous description of the low and high p T observables over a wide range of multiplicity remains challenging within any single theoretical framework. Of course, several outstanding issues remain on the theory side and the development of a unified framework that combines 1) dynamics of multiple-soft interactions with 2) a first principle calculation of jet production and 3) a state-ofthe-art fragmentation scheme is essential for a wide range of phenomenological applications. So far experimental efforts have focused on soft observables such as e.g. the v n s, however complimentary information on the dynamics can be obtained by considering higher-momentum probes. While in heavy-ion collisions, the observation of strong jet-quenching phenomena provides an important indication for the formation of a strongly interacting Quark-Gluon plasma, no such features have been reported so far in small systems. Even though such an analysis is experimentally challenging, it would be important to establish to what extent mini-jets and actual jets are modified in high-multiplicity events to further constrain the relative importance of initial state and final state effects in small systems.