Light-Heavy Ion Collisions: A window into pre-equilibrium QCD dynamics?

Relativistic collisions of light on heavy ions (p+Au at sqrt(s)=7.7 GeV, p+Au, d+Au,3He+Au at sqrt(s)=62.4 GeV and 200 GeV and p+Pb, 3He+Pb at sqrt(s)=5.02 TeV) are simulated using"superSONIC", a model that includes pre-equilibrium flow, viscous hydrodynamics and a hadronic cascade afterburner. Even though these systems have strong gradients and only consist of at most a few tens of charged particles per unit rapidity, one finds evidence that a hydrodynamic description applies to these systems. Based on these simulations, the presence of a triangular flow component in d+Au collisions at sqrt(s)=200 GeV is predicted to be similar in magnitude to that found in 3He+Au collisions. Furthermore, the v3(p_T) ratio of 3He+Au to d+Au is found to be sensitive to the presence of pre-equilibrium flow. This would imply an experimentally accessible window into pre-equilibrium QCD dynamics using light-heavy ion collisions.


I. INTRODUCTION
Experimental data on d+Au , p+Pb and 3 He+Au collisions seems to indicate the presence of collective flow in an (at least partially) equilibrated system [1][2][3][4][5]. While alternative explanation have been put forward [6], the result from measurements of multi-particle cumulants are strongly suggestive of a hydrodynamic origin of this collectivity [7]. This study is motivated by the following, interrelated, questions: 1. Does hydrodynamics apply to systems created in light-heavy ion collisions?
2. At which system size or collision energy does hydrodynamics break down?
3. What experimental observables would confirm or rule out a theorist's answer to the above questions?
In order to obtain answers for these questions, I use the new simulation package "superSONIC", which is an eventby-event generalization of the SONIC model [8], including pre-equilibrium flow, relativistic viscous hydrodynamics as well as a hadronic cascade afterburner (details are given below). Using a Monte-Carlo Glauber model for generating event-by-event initial conditions, superSONIC can be used to simulate a range of different collision systems at different collision energies, such that results are directly comparable. Comparison to data will be performed wherever experimental information is available.
Quite a number of theoretical studies on light-heavy ion collisions have appeared in the recent literature, and in the following similarities and notable differences of these works with respect to the present study are highlighted.
In Ref. [9], Piotr Bozek studied p+Pb and d+Pb collisions at √ s = 4.4 TeV and √ s = 3.11 TeV, respectively, using Monte-Carlo Glauber initial conditions followed by a 3+1d viscous hydrodynamics evolution and a hadronic cascade afterburner. He successfully predicted the large flow signal in p + Pb collisions at √ s = 5.02 TeV.
In Ref. [10], Nagle et al. studied p+Au , d+Au and 3 He+Au collisions at √ s = 200 GeV and p+Pb collisions at √ s = 5.02 TeV using Monte-Carlo Glauber initial conditions followed by a 2+1 viscous hydrodynamics evolution and a hadronic cascade afterburner ("SONIC", the predecessor of the package used in the present study). Nagle et al. found results consistent with Ref. [9], and first proposed 3 He+Au collisions as an interesting handle on QCD transport properties. Elliptic and triangular flow in 3 He+Au at √ s have since been measured experimentally [11]. In Refs. [12,13], Schenke and Venugopalan studied light-heavy ion collisions using IP-Glasma initial conditions, followed by a 2+1d viscous hydrodynamic evolution, but no hadronic cascade afterburner. Their study is similar to the present one in that the sudden matching from the Glasma evolution to hydrodynamics includes non-vanishing pre-equilibrium flow. However, given that the Glasma evolution never drives the system to equilibrium, one can expect the results to be strongly dependent on the Glasma-hydro switching time. In their study, Schenke and Venugopalan found sizable flow components v 2 , v 3 for d+Au , 3 He+Au collisions (in agreement with experiment), but very little flow for p+Pb collisions (in disagreement with experiment).
In Ref. [14], Kozlov et al. studied p+Pb collisions at √ s = 5.02 TeV using 3+1d viscous hydrodynamics, but without any pre-equilibrium flow or hadronic afterburner. They implemented initial conditions based on a Monte-Carlo Glauber model supplemented with negative binomial energy fluctuations, which turn out to be very important for describing rare "high-multiplicity" events. Kozlov et al. found that they could successfully describe experimental flow data in p+Pb using their model. Compared to some of the studies listed above, the present study suffers from the drawback of only simulating boostinvariant (2+1d) dynamics, which clearly is not a good approximation to the longitudinal dynamics for light-heavy ion collisions. However, many of the conclusions in the present article will be based on ratios of flow observables in the transverse plane, so that one can expect conclusions to not be dominated by longitudinal artifacts. Nevertheless, it would be good if the present results could be checked by a full 3+1d calculation in the future. Also, the present study does not include negative binomial energy fluctuations, on the basis that these fluctuations are most relevant only for rare high-multiplicity events. In this article, the emphasis is on central (impact parameter b < 2 fm/c), but not high-multiplicity selected events. Finally, the present study is based on a geometric Monte-Carlo Glauber model [15] rather than the IP-Glasma model [16], mostly to offer a baseline result using transparent ingredients. If it would turn out that some experimental result could not be described by the present approach based on Monte-Carlo Glauber, but can be described using IP-Glasma initial conditions, this could be regarded as experimental evidence for Color-Glass-Condensate dynamics in QCD.
On the other hand, the present simulation package for the first time combines pre-equilibrium flow, viscous hydrodynamics and hadron cascade dynamics in an event-by-event study that is applicable to different collision geometries. In this sense, it is the most realistic description currently available, and the direct comparison to experimental data may therefore offer interesting clues about the nature of hydrodynamics and transport in strongly coupled QCD.
The remainder of the article is structured as follows: in Sec.II, a brief description of the components of superSONIC is given, followed by results for p+Au ,d+Au , 3 He+Au and p+Pb collisions at various energies in Sec.III, and the conclusions in Sec.IV.

