Angular correlations in pA collisions from CGC: multiplicity and mean transverse momentum dependence of $v_2$

Within the dense-dilute Color Glass Condensate approach, and using the Golec-Biernat-Wuesthoff model for the dipole scattering amplitude, we calculate $v_2^2$ as well as the correlations between $v_2^2$ and both the total multiplicity and the mean transverse momentum of produced particles. We find that the correlations are generally very small consistent with the observations. We note an interesting sharp change in the value of $v^2_2$ as well as of its correlations as a function of the width of the transverse momentum bin. This crossover is associated with the change from Bose enhancement dominance of the correlation for narrow bin to HBT dominated correlations for larger bin width.


I. INTRODUCTION
The observation at the Large Hadron Collider (LHC) of strong long range rapidity correlations with a characteristic structure in azimuth -"the ridge", in small size collision systems, pp and pA [1], is a very interesting result. The long range in rapidity implies, by causality arguments, that the correlations must originate in the initial stages of the collision, where the collectivity must emerge from the underlying microscopic dynamics. Two very different approaches have been pursued to provide a comprehensive explanation of the observed azimuthal correlations. The first one relies on strong final state interactions between the many produced particles, that lead to a relativistic hydrodynamic description of the evolution of the final state fireball [2]. This approach is quite successful in describing experimental data. It is well justified in nucleus-nucleus collisions where the ridge is also observed [3], but the theoretical rationale for its application to small systems (such as the ones produced in pp collisions) is still wanting. An alternative approach stresses the initial state correlations encoded in the light cone wave-functions of the incoming hadrons [4], and, in first approximation, ignores any final state effects. The present paper is devoted to further exploration of correlations within the latter framework.
There is a vast literature devoted to initial state correlations, see e.g. [5] for the most recent review and references therein. The main theoretical framework to compute multiparticle production at high energies is the Color Glass Condensate (CGC) effective theory, see e.g. [6,7]. Building on earlier calculations [8][9][10][11], most previous work on the subject has been focused on understanding angular Fourier harmonics v n (e.g. [12]). Two distinct quantum interference effects have been identified as giving rise to the emergence of even harmonics: the Bose enhancement of gluons in the incoming projectile wave function and the Hanbury Brown-Twiss (HBT) interference in the emission of gluons in the final state [15][16][17] 1 . The quantum interference effects between quarks and antiquarks were subsequently studied [18][19][20]25], and the formalism has been extended to calculate inclusive production of more than two particles in pp [21,22] and pA [23][24][25]. Further extensions allowed for the inclusion of odd harmonics via density [26] and non-eikonal [27,28] corrections, and considering anisotropies in the target produced via fluctuations [29].
In the famous pioneering publication of CMS on the ridge in pp collisions (first paper in [1]), the striking feature was the prominence of the ridge signal in rare high multiplicity events. The ridge signal was not observed in events approaching mean multiplicity. More recent analysis (utilizing modified background subtraction procedures) revealed azimuthal correlations also in minimal bias events, albeit with a somewhat weaker strength. It has also been observed that the momentum integrated v 2 , once it can be measured, shows very little dependence on multiplicity in a rather wide multiplicity range (see e.g. the ATLAS analysis in [1]). It is clearly important to explore further the multiplicity dependence of the ridge, both experimentally and theoretically.
Our work is devoted to this theoretical problem. Our calculations are performed for large values of transverse momenta, k 2 Q 2 s and keeping only leading contributions at large N c . We calculate the correlation between v 2 and the multiplicity within the CGC framework. We also study the correlation of v 2 with the mean transverse momentum of the particles produced in the collisions, motivated by experimental data [30] and by phenomenological works [31] where such correlations have been proposed as a tool to constrain the initial conditions for hydrodynamic evolution and to discriminate initial from final state correlations. Despite being motivated by experimental findings, we do not consider our approach as sufficiently rigorous from the phenomenological perspective and, hence, we are not making any attempt to compare with data. The main take home conclusion of our work is rather the qualitative features of the above correlations, which show an interesting characteristic dependence on the transverse momentum bin size used to calculate v 2 . We discuss these features in detail in the text as well as in the discussion section.
The manuscript is structured as follows. In Section II, we present the details of dense-dilute CGC approach as well as the model assumptions that we are using to compute the multiparticle production. Section III is devoted to the calculation of the second flow harmonic v 2 2 as well as correlations of v 2 2 with multiplicity, and the correlation of v 2 2 with the average transverse momentum of produced particles. Numerical results for these quantities are presented in Section IV. We conclude with a discussion in Section V. Technical details are provided in the Appendices.

