Azimuthal anisotropies in p+Pb collisions from classical Yang-Mills dynamics

We compute single and double inclusive gluon distributions in classical Yang-Mills simulations of proton-lead collisions and extract the associated transverse momentum dependent Fourier harmonics $v_2(p_T)$ and $v_3(p_T)$. Gluons have a large $v_2$ in the initial state, while odd harmonics such as $v_3$ vanish identically at the initial time $\tau=0^{+}$. By the time $\tau \lesssim 0.4\,{\rm fm/c}$ final state effects in the classical Yang-Mills evolution generate a non-zero $v_{3}$ and only mildly modify the gluon $v_{2}$. Unlike hydrodynamic flow, these momentum space anisotropies are uncorrelated with the global spatial anisotropy of the collision. A principal ingredient for the generation of $v_2$ and $v_3$ in this framework is the event-by-event breaking of rotational invariance in domains the size of the inverse of the saturation scale $Q_s$. In contrast to our findings in p+Pb collisions Yang-Mills simulations of lead-lead collisions generate much smaller values of $v_{2,3} (p_T)$ and additional collective flow effects are needed to explain experimental data. This is because the locally generated anisotropy due to the breaking of rotational invariance is depleted with the increase in the number of uncorrelated domains.


I. INTRODUCTION
A striking result from high-multiplicity proton-proton (p+p) and proton-lead (p+Pb) collisions at the LHC was the discovery of a "double-ridge" structure in twoparticle correlations that is long range in their rapidity difference ∆η and includes a dominant cos(2∆φ) modulation in their relative azimuthal angle [1][2][3][4][5][6][7]. Key aspects of the structure of the observed correlations in protonnucleus collisions bear a striking similarity to those observed in heavy-ion collisions [8][9][10] and may point to a form of collective behavior where many particles are correlated with each other.
In heavy-ion collisions the azimuthal structure of multi-particle correlations is quantitatively well described by viscous fluid dynamic calculations with fluctuating initial state geometries [11]. This naturally leads to the assumption that the physics responsible for the ridge structure in high multiplicity p+p and p+Pb events may also be driven by the final state collective flow of the system. Indeed, calculations using hydrodynamics (or microscopic models of final state interactions) are able to reproduce features of the azimuthal structure in p+Pb collisions [12][13][14][15][16]. Specifically, the Fourier coefficients v 2 (p T ) and v 3 (p T ) in the expansion of the particle distribution are described by several of these models, albeit with different parameters. In Eq. (1), ψ n = 1 n arctan sin(nφ) cos(nφ) is the event plane angle associated with the n th harmonic, and · denotes the average with respect to the single inclusive particle distribution.
While some detailed features of the data -such as the dependence of v 2 (p T ) on the particle mass -can be explained quite naturally by final state collective effects, a number of conceptual problems can be identified within this theoretical approach. One concerns the applicability of viscous hydrodynamics due to large pressure gradients [17,18] that are present for a significant fraction of the space-time evolution. Another concerns the sensitivity to the initial state in small systems [19], which requires a better theoretical understanding of the early time dynamics. Further, the observation of pronounced azimuthal anisotropies even at high transverse momenta p T > ∼ 3 GeV challenges the hydrodynamic paradigm which is best applied to a description of low momentum excitations. It is therefore important to explore if multiparticle correlations at different transverse momentum scales can be understood in part or whole in alternative approaches.
For instance, computations of intrinsic two particle correlations in the Color Glass Condensate (CGC) framework have been shown to produce azimuthal anisotropies compatible with ridge data for even harmonics for p T > 1 GeV in p+p and p+Pb collision systems without the need for additional final state collective effects [20][21][22][23][24]. However, no odd harmonics were generated because rescattering contributions to the intrinsic correlations were not included [25]. In addition to these intrinsic twoparticle correlations, the presence of domains of directed chromo-electric fields inside the proton and the nucleus breaks rotational invariance on an event-by-event basis and thereby generates azimuthal anisotropies [26]. For quarks scattering off a colored target, these colored domains generate both v 2 and v 3 [27,28]. However, gluons scattering off this target only generate even harmonics. We note that the possibility that azimuthal anisotropies could arise from the event-by-event breaking of rotational invariance by color fields has also been considered in a related approach [29].
Thus far the studies of initial state correlations generated by color-electric domains in [27,28], and recent extensions thereof in [30], have been based on a description of the proton as a dilute projectile of valence quarks scattering off the nucleus. However, experimentally significant anisotropies are only measured in events with very high multiplicities, where it is more appropriate to describe both the proton and the lead nucleus as dense colored objects. Because of the high gluon occupancies in both the proton and the lead nucleus, they can be approximated as classical gluon fields, and their leading order space-time evolution is described by solving classical Yang-Mills equations.
Such classical Yang-Mills (CYM) simulations were performed previously for nucleus-nucleus collisions [31,32] (including an early study of elliptic flow in [33]) and proton-nucleus collisions [34,35]. Recently, these studies were extended to include more realistic initial conditions in the IP-Glasma model [36,37], providing a satisfactory description of multiplicities in high energy hadron collisions [35]. Since the IP-Glasma model includes fluctuations of color charges which break rotational invariance on an event-by-event basis, we will use this framework to study azimuthal anisotropies of gluons in the initial state and during the early time dynamics of proton-nucleus and nucleus-nucleus collisions. When studying p+Pb collisions, we will also consider additional fluctuations of the proton geometry on sub-nucleon size scales [38] and study their effect on the Fourier harmonics.
This letter is organized as follows. In Sec. II we briefly outline the theoretical framework underlying the classical Yang-Mills simulations. Our discussion follows the literature in the context of the IP-Glasma model [37] with the significant modification that we will also consider 'eccentric' proton configurations following [38]. We then discuss the measurement of azimuthal anisotropies in this framework in Sec. III and present numerical results for azimuthal Fourier harmonics of gluons in protonnucleus collisions in Sec. IV. We investigate the sensitivity of our results with respect to variations in the spatial color structure of the proton and perform a comparison of the effects in proton-nucleus and nucleus-nucleus collisions in Sec. IV. The final section summarizes our conclusions and their implications for collective dynamics in proton-nucleus and nucleus-nucleus collisions.

