When gluons go odd: how classical gluon fields generate odd azimuthal harmonics for the two-gluon correlation function in high-energy collisions

We show that, in the saturation/Color Glass Condensate framework, odd azimuthal harmonics of the two-gluon correlation function with a long-range separation in rapidity are generated by the higher-order saturation corrections in the interactions with the projectile and the target. At the very least, the odd harmonics require three scatterings in the projectile and three scatterings in the target. We derive the leading-order expression for the two-gluon production cross section which generates odd harmonics: the expression includes all-order interactions with the target and three interactions with the projectile. We evaluate the obtained expression both analytically and numerically, confirming that the odd-harmonics contribution to the two-gluon production in the saturation framework is non-zero.


I. INTRODUCTION
Over the past decade long-range rapidity correlations between the produced hadrons were observed in heavy ion (AA) and high-multiplicity proton-nucleus (pA) and proton-proton (pp) collisions at RHIC [1-4] and LHC [5][6][7][8]. The novel di-hadron correlations enhance the ∆φ = 0 and ∆φ = π regions of the azimuthal opening angle between the hadrons and stretch over several units in the rapidity interval ∆η between the hadrons. Finding an explanation of these previously unobserved correlations is important for the understanding of the strong interactions dynamics in high energy hadronic and nuclear collisions.
Since the long-range rapidity correlations were first discovered in heavy ion collisions, it is natural to ascribe their origin to the dynamics of quark-gluon plasma (QGP) produced in such collisions. A vigorous activity is under way to account for the long-range rapidity di-hadron correlations within hydrodynamic models of QGP [9][10][11][12][13][14][15][16]. However, these approaches suffer from one conceptual problem: long-range rapidity correlations cannot originate at later proper times in the collision, when plasma is produced. There exists a causality argument [17] demonstrating that the later the proper time of the interaction, the narrower in ∆η the corresponding correlation will be. Therefore, hydrodynamics, being applicable at a relatively late times, is not likely to generate these long-range rapidity correlations: rather, the rapidity correlations have to be included into the initial conditions for hydrodynamic evolution. Thus, the question about the origin of the long-range rapidity correlations remains, with the initial-state early-proper-time dynamics being the most probable suspect.
Long-range rapidity correlations were first advocated in the saturation/Color Glass Condensate (CGC) framework in [17] (see also [18]). (See [19][20][21][22][23][24][25] for reviews of saturation/CGC physics.) The weakly-coupled CGC dynamics can generate long-range correlations which at large rapidity intervals ∆η are independent of ∆η, in qualitative agreement with the experimentally observed correlations. Efforts to obtain similar correlation function behavior originating in a different initial-state dynamics scenario so far came up short: at strong coupling, calculations employing the anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [26,27] so far give a correlation function that grows with ∆η [28], in disagreement with the experiment.
The ridge correlation may result from the CGC dynamics giving long-range correlations in rapidity, with the subsequent hydrodynamic evolution generating the azimuthal collimation observed in the data [29]. This explanation of the ridge phenomenon requires a thermal medium to be generated in high-multiplicity pp and pA collisions. Thermalization of the medium produced in heavy ion collisions requires multiple interactions between the quarks and gluons: it is natural to assume that thermalization in high-multiplicity pp and pA collisions, if it does take place, would require similar multi-parton interactions. Perhaps a more conservative way of generating the ridge correlation is due to the so-called "glasma graphs" proposed in the saturation/CGC framework in [17]: they require only a double interaction in both the projectile and the target. The "glasma graphs" generate both the long-range correlation in rapidity [17,[29][30][31][32][33][34][35][36][37][38][39][40][41], along with the near-and away-side ridge correlations [42,43]. Since the "glasma graphs" had been originally put forward in [17], the corresponding two-gluon production cross section was improved by including multiple interactions (saturation effects) with one of the nuclei (the target) in [38,42]. As usual in the saturation physics, the effect of extra rescatterings appears to be mainly in screening the infrared (IR) divergences [43]. In addition, one can show that multiple rescatterings violate the k T -factorization formula one could conjecture arXiv:1802.08166v1 [hep-ph] 22 Feb 2018 for two-gluon production based on the "glasma graphs" alone [43].
The approach based on the "glasma graphs" with the infrared screening provided by the saturation scale Q s enjoyed successful phenomenology [39][40][41]. However, the two-gluon production cross section given by the "glasma graphs" [17] along with the saturation corrections in the target hadron/nucleus [38,42] turned out to be invariant under k 1 ↔ k 2 interchange and, also, is separately invariant under k 1 → −k 1 or k 2 → −k 2 replacements. Here k 1 and k 2 are the transverse momenta of the two produced gluons. The result of this symmetry is that the corresponding gluon-gluon correlation function only contains even azimuthal harmonics v 2 2n {2}, with all the odd harmonics v 2 2n+1 {2} being zero. At the same time, odd azimuthal harmonics have been observed in the di-hadron correlators measured at RHIC and at LHC [8,44,45]. Odd harmonics are somewhat smaller than the even ones, but are non-zero. This experimental result presented a conundrum for the saturation community: can the saturation dynamics account for the observed long-range rapidity correlations with the non-zero odd azimuthal harmonics?
To resolve this ambiguity, several observations have been put forward. In [46] the authors observed that the symmetry of the di-gluon correlator under k 1 ↔ k 2 , k 1 → −k 1 and k 2 → −k 2 is "accidental", and is not required by the symmetries in the problem. They have then argued that odd harmonics may arise if one includes saturation effects in the wave function of the projectile: thus to find odd harmonics one needs to augment the existing calculations of the two-gluon production cross section [38,42], in which the saturation effects were only included in the interaction with the target. The idea of saturation corrections in the projectile being responsible for the odd harmonics was developed in [47], where it was shown that such corrections indeed have the potential to generate odd harmonics by violating the k 1 → −k 1 and k 2 → −k 2 symmetries of the two-gluon correlation function. In addition, a numerical simulation of the classical gluon fields produced in heavy ion collisions in the McLerran-Venugopalan (MV) model [48][49][50] appears to generate odd harmonics as well [51,52]: since the difference between this numerical result for inclusive two-gluon production and the expressions obtained in [38,42] is due to the saturation effects in the projectile, and since the calculation in [38,42] only gave even harmonics, it is natural to conclude that the odd harmonics likely originate in the higher-order projectile interactions.
Our goal in the present paper is to construct a complete expression for the part of the two-gluon inclusive production cross section responsible for the odd harmonics by calculating the first saturation correction in the interactions with the projectile. Our goal is complicated by the fact that even for the single inclusive gluon production cross section the first saturation correction in the projectile has not been found yet. However, partial results exist in [53,54]. Let us first consider the single inclusive gluon production cross section. Suppressing the transverse momentum dependence, the cross section in the classical MV model can be written as (in the MV model power counting, using the approach to it from [55,56]) where A 1 and A 2 are the atomic numbers of the projectile and target nuclei respectively, α s is the strong coupling constant, B is the impact parameter between the nuclei, b is the transverse position of the gluon, and the angle brackets denote averaging in the wave functions of the projectile and target nuclei, which is equivalent in the MV model to averaging over their color charge densities ρ p and ρ T [48-50, 55, 56]. The function f was only studied numerically [57][58][59][60]. Analytically we only know its expansion in either one of its arguments. Assuming that the projectile is a dilute object with α 2 s A 1/3 1 < ∼ 1 one can expand the cross section in this parameter The function f 1 is known analytically from the gluon production cross section in the proton-nucleus (pA) collisions [61,62], since it comes in with the term involving only one power of A 1/3 1 , that is, only one nucleon in the projectile, making the projectile effectively a "proton" in this power counting. Functions f 2 , f 3 , . . . are not known analytically at present.
The function f 2 gives the first saturation correction in the projectile, since it comes in with two powers of A 1/3 1 , corresponding to interactions with two nucleons in the projectile nucleus. The efforts to calculate f 2 analytically was started in [53] and more recently revisited in [54]. The calculation of the order-α 2 s A 1/3 1 2 correction implies including an order-α 2 s correction to the projectile interaction as compared to the leading order-α 2 s A 1/3 1 term from [61]. This order-α 2 s correction involves interaction with the extra nucleon in the projectile, which brings in an additional A amplitude was found. No one has yet analytically calculated the order-α 2 s correction to the same amplitude to complete the efforts to determine f 2 ! The same philosophy applies to the two-gluon production. For the classical two-gluon production cross section one can write with the new unknown function h. Here k 1 and k 2 are the gluons' transverse momenta, while b 1 and b 2 are their transverse positions. Again, assuming a dilute projectile we expand in The function h 1 can be found from the results of [38,42]. As described above, this part of the two-gluon production cross section generates only even harmonics. Finding the function h 2 requires a rather lengthy calculation. However, as we will argue below (basing our argument on [47]), the part of h 2 responsible for the odd harmonics can be found using the results of [54]. The corresponding diagrams are shown below in Fig. 4. We thus find the part of the two-gluon production cross section responsible for odd harmonics at the order α 2 s A 1/3 1 3 : it is given by Eq. (17). This is the leading contribution to the odd harmonics resulting from the two-gluon correlation function.
The paper is structured as follows: in Sec. II we derive the contribution (17) to the two-gluon production cross section giving the odd harmonics. We then proceed in Sec. III by evaluating the expression (17) analytically to the lowest order in the interactions with the target using the Golec-Biernat-Wusthoff (GBW) [63,64] approximation to the full MV interaction with the target. The resulting two-gluon production cross section is given in Eq. (77). The most important conclusion we draw here is that Eq. (77) is non-zero: hence the saturation dynamics does generate non-trivial odd harmonics. Interestingly, a prominent contribution to the correlation function comes from the δfunctions resulting from the gluon Hanbury-Brown-Twiss (HBT) [65] diagrams, of the same general type as those advocated in [43,66]. The odd-harmonics correlation function resulting from Eq. (17) with all-orders interactions with the target is evaluated in Sec. IV numerically, with the resulting odd harmonic coefficients plotted in Fig. 6 and the odd part of the two-gluon correlation function shown in Fig. 7. Keeping the interaction with the target to the lowest non-trivial order in the numerical simulations we observe qualitative and even quantitative agreement with Eq. (77). We conclude in Sec. V by observing that saturation/CGC dynamics does lead to the odd harmonics in di-gluon correlation functions, which may allow for further successes of the saturation approach to the correlation function phenomenology.

II. APPEARANCE OF ODD HARMONICS IN GLUON PRODUCTION DIAGRAMS
A. General discussion Let us begin by analyzing how the two-gluon production cross section involving two classical gluon fields can in principle generate the odd azimuthal harmonics. As discussed in the Introduction, we need to violate the k 1 ↔ k 2 , k 1 → −k 1 and k 2 → −k 2 symmetry of the the two-gluon production cross section. For two gluons originating from two identical classical fields the k 1 ↔ k 2 symmetry appears impossible to break: we thus conclude that the classical two-gluon production cross section should always be k 1 ↔ k 2 symmetric. This implies that the dependence on the azimuthal angle ∆φ between the two produced gluon enters the cross section via terms proportional to cos(n∆φ) with the integer values of n, while terms proportional to sin(n∆φ) do not enter the cross section.
We thus see that the only way to generate odd harmonics is to violate the k 1 → −k 1 and/or k 2 → −k 2 symmetries. Let us describe how this violation may happen in general. Imagine a production cross section for a particle with the transverse momentum k = (k 1 , k 2 ), possibly alongside with a number of other particles whose momenta we do not explicitly display below for brevity. The production cross section is proportional to Here M (x) is the Fourier transform of the scattering amplitude M (k) into transverse coordinate space and the asterisk denotes complex conjugation. We want to find a condition under which this cross section has a contribution that changes sign under k → −k. To do this, imagine that the scattering amplitude in the transverse coordinate space can be written as for instance due to an expansion in the coupling constant (that is, M 1 is the leading contribution, and M 3 is one of the higher-order corrections with the ellipsis denoting other higher-order corrections). Using Eq. (6) in Eq. (5) (while keeping only M 1 and M 3 ) one can easily see that only the interference terms between M 1 and M 3 may lead to a contribution odd under k → −k. Keeping only the interference terms we write their contribution to Eq. (5) as Requiring that the expression (7) changes sign under k → −k results in the following condition on M 1 and M 3 : Since M 1 and M 3 , in general, are very different functions of their arguments, Eq. (8) is most easily satisfied by requiring that which, in turn, means that M 1 (x) M * 3 (y) is imaginary. Therefore, we conclude that the phases of M 1 and M 3 (in transverse coordinate space) should be off by π. More generally, for the higher-order corrections to the leading-order amplitude M 1 to generate k → −k odd contribution one needs a phase difference between M 1 and the full higher-order amplitude M 3 + . . ., where the ellipsis denote other higher-order corrections, some of them possibly of the same order as M 3 . Hence, to find the odd harmonics we need to find a higher-order correction to the results of [38,42] that generates a phase difference between the higher-and leading-order amplitudes (see also Appendix A).
The situation is not dissimilar to the single transverse spin asymmetry (STSA) which is observed in the scattering of the transversely polarized protons on the unpolarized ones and in semi-inclusive deep inelastic scattering (SIDIS) on a transversely polarized proton [67][68][69][70][71][72][73]. There, generating the asymmetry requires the amplitude dependence on the transverse proton polarization and a phase difference between the leading-order amplitude and the higher-order correction to it: the asymmetry is then given by the interference between the two amplitude contributions [74][75][76][77][78]. The difference between the STSA case and the odd harmonics problem at hand is that the phase difference in the former case is between the momentum-space amplitudes, while for the latter, Eq. (9) implies a phase difference between the Fourier transforms of the amplitudes into the transverse coordinate space.

B. Order g 3 amplitudes
Let us now apply the insight we obtained in the previous section to the calculation of gluon production in the saturation framework. As discussed in the Introduction, we will assume that our calculation resums all orders in the interaction with the target via the Wilson lines for quarks and gluons traversing the target. We will perform expansion of the gluon production amplitude to the first two non-trivial orders in the interaction with the projectile.
The leading-order gluon production is described by the amplitude * λ · M 1 which is given by the diagrams shown in Fig. 1. They yield [61] with the adjoint and fundamental Wilson lines directed along the "+" light cone, defined by the direction of the projectile motion. All the notation is explained in Fig. 1: the outgoing gluon has polarization λ with λ = −(λ, i)/ √ 2 and color a. The horizontal straight lines in Fig. 1 represent (valence) quarks in the nucleons of the projectile nucleus. The vertical bar denotes the target, which is Lorentz-contracted to a shock wave, whose interaction with the projectile quark and emitted gluon is very fast, practically instantaneous on the time scales of the projectile wave functions: hence the gluon produced in Fig. 1 can be emitted either before or after the interaction with the target, but not during that interaction [61]. The gluon can be emitted by either one of the many projectile quarks: we show the emission by only one of the quarks, with the sum over other emissions implied implicitly. Following the power counting from [54,56] we assume that the interaction with the target is strong, , and, therefore, all the Wilson lines are also of order one, U ∼ V = O(1). Hence we see that M 1 = O(g).
The next order correction to M 1 could be an order-g 2 amplitude with two gluons emitted by the quarks (with only one of these gluons being tagged on if one wants to calculate the single inclusive gluon production cross section). However, such contribution in the counting of powers of α 2 s A 1/3 1 can be rearranged and absorbed into the complex conjugate amplitude where one uses retarded Green functions instead of Feynman propagators for the gluons (see [53,54]). This way, the first correction to M 1 is O(g 3 ). Some of (the large number of) the relevant diagrams are shown in Fig. 2: in the notation of the previous sub-section, these are the diagrams contributing to x 1⊥ The diagrams contributing to M 3 and involving two nucleons from the projectile, as shown in Fig. 2, were calculated in [54], with the result * The label M ABC 3 reflects the fact that Eq. (13) contains contributions of the diagrams A, B and C in the notation of [54]. Here Λ is the IR cutoff.
There are also order-g 3 diagrams involving interaction with only one nucleon in the projectile, labeled D and E in [54]: those are shown in Fig. 3 and their sum is given by * Both equations (13) and (14) are written in transverse coordinate space. Therefore, the formalism of Sec. II A applies. Since M 1 in Eq. (10) contains a factor of i, and we require a phase shift of π between M 1 and M 3 , we need the "real" part of M ABC to obtain the odd harmonics via the formula Eq. (7). We conclude that the odd harmonics should be given by the interference between M 1 from Eq. (10) with the "real" part of M ABC Eqs. (13) and (14). Note that correlators of Wilson lines in the MV model are all real, if one keeps the leading C-even parts: hence by the real part of M ABC we mean only the real part of the associated light-cone wave functions, such that the Re sign does not apply to the Wilson lines. With this caveat, the real parts are * and * The factor of (V b 1 ) in Eq. (16) denotes the Wilson line of the spectator quark line which simply cancels the same (but conjugate) Wilson line in the diagrams of Fig. 1 contributing to M 1 (but not shown explicitly in the expression for M 1 ) when D 3 is used in place of M 3 in the interference formula (7): due to such a trivial role, the dependence on b 1 is not explicitly included in the arguments of D 3 . According to Eq. (7), the odd azimuthal harmonics can be given by the interference between M 1 from Eq. (10) with M 3 and D 3 from Eqs. (15) and (16) respectively.
C. Odd-harmonic part of the two-gluon production cross section at the order α 4 s Equation (7) is written down for the single-gluon production, where it is impossible to have a cross section contribution which is odd under k → −k in the case of unpolarized initial and final states: indeed, if we take the interference of a single M 1 with a single M 3 contribution from above, as shown in Fig. 8 below, the color averaging would give zero.
However, we are interested in odd harmonics in two-gluon production. Hence we need to iterate Eq. (7) twice, once for one produced gluon, one for another. The corresponding diagrams we need to calculate are given in Fig. 4: they involve quarks from three nucleons in the projectile. Each diagram in Fig. 4 denotes a class of diagrams including all the diagrams obtained from it by reflecting either one (or both) of the connected gluon interactions with respect to the final-state cut (the vertical solid line in Fig. 4). Given the above ingredients, the diagram calculation is straightforward. The part of their contribution to the two-gluon production cross section that gives odd harmonics is dσ odd The dot-product in Eq. (18) arises after a sum over the final-state gluon polarizations λ. Each diagram of A-F in Fig. 4 contributes four terms to Eq. (18), with the terms in the latter ordered in the same way as the diagrams A-F: the first four terms come from the diagram A, the next four terms are from B, etc. FIG. 4: The contributing diagrams for the odd-harmonic part of the two-gluon production cross section.
Equation (17) is the main analytic result of this work: it gives the leading contribution to the inclusive two-gluon production cross section in the saturation framework which generates odd azimuthal harmonics in the correlation function. It involves interactions with three nucleons in the projectile nucleus and all-order interactions with the nucleons in the target nucleus. The Wilson line correlators also include the linear and non-linear small-x evolution between the target nucleus and the produced gluon [79][80][81][82][83][84][85][86][87][88].
However, we are not done yet: to convincingly demonstrate that Eq. (17) does generate odd harmonics in the correlation function, one needs to show that it is not zero. While there appear to be no symmetries requiring Eq. (17) to be zero, an actual evaluation of this expression is required to establish this with full certainty. Below we will evaluate Eq. (17) using a quasi-classical MV-based approximation for the correlators of Wilson lines, expanding those correlators to the lowest non-trivial order, which involves interaction with three nucleons in the target nucleus.
We will use the following Fourier representations of the amplitude contributions to evaluate Eq. (17): * and * III. EVALUATING THE ODD-HARMONIC PART OF THE TWO-GLUON PRODUCTION CROSS SECTION: ANALYTIC APPROACH

A. Diagram A
Our goal here is to evaluate the expression (17). We begin with the first four terms in Eq. (18), whose contribution to Eq. (17) is proportional to and corresponds to the diagram A in Fig. 4 along with all the gluon reflections with respect to the final-state cut. Let us evaluate the first line of Eq. (21) first. Disregarding the conjugation of the V s, since they cancel in the net expression Eq. (21) anyway, we write This way Eq. (21) can be written as Plugging in Eqs. (19) and (20) into the first line of Eq. (23) and taking all the fundamental color traces while averaging over colors of all three quarks in the projectile yields Here we took a step back and reintroduced the polarization vectors (e.g. we replaced M 1 · M * 3 → λ * λ · M 1 λ · M * 3 ) to simplify the derivation.
The interaction with the target from Eq. (24) (the expression in angle brackets) is depicted in Fig. 5 diagrammatically: there we have used the crossing symmetry [89,90] to reflect all the adjoint Wilson lines, represented by gluons, from the complex conjugate amplitude into the actual amplitude. We evaluate this interaction with the target in the large-N c limit, in which it reduces to a sum of products of three fundamental quadrupoles [91,92], Here "permutations" denote 63 other 3-quadrupole products. Hence the exact analytic evaluation of Eq. (24) appears to be prohibitively complicated. Instead we will expand the interaction with the target (25) to the lowest non-trivial order in gluon exchanges, or, equivalently, in the saturation scale. To do this, we employ the large-N c fundamental quadrupole amplitude evaluated in the quasi-classical MV approximation in [91] (see Eq. (14) there) with Here Q s0 is the gluon saturation scale of the target taken in the quasi-classical limit [93]. Inserting Eq. (26) into Eq. (25) with all the permutations included, and expanding in Q s0 to the lowest non-trivial order, which is order-Q 6 s0 corresponding to the 6-gluon exchange, yields with Alternatively, these results can be obtained by expanding the Wilson lines in Eq. (28) to the first nontrivial order in the target field and averaging with respect to the target ensemble; this leaves rather complicated color sums involving eight anti-symmetric structure constants, see Appendix B. This method is used to compute the diagrams B and C below.
In arriving at Eq. (29) we have neglected the logarithms of Eq. (27) by putting them equal to 1. Henceforth we will refer to this approximation as the Golec-Biernat-Wusthoff (GBW) approximation, since it was used in [63,64] for the Glauber-Mueller dipole amplitude [93] to successfully describe the small-x HERA data on structure functions. Despite this phenomenological success, this is admittedly a dangerous approximation for observables depending on transverse momenta: the expansion to the lowest order of gluon exchanges must be valid at high transverse momenta, where logarithms from Eq. (27) are most important. For a number of transverse-momentum dependent gluon distributions (gluon TMDs) at small x, evaluated in the quasi-classical approximation, neglecting such logarithms either leads to a zero result (e.g. for h (1) g , as discussed in [94]) or to an incorrect high-k T asymptotics, which becomes Gaussian in k T instead of giving the correct ∼ 1/k 2 T power law. We are encouraged, however, by the fact that such problem does not arise for the single inclusive gluon production cross section [21,61], where neglecting transverse coordinate logarithms still leads to the correct power-law high-k T asymptotics (but does not capture the factor of ln(k T /Λ)).
Substituting Eq. (28) into Eq. (24) we obtain It is more convenient to evaluate the expression (30) by pieces, which is allowed by the factorized form of each term in Eq. (29). Using for any well-behaved function f (q) we can evaluate The other term depending on the same gluon polarization λ is Using we can multiply Eq. (32) by Eq. (33) and sum over polarizations obtaining where we have defined Here ij is the Levi-Civita symbol in two dimensions. Using this result in Eq. (30) yields g 8 32 Q 6 s0 d 2 l (2π) 2 e il·b 21 −ik 1 ·b 23 d 2 l (2π) 2 e −il ·b 32 +ik 2 ·b 31 H ijn (k 1 , l) H * i j n (k 2 , l ) − δ ij δ jn δ ni + 2 δ in δ jn δ i j − δ in δ jn δ i j + 8 δ in δ nj δ ji − δ ii δ jj δ nn − δ ij δ i j δ nn − δ ij δ jn δ n i − δ in δ jj δ i n + 2 δ ij δ nj δ i n − δ ii δ jn δ j n + 2 δ in δ ji δ n j − δ ij δ i n δ n j .
Next, in Eq. (17) we approximate (cf. [42]) The rationale here is that the nuclear profile function T 1 (b) does not vary much across the size of a single nucleon. At the same time, from the standpoint of our perturbative calculation, integrating over the impact parameter over the distances comparable to the nucleon radius is approximately equivalent to integrating to infinity. Integrating Eq. (37) over b 12 and b 23 , and over l with l , we arrive at Summing over all the indices in Eq. (39) yields Finally, anti-symmetrizing Eq. (40) under k 1 → −k 1 and k 2 → −k 2 according to the prescription (23) we obtain the expression for the contribution to the two-gluon production cross section (17) of the first four M-terms corresponding to the diagram A in Fig. 4, For brevity we are omitting the factor of d 3 ] 2 which will be reinstated later.