II. MULTI GLUON PRODUCTION IN DILUTE-DENSE SCATTERING
A. The basic setup Our goal in this paper is to study correlations between the second flow harmonic coefficient v 2 2 and the particle multiplicity and average transverse momentum. We will calculate these correlations in the so called dense-dilute approximation. Below, we closely follow Ref. [23], in which double and triple inclusive gluon production in this approach were calculated in the CGC framework [6,7].
In this set up, the number of produced gluons for a given configuration of the projectile (proton) and the target (nucleus) is given by Here the Lipatov vertex L and its square Γ are defined as Besides, ρ p (p) is a given configuration of the color charged density in the projectile, and U (q) is the eikonal scattering matrix -the adjoint Wilson line -for scattering of a single gluon on a fixed configuration of target fields. The target Wilson lines depend on the target color sources, ρ t ; we suppress this in our notation for simplicity. The single inclusive and double inclusive gluon production in this approach are given by and where the averaging is performed over the projectile and target color charge configurations: and The normalization factors, Z p,t , are fixed so that The total multiplicity N (per unit of rapidity) is The mean transverse momentum squared per particlek 2 is calculated as Here and below we suppress the rapidity index. Note that formally both of these quantities may be divergent -the multiplicity diverges in the infrared if we treat the projectile as translationally invariant in the transverse plane, while the mean momentum diverges in the ultraviolet under fairly general assumptions about the structure of the projectile and target averages. We will regulate these divergencies by introducing physically motivated cutoffs (see below).
The azimuthal flow harmonics v 2 n are defined as where φ 1 , φ 2 are the azimuthal angles of the corresponding transverse momenta. Below, we focus on v 2 2 only. In this framework, each collision event corresponds to a fixed configuration of ρ p and ρ t . The averaging introduced in (6) and (7) is equivalent to averaging over all possible events.
An ideal way fo studying the dependence of v 2 on multiplicity would be to select from the total ensemble only events with the total multiplicity N i in some multiplicity bin (labeled by index i) and calculate v 2 by averaging only over those events. Mathematically one would have to introduce a projector P i (ρ p , ρ t ) on this subspace of configurations and then modify the averaging procedure as In practice, constructing P i is a complicated task that so far has not been accomplished (see, however, [32,33]). A simpler but nevertheless quite informative observable is a correlation between v 2 and N , i.e. v 2 (k 1 , k 2 )| ρp,ρt N | ρp,ρt p t , and similarly between v 2 and the transverse momentum per particle. The averaging in these expressions goes over the whole ensemble of events, and thus there is no need to consider particular sub ensembles. The calculation of these correlations requires us to calculate the three gluon inclusive production and Fortunately, three gluon inclusive production at mid rapidity has already been investigated in Ref. [23] and we will use many of the results of this paper. In order to be able to progress with the computations, we need to know the projectile and target averaging functionals. Here we are going to use simple models frequently used in a variety of CGC based calculations. For the projectile averaging we will use the simple Gaussian McLerran-Venugopalan (MV) model [34] specified by which corresponds to the weight functional Note that this weight functional defines a statistical ensemble that is translationally invariant in the transverse plane. This approximation is only reasonable for momenta of produced particles larger than the inverse radius of the projectile. Thus, in the following we will always assume k min > 1/R. As mentioned above, for the calculation of total multiplicity we will need to impose an infrared cutoff, which we will choose to be of the order of the transverse radius of the proton. We will further simplify our calculations by choosing µ to be independent of momentum. This last assumption means that the color charge density is completely uncorrelated in the transverse plane. Although this is clearly not the whole truth, one expects the density-density correlations to be only important at distance scales of the order of the confinement radius and above. Since we limit ourselves to consider momenta greater than the inverse radius of the proton, the presence of such correlations should not affect our results.
The averaging over the Wilson lines of the target will be performed in the approximation articulated recently in Ref. [25]. In this approximation any product of Wilson lines is factored into pairs according to the basic Wick contraction with the "adjoint dipole" amplitude defined as As explained in Ref. [25] this approximation is appropriate for the dense target regime. It collects all terms in the n-particle cross section which have the leading dependence on the area of the projectile. The approximation only includes terms which contain "small size" color singlets in the projectile propagating through the target. Any non singlet state that in the transverse plane is removed by more than 1/Q s ( where Q s is the saturation momentum of the target) from other propagating partons must have a vanishing S-matrix on the dense (black disk) target. On the other hand if the singlet state contains more than two partons, one looses a power of the area when integrating over the coordinates of the partons. Thus the leading contribution to the S-matrix in the (target) black disk limit is due to colorless projectile dipoles of the size smaller than the inverse target saturation momentum. The same approximation for the quadrupole amplitude has been used previously in Ref. [16] where its consistency with the explicit modelling of the Wilson line correlators via MV model has been verified. Note that this averaging procedure for the target is formally equivalent (disregarding subtleties related to the definition of the Haar measure) to the following form of the weight functional: Finally, we will adopt the Golec-Biernat -Wusthoff (GBW) model [35] for the dipole This model is known to be a good description for momentum transfer p of order of the saturation momentum and below. Although it does not properly account for the perturbative high momentum tail of the momentum transfer, we believe that it is quite adequate for the purposes of qualitative understanding of correlations, which is the goal of this paper.

