CGC/saturation approach for high energy soft interactions: `soft' Pomeron structure and $v_{n}$ in hadron and nucleus collisions from Bose-Einstein correlation

In the framework of our model of soft interactions at high energy based on CGC/saturation approach,we show that Bose-Einstein correlations of identical gluons lead to large values of $v_n$. We demonstrate how three dimensional scales of high energy interactions: hadron radius, typical size of the wave function in diffractive production of small masses (size of the constituent quark), and the saturation momentum, influence the values of BE correlations, and in particular, the values of $v_n$. Our calculation shows that the structure of the `dressed' Pomeron leads to values of $v_n$ which are close to experimental values for proton-proton scattering, 20\% smaller than the observed values for proton-lead collisions, and close to lead-lead collisions for 0-10\% centrality. Bearing this result in mind, we conclude that it is premature to consider, that the appearance of long range rapidity azimuthal correlations are due only to the hydrodynamical behaviour of the quark-gluon plasma.


I. INTRODUCTION
In our previous paper [1] we showed that Bose-Einstein correlations lead to strong azimuthal angle correlations, which do not depend on the difference in rapidity of the two produced hadrons (long range rapidity LRR correlations). The mechanism suggested by us, has a general origin, and thus manifests itself in hadron-hadron, hadron-nucleus and nucleus-nucleus interactions, and generates the correlation that has been observed experimentally [2][3][4][5][6][7][8][9][10][11][12]. The fact that Bose-Einstein correlations lead to strong LRR azimuthal angle correlations, was found long ago in the framework of Gribov Pomeron calculus [13], and it has been re-discovered recently in Refs. [14,15] in CGC/saturation approach [16]. In Ref. [1] it was noticed, that these correlations give rise to v n for odd and even n, while all other mechanisms in CGC/saturation approach, including the correlations observed in [14,15], generate only v n with even n.
The LRR correlations in CGC/saturation approach originate from the production of two parton showers (see Fig. 1). The double inclusive cross section is described by the Mueller diagram of Fig. 1-b, in which the production of gluons from the parton cascade, is described by the exchange of the BFKL Pomeron (wavy double line in Fig. 1-b), while, due to our poor theoretical knowledge of the confinement of quarks and gluons, the upper and lower blobs in Fig. 1-b require modeling. . FIG. 1: Production of two gluons with (y1, p T 1 ) and (y2, p T 2 ) in two parton showers ( Fig. 1-a). Fig. 1-b shows the double inclusive cross section in the Mueller diagram technique [19]. The wavy lines denote the BFKL Pomerons [20,21].
If the two produced gluons have the same quantum numbers, one can see that in addition to the Mueller diagram for different gluons (see Fig. 1-b), we need to take into account a second Mueller diagram of Fig. 2-b, in which two gluons with (y 1 , p T 2 ) and (y 2 , p T 1 ) are produced. When p T 1 → p T 2 , the two production processes become identical, leading to the cross section σ (two identical gluons) = 2σ (two different gluons), as one expects. When |p T 2 − p T 1 | 1/R, where R is the size of the emitter [17], the interference diagram becomes small and can be neglected.
The angular correlation emanates from the diagram of Fig. 2-b, in which the upper BFKL Pomerons carry momentum k − p T,12 with p T,12 = p T 1 − p T 2 , while the lower BFKL Pomerons have momenta k. The Mueller diagrams for the correlation between two gluons are shown in Fig. 2.
After integration over k T , the sum of diagrams Fig. 2-a and Fig. 2-b can be written as Eq. (1) coincides with the general formula for the Bose-Einstein correlations [17,18] d 2 σ dy 1 dy 2 d 2 p T 1 d 2 p T 2 (identical gluons) ∝ 1 + e irµQµ (2) where averaging . . . includes the integration over r µ = r 1,µ − r 2,µ . There is only one difference: Q µ = p 1.µ − p 2,µ degenerates to Q ≡ p T,12 , due to the fact that the production of two gluons from the two parton showers do not depend on rapidities. Note, that the contribution of Fig. 2-b does not depend on the rapidity difference y 1 − y 2 nor on y 1 and y 2 . For y 1 = y 2 Eq. (1) follows directly from the general Eq. (2), and the interference diagram of Fig. 3-b leads to Eq. (1), and allows us to calculate the typical correlation radius and the correlation function C (R|p T 2 − p T 1 |). On the other hand, for y 1 = y 2 but for p T 1 = p T,2 Eq. (1), gives a constant which does not depend on y 1 and y 2 . However, in general case y 1 = y 2 and p T 1 = p T,2 the diagram of Fig. 2-b looks problematic 1 , since it seems to describe the interference between two different final states. In appendix A we demonstrate that the contribution of Fig. 2-b does not vanish even in this general case. Note, that for y 1 = y 2 , the sum of two Mueller diagrams , indeed, relates to the interference between two diagrams, as is shown in Fig. 3-a and Fig. 3-b. For these kinematics, as we have mentioned For p T 1 = p T 2 , the sum of two Mueller diagrams can also be viewed as the interference of two diagrams of Fig. 3-c and Fig. 3-d, leading to The calculation of the Mueller diagram shows that this average does not depend on y 1 and y 2 .
Remembering that for two parton showers in each order of perturbative QCD, (or in other words at fixed multiplicity of the produced gluons) the amplitude can be written in the factorized form A = A L (r + , r − ) A T (r T ) leading to e irµ Qµ = e ir T ·Q T averaging over r T × e ir + Q − + ir − Q + averaging over r+,r− In our opinion, the above discussion shows that the Mueller diagram of Fig. 2-b, does not characterize the interference between two orthogonal states, but is an economical way to describe the independence of identical gluon production on rapidities, providing the smooth analytical description of the cross section from y 1 = y 2 to the general case y 1 = y 2 . Since this point is not obvious we would like to recall the main features of the leading log(1/x) approximation (LLA). In the LLA we account for the following kinematic region [20,21] for the production of two parton showers (see Fig. 4): The cross sections of double inclusive productions can be calculated in LLA for the production of two parton showers 1 We thank Alex Kovner for vigorous discussions on this subject. in the following way: where dΦ (1) n1+n2 and dΦ (2) n3+n4 are the phase spaces of the produced gluons in the first and second parton show- Fig. 4) and all other notations are shown in Fig. 4.
The transition from Eq. (7a) to Eq. (7b), occurs due to the fact that we want to obtain the log contribution ∝ (y i+1 − y i−1 ) for each dy i , These logarithms stem from the integration of the phase space of produced particles, while we can neglect the y i dependence of the production amplitude. In other words, the production amplitudes are functions only of the transverse momenta and Eq. (7b) shows that the longitudinal degrees of freedom can be factorize out [20,21].
Eq. (7a) after integrations over y i , can be re-written in a more efficacious form , viz.
integral over the longitudinal phase space Summing over n i we obtain the Mueller diagram of Fig. 1.
For identical particles we need to replace We wish to stress, that in Eq. (9) we use the Bose-Einstein symmetry for the production amplitudes, which are only functions of the transverse momenta of produced particles. Such a replacement leads to the sum of the diagrams of Fig. 2-a and Fig. 2 The goal of this paper is to calculate the function C (R|p T 2 − p T 1 |) which tends to 1 at p T 2 → p T 1 , and vanishes for R|p T 2 − p T 1 | 1. To estimate C (R|p T 2 − p T 1 |), it is sufficient to know the double inclusive cross section for y 1 = y 2 , where Fig. 2-b contributes significantly.
To obtain the double inclusive cross section, we need to add the cross section for two different gluon production which has the form FIG. 2: Production of two identical gluons with (y1, p T 1 ) and (y2, p T 2 ) in two parton showers. The diagrams in the Mueller diagram technique [19] are shown in Fig. 2-a and Fig. 2-b. The wavy lines denote the BFKL Pomerons [20,21].
In Eq. (10) we take into account, that we have N 2 c − 1 pairs of the identical gluons, where N c is the number of colours, and that the polarizations of the identical gluons should be the same. The latter leads to a suppression of 1 2 of the second term in Eq. (10). Using Eq. (10) we can find v n , since where ∆ϕ is the angle between p T 1 and Eq. (12)-1 and Eq. (12)-2 depict two methods of how the values of v n have been extracted from the experimentally denotes the momentum of the reference trigger. These two definitions are equivalent if V n∆ (p T 1 , p T 2 ) can be factorized as V n∆ (p T 1 , p T 2 ) = v n (p T 1 ) v n (p T 2 ). We will show below that in our approach this is the case for the restricted kinematic region R p T i 1. The first problem that we face in calculating C (R|p T 2 − p T 1 |), is to estimate the value of R, which increases with energy (see for example LHC data of Ref. [22]). On the other hand, the BFKL Pomeron [20,21] does not lead to the shrinkage of the diffraction peak, as it has no slope for the Pomeron trajectory. The only way to obtain a size which increases with energy, is to use the unitarity constraints , A BFKL (Y, b) ∝ e ∆BFKL Y a(b) < 1 [23], where ∆ BFKL is the intercept of the BFKL Pomeron and b is the impact factor. However, in QCD a (b) decreases as a power of b and the unitarity constraints lead to R ∝ exp (∆ BFKL Y ) [14]. Therefore, to obtain the energy behaviour of R, we need to introduce a non-perturbative correction at large b, which assures a(b) ∝ exp (−µ soft b), and also to take into account the multi Pomeron interactions which satisfy the unitarity constraints. Fortunately, the second part of the problem has been solved in the CGC/saturation approach [16], but the first needs modelling of the unknown confinement of quarks and gluons. Hence, we are doomed to build a model which includes everything that we know theoretically regarding the CGC/saturation approach, but in addition, one needs to introduce some phenomenological descriptions of the hadron structure, and the large b behaviour of the BFKL Pomeron.
Such a model for hadron-hadron interactions at high energy has been developed in Refs. [25][26][27][28], and it successfully describes the experimental data on total, inelastic and diffractive cross sections, as well as the inclusive production and LRR correlations. The goal of this paper is to show that the structure of the 'dressed' Pomeron in this model leads to strong BE correlations, and generates v n both for even and odd n, in hadron and nucleus interactions. In the next section we consider the contribution to C (R|p T 2 − p T 1 |) from the first Mueller diagram, and discuss the different sources of BE-correlations. In section 3 we give a brief review of the structure of the Pomeron in our model, in which we incorporate the solution to the CGC/saturation equations with additional non-perturbative assumptions: the large b behaviour for the saturation momentum, and the structure of hadrons. It has been known for a long time [13,29] in the framework of Gribov Pomeron calculus, and has been re-considered in CGC/saturation approach [30], that the LRR correlations stem from the production of gluon jets from two different parton showers (see Fig. 1). In section 4 we evaluate the BE correlations that result from the dressed Pomeron of our model, and show that they are able to describe the main features of the experimental data.

