Hybrid Color Glass Condensate and hydrodynamic description of the Relativistic Heavy Ion Collider small system scan

Multi-particle correlation observables in the Relativistic Heavy Ion Collider small system scan are computed in a framework that contains both initial state momentum anisotropies from the Color Glass Condensate effective theory and final state hydrodynamic evolution. The initial state is computed using the IP-Glasma model and coupled to viscous relativistic hydrodynamic simulations, which are followed by microscopic hadronic transport. All parameters of the calculation were previously constrained using experimental data on Au+Au collisions at the same center of mass energy. We find that the qualitative features of the experimental data, such as the system and centrality dependence of the charged hadron momentum anisotropy, can only be reproduced when final state interactions are present. On the other hand, we also demonstrate that the details of the initial state are crucially important for the quantitative description of observables in the studied small systems, as neglecting the initial transverse flow profile or the initial shear stress tensor, which contain information on the momentum anisotropy from the Color Glass Condensate, has dramatic effects on the produced final state anisotropy. We further show that the initial state momentum anisotropy is correlated with the observed elliptic flow in all small systems, with the effect increasing with decreasing multiplicity. We identify the precise measurement of $v_2$ in d+Au and Au+Au collisions at RHIC energy at the same multiplicity as a means to reveal effects of the initial state momentum anisotropy.