B. From one to three gluons
In this subsection, we briefly summarize the results of [23] focusing on the contributions that will be essential for the computation of the observables that we are interested in.
Let us start with the single inclusive gluon production with transverse momentum k 1 : where q 1 is the momentum transfer from the target during the interaction and (k 1 − q 1 ) is the momentum of the gluon in the incoming projectile wave function. The expression for the double inclusive gluon spectrum can be organized as Here the first term is the square of the single inclusive spectrum and is the uncorrelated contribution to the double inclusive spectrum which is unimportant for our purposes. On the other hand, the second term in Eq. (22), is a quadrupole contribution which encodes the quantum interference effects. In the approximation of Eq. (17) its explicit expression reads where Our calculations in this paper will be always to leading nontrivial order in 1/N 2 c , and thus we do not specify the subleading term in Eq. (24).
The physical meaning of the two terms in Eq. (24) has been extensively discussed in [23]. The first term, I Q,1 , in the translationally invariant approximation contains two contributions. One is proportional to δ (2) Given that k i − q i is the momentum of the i-th gluon in the incoming wave function, this contribution clearly is due to the standard Bose enhancement of gluons in the incoming projectile state. The second contribution to I Q,1 is proportional to δ (2) (k 2 − q 2 + (k 3 − q 3 )). This contribution is due to the gluonic condensate in the projectile wave function, which is equal in magnitude to the Bose enhanced contribution. For simplicity we will refer to it as "backward" Bose enhancement, although one should keep in mind that the physics of this term is distinct from the physics of Bose enhancement. The second term, I Q,2 contains a delta-function of the final state momenta δ (2) (k 2 ±k 3 ). These correspond to the Hanbury-Brown Twiss correlations between the emitted gluons. The HBT correlations here exist for collinear as well as anti collinear momenta, since the gluons do not carry any physical charges. We will refer to these contributions as forward and backward HBT correlations respectively.
The triple inclusive gluon spectrum can be written as The first two terms in Eq. (27) correspond to the case where at least one of the gluons is uncorrelated with the other two. It is easy to see that in order to have a nontrivial correlation of v 2 with either multiplicity or average momentum, all three gluons have to be correlated. Therefore the first two terms in Eq. (27) do not contribute to the observables of interest. The only nontrivial contributions to the observables of interest arise from the last term in Eq. (27) which corresponds to fully correlated part of the triple inclusive gluon spectrum at leading N c . Its explicit expression reads The various expressions that enter Eq. (28) at leading order in 1/N 2 c are given and discussed below. The first term reads withĨ X,1 andĨ X,1 defined as Assuming translational invariance (15) one can clearly identify the various quantum interference effects contributing to three gluon correlations. Below we briefly review them following [23]: arises due to forward HBT correlation between the gluons 1 and 2, and forward Bose enhancement between the gluons 1 and 3.
• The second term in ) arises due to the backward HBT between gluons 1 and 2, and forward Bose correlation between the gluons 1 and 3.
• The third term in ) arises due to the forward HBT of the gluons 1 and 2, and backward Bose correlation between the gluons 1 and 3.
• The forth term in ) arises due to the backward HBT correlation between the gluons 1 and 2, and backward Bose correlation between the gluons 1 and 3. The terms I X,2 and I X,3 are obtained from I X,1 by permutation of the momenta of produced gluons The identification of various quantum interference effects in I X,2 and I X,3 follows straightforwardly from the earlier discussion using the very same permutation transformation. Note that the terms I X,2 contain an explicit contribution to forward/backward HBT of the gluons 2 and 3. These terms will be important later when we calculate v 2 .
Next is the explicit expressions for I X,4 : Again, all the physical effects behind each term can be identified: Bose correlation of gluons 1 and 1, forward Bose correlation of gluons 1 and 1 and forward Bose correlation of gluons 2 and 3.
Bose correlation of gluons 2 and 3, backward Bose correlation of gluons 1 and 3 and backward Bose correlation of gluons 2 and 1.
Bose correlation of gluons 1 and 2, forward Bose correlation of gluons 1 and 3 and backward Bose correlation of gluons 2 and 3.
• The fourth term in Finally, for I X,5 we have with the corresponding identification: -the forward HBT correlation of gluons 1 and 2, forward HBT of gluons 1 and 3, and forward HBT of gluons 2 and 3.
• The second term in -the forward HBT of gluons 2 and 3, backward HBT of gluons 1 and 3, backward HBT of gluons 1 and 2.
• Third term in -the backward HBT of gluons 1 and 2, forward HBT of gluons 1 and 3, backward HBT of gluons 2 and 3.
• The fourth term in -the forward HBT of gluons 1 and 2, backward HBT of gluons 1 and 3, backward HBT of gluons 2 and 3.