II. CALCULATION OF THE FIRST DIAGRAM
A. Proton-proton scattering The first Mueller diagram which contributes to C (R|p T 2 − p T 1 |) and which we need to calculate, is shown in Fig. 2-e and can be written in the form [31]: In Eq. (14) φ BFKL denotes the parton density of the BFKL Pomeron, with momentum transferred by the Pomeron k T or k T + p T,12 . The Lipatov vertex Γ µ , as well as the equations for φ BFKL will be discussed in the appendix A. Generally speaking, N I P h has a structure which is shown in Fig. 5: where M n denotes the mass of resonances, ∆ BFKL the intercept of the BFKL Pomeron, and G 3I P the triple Pomeron vertex. Considering the contribution of the first term to N I P h , we can neglect, in the first approximation, the dependence of φ BFKL on the momentum transferred, since Q T turns out to be of the order of the saturation momentum This is not the case for the second term in Eq. (15), which has Q T ∼ Q s . It leads to the BFKL Pomeron calculus which takes the Pomeron interactions into account. We will discuss this contribution in sections 3 and 4. In this section we restrict ourselves to the first term in the sum in Eq. (5). Collecting all formulae, we obtain that in the first diagram To obtain the first estimates for the vertices of the soft Pomeron interaction with the projectile and target, we use the following parameterizations: For proton-proton collisions we take B pr = B tr = B.
In this case [1,13] with B R = B pr B tr / (B pr + B tr ). B R = 1 2 B for proton-proton scattering. In Ref. [1] it is shown that Eq. (18) leads to V ∆n of Eq. (11) which is equal to where I n is the modified Bessel function of the first kind. In Fig. 6 taking B R = 5 GeV −2 , we plot the prediction for v n using Eq. (19) and Eq. (12). This value of B R , corresponds to the slope of the elastic cross section for proton-proton scattering at W = 13 GeV . One can see that Eq. (12)-1 and Eq. (12)-2 give different predictions, demonstrating that we do not have factorization with p min T = 0.5 GeV and p max T = 5 GeV , as done in Ref. [10]. To calculate such a v n , we need to know the dependence of the cross section on p T 2 . Indeed, we need to take Eq. (10) and integrate it over p T 2 : viz.
For Fig. 6-c, we need to know the behaviour of the double inclusive cross section on p T 1 and p T 2 . We assume that   Fig. 6-c describes the same vn as in Fig. 6 We took the energy dependence into account, by calculating the B R from the slope of the elastic scattering at given energy W , which was taken from Ref. [25]. which is taken in the interval 0.5 to 5 GeV, as is measured in Ref. [10]. Fig. 7-b exhibits the experimental data taken from Ref. [10].
One can see that the calculated values, as well as energy dependence are close to the experimental data of Ref. [10]. The main difference is in the p T dependence, which suggests the necessity to include the diffractive dissociation process or, in other words, the entire sum in Eq. (15), as well as the enhanced diagrams that are generated by the BFKL Pomeron calculus (see Fig. 5).
We can estimate the sum over resonances or, in other words, the diffraction production of states with low mass, by using our model (see appendix B for necessary formulae). In Fig. 8 we plot the correlation function C (p T,12 ) as defined in Eq. (10) for |p T 1 | = |p T 2 |, which is the result of these calculations. One can see that the effective p T,12 dependence of the slope, turns out to be much smaller than our estimates from the first diagrams that we obtained above. The slope that we used for the calculation shown in Fig. 6 was estimated as 1 4 B el , where B el = 20 GeV −2 is the slope of the elastic cross section at W = 13 T eV . We see two reasons for such a drastic change in the p T,12 dependence: first, we took into account the diffractive production processes which were neglected in Fig. 6; and second, in our model the effective shrinkage of the diffraction peak originates from the shadowing corrections, as the BFKL Pomeron has no inherent shrinkage. Such corrections are stronger in net-diagrams of Fig. 14-b that are responsible for elastic scattering, than for the fan diagrams of Fig. 24 that contribute to inclusive production. Recall that B shr ≈ 10 GeV −2 at W = 13 T eV , comes from the shrinkage of the diffraction peak. The calculation of v n are shown in Fig. 9. One can see that we obtain large v n for both odd and even n. The value of v 2 from Fig. 9-c is about 10% larger, than the experimental one from Ref. [10] (see Fig. 7-b). However, our calculations lead to narrower distributions in p T than the experimental one. The factorization is strongly violated, as in the case of estimates of the first diagram. is taken in the interval 0.5 to 5 GeV, as is measured in Ref. [10]. Fig. 7 illustrates the energy dependence of v n for proton-proton scattering, showing v 2 for two energies W = 2.56 T eV and W = 13 T eV . Note, that v 2 does not depend on energy, in accord with the experimental data of Ref. [10].
Therefore, we can conclude that the first term in Eq. (15) leads to a value of v n , which is large and of the order of the experimental one; the inclusion of diffraction in the region of small mass (sum over resonances in Eq. (15)) leads to a decrease of the interaction volume, but cannot reproduce the experimental p T distributions of v n , and BE correlations show the experimentally observed independence on energy.