II. THEORETICAL FRAMEWORK
Within the CGC framework [39][40][41], the dynamics of a high-energy collision is -to leading order in α s -described by solutions of the classical Yang-Mills equations, in the presence of an eikonal color current J ν . Here, D µ is the covariant derivative in the presence of the field A µ and F µν is the gluon field strength tensor. For a right moving (projectile) proton and left moving (target) nucleus one has and each event is characterized by a different color neutral distribution of random color charges ρ p/Pb (x ∓ , x) inside the proton and nucleus.
Before the collision, the small-x gluon fields inside the target and projectile nucleus are determined by the solution of (2) and can be compactly expressed 1 as [42][43][44] in terms of the fundamental Wilson lines V p/Pb (x ⊥ ) of the projectile and target nucleus. By dividing the longitudinal direction into N Y discrete rapidity intervals, these can be computed as [45] for a given configuration of color charges. 2 Here, ∇ 2 ⊥ = ∂ i ∂ i . We employ a McLerran-Venugopalan type model [42] for the color charge densities, which follow local Gaussian distributions with variance where b = (x + y)/2. The spatial distributions of color charge inside the proton and the nucleus are described by S p/Pb (b) and the respective models are outlined below.

A. The proton
We will consider two different models for the distribution of color charge S p (b) inside the proton to study the effect of the spatial sub-structure of the proton on the observed correlations.
The spherical proton model is a variant of the IP-Sat model [46] where the color charge distribution inside the proton is spherically symmetric in impact parameter space. In this case, the color charge density is proportional to the saturation scale, i.e., Q s itself depends on the transverse position b via the nucleon thickness function The Gaussian width is related to the (two dimensional) proton radius relevant to strong interactions as R p = √ 2 · B G [47]. Its value B G = 4 ± 0.4 GeV −2 as well as the dependence of Q 2 s on the nucleon thickness t(b) and the center of mass energy √ s have been extracted from fits to deep-inelastic scattering (DIS) data [48]. 3 We set B G = 4 GeV −2 and the only free parameter in this model is the proportionality factor c in Eq. (7), defined as g 4 µ 2 /Q 2 s , where g 2 µ 2 is the color charge square per unit transverse area. We choose to use a fixed value of c = 2 and shall comment later on the sensitivity of our results under variation of c.
The constituent quark proton model was previously outlined in [38], where the distribution of the color charge density is concentrated around the (transverse) positions x CQ of three constituent quarks which fluctuate from event to event according to a Gaussian distribution with expectation value x 2 CQ = B G . The gluon distribution around each constituent quark is spherically symmetric with a radius denoted by R CQ . We will use R CQ = √ B G /2 and adjust the overall strength Q = 3 GeV, which yields similar results for the dipole scattering amplitude (relevant to DIS) as the spherical proton model.