III. THE v2 AND THE CORRELATIONS
A. Total multiplicity, mean transverse momentum and v2 We start with calculating the total multiplicity and the second flow harmonic coefficient v 2 . Starting from the expression for the single inclusive spectrum (21), and carrying out the q 1 integral we obtain where S ⊥ is the transverse area of the projectile and Q 2 s is the saturation momentum of the target as defined in Eq. (20). Here we have introduced the infrared cutoff λ 1 by regulating the integration over the Schwinger parameter, see Appendix A. In terms of physical quantities the value of the IR cutoff is of order λ ∼ 1/(S ⊥ Q 2 s ). As is well known, the spectrum is divergent in the infrared, in the limit k 2 1 /Q 2 s → 0: For large momenta k 2 1 /Q 2 s → ∞, the spectrum reduces to the usual perturbative one: At finite λ the IR asymptotics of the expression in Eq.
Interestingly, the range of momenta in which this asymptotic behavior holds at very small λ is narrow, k 2 < Q 2 s /|2 + ln λ|. For very small values of λ the total multiplicity is dominated by momenta in the range Q 2 s /|2 + ln λ| < k 2 < Q 2 s where the spectrum is actually, flat dN (1) d 2 k1 ≈ µ 2 S ⊥ 2 + ln(λ) . This is an interesting feature of the spectrum, but it is only apparent at very small values of λ. In pA scattering for reasonable values of parameters (Q s ∼ 1 GeV, S ⊥ ∼ 1/Λ 2 QCD ) we use λ ∼ 1/S ⊥ Q 2 s ∼ 1/25. At this value of λ the interval of momenta where the spectrum is flat shrinks almost completely, and the spectrum exhibits a rather sharp crossover from a 1/k 2 behavior at k 2 < Q 2 s to 1/k 4 at k 2 > Q 2 s . In the next section we evaluate the k 2 1 integral of (36) numerically taking λ = 1/25. Since the dependence on λ is only logarithmic, the result is quite insensitive to the precise value of the IR cutoff.
The mean transverse momentumk 2 is formally defined as the average of k 2 1 with the weight Eq. (36). From the expression in Eq. (36), The integral over k 1 diverges logarithmically in the UV Being divergent, the average momentum defined this way is not a very useful quantity. Instead, whenever we need a quantity that represents a typical momentum of produced particles we will usē The second flow coefficient is evaluated using our expressions for the double inclusive gluon spectrum introduced in (24).
where in the large N c limit and in the approximation of translationally invariant projectile (see Appendix A) These expressions have been obtained assuming large values of transverse momenta k 2 One usually also averages the numerator and the denominator in Eq. (45) separately over momentum bins of finite width.
The only contribution to the numerator in Eq. (45) comes from the correlated term Eq. (24) since the uncorrelated term vanishes upon angular integration. The denominator, on the other hand, is dominated by the uncorrelated piece which at large N c is given by Eq. (23).
Although the general expressions for the two gluon inclusive spectrum have been known for a while [23,24], we are not aware of the actual calculation of v 2 in this simple dense-dilute approach. Here we evaluate Eq.(45) numerically, and present the results in the next section.

B. v2 vs total multiplicity
We now turn to our observables of interest. We first aim to study correlations between v 2 and multiplicity. The standard measure of correlation between two observables X and Y is the Pearson coefficient R which measures the strength of the correlation between X and Y relative to their autocorrelations. This type of observable was studied recently in [31] in order to flesh out the effects of initial state momentum anisotropies.
In our case however the calculation of the Pearson coefficient would involve the calculation of the four gluon inclusive production (i.e. (v 2 2 ) ) which is relatively complicated. In addition, we are not interested to compare the correlation with the autocorrelations of v 2 and N , but rather in comparing to to the average value of the observables themselves. We will therefore not calculate the Pearson coefficient, but rather define the normalized correlator as where the second equality follows since only the fully correlated part of three (two) gluon inclusive production contributes to the numerator (denominator). The numerator in this definition is precisely the same as the numerator in Eq. (46) with X = v 2 2 and Y = N , but it is normalized to the product X Y rather than to the square root of the product of variances of X and Y .
The correlation between v 2 and the total multiplicity of produced particles (per unit rapidity) is related to the inclusive three gluon production cross section (13). Starting from Eq. (28), and integrating over k 1 , the result can be split similarly as in Eq. (28): We are able to perform the k 1 integration analytically, while the remaining angular integrations are performed numerically.
Recall that we are only considering large transverse momenta of the observed particles, |k 2(3) | Q s . This large transverse momentum can be achieved in two distinct ways: either A) the incoming projectile gluons already have large transverse momentum and the momentum transfer in the scattering is relatively small, or B) most of the final state momentum is transferred to a projectile gluon in the scattering. The two contributions have very different behaviors. On the one hand, large transfer momentum is exponentially suppressed in the GBW model as exp{−k 2 /Q 2 s }, which favors contribution A. On the other hand, the number of gluons in the projectile wave function is strongly peaked at small momentum, so that N p (p)/N p (q) ∼ q 2 /p 2 . Thus the number of incoming gluons at high transverse momentum is suppressed roughly by a factor 1/(S ⊥ k 2 ). For very large transverse area this suppression may be significant enough so that contribution B can become comparable or even larger than contribution A. However for a proton projectile this factor is very unlikely to compete with the exponential suppression due to high momentum transfer. In our calculations, therefore, we only keep the contribution due to small (∼ O(Q s )) momentum transfer from the target.
The calculation is fairly lengthy and the details are given in the Appendix A. The results are presented below. In general we find two types of terms. The one type gives a correlation which in momentum space has width of order Q s . This arises from Bose correlations between the incoming gluons 2 and 3 in conjunction with either HBT or Bose correlations of any one of these gluons with gluon 1. These terms are: • X 1 : As explained in the previous section, X 1 contributes largely to the forward correlation of the produced gluons. Indeed, as seen from its final expression, X 1 is enhanced in the forward region k 2 = k 3 , where the exponential pre factor is equal to unity, The width of the forward region is clearly |k 2 − k 3 | ∼ Q s , and away from this region this expression is exponentially suppressed.
• X 3 : Comparing Eqs. (75) and (116), one notes that X 3 = X 1 (k 2 ↔ k 3 ): • X 4 : We again notice that X 4 is enhanced In the limit k 2 = k 3 : The second type of terms is due to HBT correlations between gluons 2 and 3. These correlations in the translationally invariant approximation lead to δ functional terms, contributing when k 2 = ±k 3 . Accounting for a finite projectile area would regulate the delta functions smearing them on the scale of order 1/S ⊥ . Nevertheless, the correlation due to these terms is very narrow. We will come back to this point in the next section when analyzing our numerical results. The terms of this type are: • X 2 : • X 5 : In the next section we present the results of the numerical evaluation of the angular integral of these expressions as defined in Eq. (47).

C. v2 vs mean transverse momentum
The second observable we consider is the correlation between mean transverse momentum and v 2 defined as In accordance to our discussion earlier, we have substituted Q 2 s for the average transverse momentum in the "normalization" in the denominator.
The computation of this observable proceeds very similarly to the one considered in the previous subsection. Details are given in Appendix B. Here we present the results: with (64) In the next section we present the results of the numerical evaluation.