B. Diagram B
Now let us evaluate the diagram B from Fig. 4, along with all of its reflections with respect to the final-state cut, in the same GBW approximation. The first of the next four M-terms from Eq. (18) gives Color structure of the interaction with the target is different from the case of diagram A, with again in the GBW approximation. Employing the technique from the previous Subsection to evaluate equation (42) yields where H ijn and H * i j n are defined by Eq. (36). Integrating Eq. (45) over b 12 and b 23 has to be done with care: naive integration over these variables up to infinity puts l = k 1 and l = 0, which are ill-defined limits in H ijn and H * i j n . To deal with this issue we integrate over b 12 and b 23 in the range 0 < b ij < R with R the projectile radius obtaining 1 Note that the rest of the four terms in Eq. (18) that correspond to diagram B in Fig. 4 anti-symmetrize Eq. (46) under k 1 → −k 1 and k 2 → −k 2 , giving Employing we can write an essential part of Eq. (47) as We will first integrate over the angles of l and then take the R → ∞ limit. To integrate Eq. (49) over the angles we note that where φ l is the azimuthal angle of the transverse momentum vector l, while q > = max{k, l} and q < = min{k, l}. In arriving at Eq. (50) we have used the integrals listed in Appendix C and employed the identity Write Eq. (50) as Using this along with in Eq. (47) yields Performing the summation over all repeated indices we arrive at the contribution of the diagram B as follows from comparing Eqs. (50) and (52). Equation (55) contributes only to the first azimuthal harmonic coefficient in the two-gluon correlation function. The leading IR divergent part of the |l − k 1 | and l integrals can be readily extracted from the above results, yielding with the constant where v = l R or v = |l − k 1 |R depending on the integral and R = 1/Λ. In arriving at Eq. (57) we have assumed that the IR divergence in the |l − k 1 | and l integrals is regulated by Λ = 1/R. Similar power-law divergences arise in the leading-order calculations of two-gluon production in the saturation framework [42,43]. In [43] it was argued that they are regulated by the saturation effects in the projectile. If this is the case in our calculation as well, then we should go back to Eq. (45) and subtract from it the k 1 → −k 1 and k 2 → −k 2 terms from it while adding the k 1 → −k 1 , k 2 → −k 2 term to fully account for the diagram B and its mirror reflections, obtaining × [δ in δ jj δ ni − δ ij δ jn δ ni + 2 δ in δ jn δ i j + 2 δ in δ ji δ j n + δ ii δ jn δ nj − 2 δ ij δ ji δ nn − δ ii δ jj δ nn − 2 δ ij δ i j δ nn − 2 δ ij δ jn δ i n + 2 δ ij δ nj δ i n ].
Next we regulate all the IR momentum singularities in H-factors of Eq. (59) by This is indeed not an exact way to account for the saturation effects in the projectile and should be understood in a qualitative way. After such regularization, integrating Eq. (59) over b 21 and b 23 to infinity and one gets zero, Hence the contribution of diagram B is either given by Eq. (57) or is zero depending on whether the power-law IR divergences are regulated by the IR cutoff Λ or by the saturation scale Q s due to higher-order interactions with the projectile. A more careful analysis of our main result (17) along with the inclusion of even higher-order corrections in the interaction with the projectile are needed to resolve this ambiguity. While the former can be accomplished with sufficient amount of hard work applied to Eq. (17), the latter would require diagram calculations beyond those done in [54], which is a significantly larger effort. We leave this investigation for further work, noting here that the diagram B, even if it is not zero and is given by Eq. (57), would only contribute to the first azimuthal harmonic, and is not going to cancel the contribution (41) of the diagram A, and, as we will shortly see, the contributions of other diagrams in Fig. 4.