B. Hadron-nucleus and nucleus nucleus interaction
For a nucleus we can simplify the calculation, considering cylindrical nuclei which have a form factor where J 1 is the Bessel function. Taking Eq. (21) into account one can see that We expect that S A (k T ) leads to small k T ∼ 1/R A , since the radius of nucleus is large. In Fig. 10-a we compare Eq. (22) with exp −Bp 2 T,12 which follows from Eq. (22), replacing S 2 A k 2 T by δ (k T ). The agreement is impressive. In Fig. 11 we plot the prediction for proton-gold scattering. One can see that the Bose-Einstein correlations generate large v n for n ≥ 3. Actually, we have several mechanisms (see, for example, review of Ref. [32]) for v n with even n, therefore, it is instructive to note that the simple estimates in this section lead to large v 2n−1 , larger than has been measured [11]. It should be stressed that using a more general approach which includes the diffractive production of small masses, as well as the shadowing corrections that lead to the shrinkage of diffractive peak, we obtain the predictions (see formulae in appendix C) which repeat the main features of our estimates in the simple model of Eq. (22). These calculations are plotted in Fig. 11-d -Fig. 11-f. In Fig. 10-b estimates for C (p T,12 ) in our model (see appendix C) are shown. One can see that C (p T,12 ) are different, and the model gives a smaller interaction volume. However, all qualitative features turn out to be the same: larger interaction volume than for proton-proton scattering, v 2 is much smaller than the experimental value (see Fig. 12); v 3 , v 4 and even v 5 are close to the experimental values; and the value of the typical p T is about 1 GeV instead of p T = 3 − 4 GeV in the experimental data.
For nucleus-nucleus interaction C AA (R|p T 2 − p T 1 |) takes the form Fig. 13 shows v n for gold-gold scattering. One can see three major differences: v n values turns out to be smaller than for proton-nucleus scattering, especially when p T 2 = p Ref T differs from p T 1 ; the momentum distribution is much narrower than for pA scattering, and v n are the same for all n. Comparing Fig. 6, Fig. 11 and Fig. 13 we can conclude that the simplest estimates lead to sufficiently large v n for both even and odd n, which are similar to those obtained in proton-proton and proton-nucleus collisions, but they are considerably smaller for the nucleus-nucleus case. Comparing these predictions with the experimental data of Refs. [2][3][4][5][6][7][8][9][10][11][12] we see that the BE correlations should be taken into account in all three reactions, since they give sizable contributions.  (19) and Eq. (12). Fig. 13-a shows vn that stem from Eq. (12)-1. In Fig. 13-b the estimates from Eq. (12)-2 for p Ref GeV is plotted. Fig. 13c and Fig. 13-d describe the difference between proton-gold and gold-gold interactions.