B. The nucleus
Since we expect a smaller sensitivity of our results to the impact parameter dependent structure of the nucleus, we limit ourselves to a single model for the color charge distribution inside the nucleus. We first sample the positions x i of A = 208 (for Pb) individual nucleons according to a Wood-Saxon distribution with radius and surface parameters appropriate for a Pb nucleus. We follow the IP-Sat model and set where the thickness function T (b) of the nucleus is the sum of thickness functions t(b − x i ) of individual nucleons. We note that for the case of a single nucleon (A = 1), this reduces to the spherical proton model and we employ precisely the same parametrization in both cases. 3 In practice one extracts the saturation scale Qs( √ s, t(b)) from the IP-Sat parametrization of the dipole scattering amplitude. More details on this procedure and the general features of this model can be found in [35].

C. Early-time dynamics and gluon distribution
Solving the classical Yang-Mills equations outside the forward light-cone leads to the initial state immediately after the collision (τ = 0 + ). The initial gauge fields in Fock-Schwinger gauge A τ = (x + A − + x − A + )/τ = 0 are given in terms of the projectile and target fields as [49,50] and correspond to longitudinal chromo-electric and chromo magnetic fields [51] Starting from these field configurations, the early time dynamics in each event can be determined by numerically solving the classical Yang Mills equations in the forward light cone. Our numerical implementation is based on standard lattice gauge theory techniques and we refer the reader to [31,32,37] for more details of this procedure. We can then extract the gluon distribution at different times of the evolution by measuring equal-time correlation functions of the gauge fields. We impose the Coulomb gauge condition ∂ i A i τ at the time of each measurement and follow previous works [52,53] to compute the single particle spectrum by a projection on to transversely polarized gluon modes. The single particle distribution is then given by where g µν = (1, −1, −1, −τ −2 ) denotes the Bjorken metric, λ = 1, 2 labels the two transverse polarizations and a = 1, · · · , N 2 c − 1 is the color index. In Coulomb gauge the mode functions take the form where k ⊥ = (k x , k y ) and H (2) α denote the Hankel functions of the second type and order α (see [52] for details).
We note that the above definition of the gluon distribution is such that dN/d 2 k ⊥ is exactly conserved for a non-interacting system. This property is important, as it will enable us to clearly distinguish between the properties of the initial state at τ = 0 + and the effect of final state interactions at later times.

III. AZIMUTHAL ANISOTROPIES
We will now discuss the measurement of azimuthal anisotropies, in particular the extraction of v 2 (p T ) and v 3 (p T ) for gluons, using two different methods. One is based on the traditional event plane method and the other follows more closely experimental measurements based on two-particle correlations.

A. Event-plane method
The event plane method determines the Fourier coefficients v 2 (p T ) and v 3 (p T ) from the single particle spectra. Since the lattice simulation yields the single particle spectrum at discrete values of the transverse momenta k x and k y , we first perform a bi-linear interpolation and divide the data into transverse momentum |k ⊥ | and azimuthal angle φ bins. We then compute the second and third order event plane angles Ψ 2 and Ψ 3 in each event, for gluons in the reference momentum region chosen to be 1 GeV < p ref T < 6 GeV as discussed below. This allows us thereby, event-by-event, to determine the coefficients v 2 (EP )(p T ) and v 3 (EP )(p T ) in the Fourier expansion of the single particle spectrum in Eq. (1). The final result is obtained by performing an average of these quantities over all events.