A. Early Stage: Initial Conditions and Pre-equlibrium
In principle, obtaining information about real-time evolution in relativistic ion collisions is a problem that can be well defined in QCD. However, currently there are no known techniques to solve this problem, and thus obtain information about the pre-equilibrium stage in these systems. Nevertheless, a number of approximate approaches exits. Among these, there are weak-coupling inspired methods that rigorously apply to QCD in the asymptotically high energy limit (cf. [17][18][19][20][21][22][23][24][25]), and strong-coupling methods that rigorously apply to a certain class of non-abelian gauge theories (but not QCD) in the limit of large number of colors and large t'Hooft coupling (cf. [26][27][28][29][30][31][32][33][34][35]).
From a modeling perspective, what is needed from any of these methods is information about when the system behaves approximately hydro-dynamically (equilibration time) as well as the values of the hydrodynamic degrees of freedom (hydrodynamic initial conditions). Other than the hydrodynamic starting time [22,25], results for hydrodynamic initial conditions are presently available only from the strong-coupling method [32,34,35], so this is what will be used in the present study.
In more detail, from numerical relativity simulations of space-times modeling the relativistic collision of "ions" in N = 4 SYM it has been found that the radial fluid velocity is proportional to the gradient of the initial density distribution [32,35,36]. In this study, this finding is promoted to the fluid velocity in all directions, such that where τ ≡ √ t 2 − z 2 and R 2 (x) is the product of the particle densities at the time of collision (τ = 0). Note that this result is consistent with Ref. [37], but the pre-factor is non-trivial (cf. [8]). The result (1) has the added benefit that final particle momenta were found to be almost insensitive to the choice of the hydrodynamic starting time τ 0 in heavy-ion collisions [32].
For the energy-density ǫ(τ 0 , x), initial conditions for each event are constructed as follows. Using Woods-Saxon distribution functions for the heavy ions such as Au, Pb [38,39], the Hulthen wavefunction for the deuteron (cf. [40]) and realistic calculations for the 3 He wavefunction [41], probability distributions of the nucleons within the nuclei of interest (cf. [10]) are obtained. Using a Monte-Carlo Glauber [15], these probability distributions are mapped to positions of individual nucleons in the transverse (x, y) plane on an event-by-event basis implementing a hardcore repulsive potential of radius 0.4 fm between nucleons. The positions of nucleons undergoing at least one inelastic collisions are recorded ("participants") and converted into a density function R 2 (x) by assuming that each participant contributes equally as a Gaussian with a width of w = 0.4 fm (to match the RMS radius of a single nucleon). This fixes the initial fluid velocities via Eq. (1), and the energy density for hydrodynamics is assumed to be given as with E 0 an overall constant (dependent on τ 0 , collision energy and collision system) that is related to the total multiplicity of the event. Typically 100 event initial conditions are generated for each collision system. While the procedure to obtain the energy-density for hydrodynamics in Eq. (2) is fairly standard in relativistic ion collision literature, the presence of the pre-equilibrium flow in Eq. (1) is relatively new. For this reason, also results without pre-equilibrium flow will be presented, and potential experimental ways to determine if pre-equilibrium flow is present in light-heavy ion collisions will be discussed. Other than the energy density and flow velocity field, relativistic hydrodynamic formulations also require the initial value for the shear and bulk stress tensors, which will be set to zero for simplicity. In heavy ion collisions, this assumption is harmless, as it has very little influence on final results cf. [32,42]. However, given the short life-time of systems created in light-heavy ion collisions compared to heavy ion collisions, this assumption should be carefully revisited in follow-up studies.