C. Diagram C
Moving on to the diagram C along with the mirror reflections of its gluon lines with respect to the final-state cut we see that the first of the next four M-terms from Eq. (18) gives Evaluating the interaction with the target in the GBW approximation and performing a number of integrals along the lines used for other diagram above, we rewrite this as s0 32 d 2 l (2π) 2 e il·b 32 +ik 1 ·b 13 d 2 l (2π) 2 e −il ·b 32 −ik 2 ·b 13 H ijn (k 1 , l) H * i j n (k 2 , l ) [2 δ i,n δ i ,n δ j,j + 8 δ i,i δ j,j δ n,n − δ i,n δ i ,j δ j,n − δ i,j δ i ,n δ j,n − δ i,n δ i ,j δ j,n − δ i,j δ i ,n δ j,n − δ i,n δ i ,j δ j ,n − δ i,j δ i ,n δ j ,n − δ i,n δ i ,j δ j ,n − δ i,j δ i ,n δ j ,n + 2δ i,i δ j,n δ j ,n + 2δ i,j δ i ,j δ n,n ].
The delta-functions indicate that these are indeed gluon HBT and anti-HBT diagrams in the terminology of [42]: such contributions were previously observed in [17,42] in the leading-order (even-harmonics) two-gluon production calculations. Summing over all the indices we arrive at where the integral is dominated by the IR divergences, which we regulated by the IR cutoff Λ. (Once again, it is likely that Λ should be replaced by the saturation scale of the projectile [43].) We thus have for the diagram C, In the GBW approximation the diagram gives only HBT (and anti-HBT)-type correlations. It appears likely that in the full MV model, without the GBW simplification, diagram C would also lead to non-HBT types of terms (that is, the contribution of the full diagram C is probably not limited to the delta-functions of Eq. (67)).

D. Diagram D
Let us evaluate diagram D next. The first of the next four M-terms from Eq. (18) gives Evaluating the interaction with the target and performing a number of integrals we rewrite this as s0 128 e ik 1 ·b 21 +ik 2 ·b 32 with the exact structure of the Kronecker delta functions not being important for what follows. Therefore, we do not show it explicitly. Integrating over the impact parameters b 21 and b 32 we arrive at Subtracting from Eq. (70) the k 1 → −k 1 and k 2 → −k 2 terms and adding the k 1 → −k 1 , k 2 → −k 2 term gives zero: (Here D 1 , . . . , D 4 denote the four contributions one could get from the diagram D in Fig. 4 by reflecting the gluons with respect to the final-state cut.)

E. Diagram E and F
We conclude by evaluating diagrams E and F in Fig. 4. The next four M-terms in Eq. (18) give Evaluating the interaction with the target and performing a number of integrals we rewrite this as Integrating over the impact parameters b 31 and b 32 we arrive at While the H * i j n (k 2 , k 2 ) in Eq. (74) requires proper regularization, like it was done above in evaluating the diagram B, it is clear that anti-symmetrization of (74) under k 1 → −k 1 gives zero, Diagram F is obtained from the diagram E by interchanging k 1 ↔ k 2 . Therefore, the sum of the four F-graphs is also zero, F. Sum of all diagrams Adding all the above results for diagrams A, B, . . . , F together we get The most important conclusion of our approximate analytical calculation that led to Eq. (77) is that we get a nonzero contribution, which would yield odd azimuthal harmonics in the two-gluon (and, hence, di-hadron) correlation function. Hence, our main exact result (17) is not zero either. We conclude that we have identified a source of odd harmonics in the two-gluon correlation function in the saturation framework. Interestingly, the HBT term in Eq. (77) has the largest IR divergence and, therefore, dominates Eq. (77) and the corresponding correlation function. As was already suggested in [17], fragmentation may "wash out" the deltafunctions in the HBT term to some degree, such that the hadronic correlation function would not contain the literal delta-function correlations from Eq. (77). While the delta-function shape may not survive fragmentation, the HBTtype gluon correlations should still manifest itself in the hadronic correlation function, and may also dominate in it just like the HBT correlations dominate the two-gluon correlator.
Comparing Eq. (77) to the leading-order two-gluon production cross section in the classical formalism, e.g., to Eqs. (49), (58) and (59) in [42], we observe the following: Eq. (77) has an extra power of Q 2 s0 as compared to the leading-order cross section, corresponding to an extra rescattering in the target. Similarly, Eq. (77) includes an extra factor of α 2 s as compared to the leading-order cross section: this factor arises through the extra interaction with the projectile. In addition, the IR divergence in Eq. (77) is ∼ 1/Λ 4 , which is a higher degree of divergence than ∼ 1/Λ 2 observed in the leading-order result [43]. This is an indication that the higher-order rescatterings in the target and in the projectile may screen the IR divergence in the leading-order expression for two-gluon production cross section, effectively replacing 1/Λ 2 by 1/Q 2 s . Further work is needed to firmly establish this conclusion.