B. Two particle correlations
The measurement of two-particle correlations in experiments is based on the per trigger yield defined to be [3,6] 1 N trig Here N trig is the number of trigger particles in the momentum bin under consideration, is the signal and is the background contribution estimated from mixed events.
Since there are no acceptance and efficiency corrections in the theory calculation, we can directly compute the quantity [54,55] 2π where N assoc is the number of associated particles in the momentum bin considered and · denotes an average over events. Note that while in experiments a rapidity gap is introduced to suppress jet-like correlations around ∆η = 0, jet-like correlations are not present in our calculation at this order [22] and ∆η gaps are therefore unnecessary. The Fourier expansion of Eq. (21), yields the coefficients V n∆ , from which we define the two particle gluon v 2 (p T ) and v 3 (p T ) to be as is done by the experimental collaborations to measure v n (p T ) for inclusive hadrons. We choose the reference momentum range for gluons to be 1 GeV < p ref T < 6 GeV. The upper limit in p T of this range extends to somewhat larger p T than that of experimental measurements for inclusive hadrons [3,6]. This is to account very roughly for the fragmentation of higher momentum gluons into hadrons in the p T bin of interest. We will compute V n∆ (p T , p ref T ) in p T bins of width 0.25 GeV defined symmetrically around the quoted p T values. We will comment in Sec. V on the sensitivity of our results to the reference momentum range.