IV. NUMERICAL RESULTS
We now turn to numerical evaluation of the correlators discussed above. Here we mainly present the results, keeping their discussion for the next Section.
Note that in all the figures we plot momentum in units of Q s and the quantities of interest multiplied by the factor (N 2 c − 1)S ⊥ Q 2 s in order to exhibit the universal features of the result applicable to any target (any value of Q s ) and projectile (any value of S ⊥ ). The ratios we calculate also do not depend on the projectile scale µ 2 . To extract a number relevant for p-Pb or p-Au scattering one should take the realistic value (N 2 c − 1)S ⊥ Q 2 s ∼ 200. For the normalization in Eqs. (47) and (56), the value of the cutoff λ has to be specified in the integration Eq. (36). While λ = 1/25 was selected, we have checked that varying λ in reasonable limits does not appreciably change the results.
We start with calculating v 2 , Eq. (45). In addition to the angular integration we also integrate the absolute values of transverse momenta within finite width bins. Thus we calculate We take k ∆, k ∆ and ∆ ∼ Q s . We find it interesting to explore the interplay between the relative position of the centers of the two bins, k and k and the width of a bin ∆. As discussed above, v 2 2 receives contributions form two types of correlations: the Bose and the HBT correlations. While the width of the Bose correlation in momentum space is naturally of order Q s , the HBT correlations have much shorter range (in our expressions they are formally represented by a delta function). Thus we expect that when |k − k | < ∆ both, the HBT and Bose effects will contribute to v 2 2 , however when there is no overlap between the two bins, the HBT correlation should disappear. We thus expect a characteristic dependence of v 2 2 on ∆ (at fixed k − k ) such that v 2 2 should vary steeply when k − k ≈ ∆. Fig. 1 shows our results for v 2 2 . In the left panel we see that the dependence of v 2 2 on the transverse momentum is rather different for overlapping and non overlapping momentum bins. In the right panel we observe, as expected, a sharp change in v 2 2 at the point when the width of the interval equals the distance between the interval midpoints. Interestingly we learn from Fig. 1 that the contribution of the HBT correlations to v 2 2 is overwhelmingly large: it is by about a factor of ∼ 50 dominates over the contribution of Bose enhancement (right panel of Fig. 1).
Next up is the correlation of v 2 2 with multiplicity, Eq. (47). Again we integrate over bins of width ∆ for the two momenta, Our numerical results for the correlation function between v 2 2 and the total multiplicity are presented in Fig. 2. We first take coinciding bins, that is k = k and the bin width ∆ = Q s /2. In this kinematics v 2 2 is dominated by HBT. The result is the solid (blue) curve in Fig. 2. The dashed curve in Fig. 2 displays the situation when the momenta are offset by Q s , that is k = k + Q s . This choice eliminates the HBT contribution to the azimuthal anisotropy v 2 2 . Fig. 2 shows that the normalized correlation function is strongly suppressed for values of bin width for which v 2 2 is sizable, which is when the HBT effect in v 2 2 is dominant. The same effect is also demonstrated in Fig. 2, where we show the correlation function as a function of the bin width ∆. For illustration, we chose the centers of the bins at k = 4.5Q s and k = 5Q s . When ∆/Q s is small, the bins are not overlapping and no HBT contribution is present in v 2 2 . At these values of bin width the correlation between v 2 2 and multiplicity is sizable. However for ∆ > 1 2 Q s = |k − k | there is a steep decrease of the correlation and it very sharply drops to negligible values.
We observe a similar behavior for the correlation of v 2 2 with transverse momentum. Fig. 3 shows this correlation as a function of transverse momentum and the same quantity as a function of the bin width.
Finally, Fig. 4 shows the ratio R ≡ O k,v2 /O N,v2 as a function of transverse momentum. The correlation with transverse momentum clearly drops with k slower than the correlation with multiplicity.

