Multiparticle collectivity from initial state correlations in high energy proton-nucleus collisions

Qualitative features of multiparticle correlations in light-heavy ion ($p+A$) collisions at RHIC and LHC are reproduced in a simple initial state model of partons in the projectile coherently scattering off localized domains of color charge in the heavy nuclear target. These include i) the ordering of the magnitudes of the azimuthal angle $n$th Fourier harmonics of two-particle correlations $v_n\{2\}$, ii) the energy and transverse momentum dependence of the four-particle Fourier harmonic $v_2\{4\}$, and iii) the energy dependence of four-particle symmetric cumulants measuring correlations between different Fourier harmonics. Similar patterns are seen in an Abelian version of the model, where we observe $v_2\{2\}>v_2\{4\}\approx v_2\{6\}\approx v_2\{8\}$ of two, four, six, and eight particle correlations. While such patterns are often interpreted as signatures of collectivity arising from hydrodynamic flow, our results provide an alternative description of the multiparticle correlations seen in $p+A$ collisions.

A remarkable series of recent experiments at CERN's Large Hadron Collider (LHC) and at the Relativistic Heavy Ion Collider (RHIC) at BNL have demonstrated the existence of collective multiparticle dynamics in proton-proton (p + p) and light-heavy ion (h + A) collisions. Collectivity is represented by the behavior of n-th Fourier moments of the cumulants c n {m} of m-particle (m ≥ 4) azimuthal angular anisotropy correlations; it is observed that corresponding real valued m-th roots of c n {m}, the anisotropy coefficients v n {m}, have nearly identical values for high multiplicity events. These results are similar to those obtained in peripherally overlapping collisions of heavy nuclei and even exhibit some of the systematics observed in fully overlapping central heavy-ion collisions. The collective dynamics of the latter is well described by sophisticated hydrodynamic models which presume the creation of a thermalized strongly interacting Quark-Gluon Plasma (QGP).
Hydrodynamic models have also been employed to describe the experimental results on multiparticle correlations in the smaller systems. Their agreement with data is however sensitive to the initial conditions for hydrodynamic flow [1], specifically to the gluon "shape" fluctuations within the proton [2,3]. Since these shape fluctuations are themselves a consequence of strong initial state correlations, it is interesting to ask whether the collective properties of quark-gluon matter in the aforementioned small systems are those of nature's smallest fluids or whether there are alternative explanations for this collective behavior from initial state correlations alone.
An initial state correlation scenario to explain the "ridge-like" structure of two-particle correlations in small systems was advocated in [4] based on the Color Glass Condensate (CGC) effective theory [5,6]. (Refs. [7,8,35] review these and related frameworks.) For hadron transverse momenta p ⊥ , q ⊥ ≥ Q s,T , where Q s,T ∼ 1 − 2 GeV is the saturation scale of strongly correlated gluons in the nuclear target, a so-called Glasma graph scenario describes both ridge-like and jet-like correlations in the data on p + p and h + A collisions [9][10][11][12][13][14]. For lower momenta p ⊥ ≤ Q s,T , corrections of order Q s,T /p ⊥ are large. For two-particle correlations, the contribution of these corrections were quantified in [15] and the qualitative behavior of the anisotropy coefficients v n {2} were reproduced 1 . Collectivity from four-particle initial state correlations has remained elusive. Prior discussions were within a "color domain" model [20][21][22][23][24] whose theoretical foundations are unclear [15]. In this Letter, we will compute v 2 {m} for m ≥ 4 systematically for the first time in an initial state framework. An important ingredient is a first computation of the average of the product of four lightlike "dipole" Wilson-line correlators. We will explore the systematics of v 2 {4} as a function of Q 2 s,T , both integrated and differential in p ⊥ . In addition, we will study so-called four-particle symmetric cumulants which have recently been measured in light-heavy ion collisions. We observe strikingly that qualitative features of v n {2} and v 2 {4} measured in small systems are reproduced in this initial state framework. Since the computational effort for v 2 {m} for m > 4 increases rapidly with m, the same computation can be carried out in an Abelian variant of the model. Remarkably, the behavior of these higher cumulants are consistent with those observed in the LHC data. 1 These computations are within a dilute-dense approximation Q s,P Q s,T in the CGC framework, where Q s,P is the projectile saturation scale. For high multiplicity events, corrections of order Q s,P /p ⊥ are important and can be computed by numerically solving classical Yang-Mills equations [16,17] for two dense sources.

arXiv:1705.00745v2 [hep-ph] 29 Jan 2018
For our proof of principle computation 2 , we model the proton-nucleus collision very simply as the eikonal scattering of nearly collinear quarks in the projectile scattering off color domains of size 1/Q s,T inside the nuclear target [15,26,27]. The m-particle correlation can be expressed as Here we have made the simplifying assumption that the m-particle Wigner function representing quark distributions in the incoming proton factorizes as where N c is the number of colors and U (x) (U † (y)) are light-like Wilson lines appearing in the amplitude (complex conjugate amplitude) for quarks multiple scattering off gluons in the target. In the McLerran-Venugopalan (MV) model [28][29][30], these Wilson lines are path ordered exponentials of color charges in the target, and the average · · · in Eq.(1) is performed over a Gaussian distribution of color charges with a weight proportional to Q 2 s,T [5,31]. We will assume further that the Wigner distributions of the nearly collinear quarks have the Gaussian form where B p = 4 GeV −2 , a scale controlling the quark transverse momentum and spatial resolution, is fixed 3 using dipole model fits to HERA data [32,33]. We can perform the k i integrals in Eq.(1) explicitly, which gives Before we proceed to the computation of multiparticle cumulants, we will address some of the features and limitations of this simple model. First, we note that even though rapidity is not explicit in this model, particle correlations are ridge-like and long range in rapidity. As we demonstrate explicitly in [25], these correlations can be obtained in our model by convoluting the longitudinal momentum distributions of the quarks with their parton distributions in the incoming proton. Our model shares these features with the hybrid framework of multiparticle correlations discussed in [26,34,35]. This hybrid scenario will receive significant modifications when high parton density effects in the projectile become important. These effects quantitatively go as Q 2 s,P /p 2 ⊥ ; saturation models fit to HERA data conservatively suggest that these effects become non-negligible around x = 0.01 [33]. However, depending on the transverse momentum range studied, the qualitative features we observe could persist to smaller values of x. Parametrically the rapidity range where corrections to the hybrid model occur is Second, an obvious limitation of our model is that it only includes quarks. This is clearly not sufficient at the highest RHIC energies and at the LHC, though it may suffice to explain the ridge like correlations now seen at fairly low energies in deuteron-gold collisions at RHIC. Our model can be extended to include gluon degrees of freedom from the projectile; the only modification is that they will be color rotated by adjoint Wilson lines from the target and one has to compute color traces of products of these Wilson lines instead. However such computations alone are insufficient because they do not generate odd moments of the azimuthal distributions, a consequence of the generators of the adjoint representation being real [13,34,36]. The odd moments can, however, be recovered by going beyond the strict dilute-dense limit and including gluon exchanges between spectator partons and the scattered gluons in the projectile [17][18][19]. These considerations are at present beyond the scope of this Letter.
The two-and four-particle cumulants are defined as [37] and with where the integration over the two and four-particle phase space is implicit. The computation of two-particle cumulants is straightforward. The corresponding anisotropy coefficients are defined to be [37] v n {2} = (c n {2}) 1/2 , and were computed previously for the MV model in [15].
In computing these, we first fix the value of Q s,T and integrate over the momenta of all m particles in the range p ⊥ ∈ [0, p max ⊥ ]. We then study the variation of v n {2} with increasing Q 2 s,T ; in our simple model, this corresponds to increasing the center-of-mass energy or the centrality of the collision. By construction, our results are independent of the number of charged particles, N ch , produced in the collision. In Fig.(1), we plot the two-particle Fourier harmonics for n = 2, .., 5 as a function of Q 2 s,T . The upper limit of the transverse momentum integration is taken for multiple values to ensure convergence; by p max ⊥ = 3 GeV the results are no longer dependent on p max ⊥ . We observe a clear ordering of the n harmonics. These observations are in qualitative agreement with experiment. We should caution, however, that there is not a simple one-to-one map between Q s,T and the energy or centrality. Further, as noted in [15], the QCD evolution of the MV model with energy will lead to lower values of v n {2}. Fragmentation of gluons into hadrons will further soften the signal [16]. Our results therefore represent maximal values for azimuthal correlations in this initial state framework.
We will now go beyond the study in [15] and discuss the behavior of the four-particle flow coefficient defined as [37] v n {4} = (−c n {4}) 1/4 .
The computation of four-particle cumulants is significantly more complex for two reasons. First, as is clear from Eqs.(3)-(6), one has 24 integrals to perform relative to 12 previously for two-particle cumulants. However, more importantly, computing the expectation value of the product of four dipole correlators is non-trivial. While an analytical expression exists for the expectation value of the product of two dipoles [15,38,39], such an expression is not known for the product of four dipoles. For work in this direction, see [40]. Our strategy for computing n-dipole correlators follows the framework introduced in [38]. For four dipoles, the building blocks are four light-like Wilson lines, localized at distinct transverse positions, in the scattering amplitude for four quarks along the light cone from x + = −∞ to x + = +∞ and their counterparts in the complex conjugate amplitude. The terms in the expansion of the Wilson lines correspond to multiple gluon exchanges, ordered in the x + direction, between the dipoles. These exchanges generate quadrupole, sextupole, and octupole configurations that are, respectively, the traces over the product of four, six, and eight light-like Wilson lines 4 . Permutations of the coordinates for each of these topologies results in 24 distinct basis elements. Subsequent gluon exchanges generate transitions between elements of this basis, corresponding to a 24 × 24 matrix that can be exponentiated numerically.
All the basis elements are known in the MV model and the result of the computation can be expressed in terms of the saturation scale Q s,T and a cutoff Λ regulating the infrared behavior of the two-dimensional gluon exchange propagator. We choose Λ = 0.241 GeV [41]; our results are insensitive to variations in this scale 5 . The procedure can be extended to the average of m > 4 dipole correlators and is discussed at length in [25].
Computing Eq.(1) as outlined, we can evaluate v n {4} using Eqs. (5), (6), and (8). The results are shown in Fig.(2). For reference, we also plot the v 2 {2} values shown in Fig.(1) for p max ⊥ = 2, 3 GeV. First, we see striking evidence of collectivity as defined: the value of c 2 {4} is negative allowing us to extract a real valued v 2 {4}. It is smaller than the value of v 2 {2} and both of these are relatively flat as a function of Q 2 s,T . Since the increase in Q s,T corresponds to an increase in the center-of-mass energy, the two and four-particle elliptic anisotropy coefficients are independent of energy in our initial state model. Experimental results from light-heavy ion collisions at RHIC and LHC similarly show a weak variation of these quantities across a very wide window of centerof-mass energies [43][44][45][46][47]. As noted previously, our results are independent of N ch , as is also approximately the case in experiment. 4 For a number of dilute-dense multiparticle processes only dipoles and quadrupoles contribute in the high energy limit [42]. However because the leading contributions to the correlation observables here are themselves Nc suppressed, sextupoles and octupoles are of equal importance. 5 For values Q 2 s,T ≈ (p max ⊥ ) 2 Λ 2 , the results are insensitive to Q 2 s,T Bp, the number of color domains in the target [25]. To compute the two and four-particle Fourier harmonics for a fixed transverse momentum, we define and Theκ n are the differential form of κ where one of the momenta are not integrated over. The v n anisotropies as a function of p ⊥ are given by [47] v We plot these for n = 2 and m = 2, 4 in Fig.(3) for a fixed Q 2 s,T = 2 GeV 2 . The error bands represent the systematic uncertainty in the integrated four-particle cumulant. We again see that v 2 {m}(p ⊥ ) has the same qualitative behavior as data in light-heavy ion experiments at RHIC and LHC. Energy evolution of parton distributions and parton to hadron fragmentation will decrease the values shown. Increasing Q 2 s,T has the effect of flattening out v n {m}(p ⊥ ) at higher p ⊥ [25].
Symmetric cumulants defined as have recently been measured at the LHC. Theκ n,n are analogous to κ n in Eq. even angles have the harmonic n . The SC(n, n ) cumulants directly measure correlations between the different flow harmonics [48]. In heavy-ion collisions, the data show that SC(2, 3) are increasingly anti-correlated with increasing centrality percentile, while the SC(2, 4) cumulant are increasingly correlated; these systematics are also seen in hydrodynamic models [49]. The SC(2, 3) and SC (2,4) cumulants have also been measured in light-heavy ion collisions and show the same pattern of correlations as for heavier systems [50]. This is perhaps not too surprising because the correlations/anticorrelations are most significant in peripheral heavy-ion collisions. Our results for SC (2,3) and SC (2,4) are shown in Fig.(4) as a function of Q 2 s,T . We see that SC(2, 3) is negative by Q s,T = 1 GeV while SC (2,4) is positive for all Q s,T . Our results demonstrate clearly that such patterns are not unique to an interpretation requiring hydrodynamic flow. Results from hydrodynamic computations for these cumulants in light-heavy ion collisions are not yet available.
To gain further insight into our results, it is useful to ask whether coherent multiple scattering off the target is crucial. One way to test this within our framework is to employ the Glasma graph approximation [4,[9][10][11][12][13][14], valid for p ⊥ > Q s,T . For two partons scattering off the target color fields, the Glasma graph approximation corresponds to two gluon exchange in the scattering amplitude [15]. Non-linearities, that are large for p ⊥ ≤ Q s,T in the MV model, arise from multiple gluon exchanges between the projectile and the target. One can similarly implement the "linear" Glasma graph approximation for four partons scattering off the target; we find c 2 {4} for the Glasma graphs is positive [25]-this confirms the importance of coherent multiple scattering. It is also interesting to consider coherent multiple scattering in the Abelian limit of this model. In this case, the Wilson lines are not matrices in color space, but simply path ordered exponentials [51]. The product of dipoles in Eq.(1) is significantly simpler to compute [25], enabling one to extract v 2 {6} and v 2 {8} from the corresponding cumulants [37,52]. Our results, shown in Fig.(5), demon- , as also seen in the LHC data on multiparticle harmonics [46,52].
The fact that this behavior is reproduced in a simple initial state model is a proof of principle that it is not unique to interpretations of collectivity arising from the hydrodynamic response of the system to the n-th moments of m particle spatial eccentricities [53][54][55][56]. For a recent review on hydrodynamic collectivity and relevant references, see [57]. Our results do not necessarily mean that an initial state interpretation of the data is favored. We instead conclude that the v n {m} measurements alone are insufficient to unambiguously distinguish between initial and final state approaches.
While it is remarkable that our results qualitatively explain observed multiparticle correlations, it is also clear that the model is missing key features of QCD dynamics that should be important at high energies. In this regard, the initial state framework in [16,17] includes a more systematic treatment, albeit at an enormously greater computational effort. Nevertheless, since multiparticle correlations display similar features in light-heavy ion collisions spanning two orders of magnitude in centerof-mass energies, where QCD degrees of freedom evolve significantly, it is worth thinking further why this simple model appears to capture the underlying dynamics.