IV. NUMERICAL RESULTS FOR p+PB COLLISIONS
We will now present numerical results for the azimuthal correlations of gluons in p+Pb collisions. We first study the behavior for collisions of a constituent quark proton with a lead nucleus and focus on collisions with zero impact parameter b = 0 fm without any further multiplicity selections applied. Our results for v 2 (p T ) and v 3 (p T ) of gluons -obtained from an average over N evt = 128 events -are shown in the left and right panels of Fig. 1. We also show experimental results for inclusive hadron v 2 (p T ) and v 3 (p T ) from the ATLAS [6] and CMS [3] collaborations to provide an estimate of the relative magnitude and momentum dependence of the observed correlations. We emphasize that since the gluon to hadron conversion is by no means straightforward, we do not expect a quantitative description of the data within the present framework. We find that already at the initial time τ = 0 + the v 2 (p T ) of gluons is quite large and extends up to relatively high transverse momenta, which can be seen in both the two particle correlation and event plane measurements. The fact that both methods give rise to very similar results points to the fact that the origin of the observed v 2 (p T ) is due to a breaking of rotational symmetry of the single particle spectrum. In other words, we conclude that gluons are produced with a preferred azimuthal direction in each event.
While the emergence of a preferred azimuthal direction is well understood in the context of a collective expansion of the system, we emphasize that the nonzero gluon v 2 (p T ) at the initial time can not be attributed to collective flow effects. As a consequence of the initial color fields in Eq. (14) the energy momentum tensor at τ = 0 + is strictly diagonal ) and the Poynting vector E × B characterizing the energy momentum flow vanishes identically. Instead, the observed v 2 (p T ) should be attributed to anisotropic gluon production from the fluctuating color fields inside the projectile and target.
Our result for v 2 and v 3 from two-particle correlations can be separated into three contributions. First, the Glasma graph contribution, which corresponds to the connected graphs analyzed in [20][21][22][23][24][54][55][56]. Second, a contribution from disconnected graphs, which takes into account event-by-event breaking of rotational symmetry similar to the effect discussed in [26-28, 30, 57], and finally, a contribution from final state interactions repre-sented by the Yang-Mills dynamics of the produced gluon fields.
The contribution from disconnected graphs can be understood in the following simplified picture: One can imagine the gluon production process as the partons inside the projectile scattering off a high-energy nucleus. Each parton experiences the color electric field inside the nucleus whereby it receives a momentum kick in a certain direction. In each event the color electric field of the nucleus is characterized by several domains with different orientations in azimuthal direction and color space. While partons scattering off different domains receive a different kick, partons (with the same color charge) scattering off the same domain will receive a kick in the same azimuthal direction. Hence, depending on the number of uncorrelated domains probed by the projectile, the initial state particle production can be significantly anisotropic.
The Glasma graph contribution directly probes the color structure within each domain and corresponds to genuine non-factorizable two-particle correlations. Depending on the degree of polarization of a single domain, the two contributions can be of comparable size [57].
With regard to the effect of final state interactions on the observed gluon v 2 (p T ), we find that the classical Yang-Mills evolution leads to slight decrease of the observed v 2 (p T ) with time. This decrease may be attributed to gluon rescattering after the collision. The overall effect of the classical Yang-Mills evolution on v 2 is very small and the initial state value provides a very good estimate for the correlation at later times. We note again that this is conceptually quite different from a hydrodynamic mechanism, where v 2 (p T ) is gradually build up as a function of time.
We now turn to a discussion of the gluon v 3 (p T ) shown in the right panel of Fig. 1 for which the interpretation is drastically different. We find that the initial state v 3 (p T ) vanishes identically at τ = 0 + . This is a consequence of the vanishing transverse color electric (E ⊥ ) and color magnetic fields (B ⊥ ) at the initial time. With only the longitudinal color electric (E η ) and color magnetic (B η ) fields being non-zero at τ = 0 + , the gluon spectrum in Eq. (15) takes the form The above spectrum is manifestly symmetric under k T ↔ −k T -hence odd Fourier harmonics vanish at the initial time. As a result, we can attribute the observation of a non-zero v 3 (p T ) in our framework at later times exclusively to coherent final state effects included via the classical Yang-Mills evolution. Quantitatively, we find that already at very early times τ ∼ 0.2 fm/c the classical Yang-Mills evolution leads to the build up of a sizable gluon v 3 (p T ) extending to relatively large transverse momenta. Beyond τ = 0.2 fm/c, the v 3 of hard gluons remains approximately constant while the low momentum v 3 continues to show an increase with time. When we follow the classical Yang-Mills dynamics to even later times the system becomes more and more dilute and approaches a free streaming behavior around τ = 0.4 fm/c as previously reported in [37].
We note the agreement between the two particle correlation and event plane measurement of v 3 (p T ) again points to a correlation between many particles in each event. While it may appear suggestive that the build up of energy-momentum flow may cause a non-zero v 3 (p T ) at later times, to our surprise, we did not observe a significant correlation between the global initial state eccentric-ity ε 3 and the p T integrated momentum space anisotropy v 3 on an event-by-event basis. However, it is possible that there is a correlation between the final v 3 and triangular anisotropies on shorter length scales than the ones probed by the global ε 3 -such geometrical features are however difficult to extract. Generally, a simple description of the non-linear dynamics underlying the emergence of a nonzero v 3 (p T ) remains elusive -a deeper understanding of this effect is clearly desirable.

V. SENSITIVITY TO PROTON STRUCTURE AND COLLISION GEOMETRY
We will now study the effect of the collision geometry and system size on the azimuthal correlations. Our results are compactly summarized in Fig. 2 which shows a comparison of the gluon v 2 (p T ) and v 3 (p T ) at time τ = 0.4 fm/c after the collision.
We first analyze the effect of a more detailed substructure of the proton projectile on the correlations in p+Pb collisions by comparing a spherical proton with one composed of three valence quarks. Generally, a finer substructure leads to larger v 2 and v 3 at transverse momenta p T > ∼ 2 GeV -corresponding to wave lengths on size scales much smaller than the nucleon size. However the overall effect of the proton's geometry for the observed azimuthal anisotropy of gluons is far less significant than one would expect in a mechanism that generates azimuthal anisotropies via final state collective effects [19]. Since the origin of the observed correlations is due to the microscopic structure of color fields, one instead expects the correlations to be approximately independent of the global event geometry. Our result in Fig. 2 confirms this picture.
We have also considered variations of the nonperturbative mass scale m and the coefficient c in Eq. (7) by a factor of two. While in both cases we did not observe a significant effect on the overall magnitude of the observed correlations, we found that for smaller (larger) values of m the correlations extend over a slightly larger (smaller) transverse momentum range. We note further that changing the reference momentum range p ref T can also have a significant effect on the transverse momentum dependence of the signal which is somewhat more pronounced for v 2 (p T ) as compared to v 3 (p T ). While the gluon spectrum is generally anisotropic up to high momenta, the event plane angles at different transverse momenta Ψ n (p T ) are not strongly correlated across the entire range of momenta. Considering a larger (smaller) reference momentum range can therefore lead to a slower (faster) decrease at high momenta. Similar observations have recently been reported in a related study [30] in which the proton was treated as a dilute object.
We conclude our study with a comparison between proton-nucleus (p+Pb) and nucleus-nucleus (Pb+Pb) collisions. In the latter case, we have analyzed collisions at two different impact parameters corresponding to central (b = 0 fm) and peripheral (b = 11 fm) collision. 4 We find that in both central and peripheral nucleusnucleus collisions the correlations between gluons are much smaller compared to p+Pb collisions. Qualitatively this difference can be understood when considering the different number of localized domains of fluctuating color fields responsible for the production of gluons in the different collision systems. While in proton-lead collisions, particles are produced from a small number of different domains inside the nucleus, the larger overlap area in lead-lead collisions gives rise to particle production from a much larger number of different domains. Since different domains are uncorrelated with each other the azimuthal anisotropy of the gluon spectrum decreases with the number of domains (see also [27]). Consequently, the initial state momentum space anisotropy is much smaller in Pb+Pb collisions as compared to p+Pb collisions.

VI. CONCLUSIONS
We have presented results for the azimuthal anisotropy of the single and double inclusive gluon distributions in p+Pb and Pb+Pb collisions obtained from classical Yang-Mills simulations. Both the proton and the nucleus have been treated as dense QCD objects with high gluon occupancy. This description is appropriate for the early time space-time evolution of high multiplicity p+Pb as well as heavy ion collisions at high energies.
Gluons produced in p+Pb collisions show a significant v 2 (p T ) already at the initial time immediately after the collision. Further evolution governed by the Yang-Mills equations modifies this v 2 (p T ) only slightly. In contrast, odd harmonics of gluons are initially exactly zero, but significant values of v 3 (p T ) are built up within times τ < ∼ 0.4 fm of the classical Yang-Mills evolution. These momentum space anisotropies at early times are uncorrelated with the global spatial anisotropy, in contrast to anisotropies generated by collective flow.
Our results indicate that in p+Pb collisions there are significant contributions to v 2 and v 3 from the initial production (in the case of v 2 alone) and the early time nonequilibrium dynamics within the first half fermi of evolution. These effects cannot be neglected and any calculation based merely on final state effects is thus incomplete.
A similar analysis of Pb+Pb collisions reveals a different picture. Initial and early-time contributions to v 2 and v 3 are small, indicating that for larger systems final state collective effects are indeed the dominant mechanism for generating the observed anisotropies -at least at momenta p T ≤ 2 GeV, where the presence of such effects is very plausible. The difference between p+Pb and Pb+Pb collisions can be understood as a consequence, in this framework, of anisotropies being generated due to localized domains of color fields. A large number of mutually uncorrelated domains probed in Pb+Pb collisions leads to a nearly isotropic gluon spectrum. Hence initial state contributions to v 2 and v 3 are small for large collision systems.
While our study provides a first attempt to quantify the importance of initial state effects in high-multiplicity proton-nucleus collisions, we expect systematic comparisons of p+A collisions with deuteron-nucleus (d+A) [58][59][60][61] and 3 He+Au [62] collisions at RHIC to provide further insight into the relative significance of initial state and final state effects in small systems. Forthcoming p+A collisions at RHIC will also help clarify the role of nucleon fluctuations relative to sub-nucleon scale effects in small systems besides providing a comparison of results for identical systems at vastly different center of mass energies.