III. A BRIEF REVIEW OF OUR MODEL
In this section we will give a brief review of our model which has been developed in our papers [25,26]. The advantage of the model is that it describes the experimental data on diffractive and elastic production [25]; the inclusive production [27] and large rapidity range (LRR) correlations [28].
As has been mentioned we need to build a model which incorporates at least two non-perturbative phenomena: the correct large b behaviour of the amplitude and the hadron structure. These need to be incorporated so as to reproduce in the framework of one approach, the main features of the experimental data, such as the increase of the interaction radius with energy, a sufficiently large cross section of diffraction production, as well as energy and multiplicity dependence of inclusive cross sections and two particle correlations. On the other hand, we wish to include as much information as possible from a theoretical approach based on QCD.

A. Theoretical input and 'dressed' Pomeron Green function
At the moment, the effective theory for QCD at high energies exists in two different formulations: the CGC/saturation approach [33][34][35][36], and the BFKL Pomeron calculus [20,[37][38][39][40][41][42][43][44][45][46][47]. In building our model we rely on the BFKL Pomeron calculus, since the relation to diffractive physics is more evident in this approach. However, we are aware that the CGC/saturation approach gives a more general pattern [45,46]. In Ref. [46] it was proven that these two approaches are equivalent for where ∆ BFKL denotes the intercept of the BFKL Pomeron. As we will see, in our model ∆ BFKL ≈ 0.2 − 0.25 leading to Y max = 20 − 30, which covers all accessible energies. In addition in Ref. [46] it is shown that for such Y , we can safely use the Mueller-Patel-Salam-Iancu (MPSI) approach [48], which allows us to calculate the contribution to the resulting BFKL Pomeron Green function ( see Fig. 14-a): where A BA dipole-dipole is the dipole-dipole scattering amplitude in the Born approximation of perturbative QCD, and is shown in Fig. 14-a by the red circles.
We need to find the amplitude for the production of dipoles of size r i at impact parameters b i . This amplitude can be written as (see Fig. 14 The solution to the non-linear equation is of the following general form Comparing Eq. (26) with Eq. (27) we see Coefficients C n can be found from the solution to the Balitsky-Kovchegov equation [35] in the saturation region (see Ref. [47]).
with a = 0.65. Eq. (29) is a convenient parameterization of the numerical solution within accuracy better than 5%.
Having C n we can calculate the Green function of the dressed BFKL Pomeron using Eq. (25) and the property of the BFKL Pomeron exchange: Carrying out the integrations in Eq. (25), we obtain the Green function of the dressed Pomeron in the following form: where Γ (s, z) is the upper incomplete gamma function (see Ref. [60] formula 8.35) and T is the BFKL Pomeron in the vicinity of the saturation scale

B. Phenomenological assumptions and phenomenological parameters
The first phenomenological idea, is to fix the large impact parameter behaviour by assuming that the saturation momentum depends on b in the following way: with We have introduced a new phenomenological parameter m to describe the large b behaviour. The Y dependence as well as r 2 dependence, can be found from CGC/saturation approach [16], since φ 0 and λ can be calculated in the leading order of perturbative QCD. However, since the higher order corrections turn out to be large [49] we treat them as parameters to be fitted. m is non-perturbative parameter which determines the typical sizes of dipoles inside hadrons. As one can see from Table 1 from the fit m = 5.25 GeV, supporting our main assumption that we can apply the BFKL Pomeron calculus, based on perturbative QCD, to the soft interaction since m µ sof t where µ sof t is the scale of soft interaction, which is of the order of the mass of pion or Λ QCD .
Unfortunately, since the confinement problem is far from being solved, we have to assume a phenomenological approach for the structure of the colliding hadrons. We use a two channel model, which allows us to calculate the diffractive production in the region of small masses. In this model, we replace the rich structure of the diffractively produced states, by a single state with the wave function ψ D , a la Good-Walker [50]. The observed physical hadronic and diffractive states are written in the form Functions ψ 1 and ψ 2 form a complete set of orthogonal functions {ψ i } which diagonalize the interaction matrix T The unitarity constraints take the form where G in i,k denotes the contribution of all non diffractive inelastic processes, i.e. it is the summed probability for these final states to be produced in the scattering of a state i off a state k. In Eq. (37) √ s = W denotes the energy of the colliding hadrons, and b the impact parameter. A simple solution to Eq. (37) at high energies, has the eikonal form with an arbitrary opacity Ω ik , where the real part of the amplitude is much smaller than the imaginary part.
Eq. (39) implies that P S i,k = exp (−2 Ω i,k (s, b)), is the probability that the initial projectiles (i, k) reach the final state interaction unchanged, regardless of the initial state re-scatterings.

C. Small parameters from the fit and the scattering amplitude
The first approach is to use the eikonal approximation for Ω in which We propose a more general approach, which takes into account new small parameters, that come from the fit to the experimental data (see Table 1 and Fig. 14 for notations): The second equation in Eq. (41) leads to the fact that b in Eq. (40) is much smaller than b and b , therefore, Eq. (40) can be re-written in a simpler form Using the first small parameter of Eq. (41), we can see that the main contribution stems from the net diagrams shown in Fig. 14-b. The sum of these diagrams [25] leads to the following expression for Ω i,k (s, b) where T (r ⊥ , Y − Y 0 , b) is given by Eq. (32).
In all previous formulae, the value of the triple BFKL Pomeron vertex is known:  To simplify further discussion, we introduce the notation with a = 0.65 . Eq. (47) is an analytical approximation to the numerical solution for the BK equation [47].
. We recall that the BK equation sums the 'fan' diagrams.
For the elastic amplitude we have We will discuss the inclusive production as well as LRR correlations in appendix B.

IV. AZIMUTHAL ANGLE CORRELATION AND THE STRUCTURE OF THE 'DRESSED' POMERON
As has been discussed, our model includes three dimensional scales: m,m 1 and m 2 . m 1 and m 2 describe two typical sizes in the proton wave function, which could be associated with the distance between constituent quarks (size of proton) R p ∼ 1/m 1 and the size of the constituent quark R q ∼ 1/m 2 in the framework of the constituent quark model [55]. The third scale: m, characterizes the impact parameter behaviour of the saturation scale, and is intimately related to the structure of the dressed Pomeron in our model. In section 2 we discussed how two scales in the proton wave function arise in the BE correlations. Here, we would like to show that the third scale leads to the BE correlations which can explain the values of v n observed experimentally.
As we have discussed in section 3-A , the dressed Pomeron is the sum of enhanced diagrams (see Fig. 14-a) which is given by Eq. (31). Therefore, the exchange of the dressed Pomeron generates the production of an infinite number of the parton showers and, in particular, two parton showers which generate the BE correlations as is shown in Fig. 15. Integration over rapidities of triple Pomeron vertices [44] reduces the diagrams of Fig. 15-a and Fig. 15-b to the diagrams of Fig. 15-c and Fig. 15-d. We can calculate the probability to find two parton showers (P 2 ) inside of the dressed Pomeron expanding Eq. (31): 49) and the contribution of two parton showers production to the double inclusive cross section for the diagrams of Fig. 15-a, is equal to where a I P I P denotes the Mueller vertex of gluon emission (see Fig. 15). In our estimates for the calculation of v n , we do not need to know the probability P 2 , as well as the vertex a I P I P , assuming that a I P I P is the same in Fig. 15-a and in Fig. 15  all rapidities are in the laboratory frame. T (k T , y) is the Fourier image of T (b, y) defined in Eq. (32)-Eq. (34) and it takes the form  (50)).
For the interaction with nuclei, we need to take into account the interaction of the Pomeron with the nucleons inside the nucleus, as it is shown in Fig. 16. The equation for the resulting T A (y, k T ) takes the form (see Fig. 17-a) The triple Pomeron vertex Γ 3I P will be calculated in our model below.
The typical |k − k | ∼ 1/R A 1/m and, therefore, we can replace G A y , k − k by G A (y ) δ (2) k − k . Note that the normalization is such that the first diagram for G A = S A (b = 0) T (y, k T = 0), where S A (b) is defined in Eq. (C4). After integration over k T , Eq. (52) reduces to the following equation For G A we have the equation of Fig. 17-b, which has the following analytical form: The solution to these two equations (Eq. (53) and Eq. (54) can be written as follows whereΓ 3I P = Γ 3I P / (λγ) = P 2 .
T (y, k T ) has a physical meaning, of the BFKL amplitude in the vicinity of the saturation scale, where it has a geometric scaling behaviour [56], and it depends on one variable z = ln r 2 Q 2 s (Y ) . For diagrams of Fig. 15 the typical r ∼ 1/m i and z → λY . It is well known that the main contribution to the inclusive cross section stems from vicinity of the saturation scale, since this cross section is proportional to ∇ 2 r N (r, b; Y ), which tends to zero inside the saturation FIG. 16: The Mueller diagrams for the BE correlation for the 'dressed' Pomeron for proton-nucleus scattering. Black blob denotes the vertex for gluon emission aIP I P (see Eq. (41), the gray blob stands for the triple Pomeron vertex. domain (see Eq. (A9)). N is the scattering amplitude of the dipole with size r. The fact that we are dealing with the amplitude in the region where it has geometric scaling behaviour, is the reason why a non-linear equation of the BK type [35] is degenerate to one dimensional equations (see Eq. (53) -Eq. (55)). Using Eq. (50) we can calculate C (|p T 1 − p T 2 |) for proton-proton scattering, given by Eq. (1) which is equal to For proton-nucleus we have and for nucleus -nucleus C AA has the form: The results of calculations for C (R cor p T,12 ) using Eq. (56) -Eq. (58) are plotted in Fig. 18, One can see that the radius of correlations (R 2 cor = B) turns out to be very small in comparison with the same radius in Fig. 8 and Fig. 10. From C (R cor p T,12 ) we can calculate v n using Eq. (11) and Eq. (12)-1. However, C (R cor p T,12 ) shown in Fig. 18, are calculated for the production of gluon jets, while experimentally v n are measured for a hadron. Following Ref. [57] we explore the local parton-hadron duality(LPHD) suggested in Ref. [58]. In our approach the hadrons originate from the decay of a gluon jet, and their transverse momenta are where z is the fraction of energy of the jet, carried by the hadron. p intristic, T is the transverse momentum of the hadron in the mini-jet that has only longitudinal momentum. From Eq. (59) we obtain that the average p T of hadrons is equal to In Ref. [57] we found that we need to take z = 0.5 and p intristic, T = m π , to describe the inclusive spectra of hadron at the LHC. Using Eq. (60) we recalculate v n for a gluon jet to v n for hadrons, which are shown in Fig. 19. Comparing with the experimental data [2][3][4][5][6][7][8][9][10][11][12], and Fig. 7-b and Fig. 12, one can see that we describe the proton-proton scattering rather well, while for proton-nucleus we obtain v 2 which is smaller by 15-20%. For lead-lead collisions v 2 turns out to be two times smaller than the experimental value [12]. However, for central events with centrality 0 -10% measured v 2 is very close to our estimates. v n with n ≥ 3 are larger than the experimental values. In general the p T distribution is wider than the experimental one. The LPHD approach and Eq. (60) are very approximate, and we need to use a more advanced jet fragmentation function. Second, we need to add together the two mechanisms: one discussed in this section and one discussed in section 2. We need to include a more advanced fragmentation function, together with more careful accounting of the emission vertex in QCD (see appendix A). We will consider these in a future publication.
The estimates from our model show that the mechanism that has been discussed in section 2, yields about 10 -20% of the contribution which we now consider. In Fig. 20 one can see how the sum of two mechanism occur in v 2 . One can see that the sum has a wider p T distribution and a smaller maximal value. For proton-proton collisions both effects make predictions closer to the experimentally observed values of v 2 [10].
One of the properties that has been violated in the estimates in section 2, was the factorization r n = 1 where Fig. 21 shows that Eq. (61) holds at least for p T ≤ 4 GeV in accordance with the experimental data (see Ref. [11]).

V. CONCLUSIONS
In this paper we showed how three different dimensional scales in high energy scattering, arise in the Bose-Einstaein correlations that generates v n , for even and odd n. The first two scales are intimately related to the structure of the Fig. 19-a Fig. 19  vn versus pT at W = 13 T eV for proton-proton ( Fig. 19-a), proton-lead ( Fig. 19-b) and lead-lead (Fig. 19-c) scatterings, using Eq. (11) and Eq. (12)-1. wave function of the hadron, and have an interpretation in the constituent quark model, as the distance between the constituent quarks and the size of the quark. In a more formal way they characterize the size of the vertex of the BFKL Pomeron interaction with the hadron, and the typical size of the same vertex for the diffraction production, in the region of small mass. We demonstrated that these sizes lead to BE correlations which are large, but narrowly distributed in p T . The third size is the value of the saturation momentum in the CGC/saturation approach, and has been used in the construction of our model for the high energy soft interactions. This size is incorporated in the structure of the 'dressed' Pomeron in our model. It turns out that this size leads to values of v n which are close to the experimental values both for even and odd n, and they are broadly distributed in p T . In proton-proton scattering this mechanism is able to describe the experimental data both for even and odd v n , while for proton-nucleus and nucleus-nucleus collisions we obtain smaller values of v 2 : 20-30% smaller for proton-lead scattering, and two times smaller for leadlead collisions. However, we would like to stress that for centrality 0-10%, the structure of the Pomeron gives values of v n which are very close to the experimentally observed ones.
All estimates were made in the framework of our model for soft interactions which is based on CGC/saturation approach, but introduces non-perturbative parameters which describe the wave function of the hadron, and the large impact parameter behaviour of the saturation momentum. We describe in this model the total, elastic and diffractive cross sections as well as the inclusive production and long range rapidity correlations, and therefore, we trust that we can rely on the model when discussing the azimuthal angle correlations.
We demonstrated in this paper that BE correlations in the framework of CGC/saturation approach are able to explain a substantial part if not the entire, experimental values of v n for both even and odd n. Therefore, we believe that is premature to conclude that the origin of the observed long range rapidity correlations are only due to elliptic flow. where Eq. (A4) is the BFKL equation in the momentum representation, which has the following form in the coordinate representation [21,34]: where [16,51] 1  (14) we obtain that Eq. (A8) in the limit k T → 0, degenerates to the expression for the inclusive cross section which has the elegant form derived in Ref. [51] dσ The interesting feature of Eq. (A8) and Eq. (A9) is, that they remain correct, if we replace 2N BFKL by N G = 2 N − N 2 , where N is the solution of the Balitsky-Kovchegov equation [35]. Inside the saturation domain where N → 1, both equations lead to negligible contributions. In other words, in both equations the main contributions stem from the vicinity of the saturation scale, where x 2 12 Q 2 s ≈ 1. The solution for the scattering amplitude of two dipoles r 1 and r 2 to Eq. (A6) is known [21] where the functions φ (n) in (γ; r 2 ) are determined by the initial conditions at low energies and ω(γ, n) =ᾱ S χ(γ, n) =ᾱ S (2ψ (1) − ψ (γ + |n|/2) − ψ (1 − γ + |n|/2)) ; where ψ (γ) = d ln Γ (γ) /dγ and Γ (γ) is Euler gamma function. Functions E n,γ (ρ 1a , ρ 2a ) are given by the following equations.
In Eq. (A12) we use complex numbers to characterize the point on the plane where the indices 1 and 2 denote two transverse axes. Note that At large values of Y , the main contribution stems from the first term with n = 0. For this term, Eq. (A12) can be re-written in the form The integrals over R 1 and R 2 were taken in Refs. [21,59] and at n = 0 we have where F is hypergeometric function [60]. In Eq. (A16) w w * and b γ are equal Therefore, at large b, N BFKL decreases as a power of b which violates the Froissart theorem [14]. At present, as has been mentioned above, we cannot suggest a modification of the equation of the CGC/saturation approach in which the correct [23] exponential behaviour at large b would be incorporated. So we doomed to build a model. We discussed our model in section 3.

Born diagrams
The spirited discussions with our colleagues, showed us that it would be benificial to add a general discussion of the BFKL contribution, by calculating of the first Born diagrams for the production of two identical gluons that have rapidities y 1 and y 2 , and carry momenta p T 1 and p T 2 . These diagrams are shown in Fig. 23 for the scattering of the bound states of two oniums (two dipoles). Such a model for the scattering systems allows us to use the perturbative QCD approach, and has the analogy in the simplest bound system: deuteron.
The two onium bound state is described by the wave function Ψ (R 1 − R 2 ), where R i is the coordinate of the onium which is equal R i = 1 2 (x i + y i ) where x i and y i are coordinates of quark and antiquark in the onium (see Fig. 23). We introduce two new functions that describe the form factor of our bound state (G (q)), and the interaction of two gluons with the onium: The contribution of the diagram of Fig. 23 can be written as 2 where Γ µ (q T , p T 1 ) Γ µ (q T 1 , p T 2 ) is given by Eq. (A2) and Eq. (A3). One can see from Eq. (A19) and Eq. (3) that the typical q ≈ 1/r, where r is he size of the onium , while the typical values of k ∝ 1/R, where R is the size of the bound state. Assuming that R r, we see that k q. Anticipating p T,12 ∝ 1/R, we can reduce the contribution of the interference diagram to the following form: In Eq. (A21) we assume that p T 1 ≈ p T 2 and one can see that p T,12 from this equation is indeed of the order of 1/R, being much smaller than p T i if they are of the order of 1/r. For 1/R p T i 1/r we need to take Γ µ (q T , p T 1 ) Γ µ (q T 1 , p T 2 ) =  The Born interference diagram for production of two identical gluons with rapidities:y1 and y2 and transverse momenta p T 1 and p T 2 .Ri = 1 2 (xi + y i ), p T,12 = p T 1 − p T 2 . Red rectangle shows function Φ (k, p T 1 , p T 2 )(see text).
Appendix B: BE correlations in the model: diffractive production in the small mass region.

a. Inclusive production
The inclusive production in the framework of the CGC/saturation approach comprises two stages: the gluon minijet productions and the decay of this mini-jet into hadrons. For mini-jet production, we use the k T factorization formula, which has been proven in Ref. [51] in the framework of the CGC/saturation approach (see appendix for details).
where φ hi G denotes the probability to find a gluon that carries the fraction x i of energy with k T transverse momentum, andᾱ S = α S N c /π, with the number of colours equal to N c . 1 2 Y + y = ln(1/x 1 ) and 1 2 Y − y = ln(1/x 2 ). φ hi G is the solution of the Balitsky-Kovchegov(BK) [35] non-linear evolution equation, and can be viewed as the sum of 'fan' diagrams of the BFKL Pomeron interactions, shown in Fig. 24.  Fig. 24-a).For the sake of simplicity all other indices in φ (x1, pT − kT ) and φ (x2, kT ) are omitted. The wavy lines denote the BFKL Pomerons, while the helical lines illustrate the gluons. In Fig. 24-b the Mueller diagram for inclusive production is shown.
In our model the sum of 'fan' diagrams is given by Eq. (29). Assuming that the main contribution to dσ dy = d 2 p T dσ dy d 2 p T stems from p T ≤ Q s , we obtain the following formula: dσ dy = d 2 p T dσ dy d 2 p T = a I P I P ln (W/W 0 ) α 2 In (1) 1 2 Y + y + β 2 In (2) where whereG I P (y) and N BK have been defined in Eq. (46) and in Eq. (29), respectively. Regarding the factor in front of Eq. (B2) i.e. ln (W/W 0 ), where W = √ s is the energy of collision in c.m. frame, and W 0 is the value of energy from which we start our approach. One can see that Eq. (B1) is divergent in the region of small p T < Q s . Indeed, in this region φ's in Eq. (B1) do not depend on p T , since k T ≈ Q s > p T , and the integration over p T leads to ln Q 2 s /m 2 sof t , where m sof t is the non-perturbative scale, that includes the confinement of quarks and gluons (m sof t ∼ Λ QCD ).

b. LRR correlations
In our previous paper [28], we showed that in the framework of our model that has been described above, the main source of the long range rapidity correlation, is the correlation between two parton showers. In other words, it was shown that the contribution to the correlation function from enhanced and semi-enhanced diagrams, turns out to be negligibly small.
The appropriate Mueller diagrams are shown in Fig. 25. Examining this diagram, we see that the contribution to the double inclusive cross section, differs from the product of two single inclusive cross sections. This difference generates the rapidity correlation function, which is defined as R (y 1 , y 2 ) =