B. Thermal Stage: Viscous Hydrodynamics
In modern definitions, hydrodynamics is understood to be an effective theory of energy conservation at long wavelength [43,44]. Hydrodynamics is a good approximation of the bulk dynamics as long as higher order gradient corrections do not (strongly) change the leading order results. In standard nomenclature, viscous effects come into the energy-momentum tensor at first order in gradients (zeroth order would be ideal hydrodynamics). However, since all known consistent formulations for relativistic dissipative fluid dynamics are second order in gradients, this suggests a convenient handle on the effects of higher order gradients: the value of second-order transport coefficients, foremost the shear viscous relaxation time τ π .
The view that is adopted here is the following: while first-order gradient effects in a system may be large (viscous effects sizable), a (viscous) hydrodynamic description of the system may still be quantitatively reliable as long as higher-order gradient corrections are small compared to first order gradients. This view is informed by direct simulations of strongly coupled quantum field theories out of equilibrium where it has been shown that viscous hydrodynamics offers a reliable description even in regions where first-order (viscous) corrections to the ideal energy-stress tensor are approaching 100 percent [32]! To test if second-order corrections are small compared to first-order gradients, the value of τ π (more precisely the ratio C η = τπ(ǫ+P ) 2η ) is varied by 50 percent around the reference value C η = 3, thereby generating a "systematic" error estimate of the applicability of viscous hydrodynamics. If the hydrodynamic gradient approximation was breaking down, one would expect second-order gradient terms to be as important as first-order gradient terms. Thus, final results for particle spectra should vary considerably when changing the strength of second order terms (via the value of C η ) by 50 percent. Conversely, if final results showed very little dependence on the value for C η (as is the case for heavy ion collisions, cf. [42]), this could be considered evidence that a hydrodynamics was well applicable to such systems. In this sense, the reliability of hydrodynamics as a approximation to the system dynamics can be quantified and expressed as an error band generated by varying C η , which is the strategy adopted in the following.
Besides the second-order transport coefficients, the hydrodynamic evolution will also depend on choices for the (temperature-dependent) shear and bulk viscosity coefficients as well as the speed of sound (via the equation of state). Again, for simplicity constant values are used for the ratios of shear viscosity over entropy density η/s and the bulk viscosity is set to identically zero (ζ = 0). For the equation of state, a parametrization of lattice QCD data, given in Ref. [45], is employed. All of these choices should be revisited in a more detailed study.
The hydrodynamic evolution is performed using the open-source code VH2+1 (version 1.9) [42], adapted by smearing to prohibit code instabilities in the strong-gradient regions frequently encountered in event-by-event viscous hydrodynamics. The smearing (optimized from a version used before in Ref. [10]) is performed by replacing low-energy density values by an average over nearest-neighbor cells, and it is only implemented at temperatures below the QCD phase transition (typically below 150 MeV). I have tested that final results are not sensitive to the details of the smearing implementation. The hydrodynamic evolution solves the equations of motion on a square lattice with area 20 2 fm 2 and 200 2 grid points (lattice spacing 0.1 fm) and a time step of one hundredth of the lattice spacing. I have tested that the results are unchanged when simulating the same volume with a lattice spacing of 0.05 fm.
During the hydrodynamic evolution, cells that cool below a certain switching temperature T SW = 170 MeV are monitored. This condition defines a switching hyper-surface on which information about temperature, flow velocity, dissipative stress tensors as well as the normal vector of the hyper-surface are recorded for each cell. This information will provided the initial condition for the late-stage hadron cascade simulation.

C. Late Stage: Hadron Interactions in a Hadron Cascade
For the late-stage hadron interactions,the hadronic cascade code B3D [46] is used. Using the hyper-surface information to boost to the rest frame of each cell, the cascade is initialized with particles in the rest frame drawn from a Boltzmann distribution at a temperature T SW with modifications of the momentum distribution to include deformations from viscous stress tensors (see [47] for details). B3D includes hadron resonances in the particle data book up to masses of 2.2 GeV, which interact via simple s-wave scattering with a constant cross-section of 10 mb as well as scattering through resonances (modeled as a Breit-Wigner form). Once the resonances have stopped interacting, one can obtain final charged hadron multiplicities dN ch dY , mean charged particle momentum p T and flow coefficients v n (p T ) for n ≥ 1 from summing over individual particles with momenta p. Specifically, where ∆ T = 80 MeV, ∆ Y = 2 are the width of bins for particle p T and rapidity Y , respectively. Note that since the cascade is applied to a boost-invariance case, the large ∆ Y value is of no significance. In practice, a sum over both particles and anti-particles and division of the spectra by two is performed, in order to increase statistics. For every hydrodynamic evolution event, 100,000 B3D events are run to increase statistics. In doing so, the sums in the definition of v n above are extended over all B3D events, thereby explicitly ignoring fluctuations arising from hadronic decays. After thus obtaining results dN ch 2πpT dY dpT and v n (p T ) for each hydrodynamic event, an event average to obtain   the event-by-event mean and event-by-event fluctuation is performed, the latter of which is recast into a statistical error bar on the mean. The results from this procedure are reported on in the following. The combined simulation package of the early-stage, thermal stage and late stage evolution thus described above will be referred to as "superSONIC" in the remainder of this article.