Both contributions have been studied in a combined framework of the IP-Glasma initial state [38,39] and a microscopic description of the gluonic final state interactions via solutions of the Boltzmann equation in [40]. It was found that initial state correlations survive for transverse momenta p T > ∼ 1.5 GeV even in high multiplicity p+Pb events, with their effect and affected p T range increasing with decreasing multiplicity. We go beyond this study in two essential ways. First, we allow for the comparison to experimental data by computing hadronic final states instead of studying parton spectra, and second, we quantify how strongly the CGC initial state momentum anisotropy can affect experimental observables if final state interactions are described by realistic hydrodynamic simulations that were constrained by heavy ion data.
Previous calculations involving the IP-Glasma initial state coupled to a hydrodynamic description of the final state evolution [41][42][43] in principle also contain informa-tion on both CGC initial state and geometry driven final state sources of correlations. However, only the work on p+Pb collisions at LHC in [37] included the full initial energy momentum tensor from the IP-Glasma calculation to initialize hydrodynamics, which is necessary in order to keep the full initial state information, and to conserve energy and momentum when switching descriptions.
The striking results from the recent small system scan at RHIC [44] provide a strong motivation to extend our work to lower energies and to analyze in more detail the features of initial and final state effects on the observed momentum anisotropies. We study different collision systems (p+p, p+Au, d+Au, and 3 He+Au) at the top RHIC energy, using the hybrid framework of IP-Glasma + Music hydrodynamics + UrQMD [45,46] with parameters constrained by heavy ion collisions [47]. Apart from presenting results for two-and four-particle correlations and comparison to experimental data, we show the multiplicity and system dependence of the initial elliptic momentum anisotropy, and its correlation with the observable elliptic flow.
The energy momentum tensor from the IP-Glasma framework, used to initialize the hydrodynamic simulation, includes the single particle CGC momentum anisotropy. We study individually the effects of the initial transverse flow profile and the viscous components of the energy momentum tensor. This allows us to determine the importance of different features of the initial state description, many of which are neglected in recently explored approximations to the IP-Glasma model [48][49][50] (see also [51]).
Framework. We employ the hybrid framework described in [37] and [47], consisting of the boost-invariant IP-Glasma initial state, which provides the spatial (x ⊥ ) distribution of (locally anisotropic) energy-momentum tensors as input for the hydrodynamic simulation Music, which in turn is coupled to the hadronic transport simulation UrQMD, which describes the low energy density regime.
The IP-Glasma model samples nucleon positions from nuclear density distributions, including a three-hotspot substructure for each nucleon to generate the nuclear geometry of a single configuration. For the deuteron, we sample a Hulthen wave function [52,53], for 3 He we use the same configurations as in [54], obtained using Green's function Monte Carlo calculations [55]. For gold nuclei we use Woods-Saxon distributions as described in [47]. Given that geometry, the IPSat model [56], constrained by HERA data [57] provides the distribution of the color charge density, which determines the variance of the Gaussian distributed color charges, which are then sampled. These color charges are the sources for the gluon fields, which are found by solving the SU(3) Yang-Mills equations numerically [58]. From the two incoming sheets of gluon fields one determines the gluon fields after the collision [38,39,58], which are then evolved in time. At time τ init we switch from the Yang-Mills to the hydrodynamic description by using the gluon fields' energy momentum tensor T µν CYM (x ⊥ ) as initial condition for the hydrodynamic simulations. We use τ init = 0.4 fm/c and will comment on the effect of varying τ init between 0.2 and 0.6 fm/c.
As in [37] and [47], we include the full energy momentum tensor, which has the shear stress contribution The initial energy density ε and flow velocity u µ can be extracted by solving u µ T µν CYM = εu ν . On the Yang-Mills side, the pressure is given by P CYM = ε/3, while in the hydrodynamic simulation we employ a lattice QCD equation of state [59,60] with pressure P lat . There are various ways of dealing with this mismatch. One way is to absorb the difference in an effective initial bulk viscous term Π = ε/3 − P lat , which vanishes on a time scale approximately given by the bulk relaxation time. The initial Π constructed this way leads to an additional outward push [47]. Alternatively, one can accept a jump in the pressure at the time of matching, by setting the initial Π = 0. The former option is our preferred one, as it was used in [47]. We will also study the latter option, which will determine one contribution to our systematic uncertainties.
Once initialized at time τ init , the energy momentum tensor follows ∂ µ T µν = 0, where Gluon multiplicity Ng distributions from the IP-Glasma model in p+p, p+Au, and d+Au collisions and charged particle multiplicity N ch distribution after hydrodynamic evolution in d+Au collisions, scaled by the mean multiplicity and compared to scaled N ch distributions from STAR [67] for d+Au (uncorrected) and UA5 [68] for p+p (corrected).
with ∆ µν = g µν − u µ u ν . We use the same second order constitutive relations for the shear and bulk viscous parts as derived in [61,62] and used in [47], and solve them within the simulation Music [41][42][43]. We use the same shear and bulk viscosities as in [47]. Because the second order transport coefficients are not well known, we will vary all (except the relaxation times) between what was used in [47] (values determined from a Boltzmann gas in the small mass limit [61,62]) and zero, which constitutes the other contribution to our systematic errors. When the medium temperature drops to the switching temperature T sw = 145 MeV, the fluid is converted to particles by first computing the particle spectra according to the Cooper-Frye formula [63], using equilibrium distributions f 0 with viscous corrections δf , given in [64][65][66] for shear and bulk viscous terms. From these non-equilibrium distribution functions f = f 0 + δf , we sample particles on the switching surface 1 that then undergo the microscopic transport processes of UrQMD [45,46]. To ensure enough statistics, each individual hydrodynamic hyper-surface is oversampled until we reach at least 100,000 particles per unit of rapidity.
Multiplicity distributions. Multi-particle correlations in small systems are very sensitive to event-by-event fluctuations. Consequently, basic observables like the charged particle multiplicity distributions must be described correctly. We present both the gluon multiplicity distributions (folded with a Poisson distribution to estimate the effect of (grand canonical) sampling on the The pT of π + , K + , and p as a function of charged particle multiplicity in d+Au collisions compared to experimental data from the STAR Collaboration [67]. Dashed lines show pT in p+Au collisions. See text for details.
switching surface), which is used to determine our centrality classes, and for d+Au collisions also the final state charged hadron distribution. Results for multiplicity distributions scaled by the mean multiplicity in p+p, p+Au, and d+Au systems at √ s = 200 GeV are shown in Fig. 1, comparing to experimental data (of the uncorrected charged particle multiplicity on the Au going side) in d+Au collisions from the STAR Collaboration [67] and (to the corrected charged particle multiplicity data within |η| < 0.5) in p+p collisions from the UA5 Collaboration [68]. While the statistics for the computed charged hadron multiplicity distribution is limited, one can see that the (scaled) gluon distribution is a good proxy for the final state charged hadron distribution. Furthermore, the shape of the experimental distributions in p+p and d+Au collisions is well reproduced in our framework.
Mean transverse momentum. Being sensitive to the initial system size, equation of state, and the bulk viscosity, the mean transverse momentum p T is an important observable to constrain models of heavy ion (and smaller system) collisions. In Fig. 2 we compare our results for p T of pions, kaons, and protons in d+Au collisions as a function of charged particle multiplicity to experimental data from the STAR Collaboration [67]. The line corresponds to the same parameter set used to describe Au+Au collisions in [47]. The lower edge of the band results from instead setting the initial Π and second order transport coefficients (except for the relaxation times) to zero. Agreement with the experimental data is good for kaons and protons, while pions are overestimated, especially when taking into account initial Π and second order coefficients. We show results for p T obtained in p+Au collisions for comparison. The generally more compact p+Au collisions produce a larger p T at the same multiplicity.
Azimuthal anisotropies. We determine the transverse momentum dependent elliptic and triangular anisotropy from two-particle correlations and compare to experimental data from the PHENIX Collaboration for p+Au, d+Au, and 3 He+Au collisions at √ s = 200 GeV [44]. The experimental v 2 (p T ) is obtained using an event-plane method with the second and third order event-planes determined at forward rapidities in the Au going direction, and using charged hadrons in the mid-rapidity region. Our calculation is boost-invariant and uses the scalar product method, which should best approximate the event-plane method used in the experiment when the event-plane resolution is small [70].
We show our results for all three systems in their respective 0-5% centrality bins and compare to data from the PHENIX Collaboration in Fig. 3. The bands again result from varying the initial Π and second order transport coefficients, with the trends being the same as in the case of p T for p T < 2 GeV. The line is for the parameters constrained in Au+Au collisions [47]. Because of the different spatial geometries of the three systems ε  [44,54] one expects a corresponding ordering of the final state azimuthal momentum anisotropies, if final state response to the geometry is sizable. Indeed, we observe some of the expected ordering of v 2 (p T ) between systems, especially at larger p T . The difference in multiplicity in the three 0-5% centrality classes also plays a role for this result. The computed v 3 is generally greater than the experimental data, except at the lowest transverse momenta. The difference between the p/d+Au v 3 and the one in 3 He+Au collisions is significantly smaller than in the experimental result.
In Fig. 4 we demonstrate that v 2 (p T ) also shows the ordering expected from the geometry when selecting events with the same multiplicity. We also found that when using smoother nucleons (than the ones used and constrained in [71,72]) the difference between d+Au v 2 and p+Au v 2 at the same multiplicity is increased, highlighting a sensitivity to sub-nucleonic fluctuations. Fig. 4 also shows a comparison to PHENIX data, indicating that the relative difference between v 2 (p T ) in p+Au and d+Au collisions is smaller in our calculation than in the experimental data, but comparable given the experimental errors. Results for Au+Au collisions show a smaller v 2 than in d+Au collisions for most values of p T .
We next compute the four-particle cumulants c 2 {4} that reduce non-flow contributions from intrinsic twoand three-particle correlations without the need for a large rapidity gap. In order to facilitate a comparison with the PHENIX data we consider momentum integrated quantities plotted against the charged particle  multiplicity measured at forward rapidity, N FVTX tracks . 2 Our results for all small systems are shown in Fig. 5. We find that the c 2 {4}, that measures the momentum integrated anisotropy, is equal (within statistical errors) in p+Au and d+Au collisions for the same multiplicity. This shows that in our approach the difference observed between p+Au and d+Au for differential v 2 (p T ) seen in Fig. 4 is compensated by the difference in p T shown in Fig. 2. Within errors, c 2 {4} in 3 He+Au collisions also agrees with those in p+Au and d+Au collisions.
While in d+Au collisions c 2 {4} is negative both in the calculation (except for the lowest multiplicities) and the experimental data [73], in p+Au collisions it is positive in the experimental data, but negative in our calculation. The computed c 2 {4} in d+Au has a stronger multiplicity dependence than the experimental data, and is comparable to the AMPT result shown in [69]. The mean value of our result for 0-5% central p+p collisions is positive, but within our statistics it is consistent with zero.
Role of geometry and initial momentum anisotropy. We compute the initial momentum anisotropy defined by which is a proxy for the purely initial state v 2 (of the single-particle distribution) [74]. Note that · is defined without energy density weight. We present results for ε p as a function of multiplicity for different systems together with final v 2 {2} values, which are compared to STAR data [76,77], in Fig. 6 (a). The trend with multiplicity of the initial ε p is opposite to that of v 2 . This is in line with the color domain interpretation of the emergence of initial state anisotropies  . 6. a) Initial momentum anisotropy and final v2 as functions of multiplicity compared to STAR data [75]. STAR determined dN ch /dη using a specific quark Monte-Carlo Glauber model. b) Correlation between the initial ellipticity and final v2. c) Correlation between the magnitude of the initial momentum anisotropy εp and final v2. d) Correlation between the direction of the initial momentum anisotropy (given by ψ p 2 ) and the final v2 (with ψ v 2 2 obtained from the flow vector). [11,12,[78][79][80], where the anisotropy is suppressed in larger systems. The effect of the hydrodynamic final state interactions reverses this trend, leading to increasing v 2 with increasing multiplicity. The initial anisotropy ε p for p+Au collisions is larger than that for d+Au and Au+Au collisions at the same multiplicity, which is also obtained in other CGC calculations [18] (for p+Au and d+Au) and is expected from the color domain interpretation. We note that the direct CGC calculation of two-particle correlations also includes contributions to the anisotropy from quantum interference effects, which are lost when taking T µν and inserting it into hydrodynamics.
For Au+Au collisions we find smaller v 2 than for the "smaller" systems at the same multiplicity, in line with the result in Fig.4. Given that ε p is larger in d+Au than in Au+Au collisions at the same multiplicity (and the eccentricities are very similar in the two systems, with ε 2 being slightly larger in Au+Au collisions), the fact that v 2 (p T ) is larger in d+Au collisions than in Au+Au collisions can be interpreted as an effect of the initial momentum anisotropy. We will show below that the final v 2 can in fact retain information on the initial momentum anisotropy for this multiplicity (Fig. 6 (c) and (d)).
The STAR data [75] have large error bars, but the trend of the ordering between systems is opposite to our result. Because the STAR measurements were performed with a limited pseudorapidity separation of |∆η| > 0.7, it is expected that there are still non-flow contributions in the data. Because the difference between v 2 {2} in d+Au and Au+Au collisions at the same multiplicity could be a signal of the CGC momentum anisotropy, precise data on v n in d/p+Au and peripheral Au+Au collisions for the same kinematic cuts at RHIC are highly desirable.
At LHC, the v 2 in √ s = 5.02 TeV p+Pb collisions was found to be smaller than that in Pb+Pb collisions at the same multiplicity [81][82][83]. This is not inconsistent with our results as at the higher LHC energy hydrodynamics gains in relative importance because higher en-ergy densities lead to longer fireball lifetimes. We will explore this collision energy dependence quantitatively within our framework in the future. Fig. 6 (b) shows the correlation between the spatial ellipticity and the v 2 , demonstrating the expected increasing trend with multiplicity. For p+Au collisions the correlation is stronger than for the other systems at the same multiplicity, probably because the measure of eccentricity works best as a predictor for v 2 for a system that is spatially connected -the other systems at low multiplicity often consist of separated regions that may not interact much with one another during the hydrodynamic evolution. Fig. 6 (c) shows the effect of the initial state momentum anisotropy on the observed v 2 by presenting the correlator of ε p and v 2 as a function of multiplicity. Results for different systems are very similar, all showing a significant correlation for systems with less than 10 charged hadrons per unit rapidity, with the correlation decreasing with increasing multiplicity, as the effects of hydrodynamic response to the initial spatial geometry become dominant. For Au+Au collisions more central than 50%, the magnitude of the initial momentum anisotropy is uncorrelated with the final v 2 .
The correlation between the direction of the initial state momentum anisotropy ψ p 2 and the final elliptic flow is shown in Fig. 6 (d). For small multiplicities the directions of the final flow and initial state momentum anisotropy are strongly correlated, with the correlation decreasing with increasing multiplicity as for the case of the magnitudes. This qualitatively agrees with the recent study using AMPT [84]. All small systems have a finite correlation, indicating that effects of the initial momentum anisotropy should never be neglected. For Au+Au collisions with dN ch /dη > ∼ 100, the flow direction has no correlation with the initial momentum anisotropy anymore. Another interesting feature is that the correlation is stronger in p+Au than in d+Au collisions for the same multiplicity. This can be caused by a) the larger initial momentum anisotropy in p+Au and b) the differences in the hydrodynamic evolution, as system sizes are different in p+Au and d+Au at the same multiplicity.
Effects of initial flow and viscous stress. To further explore the effect of the details in the initial state on observables, in Fig. 7 we show results for (a) dN ch /dη, (b) p T , and (c) v 2 {2} obtained with our standard calculation using the full energy momentum tensor and matching of the initial pressure to three other situations, where we i) exclude the initial Π ii) exclude initial Π and initial π µν , and iii) exclude initial Π, initial π µν and initial flow u µ . We study both 30-40% central Au+Au and 0-5% central p+Au collisions.
As expected, the effects of the changes to the initial state affect p+Au collisions more than the studied Au+Au collisions. The multiplicity varies by maximally 20% in Au+Au collisions, mainly because removing the initial flow removes some of the deposited energy. In p+Au collisions, the maximal modification of dN ch /dη is ∼ 40%. For Au+Au, p T changes by maximally 15%, in p+Au it can change by 30%. 3 Other choices for (ζ/s)(T ), e.g. with a narrower peak around T c , would lead to a reduced τ Π and a weakened effect of the initial Π.
The most significant changes we find for v 2 , which can change by ∼ 35% in Au+Au collisions and by as much as 90% in p+Au. This enormous effect in p+Au is a result of very little flow build-up in the absence of initial flow and presence of a large bulk viscosity. While none of the scenarios, other than using the full T µν , are fully physical, because they do not conserve energy and momentum at the time of switching to hydrodynamics, they give us a good estimate on how much the details of the initial state, other than the geometry, matter in small and larger systems. Our results suggest that for any initial state model, inclusion of initial flow or offequilibrium contributions will affect final observables in small systems. Merely using spatial distributions of energy densities is likely too crude an approximation, for some observables even in Au+Au collisions.
An uncertainty we have not yet discussed stems from the choice of τ init . We found that in typical (30-40%) Au+Au collisions p T , v 2 {2}, and v 3 {2} all change by less than 1% when changing τ init to 0.2 fm. The biggest change was an increase of dN ch /dη by 20%, caused by increased viscous entropy production. On the other hand, 0-5% p+Au collisions show more sensitivity to τ init . Most observables changed between 3 and 15% from the results at τ init = 0.4 fm varying it up or down by 0.2 fm. The two larger changes were an increase of dN ch /dη by 35% and a decrease of v 3 {2} by 22% when going to τ init = 0.2 fm. Conclusions. We have presented calculations for multi-particle correlation observables in systems studied in the RHIC small system scan within a hybrid framework including an initial state from the Color Glass Condensate effective theory, hydrodynamic evolution, and microscopic hadronic transport. All free parameters of the model were constrained by Au+Au collisions previously [47]. Our framework includes CGC initial state momentum anisotropies and final state response to the initial spatial geometry. We demonstrated that the two mechanisms lead to characteristically different multiplicity dependencies of the momentum anisotropy and that final state effects are necessary to describe the qualitative features of the experimentally observed anisotropies.
This work is the first to have quantified the amount that the initial state momentum anisotropy can contribute to observables, if final state interactions are described by realistic hydrodynamic simulations. We have found non-zero correlations between the initial momentum anisotropy magnitude and the charged hadron v 2 , as well as their orientations. The correlation decreases with increasing multiplicity and vanishes in large systems (Au+Au, more central than ∼ 50%). Generally, the details of the initial state, including the initial flow profile and viscous stress tensor, affect the quantitative results for experimental observables strongly.
We overestimate experimental results for v 2 (p T ) and v 3 (p T ) from PHENIX, which could have to do with the fact that our calculations are boost-invariant, and consequently cannot employ exactly the same method as the experiment. It is also possible that a better set of parameters could be found that yields good agreement with both heavy ion and small system collision data.
While the p T dependent v 2 in p+Au and d+Au collisions show the expected ordering from geometry at the same multiplicity, the fact that v 2 {2} in d+Au has a larger value than Au+Au at the same multiplicity is interpreted as a signature of the initial state momentum anisotropy. It does not result from differences in p T , because v 2 (p T ) shows the same trend, and it is not a consequence of different eccentricities, which are close and ordered in the opposite way. Consequently, precise measurements of v 2 for these two systems at RHIC energy and at equal multiplicities would be highly desirable to help unveil a possible signature of the initial state momentum anisotropy.
Further progress will require 3D simulations either using extensions of the IP-Glasma initial state [85,86] or alternative three dimensional initial state models [87], which can also be applied to lower energy collisions, like the ones in the d+Au energy scan [73], but lack the initial state momentum anisotropy emphasized here. Better constraints on QCD transport coefficients, improved descriptions of the intermediate regime between Yang-Mills and hydrodynamic stages (see e.g. [88,89]), as well as improved particlization mechanisms (see e.g. [90]), and inclusion of hydrodynamic fluctuations [91] are potentially important future steps for the description of small system collisions.
Acknowledgments We thank Ron Belmont for his help converting to PHENIX's multiplicity measure. We thank Heikki Mäntysaari and Raju Venugopalan for helpful discussions. B.P.S. and P.T. are supported under DOE Contract No. de-sc0012704. C.S. is supported under DOE Contract No. de-sc0013460. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 and resources of the high performance computing services at Wayne State University. This work is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration.