Collectivity from interference

In hadronic collisions, interference between different production channels affects momentum distributions of multi-particle final states. As this QCD interference does not depend on the strong coupling constant, it is part of the no-interaction baseline that needs to be controlled prior to searching for other manifestations of collective dynamics. Here, we introduce a model that is based on the QCD theory of multi-parton interactions and that allows one to study interference effects in the production of $m$ particles in hadronic collisions with $N$ parton-parton interactions ("sources"). In an expansion in powers of $1/(N_c^2-1)$ and to leading order in the number of sources $N$, we calculate interference effects in the $m$-particle spectra and we determine from them the second and fourth order cumulant momentum anisotropies $v_n$. Without invoking any azimuthal asymmetry and any density dependent non-linear dynamics in the incoming state, and without invoking any interaction in the final state, we find that QCD interference alone can give rise to values for $v_n\lbrace 2\rbrace$ and $v_n\lbrace 4\rbrace$, $n$ even, that persist unattenuated for increasing number of sources, that may increase with increasing multiplicity and that agree with measurements in proton-proton (pp) collisions in terms of the order of magnitude of the signal and the approximate shape of the transverse momentum dependence. We further find that the non-abelian features of QCD interference can give rise to odd harmonic anisotropies. These findings indicate that the no-interaction baseline including QCD interference effects can make a sizeable if not dominant contribution to the measured $v_n$ coefficients in pp collisions. Prospects for analyzing QCD interference contributions further and their possible relevance for proton-nucleus and nucleus-nucleus collisions are discussed shortly.


Contents
1 Introduction 2 2 A model of multi-particle production 3 2.1 Defining the model 3 2.2 Azimuthal multi-particle correlations 6 3 The dipole interference term 7 3.1 Explicit calculation of a simple example: N = 2, m = 2 7 3.2 Dipole interference term for arbitrary m > 2 gluons from N > 2 sources 10 3.2.1 The color correction factor F (2) corr (N, m)      5) corr (N, m), F (6) corr (N, m) 45 1 Introduction Multi-particle production in proton-proton (pp) collisions is typically modeled in terms of multiple parton-parton interactions without invoking explicitly density-dependent dynamics in the incoming wave functions or final state rescattering of the outgoing partons. In particular, multi-purpose event generators provide a reasonable modeling of many characteristics of the underlying event in proton-proton collisions [1][2][3][4], but the simulation of effects that relate different parton-parton interactions is largely limited to ensuring consequences of global conservation laws (energy, momentum, color). The standard picture of multi-particle production in ultra-relativistic nucleus-nucleus (AA) collisions is radically different. Here, jet quenching provides unambiguous evidence for significant final state rescattering effects [5][6][7]. Rescattering is a precursor of fluid dynamics. Partonic systems in which rescattering is operational can be described by an effective kinetic theory that is known to hydrodynamize rapidly [8]. Indeed, fluid dynamical modeling has been demonstrated to provide a phenomenologically valid basis for the simulation of soft multi-particle production in heavy ion collisions [9]. The different dynamical pictures of multi-particle production in pp, pA and AA may be mutually compatible. The transverse size of the systems produced in pp collisions may be sufficiently small for rescattering effects to be negligible, while pA and AA collisions may be sufficiently large and dense to be dominated by multiple rescattering in the final state. However, the recent observation of heavy-ion like behavior in pp (and pA) collisions at the LHC challenges this simple interpretation. On the one hand, the observation of a strong multiplicity-dependence of (multi-)strange hadron production in pp collisions [10] and of momentum anisotropies in pp and pA collisions [11][12][13][14] seems incompatible with modeling such collisions as an essentially incoherent superposition of multiple partonic interactions supplemented by global constraints (see, e.g., Refs. [15][16][17][18] for attempts to model these phenomena). On the other hand, the apparent absence of rescattering effects in inclusive jet and hadron production (above p T ∼ O(1 GeV)) in pp and pA [19][20][21] raises the question whether final state rescattering is sufficiently effective in the smaller collision systems to give rise to measurable signs of collectivity.
This prompts us to ask whether physical phenomena could be at work that contribute to the recent observations of heavy-ion like behavior in pp collisions without invoking final state rescattering or density dependent dynamics in the incoming state. Our focus will be on QCD interference effects, as these do not depend on the coupling constant α s or on interaction probability while they are known to affect multi-particle distributions in the final state. We are mainly interested in understanding their contribution to the anisotropy coefficients v m {2n} that are measured in pp, pPb and PbPb collisions from connected (2n)-particle correlation functions via the so-called cumulant technique [22][23][24]. QCD interference is known to lead to momentum anisotropies in the 2-particle cumulant v 2 {2} (for a rederivation, see section 3 below). However, also the measured higher order cumulants v m {2n} show the same p T -, rapidity-and multiplicity-dependent sizeable values (that are typically 20% smaller than v m {2}). This is commonly refered to as a signature of collectivity [11,12], since it is consistent with a correlation amongst all particles in the event.
The technical question that we shall address with explicit calculations in this manuscript is whether QCD interference can give rise to non-vanishing higher order cumulants v m {2n} and how these are expected to scale with system size.
Under the assumption that hadronic wave functions at ultra-relativistic energies carry saturated gluon distributions, multi-particle correlations have been calculated in the socalled CGC-formalism [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]. This formalism combines effects that are non-negotiably at work (i.e., QCD interference) with effects of a saturated gluon distribution that are searched for as signatures of QCD in a novel high-density regime. In contrast, we work in a simplified model that treats QCD interference exactly but that does not invoke parton saturation effects. This may ultimately help to disentangle both classes of effects. In section 2, we define a QCD-inspired model of multi-particle production that is sufficiently simple to allow for explicit calculations of higher order cumulants. In sections 3 and 4, we calculate v 2 from 2-, 4-and 6-particle cumulants, before summarizing our preliminary analysis of odd harmonics in section 5. Section 6 discusses how the model defined in section 2 is related to the theory of multi-parton interactions. This allows us to estimate the value of the only model parameter, which we use in section 7 to obtain some numerical results. We conclude by summarizing the main conclusions as well as important open questions.
2 A model of multi-particle production

Defining the model
The model for multi-particle production introduced here views a hadronic collision as an event consisting of N parton-parton interactions occurring at positions y i , i ∈ [1, N ], in the transverse plane. To each of the transverse positions y i , the model associates a partonic line source which may be thought of as starting with initial color b i at the rapidity of the first colliding hadron, emitting gluons in the intermediate rapidity window and ending at the rapidity of the second hadron with final color c i . Each multi-particle production amplitude is therefore of the type given in Fig. 1. The model can be summarized as follows: 1. Each hadron collision is characterized by a set {y i , b i }, i ∈ [1, N ], of N particle emitting sources distributed at transverse positions y i with initial colors b i in the adjoint representation.
2. Gluon emission from a source at position y j and color b j is described by an eikonal vertex, where the integration variable x is two-dimensional transverse, and the color structure of the vertex is defined by the adjoint generators T a of SU (N c ). The vertex function f is a two-dimensional vector in the transverse plane, that in cross sections will appear dotted into another vertex function. For instance, for gluons in the non-abelian Coulomb field of an incoming source, one may write f (k) ∝ gk/k 2 . In the following calculations, however, we do not assume a specific functional shape of f (k). The vector f (k) parametrizes then the k-dependent microscopic dynamics that gives rise to gluon emission.
3. When calculating cross sections of event samples, the initial data are weighted with a classical probability distribution ρ ({y i }). Denoting coordinates in the complex conjugate amplitude with primes, this means that initial data According to the model defined above, the spectrum for emission of m particles of transverse momenta k 1 , ..., k m from N sources takes the form Here,σ({k j }, {y i }) denotes the spectrum for the production of m particles of momentum Phenomenologically relevant values for m and N may be fixed by noting that high-multiplicity proton-proton collisions at the LHC can contain m ∼ O(100) particles, and events of this multiplicity are modeled in Monte Carlo event generators typically with N ∼ O(10) parton-parton interactions. However, the main focus of the present work is not on this phenomenologically relevant parameter range but on the qualitative question of whether QCD interference can give rise to momentum anisotropies v n that persist in higher order cumulants. To this end, the main aim of this manuscript is to calculateσ ({k j }, {y i }) for arbitrary values of m and N , and to analyze in particular the limit of large N in which possible asymmetries due to fluctuations in the number of sources are absent. We do this with the following simplifications: where the indices ± denote components of light-cone coordinates and momenta. For high collision energy, however, when both the emitting sources and the emitted gluons propagate close to the light cone, one has k − ≈ 0. This implies that e i k − x + ≈ 1. If the remaining phase e i k + x − were included in the following calculations of gluon production cross sections, it would result in an additional multiplicative factor e i k + (y i −y j ) − in those terms in which a factor e i k.(y i −y j ) occurs. Here, y i , y j denote generic positions of sources from which the gluon of momentum k is emitted in the amplitude and absorbed in the complex conjugate amplitude, respectively. However, identifying the particle emitting source with an energetic parton of light cone momentum fraction p + i , it follows from the uncertainty relation that We therefore conclude that longitudinal dynamics can be neglected when discussing phase interference. After a first illustrative calculation, we shall explain at the end of section 3.1 why gluons are correlated in transverse momentum even if they are separated by a significant rapidity interval.

Emitted gluons do not cross.
In several simpler examples, it was demonstrated explicitly that contributions to the multi-gluon cross section of maximal power in ln(1/x) arise in light cone gauge from ladder diagrams in which emissions are strongly ordered in rapidity, see e.g. [40]. It is not known how these arguments extend to the more complex problem of radiation of many soft gluons from multiple sources discussed here. However, our aim is to devise a model that retains relevant features of QCD but that is simple enough to allow for the explicit calculation of soft multi-gluon interference for large m and N . Motivated by the above-mentioned results for multi-gluon cross sections in simpler systems and by the need for computational simplicity, we therefore assume that multigluon radiation is dominated by ladder-type diagrams in which gluon lines do not cross, and we think of the emitted gluons as ordered in rapidity.
3. m-particle emission cross sections will be symmetrized amongst the m emittees. We shall find that interference contributions to multi-particle emission cross sections are not always symmetric under interchange between final state momenta k i . This is so, since the color constraints on gluon emission of the first (i = 1, 2...) and last (i = .., m−1, m) gluons in the emission amplitude are different from those in between, see appendices A and B for technical details. As these differences are small and unimportant for our discussion, but since they lead to much longer expressions for higher order cumulants, we shall often randomize final results by averaging over all permutations s of the m outgoing momenta, .

No modelling of hadronization
Throughout this work, we calculate partonic spectra and momentum correlations. If hadronization would satisfy local parton-hadron duality (LPHD), then our result could be compared to measured hadron spectra and correlations. However, the simple LPHD prescription may not be phenomenologically viable for multi-particle correlations at soft transverse momentum. We regard it as a limitation of this work that we do not address uncertainties arising from the hadronization stage. We emphasize that our main focus is on addressing the qualitative question of whether QCD interference can give rise to momentum anisotropies v n that persist in higher order cumulants. Since any valid hadronization prescription conserves momentum flow in azimuth, the qualitative answer to this question should not depend on details of the hadronization model, and realistic hadronization models are expected to preserve the order of magnitude of the azimuthal asymmetries found on the partonic level. However, hadronic particle spectra and correlations are generally softened and smeared compared to their partonic parents. This is in particular a caveat for the interpretation of the transverse momentum dependencies of azimuthal anisotropy coefficient v n discussed in section 7.

Azimuthal multi-particle correlations
To define two-particle correlation functions, we average the m-particle emission spectrum in (2.2) with a phase factor e in(φ 1 −φ 2 ) , where k i , φ i denote the radial and azimuthal components of the two-dimensional transverse momenta k i . We also construct the corresponding norm where the Binomial coefficient m 2 counts the number of particle pairs in an event. The integration ρ .. ≡ N i=1 dy i ρ({y i }).. amounts to an average for a specific event sample that is defined by the source distribution ρ. The angular two-particle correlation function e in(φ 1 −φ 2 ) is then defined as [23,24] The experimentally measured (second-order cumulant) anisotropy coefficients v 2 n {2} can be identified with (2.6), We follow experimental practice by defining However, (2.7) cannot be expected to factorize, and in general e in(φ 1 −φ 2 ) (k 1 , k 2 ) = v n {2}(k 1 ) v n {2}(k 2 ) (see sections 4.1 and 7). Correlation functions for more than 2 particles can be defined analogously [23,24]. In particular, we shall calculate the normalized azimuthal 4-particle correlation functions (2.9) and from the corresponding normalization S obtained by evaluating (2.9) without phase factors. The fourth order cumulants are then defined in the standard way, which defines the fourth order cumulant anisotropy coefficient 3 The dipole interference term To illustrate the calculation of the spectrumσ in (2.2), we discuss now the case of emitting m = 2 gluons from N = 2 sources. Fig. 2 shows the 16 diagrammatic contributions. The first row of Fig. 2 shows diagonal contributions in which all gluons are emitted from the same source in the amplitude and in the complex conjugate amplitude. These contributions are free of interference effects. Denoting by a and b the colors of the two emitted gluons, we find for the top left and top right diagrams the color factor (we work in the adjoint representation) The second and third diagrams on the top row of Fig. 2 have the color factor Tr [T a T a ] Tr , so that all diagrams on the top row of Fig. 2 have the same color factor (3.1). The second row of Fig. 2 shows the four diagrammatic contributions for which both emitted gluons are off-diagonal, i.e., they are emitted from one source in the amplitude and they are absorbed by the other source in the complex conjugate amplitude. The corresponding color factor reads  Compared to the contribution (3.1) from diagonal gluon exchanges, they are O 1/(N 2 c − 1)suppressed.
The third and fourth row of Fig. 2 shows contributions with one diagonal and one off-diagonal gluon exchange. In contrast to a QED emission, these vanish in QCD since the color trace of one of the two source lines is ∝ Tr [T a ] or ∝ Tr T b .
The emission vertices (2.1) carry positive (negative) phases ∝ e ik.y i (∝ e −ik.y i ) in the amplitude (complex conjugate amplitude). The squared amplitude can then be written Here, the leading factor 4 counts the four diagrams in the first row of Fig. 2, for which all gluon emissions are diagonal and thus all phases cancel. The four O 1/(N 2 c − 1)suppressed phases in (3.3) correspond to the four diagrams in the second row of Fig. 2 (They are written in the same order in which they arise in the figure.) Here and in the following, we do not specify the normalization since it drops out of the correlation functions that we are interested in.
The two-gluon emission spectrum (2.2) can then be calculated for any given probability distribution ρ ({y i }) of sources. In particular, for a Gaussian ansatz that may be regarded as characterizing a collision at vanishing impact parameter for which the localization of sources does not have a statistically preferred azimuthal orientation, one finds after averaging over the relative distance ∆y ≡ y 1 − y 2 , This is consistent with results obtained in [25][26][27]. The simple model defined in section 2.1 thus shares important commonalities with other approaches. 1 Eq. (3.5) describes QCD dipole radiation. As the average distance ∆y ∝ √ B between the legs of the dipole increases, interference effects decrease and the second term in (3.5) becomes less important. It is a characteristic feature of the two-gluon emission spectrum (3.5) that interference effects enhance the emission equally for gluon pairs that are close in momentum space, k 1 ≈ k 2 and for those that are recoiling against each other, k 1 ≈ −k 2 , [26]. Also, it has been noted repeatedly that spectra like (3.5) are symmetric with respect to k i → −k i , so that they cannot give rise to odd harmonics [27]. The area B may be interpreted in terms of the inverse saturation scale 1/Q 2 s of a saturated parton density [25][26][27]. As we discuss in section 6, the ansatz (3.4) and a parameter range for B can also be motivated within the theory of multi-parton interactions.
We comment at this point on the physical interpretation of the two-gluon correlation in (3.5). This correlation arises from QCD interference of different production amplitudes but it should not be regarded as being the consequence of interference between the two gluons. Against the latter interpretation speaks the finding that an enhancement due to QCD interference is observed not only when the gluons sit close together in transverse momentum space (the term e −B(k 1 −k 2 ) 2 in (3.5)), but also when they are recoiling against each other (the term e −B(k 1 +k 2 ) 2 in (3.5) ). Rather, the interference pattern in (3.5) is consistent with the picture that an azimuthal asymmetry in gluon emission arises for each gluon individually from the interference of the production amplitudes in which the gluon is linked to the first and the second source, respectively. What correlates the orientation of the two emitted gluons in azimuth is not their mutual interference but the fact that both are emitted from the same source pair. Since emission from a source dipole is symmetric with respect to the plane orthogonal to the dipole orientation, each gluon has the same propensity for ending up on the left or right hand side of that plane, and the probability of both gluons ending up in the same hemisphere (the term e −B(k 1 −k 2 ) 2 in (3.5)) or in opposite ones (the term e −B(k 1 −k 2 ) 2 in (3.5)) is therefore equal.
The above argument has noteworthy consequences beyond the simple example discussed in this subsection. First, for arbitrary gluon multiplicity m and arbitrary number of sources N , the diagrams with exactly 2 off-diagonal and m − 2 diagonal gluons are of particular interest since they determine the full O(1/(N 2 c − 1)) contribution to leading order in N (see next subsection). It follows from the color traces involved that these diagrams are only non-vanishing if both off-diagonal gluons connect to the same source pair. As a consequence, the above line of argument carries over to this more general case that we discuss in the next subsection. Second, the above discussion shows explicitly that gluons need not be close to each other in transverse momentum space to be correlated, since they are correlated via common sources. The analogous argument applies to the rapidity dependence. As long as the two sources used for the calculation (3.5) are eikonal and therefore emit gluons from the same transverse positions in different rapidity windows, the two gluons will be correlated in transverse momentum due to the dipole orientation of their common source pair and irrespective of their rapidity difference. This line of argument extends to all emission patterns studied in the present paper. We therefore expect that the omission of explicit rapidity dependencies in our model calculation does not change our conclusions qualitatively.

Dipole interference term for arbitrary m > 2 gluons from N > 2 sources
The calculation of the full interference pattern for emission of a large number m of gluons from a large number N of sources is difficult. Here, we consider first the simpler problem of calculating only the dipole interference terms that include m−2 diagonal and 2 off-diagonal One specific contribution to the second term of the emission cross section (3.6) in which two off-diagonal gluons (curly lines) are supplemented by some diagonal gluons (zagged lines).
The color trace associated to this particular contribution is 2 . As argued in the text, each diagonal gluon that is sandwiched between the two off-diagonal ones and that links to the same sources as the off-diagonal ones (here, the gluons c 1 and c j ) results in a reduction of the color factor by 1/2. Therefore, diagonal gluons are not superimposed incoherently to the dipole interference pattern.
gluons. This is the leading There are N m possibilities of emitting incoherently m diagonal gluons from N sources. The color trace of an incoherent gluon emission gives one factor N c for each of the m diagonal gluons and one factor (N 2 c − 1) for each of the N sources. This explains the prefactor (N 2 c − 1) N N m c N m of the leading term in (3.6). For the subleading term in (3.6), the factor N m−2 accounts for the number of choices of connecting m−2 diagonal gluons to N sources. The sum (ij) goes over the N (N −1)/2 pairs of sources; the second term in (3.6) is therefore of the same order O (N m ) as the first one. The dipole interference of two off-diagonal gluons suppresses this term by a factor 1/(N 2 c − 1) compared to the leading one, as explained in section 3.1.

The color correction factor F
corr (N, m) To understand the factor F (2) corr (N, m) in (3.6), consider the color trace for a gluon emission diagram with one pair of off-diagonal gluons (colors a, b in Fig. 3) and with an arbitrary number of diagonal gluons. The following can be checked to be generally true: If a diagonal gluon is not connected to a source to which an off-diagonal gluon connects (color e in Fig. 3), or if it is connected to such a source but is not sandwiched between the two off-diagonal gluons (color d in Fig. 3), then the generators associated to this gluon emission stand always next to each other in some color trace and they simplify thus according to the color identity In contrast, for those diagonal gluons that are sandwiched between the off-diagonal ones and that are connected to the same sources as the off-diagonal ones (colors c 1 , c j in Fig. 3), the generators in one of the color traces can always be brought into a form where they sandwich one of the generators of an off-diagonal gluon, If one would ignore this correction factor 1/2 l , one would assume that all m−2 diagonal gluons are incoherently superimposed to the interference pattern of the two off-diagonal gluons. The number of such incoherent superpositions is Taking into account that the sum (ab) in (3.6) goes over m(m − 1)/2 choices of selecting a pair of off-diagonal gluons from the ordered list, the factor N incoh accounts for the N -and m-dependence of the prefactor of the O 1/(N 2 c − 1) -suppressed term in (3.6) if F (2) corr (N, m) equals unity.
However, each of the l diagonal gluons that are sandwiched between the off-diagonal ones comes with an extra factor 1/2 that corrects N incoh , and therefore One finds F (2) corr (N = 2, m = 2) = 1, but in general, F corr (3, 3) = 8/9 and F (2) corr (4, 4) = 27/32 are consistent with cases explicitly calculated in the appendices A and B. We note the following limiting cases: Increasing the number of sources at fixed multiplicity m favors incoherent particle production and hence lim N →∞ The same limiting value is reached for other color correction factors F ( * ) corr (N, m) that we encounter in the next section.
2. The limit m → ∞ for fixed average multiplicity per source m = m/N . This limit is consistent with analyses of LHC pp data which indicate that the multiplicity of hard processes is proportional to the soft multiplicity [41]. The color correction factors F ( * ) corr (N, m) are generally finite in this limit. In particular, For m = 1, the correction factor (3.12) is 2/e ≈ 0.73, but for m = 3, it is 0.46 and for m = 5, it equals 0.32. So, higher event multiplicity per source leads to decorrelation that reduces the interference term in (3.6).
3. The high-multiplicity limit m → ∞ for fixed N . For fixed number of sources, the color correction factor behaves asymptotically like Therefore, increasing multiplicity for a fixed number of sources leads to decorrelation.