III. RESULTS
A summary of systems that were simulated is given in Tab. I. The first column in this table gives the system configuration, the second the collision energy, the third the inelastic cross-section (from Ref. [49]) used in the Monte-Carlo Glauber. All Monte-Carlo Glauber events were generated for "central" collisions by imposing an impact parameter b < 2 fm, loosely corresponding to the 0-5% most central collisions. The fourth column in Tab. I refers to the mean number of participants obtained by averaging over 100 random Monte-Carlo Glauber events. The fifth column specifies the constant value of η/s chosen in the hydrodynamic simulations, which were all started at an initial time τ 0 = 0.25 fm/c with or without pre-equilibrium flow according to the sixth column of Tab. I. For all runs, the energy scale factor E 0 that was chosen in order to match measured or expected values ([50]) for the charged particle multiplicity, reported in column seven and eighth of Tab. I. (Note that the calculated dN ch dY is converted to experimentally measured pseudo-rapidity distribution dN ch dη by dividing by 1.1). For experimentally measured quantities, also the centrality class for the measurement as well as the corresponding reference is reported. Finally, the last two columns in Tab. I give the mean charged particle transverse momentum in the superSONIC simulation compared to experimental data where available.
As one can see from Tab. I, simulated particle multiplicities range from dN ch dη ≃ 74 down to dN ch dη ≃ 4.5. This implies that systems with very few particles are being simulated and one generally expects hydrodynamics to be less applicable to these fewer-body systems. Also, note that there is a clear change in the mean charged particle momentum for systems with compared to without pre-equilibrium flow. This is not too surprising given that systems created in lightheavy ion collisions live comparatively shorter than those created in heavy ion collisions, thus making light-heavy ion collision systems more sensitive to pre-equilibrium conditions. While other parameters (viscosity, choice of switching temperature T sw ) also affect the particle mean momentum, the strong difference between results with and without pre-equilibrium flow could serve as important discriminatory tool that is easy to implement experimentally.
In figure 2, results for the flow coefficients v n , n=2,3,4 minus v 5 are shown for p+Pb and 3 He+Pb collisions at √ s = 5.02 TeV per nucleon pair (LHC energies). Results are reported as a difference with respect to calculated v 5 (the highest harmonic calculated) because finite statistics from the 100,000 B3D runs start to pollute the results at high p T and the true v 5 (p T ) can reasonably be expected to be consistent with zero. Hence the calculated v 5 (p T ) is a good measure of the statistical error, and can be used to subtract the statistical fluctuation for the other flow harmonics. In principle, this procedure could be made unnecessary by rerunning all the simulations with at least 10 6 B3D runs per hydro event, which is left for future work.
For comparison, results with and without pre-equilibrium flow are shown. Boxes indicate uncertainty arising from both statistical fluctuations as well as systematic errors, the latter of which are quantified by performing simulations at different values of C η = 2 − 3. From these plots, the first finding is that the hydrodynamic uncertainty range for light-heavy ion collisions at √ s = 5.02 TeV is rather small compared to the overall magnitude of the flow for GeV is almost as large as for 3 He+Au , while without the presence of pre-equilibrium flow v 3 in p+Au is about half as large as in 3 He+Au . This suggests that v 3 in light-heavy ion collisions at √ s = 200 GeV could be a good experimental handle on pre-equilibrium QCD dynamics. Unlike the situation found for √ s = 5.02 TeV, simulation results predict v 4 to be consistent with zero except for p+Au ,d+Au and 3 He+Au where some small, non-vanishing v 4 is found. Since v 4 is found to be so small at √ s < 200 GeV, it is likely to be hard to measure, and hence it will be discounted as a probe for pre-equilibrium dynamics in the following.
What is striking about the results shown in Fig. 3 is the magnitude of v 2 (and even v 3 in the case of pre-equilibrium flow) predicted in p+Al collisions at √ s = 200 GeV. The system created in these collisions consists of an average multiplicity of only dN ch dη ∼ 5.4, yet superSONIC results exhibit a clear flow response much larger than the estimated uncertainty band for the applicability of hydrodynamics. According to the criterion defined above, I predict that systems created in central p+Al collisions at √ s = 200 GeV can be quantitatively described using viscous hydrodynamics, and that a v 2 flow component in these systems can be expected to be similar to that for p+Au collisions at √ s = 200 GeV. In Fig. 4  collision energies below √ s = 7.7 GeV, the zero chemical potential lattice equation of state employed in superSONIC is clearly no longer applicable, so studying proton-nucleus collisions at energies below √ s = 7.7 GeV is not feasible with the current approach.
While a complete break-down of hydrodynamics is not observed in superSONIC simulations, one does observe a gradual break-down of hydrodynamics. That is, starting with p+Pb collisions at √ s = 5.02 TeV and lowering the collision energy down to √ s = 7.7 GeV for p+Au collisions, first v 4 and then v 3 first drop and then collapse to values consistent with zero for all momenta considered. In the same fashion, one does also observe a clear drop in v 2 as a function of lowering √ s, even though at √ s = 7.7 GeV v 2 results have not (yet) collapsed to zero. Hydrodynamics is breaking down, but is has not broken down completely for p+Au collisions at √ s = 7.7 GeV. My interpretation of this finding is that the question of applicability of hydrodynamics or collectivity in small systems does not have a yes/no answer, but rather should be thought of as a gradual process similar to the confinement/deconfinement cross-over transition in QCD, where the value of the critical temperature is also dependent on the observable one considers.
In Fig. 6, the flow response v 2 , v 3 is compared for different collision systems in an attempt to quantify the presence of pre-equilibrium flow. Similar to the results shown in figures above, one finds that v 2 is fairly insensitive to preequilibrium flow, regardless of the collision energy. However, v 3 turns out to be a good indicator for the presence of pre-equilibrium flow in different collision systems at √ s ≤ 200 GeV. Specifically, while v 3 in 3 He+Au collisions has been measured at √ s = 200 GeV, a measurement of v 3 in d+Au collisions at that same energy could serve as an experimental handle on the presence of pre-equilibrium flow. Also, note that lower collision energies (such as √ s = 62.4 GeV) would be even better suited for an experimental test of pre-equilibrium flow since the systems created in these collisions live comparatively shorter, and are thus more sensitive to out-of-equilibrium QCD transport dynamics. Conversely, it would be harder to probe these pre-equilibrium transport effects at √ s = 5.04 TeV, both because 3 He+Pb collisions at the LHC are not currently planned, and because 3 He+Pb to p+Pb v 3 ratios are only sensitive to pre-equilibrium flow at p T < ∼ 0.5 GeV. Another point worth noting is that the more limited temperature range encountered in light-heavy ion collisions as compared to heavy-ion collisions, especially at RHIC energies. Fig. 1 shows the maximal temperature encountered in various collision systems from √ s = 7.7 GeV to √ s = 5.02 TeV. Also shown in Fig. 1 is the typical (event-by-event mean) starting temperature, which is a measure of the maximal temperature the average system starts out with in these collisions (dots). One finds that hydrodynamic evolution in p+Pb collisions probes and thus averages over transport properties (such as η/s and ζ/s) over a large temperature range T < ∼ 0.53 GeV (< T >= 0.4 GeV, while by contrast p+Au collisions at √ s = 7.7 GeV only probe T < 0.33 GeV (< T >= 0.25 GeV). Hence p+Au collisions at √ s = 7.7 GeV are mostly sensitive to QCD transport properties at T < 0.25 GeV, which could be a key experimental handle on probing the temperature dependence on e.g. η/s. Finally, results on HBT Radii for selected collision systems are reported in Tab. I in Figures 7,8 and 9.

IV. SUMMARY AND CONCLUSIONS
In this work, central light on heavy ion collisions from √ s = 7.7 GeV to √ s = 5.02 TeV were simulated using superSONIC, a model that combines pre-equilibrium flow, relativistic viscous hydrodynamics and a hadron cascade afterburner. By varying the strength of the second order transport coefficients, one could quantify to which extent viscous hydrodynamics offers a reliable description of these systems . It was found that even in p+Au collisions at √ s = 7.7 GeV, a sizable collective flow component v 2 much larger than the systematic hydrodynamic uncertainty is present. Thus, viscous hydrodynamics is still applicable to describe v 2 in the systems created in these small, low-energy, few-body collisions. However, there is evidence that hydrodynamics is breaking down gradually as √ s is lowered from 5.02 TeV to √ s = 7.7 GeV. Specifically, first v 4 , then v 3 and finally v 2 start to decrease and eventually v 3 and v 4 become consistent with zero as the collision energy is lowered. The question of whether hydrodynamics applies to describe light-heavy ion collisions therefore cannot be answered by a simple yes or no, but depends on the quantity in question. In future work, it would be interesting to perform simulations at collision energies below √ s = 7.7 in order to potentially observe also v 2 collapse to zero in p+Au collisions.
For d+Au collisions at √ s = 200 GeV superSONIC simulations predicted v 3 results that were of the same order of magnitude as those measured in 3 He+Au . Since similar results are observed in other theoretical models for d+Au collisions at √ s = 200 GeV, I predicted that v 3 should be observable in experiment. Furthermore, the sensitivity of flow results to the presence of pre-equilibrium flow was studied. It was found that low-energy collisions were most sensitive to the presence of pre-equilibrium flow and pointed out that the v 3 ratio of 3 He+Au to d+Au could provide an experimental handle on pre-equilibrium QCD dynamics. Moreover, simulation results clearly show that light-heavy ion collision systems probe a temperature window focused around the QCD phase transition temperature. Thus, a combination of simulation results and experimental data for these systems would offer a promising handle on the temperature dependence on QCD transport properties.
Many aspects of this study are amenable to improvements in future work, but despite the current set of limitations and approximations, the findings of this work could hopefully demonstrate the usefulness of studying light-heavy ion collisions to learn about QCD. FIG. 5. Flow harmonics vn(pT ) − v5(pT ) for unidentified charged particles (unid) for n = 2, 3, 4 from superSONIC in protonnucleus collisions, with and without pre-equilibrium flow. Boxes indicate combined statistic and estimated systematic error for hydrodynamics (latter from varying Cη = 2 − 3). Y-axis is same scale on all plots. The p+Pb system, which has the highest multiplicity (see. Tab. I), shows non-vanishing flow components up to n = 4, but as the multiplicity is decreased, first v4, then v3 and eventually also v2 start to decrease and (in the case of v4) eventually become consistent with zero.