V. DISCUSSION
Our results are quite curious. First, regarding v 2 2 we find a very characteristic sharp transition in the value of v 2 2 as a function of the bin width. Referring to Fig. 1 the value of v 2 rises sharply from v 2 ≈ 2 × 10 −3 to v 2 ≈ 1.5 × 10 −2 , i.e. almost by an order of magnitude as the bin width is increased from ∆ < |k − k | to ∆ > |k − k |. This transition is entirely due to the fact that at this value of bin width the HBT correlations start contributing to the second flow harmonic, since the HBT peak is much narrower than the Bose correlation. Although the presence of the transition is expected on these grounds, the fact that the value of v 2 rises by such a large amount is worth noting. We conclude that the contribution of the HBT correlations to v 2 completely overwhelms the contribution of Bose enhancement. The actual numerical value for v 2 that we get is quite reasonable. For bins centered at k = 4.5Q s ≈ 4.5 GeV and k = 5Q s ≈ 5 GeV and ∆ = Q s (right panel in Fig. 1) we find v 2 ≈ 1.5 × 10 −2 . This is to be compared to typical values of 0.05 − 0.08 for transverse momentum integrated v 2 in p-Pb collisions at LHC [1]. Since the integrated v 2 is dominated by the lower momenta, this discrepancy of a factor of 5 or so may be attributable to the relatively high value of transverse momentum in our calculation. The trend of v 2 rising towards low momenta is clearly seen in the left panel of Fig. 1.
We do reiterate though, that our calculation is not meant as a phenomenological fit in any way, but rather as a qualitative study of the effects of quantum statistics on the correlations. As such, we note that the sharp rise in v 2 with bin width is a very characteristic behavior, and it would be very interesting to explore such dependence experimentally.
We wish to comment on two peculiar features seen on the left panel in Fig. 1. First, it is interesting to note that in the regime dominated by Bose correlations (dashed curve), v 2 2 is only very weakly dependent on the transverse momentum. Although it does slowly decrease towards large momenta, this decrease cannot be discerned for the range of momenta on the plot. This suggests that the Bose correlated part of the two particle production scales with the same power of momentum as the square of the single particle spectrum.
Another property to note is that the ratio of the HBT to Bose contributions as seen in the left panel of Fig. 1 seems to be even greater than on the right panel. The ratio between the solid and dashed curves at k ∼ 4 − 5Q s on the left panel is around ∼ 500, rather than the factor ∼ 50 that we have inferred from the right panel. Admittedly, it is not a completely fair comparison, as the solid line on the left panel corresponds to the two momenta sampled from the same bin, while on the right panel the centers of the two bins are displaced by 0.5Q s . Still it is a little surprising that a relatively small displacement of the bin centers leads to such a dramatic effect. The reason for this is our treatment of the projectile as translationally invariant, which leads to a delta-function HBT correlation. In this situation the HBT contribution is essentially given by the overlap area of two rings corresponding to the two momentum bins. Displacing the centers of the bins relative to each other even by a small amount leads to a significant change of the overlap area, and thus the HBT contribution has a sharp peak at zero displacement. As we mentioned before, in a more realistic treatment which takes into account finite transverse size of the proton, the HBT peak should be smeared to have a width ∼ 0.2 GeV (the inverse radius of the projectile) which should significantly soften the dependence on |k − k |. We have checked numerically that smearing the HBT peak does indeed have such an effect. For that reason we limit our consideration to kinematic situations where the distance between the bin centers is greater than 0.3 − 0.4Q s .
Moving on to the correlation of v 2 2 with multiplicity as well as with the transverse momentum, we again observe a very characteristic dependence on the bin width. For small bin width where the HBT does not contribute, the correlation is sizable, ∼ 3 × 10 −3 . However for larger bin width this correlation drops by a factor of about 30 to 50, and is negligible. Note that the transition in O N,v2 and O k,v2 is the opposite to that in v 2 : whereas v 2 is smaller at small ∆, the correlations O N,v2 and O k,v2 are larger, and vice versa. This tells us that although the contribution of HBT to v 2 is much larger than that of the Bose enhancement, the HBT is much weaker correlated with total multiplicity (and transverse momentum) than the Bose enhancement is. In fact, since the magnitude of the drop in Fig. 2 is about the same as the magnitude of the rise in Fig. 1, we conclude that the numerator in Eq. (47) is a rather smooth function of ∆, and the drop in Fig. 2 is driven entirely by the sharp rise in the denominator in Eq. (47) (and the same is true for Eq. (56)).
The correlation between v 2 and N is a decreasing function of momentum. This is easy to understand, since the multiplicity is dominated by soft particles, while the correlation we track originates with gluons that already in the incoming wave function have large transverse momentum. We thus have no reason to expect a correlation at high values of k. Another aspect of this is seen in Fig. 4 which shows that the correlation of v 2 with k remains larger at high momentum than the correlation of v 2 with N , since particles with higher momentum contribute more significantly to average momentum than to the total multiplicity.
Qualitatively the smallness of the correlations between v 2 and N is consistent with the experimental data [1]. Experimentally the change in the integrated v 2 in p-Pb collisions is at most 40% over an order of magnitude change in N . Again, we note that our calculation is valid only for large k. The correlation clearly grows towards smaller values of k, but within our approximation we cannot push much below k ∼ 4Q s .
To conclude, we have calculated v 2 and correlations of v 2 with the total multiplicity and average transverse momentum per particle in the dense-dilute CGC approach using the GBW model for dipole amplitude. Our results are valid at large N c and large transverse momentum. We have not made an attempt to include any additional effects beyond multiple scattering in our calculation. We find a reasonable magnitude for v 2 and very small correlations with total multiplicity, consistent with data. An interesting observation we make is a characteristic very strong crossover in v 2 (k, k , ∆) as a function of the width of the transverse momentum bin ∆ at fixed k and k , associated with the dominance of the HBT contribution. It would be very interesting to explore such a dependence experimentally. Consider X 1 first. The k 1 integral is facilitated by Eq. (15): where X 1,a is the first term in Eq. (30), X 1,b is the second term in Eq. (30), X 1,c is the first term in Eq. (31) and X 1,d is the second term in Eq. (31), all integrated over k 1 . The explicit expression of X 1,a reads We use the first δ-function to integrate over k 1 , the second δ-function to integrate over q 1 , and the third δ-function becomes δ(0) which, as usual, is regulated by the transverse area S ⊥ of the projectile. The result reads Note that the dipole d k 2 − (k 3 − q 3 ) is enhanced when its argument is small. Hence X 1,a is a contribution to the forward (k 2 k 3 ) correlation of gluons 2 and 3.
The rest of the terms in X 1 can be computed in a similar way: . (71) X 1,b is a contribution to the backward correlation of gluons 2 and 3. By renaming q 2 → −q 2 and q 3 → −q 3 , it is easy to see that X 1,b = X 1,a (k 3 → −k 3 ).
. (72) X 1,c is a contribution to the backward correlation of gluons 2 and 3. By renaming q 3 → −q 3 , it is easy to show that X 1,c = X 1,a (k 3 → −k 3 ).
. (73) X 1,d is a contribution to the forward correlation of gluons 2 and 3. X 1,d = X 1,a . Combining all the terms, we get In X 1 , the first term is a contribution to the forward correlation of the gluons 2 and 3. The second term with (k 3 → −k 3 ) is a contribution to the backward correlation of the gluons 2 and 3. Shifting and renaming the variables the expression can be brought to a more compact form: Using the GBW model (20) for the dipoles d, the X 1 contribution can be written as Note that X 1 ∝ e −(2k 2 2 +k 2 3 )/Q 2 s which provides an exponential suppression for produced gluon momenta much larger than the target saturation scale Q s . However we still have the q 2 and q 3 integrations to perform, and those may bring a compensating exponential enhancement. Next, we will perform these integrations concentrating on such possible enhancing factors. Let us first focus on the q 2 integration (the q 3 integration is done similarly and will follow). There are three types of Gaussian integrals to consider: The first integral is trivial, The second and the third integrals are performed with the help of a Schwinger parameter t introduced for the q 2 2 factor in the denominator: Changing variables to t = 1/(1 + Q 2 s t), the result reads For the last integration, I 2,4 , a similar procedure is followed, but with two Schwinger parameters t 1 and t 2 for each 1/q 2 2 factor. This makes it possible to perform the integration over q 2 and one of the Schwinger parameters. However, the remaining integral over the second Schwinger parameter is divergent, reflecting the original IR divergence of the q 2 integration. The origin of this IR divergence is pretty clear -it is an artefact of the approximation (15) with momentum independent µ. In fact, color neutrality should be imposed at some non-perturbative IR scale Λ min below which µ must vanish, µ(p < Λ min ) = 0. Roughly Λ 2 min ∼ 1/S ⊥ . We find it more convenient introducing a cutoff λ → 0 on the Schwinger parameter than the sharp IR cutoff Λ min to regulate this divergence. The relation between the two is The final result reads where Ei is an exponential integral special function. Combining the results given in Eqs. (78), (80) and (82), the overall result of the q 2 integration reads So far the integration over q 2 was computed without any approximation. Yet, as we have mentioned above, we are interested only in terms in X 1 that are not exponentially suppressed. In other words, only exponentially enhanced terms in (83) are of interest. The λ-dependent terms are not of that type and can be neglected. For k 2 2 Q 2 s , the result can be simplified using the large argument asymptotic expansion of the exponential integral function, Thus the exponentially enhanced contribution to the total q 2 integral reads We have kept sub-leading terms of the order Q 4 s /k 4 2 since, as we will see later, those will turn out to be the first non-vanishing terms when k 2 = k 3 . After the q 2 integration has been evaluated, X 1 still contains a q 3 integral, which is our next target: where The structure of these integrals over q 3 are very similar to the ones encountered in the q 2 integrations (77). Therefore, the integration over q 3 can be performed in a similar manner. Again, keeping only the exponentially enhanced terms, the results read: T 0,2 ≈ πQ 2 s e (k2+k3) 2 /2Q 2 s 2 2 (k 2 + k 3 ) 2 T 0,4 ≈ πQ 2 s e (k2+k3) 2 /2Q 2 s 2 4 (k 2 + k 3 ) 4 T ij 2,4 ≈ πQ 2 s e (k2+k3) 2 /2Q 2 Finally, using these results in Eq. (86), X 1 can be written as which can be organized differently and rewritten as in Eq. (49). This completes our analytical analysis of X 1 .

Computation of X2
After performing the trivial integration over k 1 , X 2 reads The first term in X 2 is an explicit contribution to the forward HBT of gluons 2 and 3. The second term, which is indeed a mirror image (k 3 → −k 3 ), is an explicit contribution to the backward HBT of gluons 2 and 3. By shifting the variables and renaming them, X 2 contribution can be put into a more compact and convenient form: × d(q 1 + q 2 ) d(q 2 + k 2 ) d(q 3 + k 2 ) q k After we use the GBW model given in Eq. (20) for the dipole operators, X 2 can be organized in the following way: Let us start from the integration over q 3 . This integral has been computed and the exact result is given in Eq.
(83). Since we are interested in computing the exponentially enhanced contributions from the q 3 integration, the approximate result is given in Eq. (85). Using this result in X 2 , we get Now, let us consider the integration over q 1 which can be organized as follows: By comparing Eqs. (75) and (116), it is straight forward to realize that X 1 = X 3 (k 2 ↔ k 3 ). Therefore, one can immediately write down the result for X 3 as in Eq. (51).