IV. EVALUATING THE ODD-HARMONIC PART OF THE TWO-GLUON PRODUCTION CROSS SECTION: NUMERICAL APPROACH
In order to model the distribution of the color charges in the projectile numerically, instead of the point-charge approach used in Sec. III it is more convenient to use an alternative one based on the introduction of a continuous (light-cone) color density ρ p (x). To compute an observable, one has to averaged the corresponding operator with the weight functional W [ρ p ], similar to the way we account for the target ensemble. Nevertheless, the treatment of the projectile here is still different from the treatment of the target; the projectile charge density is considered to be dilute facilitating the expansion of the projectile Wilson lines. In [56], it was demonstrated that, in the classical approximation, the approach based on the continuous color density is completely equivalent to the one used in preceding sections.
Here we introduce the density-dependent operator describing production of a gluon with momentum k 1 , E 1 dN d 2 k1 ρ p , ρ T , see Appendix D for details. The single and double inclusive gluon multiplicities are then given by Thus, in the classical approximation, the double inclusive production factorizes on the configuration-by-configuration basis. This leads to the following interesting fact: the two-particle cumulants of azimuthal anisotropy are always positive if the magnitudes of the momenta of the produced gluons coincide, i.e., where Our goal is to compute the odd two-gluon harmonics. For this we define with the explicit expression presented in Appendix D. We thus can extract and the angular-averaged V 0 (k) contributes to the normalization of the cumulants of the azimuthal anisotropy and as such has to be computed to the leading order only. Note that in this case V 0 (k) is manifestly real. The two-gluon cumulants for odd harmonics are then The averages with respect to the ensembles . . . ρp,ρ T are performed in the Gaussian MV model. The target MV configurations are generated as described in [95,96] and are complemented by the MV configurations for the projectile, which are computed for a single slice in x − . We also assume an infinite target with the color charge density defined by a single number describing the color charge fluctuations, µ 2 = const. For the projectile we assume a finite size R p = 1/Q sp with a Gaussian profile, that is Strictly-speaking, for our analytic approach to be valid one has to have c p 1, see e.g. [42]. In our numerical simulations we fixed c p = 1/2. Our choice of the free parameters is driven by the predisposition to simplify the problem as, in the current paper, we have no ambitions to describe the experimental data quantitatively.
The odd harmonic coefficients v 3 and v 5 resulting from our numerical simulations are shown in Fig 6. We warn the reader that the calculations were performed without a high multiplicity bias; the gluon fragmentation was not accounted for; the Glauber fluctuations were neglected. Additionally, the parameters we used to model the projectile wave function are not very realistic. Nevertheless Fig. 6 can be viewed as a proof of the concept demonstrating the presence of the odd azimuthal harmonics in the saturation/CGC formalism; the magnitude of v 3 {2} is of the same order as observed experimentally in pA collisions at LHC.
Conducting phenomenologically relevant calculations would require a significant numerical effort and will be reported elsewhere [97].
By numerically computing the two-gluon correlation function, we are able to reproduce the main features of Eq. (77), including the (anti-)HBT peaks and the contribution from the diagram A. Consider the odd part of the angular correlation function C odd (|k 1 |, |k 2 |, ∆φ) defined by Performing numerical configuration-by-configuration analysis it is possible to extract C odd (|k 1 |, |k 2 |, ∆φ). The numerical results are depicted in Fig. 7  connected by the straight lines are the results of the numerical calculations, the green curve represents the fit inspired by the analytical result in Eq. (77), namely with positive a C and a A . Here the first term corresponds to the HBT peak, which, for our Gaussian projectile wave function in Eq. (86) is a Gaussian itself. This replaces the Dirac δ-function peaks, originating in Eq. (77) from the Fourier transform with respect to the impact parameter over an infinite projectile (cf. [42]). The second term in Eq. (88) corresponds to the contribution from the diagram A in which the denominator was regularized by the scale of order 1/R p . As shown in Fig. 7, the fit reproduced the numerical calculations quite well, confirming our analytical findings from (77). We believe that the disagreement between the green solid line and some of the numerical data points in the right panel of Fig. 7 is either due to the higher-order corrections originating from the MV ensemble (which are outside the precision of the classical approximation employed here) or is caused by the unavoidable discretization errors. Note that the contribution of the diagram B is not seen in our numerical analysis, which appears to be consistent with the possibility that the contribution of the diagram B is zero, as outlined above near Eq. (61).  (88) inspired by the analytical result (77). The prominent peaks at ∆φ/π = 0, 1 and 2 illustrate the HBT-type correlation from the diagram C. The peaks disappear for |k 1 | = |k 2 | (not shown).