The second order cumulant
(3.14) Here, the Bessel functions J 2 arise from the φ-integration 2π 0 dφ a e i2φa cos k a .(y i − y j ) in (2.5) and we use ∆y ij ≡ |y i − y j |. The notational shorthand ρ stands for the averaging over source distributions as defined in (2.2), ρ ... ≡ N i=1 dy i ρ({y i }).... The sum (ab) over the m(m − 1)/2 possible pairs of off-diagonal gluon momenta drops out in calculating the average e in(φ 1 −φ 2 ) in (2.7). We highlight three observations: 1. For any multiplicity m, v 2 2 {2} is finite in the limit N → ∞ of a large number of sources.
Since the sum (ij) goes over N (N − 1)/2 source pairs in (3.14), the two-particle cumulant v 2 2 {2}(k 1 , k 2 ) approaches a finite value for N → ∞. Physically, this is so since the observable v 2 2 {2}(k 1 , k 2 ) sums over all source pairs, and since each source pair contributes with a 1/(N 2 c −1)-suppressed contribution to two-particle interference terms. The signal strength v 2 2 {2}(k 1 , k 2 ) does not decrease with the number of dipoles (or the multiplicity in the event) although each dipole is oriented in a statistically independent direction.
does not factorize except for small transverse momentum. In general, due to the source average ρ , one has v 2 . For a Gaussian source distribution (3.4), however, factorization holds for soft transverse momenta to leading order in Bk 2 1 and 3. For any finite number of sources N , v 2 2 {2}(k 1 , k 2 ) vanishes in the high-multiplicity limit. This is a direct consequence of (3.13). Based on intuition from QED, one may have expected that maximal azimuthal correlation arises if all gluons are emitted from the same color dipole (i.e., N = 2). This is not the case. Emitting a large number m of gluons from a large number of sources N can yield a larger signal v 2 2 {2}(k 1 , k 2 ) than emission from a small number of sources, since the color between off-diagonal gluons is less decorrelated.
4 Beyond 2nd order cumulants: results for leading O(N ) and up to sub- We extend now the calculations of section 3 to the fourth order cumulant v 2 {4}. To this end, we have calculated the m-gluon emission cross section from N sources up to ×2 3 cos (k a .∆y lm ) cos (k b .∆y mn ) cos (k c .∆y nl ) +F (6) corr (N, m) ×2 2 cos (k c .∆y no ) cos (k d .∆y no ) 2 2 cos (k e .∆y pq ) cos (k f .∆y pq ) This expression contains color correction factors F ( * ) corr (N, m) that we determine in appendix C in close analogy to the derivation given in section 3.2.1 for F (2) corr (N, m). In equation (4.1), only those contributions are written that are O(N m ) and thus leading in the number of sources. To subleading order in the number of sources, there is a large number of additional interference diagrams. For instance, one can have four off-diagonal gluons emitted from one single source pair, and this contribution is O(N −2 ) suppressed compared to the terms given in (4.1). For the cases N = m = 3 and N = m = 4, we have determined all these contributions explicitly in appendices A and B. These appendices provide combinatorical and calculational details on how to determine (4.1). However, we have neither calculated nor classified completely the subleading contributions O(N −1 ) for an arbitrary number of sources N . For technical reasons, we therefore limit the present discussion to leading O(N ), which amounts to calculating v 2 {4} to O(N 0 ).
We are particularly interested in the question to what extent interference effects could give rise to asymmetries in the final state momentum distributions even if there are no asymmetries in the initial source distributions. We therefore specialize to a factorized ansatz of the N -source distribution in terms of a product of azimuthally symmetric singlesource probabilities, This ansatz includes factorizing Gaussian source models at vanishing impact parameter b, , that will be motivated further in section 6.
4.1 4-particle cumulant to O(N 0 ) and O 1/(N 2 c − 1) 2 Restricting our discussion to the second harmonics, we define the normalized 4-point correlation function in terms of equations (2.9) and (4.1), Here, we discuss this expression first to leading order O 1/(N 2 c − 1) 2 , when only one term of (4.1) contributes, Here, the angles φ lm , φ no denote the azimuthal orientations of the dipoles (lm) and (no). For factorizing distributions of the type (4.2), different dipoles are not correlated in angular orientation, and the source average over the phase e i4(φ lm −φno) in the first term of (4.4) vanishes. With the help of the two-point function (3.12), the fourth order cumulant (2.10) can then be written in the following compact form To order O N 0 , this expression vanishes. In more detail: 1. In the limit m = const., N → ∞, all color correction factors in (4.1) become trivial, 2. In the limit m → ∞, for fixed average multiplicity per source m = m/N , As a consequence, (4.6) vanishes in both limits and we find for the 4-particle cumulant to For the third limit discussed in section 3.2.1 (m → ∞ for constant N ), subleading terms in N would have to be kept. The present calculation therefore does not give access to this limit.
4.2 4-particle cumulant to O(N 0 ) and O 1/(N 2 c − 1) 3 We extend now the calculation of the 4-particle cumulant (4.3) to order O 1/(N 2 c − 1) 3 . All terms in (4.1) contribute to this order. In principle, the azimuthal integrations in (4.3) can be done analytically, and the result can be expressed in terms of source averages over Bessel functions. The resulting expression are straightforward to obtain but they are lengthy. They simplify significantly if one assumes Gaussian source distributions and if one limits the analysis to small transverse momenta B k 2 i 1. For part of the following discussion, we resort to this approximation in which the discussion of qualitative properties becomes more transparent.
To (2.10), we need to calculate the two-point correlation function to O 1/(N 2 c − 1) 2 , and for this we need to calculate the norm For the Gaussian source distribution (4.10) and to lowest order in transverse momenta, we find Here, the hash # denotes prefactors that are common to T and T and that will therefore drop out in the calculation of e i2(φ 1 −φ 2 ) c . (Essentially, the hash stands for the factors written in the first line of (4.1).) The numerator of (2.6) takes the form To leading order O 1/(N 2 c − 1) , the normalized azimuthal two-particle correlation function reduces then to (3.15), but there are higher order corrections The m-gluon emission cross section (4.1) includes contributions that involve interference between 3 or 4 off-diagonal gluons and that enter (4.12). For instance, the term proportional to F corr in (4.12) includes a sum (abc) over m!/(3!(m − 3)!) triplets of off-diagonal gluons. In the contribution of this term to e i2(φ 1 −φ 2 ) (k 1 , k 2 ), only those terms in (abc) survive for which two of the three gluons match the phases. As a consequence, the sum reduces to (12c) , which is a sum over (m − 2) terms. This is the reason for the factor F corr ) 2 comes from expanding the normalization 1/T to O 1/(N 2 c − 1) . Analogously, one obtains the phase factor (2.9) which determines the numerator of the 4-particle correlation function e i2(φ 1 +φ 2 −φ 3 −φ 4 ) , One finds, in close analogy to (4.11), These expressions define e i2(φ 1 +φ 2 −φ 3 −φ 4 ) according to (4.3). The resulting connected 4-particle correlation function (2.10) reads We discuss now limiting cases of this expression.

4.2.1
The limit N → ∞ for constant multiplicity m to order O(1/(N 2 c − 1) 3 ) Following eq. (4.7), all color correction factors reduce to unity in this limit, and the connected 4-point correlation function (4.16) reads According to equation (2.11), the 4-th order cumulant defines a real-valued 4-th order anisotropy coefficient v 2 {4} only if it is negative. Remarkably, this condition is satisfied for sufficiently large multiplicity m, since (4.17) turns negative for m > 3.7. The m 2 -term in (4.17) will dominate for relatively small multiplicities already (say for m > 5), and the 4-th order cumulant reads then We conclude that without any azimuthal asymmetry in the initial state [see eq. (4.10)] and without any (coupling-constant dependent) interaction in the final state, QCD interference can give rise to non-vanishing negative fourth-order cumulants that have an interpretation in terms of azimuthal harmonics v 2 {4}(k). In this sense, our calculation provides a proof of principle that the baseline of vanishing interaction in the final state and of vanishing azimuthal correlation in the initial state does not correspond to a vanishing value v 2 {4}(k). Equation (4.18) provides a proof of principle for the arguments above, but its range of validity is limited as we discuss now. We have calculated v 2 {2}(k) and v 2 {4}(k) to leading order in the number of sources and in the large-N c limit. In particular, the second order cumulant (4.13) contains a leading term ∝ F (2) corr /(N 2 c − 1) and a subleading term A similar observation can be made for the 4-particle correlator (4.16) where the term suppressed by one power 1/(N 2 c − 1) is enhanced by m 2 . This seems to indicate that the expansion in powers of 1/(N 2 c − 1) converges only as long as m 2 < (N 2 c − 1). The radius of convergence may be larger than the above estimate for the following reason: In equations (4.13) and (4.16), terms proportional to m (m 2 ) arise from integrating out the transverse momentum q of one (two) of the off-diagonal gluons involved in an interference term. In the simplest case 3 , this q-integration leads to contributions of the type For technical simplicity, we have worked in this subsection to lowest order in small transverse momentum. This amounts to the assumption that the emission vertex f (q) is dominated by transverse momenta that are much smaller than the inverse transverse size of the source. In this limit, J 0 (qz) ≈ 1 and thus a = 1. In general, however, a ≤ 1. Indeed, depending on the emission vertices f (q) and on the source density ρ(z), a may be significantly smaller than unity, and e.g. the F Integrating out transverse gluon momenta from other contributions in (4.1) can yield more complicated expressions than (4.19), but the maximal value is always obtained in the limit B q 2 1 studied here, and the value starts decreasing when the integral over the transverse momentum extends to values that start resolving the transverse distance between sources. Therefore, rather than being limited to m 2 < (N 2 c − 1), the region of validity of the calculations in section 4 is expected to extend up to higher multiplicities where a 1. Physically, a is a penalty factor for integrating out off-diagonal gluons. It depends on the source size and the emission vertex f .
One may use these relations to simplify the connected 4-particle correlation function (4.16), Since the numerator of this expression starts at O 1/(N 2 c − 1) 3 , the normalization is trivial. One finds from explicit calculation in the limit (Bk 2 i ) 1 The 6-th particle cumulant is defined as 6 c ≡ 6 − 9 2 4 + 12 2 3 (in this shorthand notation, the numbers in brackets denote the number of phases). One therefore

Higher harmonics
In sections 3 and 4, we have focussed on the calculation of the second harmonics e i2(φ 1 −φ 2 ) and e i2(φ 1 +φ 2 −φ 3 −φ 4 ) c . Here, we discuss the calculation of other even and odd harmonics.

Higher even harmonics
When calculating e in(φ 1 −φ 2 ) and e in(φ 1 +φ 2 −φ 3 −φ 4 ) c , one encounters elementary azimuthal integrals of the form dφ e inφ cos (q∆y cos(φ)) . (5.1) Here, q is the modulus of a generic transverse momentum, ∆y is the modulus of a generic transverse dipole separation, and φ is the relativ azimuthal angle between transverse momentum and dipole separation. All cosine-terms in (4.1) can be written as cos (q∆y cos(φ)) for suitable choices of q, ∆y, and φ.
The integrals (5.1) are non-vanishing for all even integers n. For n = 2, one finds for instance dφ e i2φ cos (q∆y cos(φ)) = −2πJ 2 (q∆y) This leads to the Bessel functions J 2 in the calculation of the two-particle correlation (3.14) and the four particle correlation (4.4). Explicit expressions can also be given for higher even harmonics. For instance dφ e i4φ cos (q∆y cos(φ)) = 2π J 4 (q∆y) With the help of these expressions, the analysis of section 4 could be repeated for arbitrary even harmonics. Since the small-q limit of (5.1) is ∝ q n , one concludes immediately that the small-k-behavior of v n (k) is ∝ k n . In particular, the parametric dependence of the fourth harmonic at small q is ∼ (q∆y) 4 while that of the second harmonics is ∼ (q∆y) 2 . Within the approximation of small transverse momenta, ∼ (B q 2 ) 1 explored in section 4.2, this implies that v 4 ∝ v 2 v 2 for second and fourth order cumulants. We discuss this point further in section 7.

Higher odd harmonics
In contrast to the even harmonics, the integral (5.1) vanishes for odd integers n. As a consequence, the odd harmonic flow coefficients v n vanish up to O 1/(N 2 c − 1) 4 and O (1/N ), since the spectrum (4.1) shows only cosine-terms to this order.
Based on the idea that interference patterns are momentum conjugates of spatial distributions, one may naively expect that odd harmonics make some (possibly subleading) non-vanishing contribution whenever the spatial distribution shows odd harmonic eccentricities, i.e., for N ≥ 3 sources. Motivated by this idea, we have calculated in appendix A all terms contributing to N = m = 3 in search of odd harmonics. However, the emission spectrum for N = m = 3 turned out to be free of odd harmonics.

Odd harmonics for the case N = m = 4
We have classified and calculated in appendix B all contributions to the emission spectrum for N = m = 4. Two classes of diagrams were found to lead to odd harmonics, see eqs. (B.13) and (B.15). Referring for technical details and explicit results to the appendix, we limit the discussion here in the main text to provide qualitative insight into how properties of the SU (N c ) color algebra can give rise to odd harmonics. To this end, we consider in Fig. 4 a set of diagrams with four off-diagonal gluons that are emitted from three sources l, m, n that combine to two pairs of sources (lm), (mn) .
For this contribution, one checks easily that the color trace changes depending on whether the sandwiched gluons of momentum k 2 and k 3 link from m to n or from n to m. For the four diagrams in Fig. 4, one finds

(5.4)
For each diagram that contributes with a phase e −i k j .∆y mn , there is a diagram in which the off-diagonal gluon with momentum k j links from m to n rather than from n to m while all other gluons are linked in the same way. If the prefactors of both these diagrams were the same, then these phases would add to a term e −i k j .∆y mn + e +i k j .∆y mn that is proportional to a cosine, and odd harmonics in k j would not occur. In the example above, it is only the non-abelianess of SU (N c ) that leads to different prefactors of the phases e −i k j .∆y mn and e +i k j .∆y mn , thus giving rise to terms that change sign under k j → −k j , see eq. (5.4). This is true for all odd harmonics that we have found in our calculations.

General comments on odd harmonics
In calculations of multi-particle correlations in the framework of saturation physics [28,[30][31][32][34][35][36], it has been a persistent problem to find non-vanishing values of odd harmonic anisotropy coefficients such as v 3 . Solutions to this problem include the proposal that the proton wavefunction consists of a few patches in the transverse plane in each of which color fields point in preferred direction [33], as well as the recent observation that odd harmonics arise as a high parton density effect directly from the non-linear QCD evolution [29]. The calculation in section 5.2.1 points to a third possible origin of odd harmonic contributions: without any coupling-constant dependent interaction in the initial or final state and without any asymmetry in the initial state, odd harmonic contributions to multi-particle production can arise from non-abelian properties of QCD interference. For the odd harmonic contributions analyzed in section 5.2.1, it was important that the number of off-diagonal gluons was larger than the number of sources to which they were connected. As a consequence, these contributions seem to be O(1/N )-suppressed. We note, however, that we have studied odd harmonic contributions only for N = m = 4. We did not try to determine the color correction factor when further diagonal gluons are added to the diagrams studied here (this would require combinatorical arguments that go beyond those developped in appendix C), and we do not know the full N -and m-dependence of the leading odd harmonic contribution. It remains an interesting open question to find a classification of all diagrams of O(1/N ) and to establish how odd harmonics manifest themselves in the limit of a large number of sources.

Relation to the theory of multi-parton interactions (MPIs)
Starting from ideas in the mid-80s [42,43], the treatment of multiple parton interactions (MPIs) in perturbative QCD has been developed further in recent years [44][45][46][47][48][49][50][51][52]. Here, we give simple arguments that relate this theory to the model defined in section 2 and that motivate values for the source parameter B in pp collisions.
Hadronic cross sections involving N partonic interactions are customarily parametrized as products of N independent parton-parton interactions σ i The physics of MPIs enters here in the coefficient K N that parametrizes deviations from an incoherent superposition of N interactions. K N is dimensionfull ∝ (area) N −1 . For the case of double parton interactions (N = 2), K reduces to the effective cross section σ eff that is constrained experimentally. Under certain mild assumptions (such as neglecting contributions to MPIs from 1 → 2 splittings), K N can be expressed through N -particle Generalized Parton Distributions (GPDs) G N [45] Here f (x, Q 2 ) denotes the standard single parton distribution function. The x i , ∆ i denote the longitudinal momentum fractions and the initial transverse momenta of the i-th parton in the incoming hadronic wave function. Neglecting the weak dependence of this expression on x i and Q i , K N can be expressed in terms of the function For a process involving N MPIs, an m-parton production cross section can then be written formally as where M 2 is the squared amplitude for the production of m gluons from N partons in the nucleon wave function ("N sources"). The corresponding m-particle spectrum is obtained by normalizing this expression with the cross section To arrive at eqs. (6.4) and (6.5), one assumes that m-parton production can be formulated in a pQCD factorized formalism. Here, we do not try to quantify corrections to this assumption. Rather, we treat a bold extrapolation of this perturbative approach to soft momenta as one way of getting insight into soft multi-particle production. To see the commonalities between the model in section 2 and this formalism deduced from MPI theory, we introduce one further approximation by writing the generalized Nparton distribution functions G N in a mean-field approximation [45] as products of generalized one-particle distribution functions G 1 (x, Q, ∆) = f (x, Q) F 2g (∆) [53][54][55] with a two-gluon form factor that parametrizes the transverse momentum distribution Choosing F 2 2g (∆) = exp(−B∆ 2 i ) for simplicity to be of Gaussian form (other functional dependencies could be explored), one has F N (∆ 1 , . . . , ∆ N ) = N i=1 exp(−B∆ 2 i ), and it is straightforward to switch from (6.5) to coordinate space representation 4 . (6.7) 4 Here, the density ρ ({yi}, b) is a convolution of the normalized densities of colliding partons in the two incoming hadrons. For a Gaussian ansatz In other sections, we have used this distribution for vanishing impact parameter b = 0 and normalized to unity, see eq. (4.2). This is the form of the m-particle emission spectrum (2.2) with a source probability distribution ρ of the form (4.2). In section 2, we have supplemented the structure of (6.7) with a particularly simple model for calculatingσ = |M 2 |. Here, we see that this structure arises from MPI theory once one treats N -parton GPDs in mean field approximation. We note that this framework allows one to constrain the parameter B in the source density by data. Namely, 1/K N = db dy i ρ({y i }, b) = 1/N (4πB) N −1 , and therefore σ eff = K 2 = 8πB, B = 2 GeV −2 ←→ σ eff ≈ 20 mb (6.8) Experimentally favored values for σ eff lie in the range of 15 ± 5 mb for pp collisions at the LHC [56][57][58] and for pp collisions at Tevatron [59]. Somewhat higher values σ eff of order 35-40 mb have been obtained in a mean field approach that does not include other mechanisms for MPI enhancement [47,60]. This will prompt us in the next section to scan values in the range 1 GeV −2 < B < 4 GeV −2 .

Numerical results
Whenever there are multiple partons in the final state, QCD interference contributes to the azimuthal anisotropies v n for both even (see sections 3 and 4) and odd (see section 5) harmonics. This raises the question how contributions from QCD interference compare in size and signature to those of other physically conceivable mechanisms. What is the typical signal size with which QCD interference can contribute to v n ? And what is the expected p T -, rapidity-and multiplicity-dependence of the effects discussed here? The present study is not sufficient to provide complete answers to these questions, but we summarize in this section what can be said from parametric considerations and from first numerical results. Parametrically, leading contributions to v n , n even, are O 1/(N 2 c − 1) 1/2 for the second order cumulants and O 1/(N 2 c − 1) 3/4 for the fourth-order cumulants. This N cdependence is multiplied by functions that grow ∝ |k| n for very small transverse momenta but that reach magnitudes of O(1) for sufficiently large transverse momenta, see, e.g., eq. (3.14). For N c = 3, signal sizes v n ∼ 0.1 − 0.3 seem therefore conceivable. The signal size is roughly independent of the number of sources N . While our analysis of the multiplicity dependence of v n in section 4 does not allow us to draw conclusions in the limit of large multiplicity, it points to the possibility that v n rises with multiplicity. Also, the model described in section 2 leads naturally to an approximately flat rapidity dependence. Albeit not established on the quantitative level needed for decisive tests, the above-mentioned features agree at least qualitatively with trends in the data.
It is less clear whether also the transverse momentum dependence of the interference effects studied here can account for the qualitative trends in the data. Amplitudes for the emission from different sources are expected to interfer only for sufficiently small transverse momenta k that do not resolve the separation ∆y between the emitters, k < 1/∆y. It is therefore unclear whether QCD interference can contribute significantly to v n (k) in the multiple GeV range where significant signal strength is observed. To address this question, we plot in Fig. 5 the transverse momentum dependence of the second order cumulant for values of the source parameter B favored by data on multi-parton interactions, see eq. (6.8). The main qualitative features of this second order cumulant are indeed consistent with main trends in the data, the signal strength increases up to transverse momenta of 2 GeV and a sizeable signal persists in the multi-GeV range.
As seen from Fig. 6, the signal size obtained in the present formalism varies linearly with the value of the color correction factor. The color correction factors used in Fig. 6 are for an average ofm = 3 or 5 emitted particles per source, but a more extreme choicem = 20 that may be at the upper end of what is phenomenologically viable would reduce the signal by F (2) corr = 0.1. In comparison to these uncertainties, the O(1/(N 2 c − 1))-corrections in the denominator of (3.14) turn out to be negligible.
We next point to the fact that v 2 2 {2}(k 1 , k 2 ) does not factorize. To this end, we show in the upper panel of Fig. 7 . One sees that in the transverse momentum range up to ∼ 1.5 GeV, deviations from factorization are at most ∼ 20%. For higher transverse momenta, however, particle emission starts to decorrelate as soon as the difference between the transverse momenta of the two emittees is larger than 1 -2 GeV. For the experimentally measured hadronic momentum correlations, such . Lower panel: The corresponding factorization ratio r 2 (k 1 , k 2 ) as defined in equation (7.1), plotted for selected 'trigger' momenta k 1 against the relative momentum difference k 1 − k 2 . deviations from factorization are often characterized in terms of the factorization ratio where one of the two momenta, say k 1 , is regarded as trigger, and the other is the associate particle satisfying k 2 < k 1 . The corresponding quantity is plotted in the lower panel of Fig. 7 as a function of k 1 −k 2 . For pPb and PbPb collisions, the factorization ratio r 2 (k 1 , k 2 ) has been measured at LHC [61] and it is in agreement with fluid dynamic simulations [62], see also Ref. [63] for earlier discussions of the relation of r 2 (k 1 , k 2 ) to fluid dynamics. Fig. 7 shows features qualitatively similar to those reported in Refs. [61,62] in that r 2 decreases with increasing k 1 − k 2 and that this decrease is more pronounced for increasing trigger k 1 . Clearly, the present calculation is for pp collisions, and it shows a correlation on parton level. Hadronization may be expected to affect the correlation r 2 significantly as it tends to smear and soften transverse momentum distributions. Given these substantial differences between the results of Refs. [61,62] and the present calculation, we refrain from a quantitative comparison.
In section 4, we have limited one part of our discussion to the small-k approximation in which many expressions simplify. The left hand side of Fig. 8 indicates the range of validity of this approximation. As discussed in section 5.1, it is a structural property of the small-k approximation that Similar relations between different harmonic coefficients are known to arise in fluid dynamic models as a consequence of mode-mode coupling [64][65][66]. For the specific relation (7.2), it was noted already in Ref. [67] that one expects for large transverse momentum a proportionality constant 1/2. It is therefore interesting to note that non-trivial relations between v 4 and v 2 may also arise from mechanisms that do not invoke interactions in the final state. The right hand side of Fig. 8 illustrates this point further. For small k, v 4 is seen to be approximately proportional to (v 2 ) 2 with a proportionality factor of order unity. So far, the figures shown in this section were calculated to leading order in 1/(N 2 c − 1) from the second order cumulant (3.14). To the next order in 1/(N 2 c − 1), the second order cumulant v 2 {2}(k) receives an additive correction, see eq. (4.13). However, as explained in section 4, this is an expansion in powers of m 2 /(N 2 c − 1) that does not converge for large multiplicities. An analogous statement applies to the 4-th order cumulant v 2 {4}(k) in equation (4.16) for which the contribution subleading in 1/(N 2 c − 1) grows with m 2 . Despite this caveat, we proceed here with a numerical exploration of v 2 {4}(k). To this end, we focus on the term ∝ m 2 in (4.16) that dominates for large multiplicity. Rather than going to the small-k-limit (4.18), we keep the color correction factor in (4.16), and we re-establish the large-k-behavior by replacing Bk 2 with the integral over Bessel functions from which it was obtained in a small-k-approximation, Here, d denotes a suppression factor that arises from integrating out one of the three dipoles in the term ∝ F (6) corr = F (2) corr 3 in (4.1). The factor d equals unity for B k 2 1 but it depends on the vertex function and can be smaller than unity d < 1. For any given vertex function f (k), it can be calculated in close analogy to the suppression factor a in (4.19) and one has d ∼ O(a 2 ). Fig. 9 shows the fourth root of (7.3) for multiplicities m reached in high-multiplicity proton-proton collisions at the LHC, and for a number of sources N comparable to the number of MPIs invoked in MC simulations of the underlying event of such pp collisions. We also vary the suppression factor d over ranges that are easily obtained from (4.19) and d a 2 . While the results in Fig. 9 fall short of a quantitative determination of v 4 2 {4}(k), they support the qualitative statement that the size and shape of the 4-th order cumulant v 2 {4}(k) resulting from QCD interference may be comparable to the size and shape of the second order cumulants v 2 {2}(k) within the parameter range realized in high-multiplicity proton-proton collisions.

Conclusion
Whenever multiple partons are produced in a hadronic collision, there are interference effects whose size varies with the momentum distribution and the colors of the outgoing partons. These QCD interference effects do not depend on α s . In this sense, they are part of a "no-interaction" baseline for v n that needs to be controlled prior to discussing non-linear dynamics in the incoming hadronic wavefunctions or rescattering and fluid dynamization in the final state. 5 Here we have used a simple model of multi-particle production (see section 2) to estimate how these QCD interference effects can contribute to the azimuthal anisotropy coefficients v n measured with second and higher order cumulants in pp collisions. We have pointed out (section 6) that this model can be realized in the mean field theory of multi-parton interactions (MPIs), with MPIs representing the "sources" of individual parton-parton collisions and radiated gluons corresponding to radiation associated with the MPIs.
Our calculations establish that the contribution of QCD interference to the anisotropy coefficients v n , n even, persists unattenuated for an increasing number of sources, that it can increase with increasing multiplicity, and that it persists in higher order cumulants. We have further shown that odd harmonic anisotropy coefficients arise due to the nonabelianness of the QCD interference pattern 6 . In section 7, we have supplemented these structural results about how the QCD interference impacts anisotropy coefficients with a first numerical exploration. Both, the order of magnitude of the calculated v n , as well as the shape of their transverse momentum dependence was found to be of the order of magnitude and shape of the signals observed in pp collisions at the LHC. Given the simple nature of the model in section 2 (schematic treatment of particle production, absence of hadronization, etc.), these qualitative commonalities between model and data must not be over-interpreted. However, they make it certainly conceivable that the no-interaction baseline including QCD interference effects can make a sizeable if not dominant contribution to the measured v n coefficients in pp collisions.
This leads naturally to the question whether and to what extent QCD interference effects could contribute also to the anisotropy coefficients measured in AA collisions. Here, however, marked differences need to be considered. First, jet quenching provides unambiguous experimental evidence for significant final state rescattering which is an α s -dependent microscopic mechanism that underlies hydrodynamization and that can translate spatial gradients into momentum anisotropies. Given the strength of jet quenching signals, it is thus inconceivable that QCD interference can account for the totality of the observed v n signals in AA. Second, final state scattering can destroy coherence and the resulting interference. While QCD interference clearly shifts the no-interaction baseline for v n , it is therefore questionable that effects of interactions contribute additively on top of this baseline. Rather than seeking a common explanation of v n across system size, it may therefore also be instructive to explore the opposite hypothesis, namely that the physics mechanisms underlying the v n -signals in pp and AA are qualitatively different, being dominated by QCD interference for pp while being dominated by fluid dynamics for AA. The similarities in the experimental signal of pp, pA and AA may be viewed as disfavouring this hypothesis, but the qualitatively different evidence for final state rescattering in both systems suggests that different mechanisms are at work in pp and AA and these may therefore also be at the origin of the measured v n . The present work adds to this discussion only by illustrating that QCD interference may account for the order of magnitude and main qualitative features of v n signals in pp. As far as pA collisions are concerned, the absence of experimental evidence for significant final state scattering supports the idea that the qualitative conclusions about QCD interference drawn here for pp carry over to pA.
We have aimed at controlling in this paper interference effects between an arbitrary number m of gluons emitted from an arbitrary number N of sources. This was achieved within an expansion in powers of 1/(N 2 c − 1) and to leading order in N by resumming the effects of an arbitrary number of diagonal gluons in color correction factors. We hope that these results can be useful also for the discussion of other models, such as models based on saturation physics 7 that encompass QCD interference effects and with which our calculation agrees to lowest order, see eq. (3.5). Within these models, the importance of QCD interference effects has been emphasized repeatedly (see e.g. Ref. [39]), but they are not studied in isolation and are difficult to disentangle from effects of finite partonic density. Indeed, due to the greater complexity of calculations in these models, correlation functions have not been analyzed to the same level of detail as in the present manuscript, and we are not aware of similar parametric statements about the leading N -independence and 1/(N 2 c − 1)-dependence of higher order cumulants. Future developments may also better relate the results reported here to recent efforts of improving MC simulations of the underlying event in TeV-scale proton-proton collisions. On the one hand, while all modern multi-purpose MC event generators model the underlying event in terms of multi-parton interactions, the QCD interference effects discussed here are not included in these simulations. On the other hand, there are efforts to go beyond an essentially incoherent superposition of MPIs supplemented with conservation laws, e.g. by modeling effects of overlapping strings and studying whether these could give rise to signatures of collectivity [69][70][71]. For earlier works, see e.g. Ref. [72] and approaches based on pomeron dynamics [73,74]. It would be interesting to understand how QCD interference can be included in MC simulations and how this compares e.g. to effects of overlapping strings or other models.
Finally, within the set-up of this manuscript, the expansion in O 1/(N 2 c − 1) analyzed here provided first qualitative insights, but we also discussed its limitations. In particular, a more systematic control over 1/N -suppressed terms may allow for more quantitative statements in the range of multiplicity and number of sources that are phenomenologically interesting, and it would give access to odd harmonic anisotropy coefficients. Moreover, to gain better control in the phenomenologically interesting range of multiplicities, one would ideally like to resum all contributions that come with powers in m 2 /(N 2 c − 1). We expect that such further advances are possible.
A Multi-gluon emission for N = m = 3 In this appendix, we calculate the cross section for producing m = 3 gluons from N = 3 sources. This will illustrate several statements that we have generalized to arbitrary N and m in the main text.
For N sources and m emitted gluons, there are N 2m different diagrams, since each of the m gluons can be attached to any of the N sources in the amplitude, and to any of the N sources in the complex conjugate amplitude. For N = m = 3, these N 2m = 3 6 diagrammatic contributions can be classified as follows: I 3 diagonal gluons (27 diagrams) 7 In the notation of section 2, the initial transverse density of sources is ∼ N/B. Since the v2{2} and v2{4} calculated in section 4 are N -independent and finite to leading order in 1/N , a high or saturated initial parton density is not a prerequisite for the effects discussed here while it is the basis of calculations in the so-called CGC-formalism. To illustrate how to enumerate these diagrams, consider, e.g., case II: there are 3 choices for the source to which each of the diagonal gluon can be attached. For the off-diagonal gluon, one has 3 choices for the source in the amplitude times 2 choices in the complex conjugate amplitude. As exactly one of the three gluons is off-diagonal, one has also 3 choices for selecting the off-diagonal gluon amongst all three gluons. Combining these factors leads to 3 × 3 × (3 × 2) × 3 = 6 × 27 different diagrams.
We turn now to the calculation of the different contributing cases: For three diagonal gluons radiated off three sources (case I), all 27 diagrams contributing to the 3-gluon radiation cross section have the same color factor (N 2 c − 1) 3 N 3 c . As phase factors cancel for diagonal gluons, the contribution toσ iŝ One also checks easily that the color factors of all diagrams with only 1 off-diagonal gluon (case II) vanish,σ Case III includes 8 × 27 diagrams for which the two off-diagonal gluons connect to three sources; these have vanishing color factors. In addition, there are 4 × 27 diagrams for which the two off-diagonal gluons connect to two sources. For those, we denote with a and b the colors of the off-diagonal gluons and with c the color of the diagonal one. If the diagonal gluon links to the source that is not touched by the off-diagonal gluons (case not shown in Fig. 10), then the resulting color factor is Tr If the diagonal gluon touches instead a source that is also touched by the off-diagonal gluons, but is not sandwiched between the off-diagonal gluons (see, e.g., top row of Fig. 10), then However, for those diagrams for which the diagonal gluon is sandwiched between the offdiagonal ones and for which the diagonal and off-diagonal ones have one source in common (see second row of Fig. 10), one finds The total contribution of the diagrams of case III to the 3-gluon emission cross section takes then the form σ (ord) Here, we use the subscript (ord) to indicate that the momenta k 1 , k 2 and k 3 are ordered from top to bottom in the emission diagrams. The color trace (A.4) appears only for 2 of the 3 diagrams for which the diagonal gluon carries k 2 , and this reduces the prefactor of the second term in (A.5) to 2 × 1/2 + 1 = 2. Following (2.3), we randomize the external momenta to obtain 4 cos (k a .∆y ij )cos (k b .∆y ij ) . (A.6) The prefactor 8/9 in this expression is consistent with the factor F (2) corr (N, m), obtained for N = m = 3 from eq. (3.10).
Case IV concerns diagrams with 3 off-diagonal gluons. Examples for such contributions are depicted in Fig. 11. If all three gluons are linked to the same pair of sources (see, e.g., Fig. 11(a)), we find for the symmetrized expression 8 cos (k 1 .∆y ij )cos (k 2 .∆y ij )cos (k 3 .∆y ij ) .
(A.7) Figure 11. Examples for diagrammatic contributions with 3 off-diagonal gluons that link between one (a), two (b) and three (c) different pairs of sources. Contributions of type (b) vanish while contributions of type (a) and (c) have color traces that differ by a factor (N 2 c − 1)/4, see text for more details.
If the 3 off-diagonal gluons are linked to exactly two pairs of sources (see, e.g., Fig. 11(b)), then one of the color traces vanishes,σ If the three off-diagonal gluons are linked to exactly three pairs of sources (see, e.g., Fig. 11(c)), one finds Parametrically, contributions to (A.9) are suppressed by an extra factor (N 2 c − 1) compared to (A.7). The prefactors of these expressions can be understood by noting that contributions to (A.7) (see Fig. 11(a)) have a color factor while contributions to (A.9) (see Fig. 11(c)) have a color trace The cases I and II are trivial: diagrams with only diagonal gluons are counted with color factor N 4 c (N 2 c −1) 4 and diagrams with exactly one off-diagonal gluon have a vanishing color trace,σ For the three other cases, however, qualitatively novel features arise as we discuss now. Case III includes diagrams with two off-diagonal gluons linked to 2, 3 or 4 sources. Only the first of these possibilities gives non-vanishing contributions. It includes 9×4 4 diagrams. Fig. 12 depicts some of them. One checks easily that these contributions have different color traces, namely if in between the off-diagonal exchanges there is no diagonal gluon exchange (see, e.g., Fig. 12(a)), if one of the two diagonal gluons is exchanged in between the off-diagonal ones (see, e.g., Fig. 12(b)), and if both diagonal gluon exchanges occur in between the off-diagonal ones (see, e.g., Fig. 12(c)). This illustrates the general statements about diagonal gluons sandwiched between offdiagonal ones that are made in section 3.2.1.

This coincides with F
corr (4, 4), see eq. (3.10). Case IV contains diagrams in which three off-diagonal gluons i) are linked to exactly one pair of sources (i, j).
ii) are linked to two pairs of sources (ij), (jl) made of three sources.
iii) are linked to two pairs of sources (ij), (lm) made of four sources. iv) are linked to three pairs of sources (ij), (jl), (li) made of three sources. v) are linked to 3 pairs made of four sources.
These contributions have i) 3 × 4 4 , ii) 9 × 4 5 , iii) 9 × 4 4 , iv) 3 × 4 5 and v) 3 × 4 6 diagramms respectively, that sum up to a total of 3 3 × 4 5 diagrams. The contributions ii), iii) and v) can be checked easily to have vanishing color traces. The contribution IVi) yields (we present symmetrized cross sections according to (2.3)) Similarly, for diagrams contributing to case IViv, we find color traces N 4 c (N 2 c − 1) 2 , as well as 1 2 N 4 c (N 2 c − 1) 3 for some diagrams. After symmetrization, the reduced color trace can be absorbed in a color correction factor 7/8 and one finds For the prefactors of these expressions, there is the following consistency check: One recalls that each diagram contributes phase factors with prefactor 1 times a color trace. Hence, if one ignores the color correction factor and if one multiplies the 4 terms in the sum (abc) times the N (N − 1)(N − 2) terms in (lm)(mn)(nl) times the explicit prefactor N 8, the result should match the number 3 × 4 5 of diagrams contributing to case IViv. Case V includes diagrams with four off-diagonal gluons that are i) emitted from two sources, i.e., one pair of sources (i, j).
ii) emitted from three sources that combine to two pairs of sources (lm), (mn).
iii) emitted from three sources that combine to three pairs of sources (lm), (mn), (n, l).
iv) emitted from four sources that combine to two pairs of sources (lm) and (no). v) emitted from four sources that combine to three pairs of sources (lm), (mn), (no).

C Color correction factors
In this appendix, we derive the color correction factors entering (4.1). The factor F (2) corr was derived in section 3.2.1 already.
We consider the case of m − 4 diagonal gluons and of four off-diagonal gluons that form two dipoles (lm), (no) with four different sources. This is case V iv1 in appendix B, which is the leading O 1/(N 2 c − 1) 2 contribution to the 4-particle correlator. There are m 4 possibilities to assign the 4 off-diagonal gluons to an ordered list of m gluons. For each assignment, there are • j 1 diagonal gluons between the first and second off-diagonal gluon.
• j 2 diagonal gluons between the second and third off-diagonal gluon.
• j 3 diagonal gluons between the third and fourth off-diagonal gluon.
• m − 4 − j 1 − j 2 − j 3 before the first or after the last off-diagonal gluon.
The number of possibilities of distributing these m−4 gluons such that an arbitrary number of j 1 + j 2 + j 3 gluons lies between the first and last off-diagonal gluon is This is a consistency check of our starting point that will be useful in the following. We also note that the total number of different diagrams (for given sources l, m, n, o) is The four off-diagonal gluons can be paired into dipoles in three different ways. Denoting the off-diagonal gluons as 1,1, 2,2 in a notation that makes their pairing clear, the three different paring are A. Ordering of off-diagonal gluons: 1,1, 2,2 B. Ordering of off-diagonal gluons: 1, 2,1,2 C. Ordering of off-diagonal gluons: 1, 2,2,1 For case A., there are j 1 diagonal gluons sandwiched between the first dipole pair 11, and j 3 gluons sandwiched between 22; j 2 gluons are not sandwiched between any dipole pair. l 1 of the j 1 gluons will be linked to the same source pair (lm) as 11, and analogously for l 3 out of the j 3 gluons sandwiched between 22. As each diagonal gluon sandwiched between an off-diagonal pair reduces the color trace and thus the spectrum by a factor 2, we have We note that without the above-mentioned factors 1/2, the terms in (C.3) would be multiplied by factors 2 l 1 +l 3 , and one would find F (4i,A) corr (N, m) = 1. For case B., there are again j 1 diagonal gluons sandwiched between the first dipole pair 11, and j 3 gluons sandwiched between 22. However, there are now j 2 gluons sandwiched between 11, and 22. In this case, there will be l 2 out of j 2 gluons that are linked to the dipole (lm) and there will bel 2 of the j 2 gluons linked to (no). The corresponding color factor reads (C.6) We consider the case of m − 3 diagonal gluons and of three off-diagonal gluons that are linked between the source pairs (lm), (mn) and (nl). In close analogy to appendix C.1, we denote by j 1 (j 2 ) the number of diagonal gluons between the first and second (the second and third) off-diagonal gluon. One checks easily for each of the j 1 diagonal gluons: their color trace is N c if they link to a source other than l, m or n. Also for exactly one of the three sources l, m, n, the color trace is N c , but it is N c /2 for the other two sources.
Paralleling the logic that lead to (C.3), we therefore find the norm and We consider the case of m − 4 diagonal gluons and of four off-diagonal gluons that are linked between the source pairs (lm), (mn), (no) and (ol). We denote by j 1 , j 2 and j 3 the number of diagonal gluons that are sandwiched between the first and second, the second and third and the third and fourth off-diagonal gluon, respectively. One checks easily that each of these diagonal gluons contributes to a color trace with factor N c if linked to any source other than l, m, n or o. Amongst the sources l, m, n and o, there are always exactly two sources where the color trace is N c , while the color trace for the two other sources is (m − 3 − (j 1 + j 2 + j 3 )) N m−4−(j 1 +j 2 +j 3 ) corr (N, m) The complexity of the combinatorics increases if one considers diagonal gluons linked into diagrams with 5 or 6 off-diagonal gluons. For the terms F (5) corr (N, m) and F