Azimuthal angle correlations at large rapidities: revisiting density variation mechanism

We discuss the angular correlation present in hadron–hadron collisions at large rapidity difference (α¯Sy12≫1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\alpha }_S\,y_{12}\gg \,1$$\end{document}). We find that in the CGC/saturation approach the largest contribution stems from the density variation mechanism. Our principal results are that the odd Fourier harmonics (v2n+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{2n+1}$$\end{document}) decrease substantially as a function of y12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{12}$$\end{document}, while the even harmonics (v2n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{2n}$$\end{document}) increase considerably with the growth of y12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{12}$$\end{document}.


Introduction
In this paper we address the problem of the azimuthal angle correlations of two hadrons with transverse momenta p T 1 and p T 2 and rapidities y 1 and y 2 , at large values of y 12 ≡ |y 1 − y 1 | 1/ᾱ S . Our main theoretical assumption is that these correlations stem from interactions in the initial state. We are aware that, unlike rapidity correlations which at large rapidities originate from the initial state interactions due to causality reasons [1], a substantial part of these correlations could be due to the interactions in the final states [2][3][4]. On the other hand, it has been demonstrated that at small rapidity differencē α S y 12 < 1 the interactions in the initial state [5][6][7][8][9][10][11][12][13][14] yield the value of the correlations, which describe the major part of the experimentally observed correlations .
In this paper we concentrate our efforts on calculating the long range rapidity part of angular correlations with large value of the rapidity difference y 12 . All previous calculations assumed thatᾱ S y 12 < 1 [5][6][7][8][9][10][11][12][13][14]. It turns out that in this kinematic region, the main source of the azimuthal angle correlations is the Bose-Einstein correlations of identical gluons, corresponding to the interference diagram in the production of two partonic showers. Intuitively, we expect that the correlations in the process, where two different gluons are produced from two different partonic showers, should not depend on the difference of rapidities (y 12 ), nor on the values of y 1 and y 2 . Using the AGK cutting rules [ 1 Mueller diagrams [39] for two partonic showers production. a The square of the production amplitudes, while b corresponds to the interference diagram which leads to the Bose-Einstein correlations. The wavy lines show the BFKL Pomeron [40][41][42][43], while the helical lines denote gluons. Figure 1c shows the example of a more complicated structure of the partonic cascades, than the exchange of the BFKL Pomeron. The color of the lines indicates the parton shower  Fig. 2 The generalization of Fig. 1b forᾱ S y 12 1. The wavy lines show the BFKL Pomeron [40][41][42][43], while the helical lines denote gluons. The color (blue and black) of the lines indicates the parton shower can prove that the two gluon correlations can be calculated using the Mueller diagrams [39] of Fig. 1.
The diagrams of Fig. 1 lead to correlations which do not depend on y 1 and y 2 , but only forᾱ S y 12 1. For large y 12 the contributions of Fig. 1 decrease. The main goal of this paper is to find the contributions which survive at large y 12 (ᾱ S y 12 1). At large y 12 , we have to take into account the emission gluons, with rapidities y 2 < y i < y 1 , which transform the Mueller diagram of Fig. 1b to the more general diagrams of Fig. 2. The general feature of Fig. 1b is that the lower Pomerons carry momenta Q T + p 12 and − Q T − p 12 with p 12 = p T 1 − p T 2 . Q T denotes the momentum along the BFKL Pomeron. After integration over Q T , we obtain p 12 ∼ 1/R h , where R h is the size of the target (projectile), which has a non-perturbative origin. Roughly speaking, the correlation function turns out to be proportional to G ( p 12 ), where G denotes the non-perturbative form factor of the target or projectile [12]. This conclusion stems from the value of the typical Q T for the BFKL Pomeron, which is determined by the size of the largest dipoles in the Pomeron. Figure 2 does not have these features. We will show that the azimuthal angle correlations originate from the integration over Q T (see Fig. 2), due to the structure of the vertices of emission of the Footnote 1 continued cuts of the triple BFKL Pomeron vertex with the produced gluon lead to different contributions. We do not consider such diagrams. gluons with p T 1 and p T 2 , which have contributions proportional to ( p T 1 · Q T ) n ( p T 2 · Q T ) n . Recall that these kinds of vertices are the only possibilities to obtain angular correlations in the classical Regge analysis [53]. This mechanism for azimuthal angular correlations was suggested in Ref. [54] (see also Refs. [10,[55][56][57]), and in the review of Ref. [57] it was called the density variation mechanism.
The paper is organized as follows. In the next section we discuss the contribution of the diagram of Fig. 2 in the momentum representation. In the remainder of the paper, we will use the mixed representation: the dipole sizes and momentum transferred (Q T ), which will be introduced in Sect. 3 and appendix A. Section 4 is devoted to the discussion of the single inclusive production in the Color Glass Condensate (CGC)/saturation approach. The double inclusive production is considered in Sect. 4, in which the rapidity dependence of the master diagram of Fig. 2 will be calculated. In Sect. 5, we estimate the angular correlation function and Fourier harmonics v n , and we present our prediction for the dependence of v n on the difference of rapidities (y 12 ). In Sect. 6 we draw our conclusions and outline problems for future investigation.

Correlations in the momentum representation
The double inclusive cross section of Fig. 2 takes the following form: as well as all other functions φ of this type, are the correlation functions, which at Q T = 0 give the probability to find a gluon with transverse momentum k T in the hadron (nucleus) of the projectile (target). φ k T , − k + Q T ; k T , − k T + Q T describes the interaction of two gluons with momenta k T and k T , which scatter at momentum transferred Q T . N Q T is a pure phenomenological form factor that describes the probability to find two Pomerons in the projectile or target, with transferred moment Q T and − Q T . C F = N 2 c − 1 /2N c where N c is the number of colors. The Lipatov vertex μ (k T , p T 1 ) has the following form: Using Eq. (2.2) we obtain We can simplify the master equation (see Eq. (2.1) by observing that the dependence on Q T and Q T is determined by the non-perturbative scale of the projectile (target) structure, which in Eq. (2.1) is absorbed in the phenomenological form factors N (Q T ) and N (Q T ). Therefore, the typical Q T and Q T turn out to be of the order of the soft scale μ soft , which is much smaller that the other typical momenta in Eq. (2.1), assuming that P T 1 and P T 2 are larger than μ soft . Introducing we can neglect Q T and Q T in the BFKL Pomeron Green functions and rewrite Eq. (2.1) in the form with Eq. (2.3), which takes the following form: At high energies the parton densities φ G H − l T , l T ; Y − y 2 in Eqs. (2.1) and (2.5) are proportional to exp ( BFKL Y − y 2 ) for the BFKL Pomeron, where BFKL = 2.8ᾱ S is the intercept of the BFKL Pomeron. Bearing this in mind, one can see that the interference diagram for the double inclusive cross section does not depend on y 1 , y 2 or on y 12 .
As Fig. 1a also does not depend on the rapidities y 1 and y 2 , and its expression has the following form: The most economical way of calculating the diagram of Fig. 2, is to use the mixed representation of the BFKL Pomeron Green function, G r , R, Q T , Y , where r and R are the sizes of two interacting dipoles, Q T denotes the momentum transferred by the Pomeron, and Y the rapidity between the two dipoles. This Green function is well known [42,43], and for the completeness of presentation we discuss it in Appendix A, referring to Refs. [42,43,[58][59][60] for all details.

Single inclusive production in a one parton shower
3.1 BFKL Pomeron: the simplest approach for a one parton shower The single inclusive cross section resulting from the one BFKL Pomeron is well known, and it is equal to The relation between the parton densities φ and the Green function of the BFKL Pomeron has been given in Ref. [46]: where N BFKL (r, r 1 ; Y ) is given by Eq. (A.1) or by Eq. (A.9), in the high energy limit. Equation (3.2) can be rewritten as follows: We have Substituting Eq. (3.2) and also Eq. (3.4) into Eq. (3.1) we obtain [46] where N pr and N tr denote the probability to find a dipole in the projectile and target, respectively. r 1 and r 2 are the typical dipoles sizes in the projectile and target. As can be seen from Eq. (2.1) we need to generalize Eq. (3.5) for the case Q T = 0. Equation (3.1) has to be replaced by Taking into account Eq. (3.2) for Q T = 0 and we rewrite Eq. (3.5) in the form (3.8)

General estimates
It should be stressed that the single inclusive production has the form of Eqs. (3.5) and (3.8) for the general structure of the single parton shower, as was shown in Ref. [46]. For example, for the process shown in Fig. 1c. We need only to substitute N tr (r, r 2 ; y, Q T ) is a solution to the non-linear evolution equation. For the case of inclusive production, we can considerably simplify the estimates noting that where Q s (y) denotes the saturation momentum.
In other words, the main contribution to inclusive production comes from the vicinity of the saturation scale, where r 2 Q 2 s ≈ 1. Fortunately, the behavior of N in this kinematic region is determined by the linear BFKL evolution equation [62][63][64][65][66][67] and has the following form [68]: where γ cr = 0.37. From Eq. (A.8) we see that, for Q T = 0, the scattering amplitude decreases at Q 4 T r 2 r 2 2 1. Therefore, we need to consider rather small values of Q T : Q 4 T r 2 r 2 2 ≤ 1. The product of vertices that determines the amplitude has two terms (see Eq. (A.5)) which are proportional to r 2 /r 2 2 iν and to Q 4 T r 2 r 2 2 I ν . Therefore, the maximum of ∇ 2 r N can be reached if r 2 /r 2 2 e κ y ∼ 1 and Q 4 T r 2 r 2 2 e κ y ∼ 1 and the amplitude then has the following form: The first term does not depend on Q T and, therefore, the upper limit of the integral over Q T , goes up to Q max T 2 ≈ 1/(r r 2 ). The second term, both for Q 2 T r r 2 < e − 1 2 κ y and for Q 2 T r r 2 > e − 1 2 κ y , turns out to be small. Indeed, in the first region the amplitude is small, while in the second region we are deep in the saturation domain where ∇ 2 r N → 0. Hence, we expect that in the integral over Q T , the first term gives a larger contribution than the second term, and we will only keep this contribution in our estimates.

The simplest diagram
In this section we calculate the simplest diagram of Fig. 2. We need to integrate the product of two BFKL Pomerons over Q T (see Eq. (2.5)): From Eq. (2.5) in the momentum representation, we see that r 1 = r 1 (r 2 = r 2 ), however, they are close to each other, being determined by the same momentum k T . We assume that Considering r 1 ≈ r 1 r 2 ≈ r 2 we will show that in the integral over Q T , the typical Q T ∼ 1/r 2 . In other words, the dependence of Q T is determined by the largest of interacting dipoles.
Assuming that both ν 1 and ν 2 are small, we see that all four terms are equal to each other, and the integral can be written as follows: The appearance of the pole ν 1 = −ν 2 indicates that the contribution from this kinematic region is large.
Closing the contour of integration on ν 2 over the pole, we obtain (4.4) Actually, the double inclusive cross section depends on ∇ 2 N as we argued in the previous section. Repeating the procedure for we obtain for small ν 1 and ν 2 Using the method of steepest descent, to integrate over ν 1 , we obtain the following contribution: where BFKL and D are defined in Eq. (A.2). Rewriting Eq. (2.5) in the coordinate representation we obtain (4.8) In Eq. (4.8) we have neglected the terms which are proportional to Q 2 T (see Eq. (2.5)), since, as we have argued, the typical Q T are small, and because these terms do not lead to additional correlations in the azimuthal angles. In Appendix B we calculate this integral and obtain the final expression for the double inclusive cross section: (4.9)

The CGC/saturation approach
The integral over k T in Eq. (B.6) has an infrared singularity with a cutoff at p T 2 , since we assume that p T 2 is the smallest momentum. This reflects the principal feature of the BFKL Pomeron parton cascade, which has diffusion, both in the region of small and large transverse momenta. On the other hand, we know that the CGC/saturation approach suppressed the diffusion in the small momenta [44], providing the natural cutoff for the infrared divergency. We expect that such a cutoff will be the value of the smallest saturation momenta: Q s (Y − y 1 ) or Q s (y 2 ), which will replace one of the p 2 T 2 in the dominator of Eq. (4.9). Therefore, we anticipate that for a realistic structure of the one parton shower cascade, (see Fig. 1c for example), the contribution for the double inclusive cross section will be different.
We need to specify the behavior of the scattering amplitude in the vicinity of the saturation scale. We have discussed the basic formulas [68] of Eq. (3.11), but for integration over the dipole sizes we need to know the size of this region. The scattering amplitude can be written in the form where ω(γ, 0) is given by Eq. (A.2), replacing 1 2 + iν ≡ γ and ξ = ln r 2 1 /r 2 2 . The saturation scale is determined by the line on which the amplitude is a constant (C), of the order one. This leads to the following equation for the saturation scale [62,68]: which results in the value of γ cr given by the equation and gives γ cr = 0.37, with the equation for the saturation momentum: Expanding the phase ω (γ, 0) Y − (1 − γ )ξ in the vicinity ξ = ξ − ξ s and γ = γ − γ cr we obtain (4.14) At first sight, Eq. (4.14) shows that the amplitude has a maximum at τ = r 2 1 Q 2 s = 1. However, this is not correct. Equation (4.14) gives the correct behavior for τ < 1, while for τ > 1 we need to take into account the interaction of the BFKL Pomerons and the non-linear evolution, generated by these interactions. The general result of this evolution is the fact that the amplitude depends on one variable [69][70][71][72] τ , i.e. N (τ ) (as it shows geometric scaling behavior). The peak at τ = 1 appears in From Eq. (4.15) we can conclude that the width of the distribution in r 2 1 is of the order of Q 2 s , but it depends crucially on the model for the Pomeron interaction. In Fig. 3a we plot this value for the behavior of the scattering amplitude deep in the saturation domain (see Ref. [72]).
This approach is not correct for τ → 1 and −∇ 2 N = 1.58 at τ = 1, but it starts to be small at τ > 2, which could be large enough to trust the formulas of Ref. [72]. At least such a conclusion can be justified considering the fit of the DIS data in the saturation model of Refs. [73,74], which is based on the idea of Ref. [75], and which has the correct behavior of the scattering amplitude, both deep in the saturation domain, and near τ = 1. Hence, we expect that ∇ 2 N decreases faster than we can see from Eq. (4.14). Bearing these conclusions in mind, we will calculate the contribution of Fig. 2, keeping all N in Eq. (4.8) in the vicinities of the saturation scales, by replacing We will show in the following that we cannot integrate over the dipole sizes, so that all six Pomerons will be in the vicinity of the saturation scale. At least two of the Pomerons occur either deep in the saturation domain, or in the perturbative QCD region. We believe that the largest contribution stems from the exchange of two Pomerons between rapidities y 1 and y 2 (see Fig. 3b), which are in the perturbative QCD region. Unfortunately, we cannot use the AGK cutting rules [38], which state that these Pomerons will not be affected by the Pomeron interaction, and the contributions of these interactions (see the red Pomeron in Fig. 3b) are canceled. Indeed, it has been proven that for the double inclusive production [47] they are not applicable in perturbative QCD. On the other hand, these Pomerons carry transverse momentum Q T , unlike the others in the diagram, which is larger than the saturation scale Q s (y 2 ); hence, their contributions are suppressed in comparison with the other Pomerons in Fig. 2. In addition our choice leads to the natural matching with the regionᾱ S y 12 < 1.
The integration over Q T will produce the same result as Eq. (4.7), as in the previous section. In Appendix C we discuss making estimates for the integrals over the dipole sizes which lead for p T 1 Q 2 s (y 1 ) to the following cross section:  Fig. 4 The ratio of Eq. (4.17) at W=13 TeV versus y 12 , assuming that the experiment has a symmetric pattern with Y − y 1 = y 2 = 1 2 (Y − y 12 ). The dotted line in a is for the estimates for the y 12 dependence of the Bose-Einstein contribution at small y 12 [11,76]. a and b The estimates in the leading order of perturbative QCD withᾱ S = 0.25. In c we take BFKL = 0.25 and Q 2 s (Y ) ∝ exp (λY ) with λ = 0.25. These numbers correspond to the BFKL phenomenology One can see that Eq. (4.17) demonstrates the additional suppression due to the infrared cutoff at Q s (y 2 ) instead of at p T 2 , as taken in the calculation of the simplest diagram. The factor exp (2 BFKL y 12 ) reflects the fact that two BFKL Pomerons between rapidities y 1 and y 2 are taken in the perturbative QCD region. It should be stressed that we can only trust our estimates for values of y 12 at which the exchange of the BFKL Pomeron with rapidity y 12 give a contribution smaller than C. This condition means that 1 (2 D y 12 ) e 2 BFKL y 12 < C. (4.18) Taking BFKL = 0.25 and Q 2 s (Y ) ∝ exp (λY ) with λ = 0.25 (these values correspond to the BFKL phenomenology) we see that the l.h.s. of Eq. (4.18) is smaller than 0.15 for y 12 ≤ 7. Therefore, we can trust our estimates shown in Fig. 4 for C > 0. 15. We take C = 0.3, which leads to the contribution of the shadowing corrections of the order of 30%.
In Fig. 4 we plot the ratio R as a function of y 12 for y 12 ≤ 7 (see Eq. (4.18). One can see that the ratio increases for large y 12 .

Azimuthal angle correlations
The azimuthal angle correlations stem from terms ( Q T · r i ) n in the vertices (see Eqs. (A.5), (A.6)). Indeed, after integrating over r i these terms transform to expressions of the following type [54]: , which lead to the term ( p T 1 · p T 2 ) m . We have illustrated in Eqs. (A.5) and (A.6) how these originate from the general form of the BFKL Pomeron vertices in the coordinate representation. From Eqs. (A.5) and (A.6) only terms proportional to Q T · r i n with even n appear in the expansion. Therefore, the azimuthal angle (φ) correlation function contains only terms cos 2n (φ), and it is invariant with respect to φ → π −φ. In other words, v n with odd n turn out to be zero. Hence, we have the first prediction: the value v n with odd n should decrease with y 12 , and their dependence should follow the dotted lines in Fig. 4a.
We return to Eq. (4.1) and integrate over Q T , collecting terms that depend on the angles between Q T and r i , which we have neglected in the previous section. As we have learned, the typical values of Q T ∝ 1/r 2 ∼ 1/r 2 where r 2 and r 2 are larger than r 1 and r 1 . In other words, we showed that the main contributions stem from the kinematic regions The typical Q T is determined by the largest dipoles and, therefore, we expect Q T ≈ 1/r 2 (1/r 2 ), as has been demonstrated above. Bearing these estimates in mind, we can replace vertices V ν 1 r 1 , Q T and V ν 2 r 1 , Q T in Eq. (4.1) by Eq. (A.6) in which we put Q T = 1/r 2 and Q T = 1/r 2 , respectively. Taking into account that r 1 /r 2 1(r 1 /r 2 1), we obtain At first sight Eq. (5.1) should enter two angles between Q T and r 1 and r 1 , respectively. However, in the integrand for integration over r i (see Eq. (B.1)) it depends only on one vector p T 1 . Therefore, after integration over all angles, we find that the angle φ in Eq. (5.1) is the angle between Q T and p T 1 . For vertices V * ν 1 r 2 , Q T and V * ν 2 r 2 , Q T in Eq. (4.1) we use Eq. (A.8). Finally, we need to evaluate the integral with better accuracy that we did in Sect. 5.1, keeping the dependence on the angle between Q T and r 2 . Note that the factor cos Q T · r 2 comes from exp i Q T · r 2 in Eq. (4.8). Taking this integral we substitute for the terms in parentheses in Eq. (5.1), |Q T | = 1/r 2 2 (1/r 2 2 ). The integral is equal to where n = Q T /Q T , and φ 2 is the angle between n and n 2 = r 2 /r 2 . In Eq. where n 1 = r 1 /r 1 , n 1 = r 1 /r 1 and n 2 = r 2 /r 2 . Deriving Eq. (5.4) we neglected the extra powers of r 2 1 /r 2 2 , which are small. Finally From Eq. (4.8) we can see that the integration over r i can be written in the form d 2 σ dy 1 d 2 p T 1 dy 2 d 2 p T 2 (Fig. 2 The tedious calculations of the integral over the dipole sizes in Eq. (5.6) are collected in Appendix D.
The values of v 2 and v 4 can be determined from the following representation of the double inclusive cross section: where ϕ is the angle between p T 1 and p T 2 . v n is calculated from v n,n ( p T 1 , p T 2 ), Introducing the angular correlation function as ; v n = √ v n,n ; (5.10) In Eq. (4.17) we have calculated the part of C ( p T , φ) which does not depend on φ, which coincides with C( p T , φ = 0) = R of Eq. (4.17) for Q s (Y − y 1 ) Q s (y 2 ). To calculate the contribution to C, which depends on φ, we need to take the separate integrals over ν 1 and ν 2 , since the terms, which are proportional to cos 2 (φ) and cos 4 (φ) do not have a pole at ν 1 = −ν 2 (see Eq. (5.5)). These integrations lead to the following extra factor in C ( p T , φ) − C ( p T , φ = 0): . We took factors proportional to p T from the expression for A (k T , p T 1 ) and A (4) (k T , p T 1 ) putting p T 1 = p T 2 = p T . To find the final correlation function and v 2,2 and v 4,4 , we need to collect all numerical factors that come from A (k T , p T 1 ), A (4) (k T , p T 1 ) and Eq. (5.6), and to integrate over φ, as given in Eq. (5.9).
Note that in the case of symmetric kinematics, where Y − y 1 = y 2 = 1 2 (Y − y 12 ), ξ = 0 and Eq. (D.7 vanishes. In this case, we have to use Eq. (A.5) instead of Eq. (A.6), keeping track of the corrections, which are proportional to ν i . As a result, we can consider ξ = 0 in Eq. (D.7), but we need to replace the factor ξ 2 by 1.
Equations (4.17) and (D.7) contain numerical uncertainties, which stem both from the values of the soft parametersμ soft and μ soft , as well as the values of the saturation scale at low energies, and from the integration in Eqs. (C.3) and (C.5), which were taken neglecting contribution from the region τ < 1. On the other hand, the contribution to the double inclusive cross sections of the diagram of Fig. 2 atᾱ S y 12 1 coincides with the contribution of Fig. 1b, (Fig. 1b) dy 1 dy 2 d 2 p T 1 d 2 p T 2 .
(5.12) Therefore, to obtain the realistic estimate we use the following procedure of matching: Fig. 1b , (5.13) where v 2 ( p T = 5 GeV) | Fig. 1b and v 4 ( p T = 5 GeV) | Fig. 1b are taken from Ref. [14] where the estimates were performed based on the model for soft interaction which describes all features of the soft interaction at high energy and provides an interface with the hard processes. Figure 5 shows the p T and y dependence of the v 2 and v 4 using Eq. (D.8) for normalization. In addition we take BFKL = 0.25 and Q 2 s (y) ∝ exp (λ y) with λ = 0.25. These values correspond to the BFKL Pomeron phenomenology. We believe that this figure illustrates the scale of rapidity dependence and will be instructive for future experimental observations.

Conclusions
In this paper we generalize the interference diagram that described the Bose-Einstein correlation for small rapidity differenceᾱ S y 12 1, to include the emission of the gluons with rapidities (y i ) between y 1 and y 2 (y 1 , y i < y 2 ). We calculate the resulting diagram in the CGC/saturation approach and make two observations, which we consider as the main result of this paper. The first one is a substantial decrease of the odd Fourier harmonics v 2n+1 as a function of the rapidity difference y 12 ( see Fig. 4c). The second result is that even Fourier harmonics v 2n have a rather strong dependence on y 12 , showing a considerable increase in the region of large y 12 (see Fig. 5). We believe that our calculations, which have been performed both for the simplest diagrams  and for the CGC/saturation approach, will be instructive for further development of the approach especially in the part that is related to the integration of the momenta transferred by the BFKL Pomerons.
We demonstrated in this paper the general origin of the density variation mechanism, whose nature does not depend on the technique that has been used. This mechanism has to be taken into account, since it leads to the values of the Fourier harmonics that are large and of the order of v n , which have been observed experimentally.
We hope that the paper will be useful in the clarification of the origin of the angular correlation, especially for hadron-hadron scattering at high energy. We firmly believe that the experimental observation of both phenomena, the sharp decrease of v n with odd n and the substantial increase of v n with even n as a function of y 12 , will be a strong argument for the CGC/saturation nature of the angular correlations.
where ω (ν, n) = 2ᾱ S Re ψ 1 2 where ψ(z) is the Euler ψ-function (see Ref. [61] formulas 8.36) and BFKL =ᾱ S 4 ln 2, D =ᾱ S 14ζ (3), ξ = ln r 2 1 /r 2 2 . Each term in Eq. (A.1) has a very simple structure, being the typical contribution of the Regge-pole exchange: the product of two vertices and Regge-pole propagator. From Eq. (A.2) one can see that at large Y the main contribution comes from the term with n = 0, and in what follows we will concentrate on this particular term. The vertices with n = 0 have been determined in Refs. [42,43,[58][59][60], and they have an elegant form in the complex number representation for the point on the two dimensional plane, viz., Using this notation the vertices have the following structure: At Q T → 0 this vertex takes the form iν × ⎛ ⎝ (ν − 2i) 8( Q T · r) 4 − 8( Q T · r ) 2 Q 2 T r 2 + 5 1 2 Q 4 T r 4 ) 2 12 ((2 + iν)(1 + iν)) 2 + 2(1 + iν) (2( Q T · r ) 2 − Q 2 T r 2 2 6 (1 + iν) 2 − 1 ⎞ ⎠ . (A.5) For small values of ν (which are related to the region of largeᾱ S Y 1), Eq. (A.5) can be simplified and reduced to the form (A.6) Using at ν 1 we obtain for Q 2 T r 2 1 The contribution of the first term in Eq. (A.1) can be reduced to the following form for the scattering amplitude of two dipoles with sizes r 1 and r 2 : 1. To integrate over ν, we use the method of steepest descent, and the expansion of ω (ν, 0) at small ν (diffusion approximation; see the second equation in Eq. (A.2)). N (r 1 , r 2 ; Y ) denotes the imaginary part of the dipoledipole sacttering amplitude at Q T = 0, which is related to the cross section. One can check that Eq. (A.9) has the correct dimension. whereγ = 1 − γ cr , we see from Eq. (C.1) that integration over r takes the form Recall that we consider Q s = Q s (Y − Y 1 ) in Eq. (C.3).
For p T 1 Q s (Y − y 1 ) we can replace e −i p T 1 · r 1 = 1 in Eq. (C.1). In this case the integral has the form Finally, collecting all numerical coefficients, we obtain d 2 σ dy 1 d 2 p T 1 dy 2 d 2 p T 2 (Fig. 3b) =