V. CONCLUSIONS AND OUTLOOK
In this paper we have demonstrated analytically that the classical gluons fields of the saturation/CGC approach to heavy ion collisions do generate odd azimuthal harmonics in the two-gluon correlation function. Since the classical fields are the leading-order contribution to the two-gluon production cross section, we conclude that the odd azimuthal harmonics are an inherent property of particle production in the saturation framework. This conclusion is consistent with the results of numerical simulations for the classical gluon fields of two colliding heavy ions carried out in [51,52]. The difficulty in identifying the odd-harmonics contribution analytically is related to the fact the analytic expressions for the classical single-and double-gluon production cross sections in heavy ion collision do not exist: instead, as explained in the Introduction, to obtain analytic results one has to assume that one of the nuclei is dilute and expand in the interactions with this dilute projectile order-by-order. As we have shown in this work, odd azimuthal harmonics appear only in the terms contributing at least three interactions with the projectile to the two-gluon production cross section. The part of this term giving the odd harmonics was found above and is given by Eq. (17).
We evaluated this odd-harmonics contribution to the two-gluon production cross section analytically in the GBW model by expanding the interaction with the target to the lowest non-trivial order (six gluon exchanges). The result is given in Eq. (77) and is non-zero: hence the classical gluon fields do generate odd harmonics. In Sec. IV we evaluate the same odd-harmonics correlation function numerically in the full MV model, taking into account the full interaction with the target. The resulting odd harmonics coefficients are plotted in Fig. 6 while the correlation function is given by Fig. 7.
Both the analytic expression (77) and the correlation function in Fig. 7 appear to be dominated by the δ-functionlike peaks at ∆φ = 0 and ∆φ = π. These peaks result from the so-called gluon HBT diagrams in the notation of [42]. We believe the dominance of these peaks is responsible for the v 3 and v 5 in Fig. 6 being so close to each other in their values. To see this imagine a toy two-particle distribution given by with some coefficients A and B. Clearly the expectation values of cos(n ∆φ) averaged with the distribution (89) are independent of n for even and odd n separately, implying that all v 2n+1 are equal to each other, and all the v 2n are equal to each other (but different from v 2n+1 ). Something similar happens in Fig. 6 due to the dominance of the δ-function contribution to the correlator. Certainly, in the actual collisions, fragmentation will turn gluons into hadrons, in the process broadening these δ-function peaks: while a detailed investigation of the fragmentation effects on v n 's is left for further work, one could hope that part of the toy model mechanism suggested here remains, contributing to the similarity of all v 2n+1 and of all v 2n coefficients. (The importance of the ∆φ = 0 and ∆φ = π peaks for v n values was also stressed in [98].) Let us stress further that the gluon correlation function in Fig. 7 should not be compared directly to the data on the di-hadron correlation function. At the very least we expect the fragmentation functions to modify the shape of the correlator, most probably broadening the δ-function peaks. Note also that even in the numerical part of this work, we expand the interaction with the projectile nucleus to the lowest order needed for odd harmonics: higher-order interactions with the projectile may need to be included for the comparison with the experimental data. Unfortunately, at this point, this appears possible to do only numerically. In addition, further theoretical work should include the small-x evolution effects [81][82][83][84][85][86][87][88] into the two-gluon production cross section. Only after all of the above effects are included can one try comparing the resulting correlation function to the experimental data. Finally, we would like to point out that while our MV-model power counting in the calculations presented in this paper assumed a collision of two nuclei (with one of them being more dilute than the other, A 1 A 2 ), the results of our calculations can also be applied to describe the data on di-hadron correlations reported in high-multiplicity pp and pA collisions (with all the caveats listed above) if the saturation scales of both the target and the projectile are perturbatively large, with the target saturation scale being larger than the projectile one.