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

TL;DR

This paper uses simulations to explore whether hydrodynamics can describe small light-heavy ion collision systems, revealing potential experimental signatures of pre-equilibrium QCD dynamics through flow measurements.

Contribution

It demonstrates that hydrodynamic models can apply to small systems and predicts observable flow signatures sensitive to pre-equilibrium dynamics in light-heavy ion collisions.

Findings

Triangular flow in d+Au at 200 GeV is comparable to 3He+Au.

The v3(p_T) ratio between 3He+Au and d+Au is sensitive to pre-equilibrium flow.

Hydrodynamics may describe systems with only a few tens of particles.

Abstract

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.

Figures

Tables

Table 1. Table 1: Comparison of superSONIC runs (“th”) with available experimental data (“exp”). For all model runs w = 0.4 𝑤 0.4 w=0.4 fm, τ 0 = 0.25 subscript 𝜏 0 0.25 \tau_{0}=0.25 fm/c and T sw = 0.17 subscript 𝑇 sw 0.17 T_{\rm sw}=0.17 GeV as well as “central collisions” with impact parameter b < 2 𝑏 2 b<2 fm/c are chosen. The effect of changing C η = 2 − 3 subscript 𝐶 𝜂 2 3 C_{\eta}=2-3 is contained in the reported theoretical error estimates. Mean particle p T subscript 𝑝 𝑇 p_{T} is for unidentified charged hadrons. For the experimental data both the centrality and the reference are reported.

System	$\sqrt{s}$ [GeV]	${\sigma }_{NN}$[mb]	$⟨N_{\mathrm{part}}⟩$	$\eta /s$	preflow?	$\frac{d{N}_{\mathrm{ch}}}{d\eta }$(th)	$\frac{{\mathrm{𝐝𝐍}}_{\mathrm{ch}}}{𝐝\eta }$(exp)	$⟨p_T⟩$ [GeV] (th)	$⟨p_T⟩$ [GeV] (exp)
p+Au	7.7	32	7.4	0.08	yes	4.5		0.550(6)
p+Au	7.7	32	7.4	0.08	no	4.6		0.497(8)
p+Au	62.4	36	8.6	0.08	yes	8		0.588(5)
p+Au	62.4	36	8.6	0.08	no	8.4		0.528(6)
p+Al	200	42	4.7	0.08	yes	5.4		0.575(6)
p+Al	200	42	4.7	0.08	no	5.5		0.515(7)
p+Au	200	42	9.6	0.08	yes	10		0.604(5)
p+Au	200	42	9.6	0.08	no	10		0.537(5)
d+Au	62.4	36	14.9	0.08	yes	15		0.570
d+Au	200	42	17.5	0.08	yes	20	$20.3\pm 1.7$ (0-5%)Huang and Mitchell (2014)	0.576(1)	0.554 (0-20%)Adare et al. (2013b)
d+Au	200	42	17.5	0.08	no	20	$20.3\pm 1.7$ (0-5%)Huang and Mitchell (2014)	0.523(5)	0.554 (0-20%)Adare et al. (2013b)
3He+Au	62.4	36	21	0.08	yes	21		0.557(1)
3He+Au	62.4	36	21	0.08	no	21		0.509(4)
3He+Au	200	42	24	0.08	yes	$27\pm 1$		0.567(2)
3He+Au	200	42	24	0.08	no	$27\pm 1$		0.520(3)
p+Pb	5020	70	15.2	0.16	yes	$39\pm 1$	$\mathrm{𝟑𝟓}\pm 0.5$ (0-20%)Abelev et al. (2013)	0.716(6)
p+Pb	5020	70	15.2	0.08	no	$38\pm 1$	$\mathrm{𝟑𝟓}\pm 0.5$ (0-20%)Abelev et al. (2013)	0.623(3)
3He+Pb	5020	70	32.4	0.16	yes	74		0.676

Equations

$$\vec{v}(\tau ,\mathbf{x})=-\frac{\tau }{3.0}\vec{\nabla }\ln R^2(\mathbf{x})$$

$$ϵ({\tau }_0,\mathbf{x})=E_0{R}^2(\mathbf{x}),$$

$$\frac{d{N}_{\mathrm{ch}}}{2\pi p_{T}dYd{p}_T}=\frac{{\sum }_{\mathrm{in}{\mathrm{p}}_{\mathrm{T}}\mathrm{bin}}^{ch.\mathrm{particles}}}{2\pi p_T{\Delta }_T{\Delta }_Y},\frac{d{N}_{\mathrm{ch}}}{dY}={\int }_0^{\infty }d{p}_T\frac{d{N}_{\mathrm{ch}}}{dYd{p}_T},⟨p_T⟩=\frac{{\int }_0^{\infty }d{p}_T{p}_T\frac{d{N}_{\mathrm{ch}}}{dYd{p}_T}}{\frac{d{N}_{\mathrm{ch}}}{dY}}$$

$$\displaystyle|v_{n}|(p_{T})=\sqrt{s_{n}(p_{T})^{2}+c_{n}(p_{T})^{2}}\,,\quad\left(\begin{array}[]{cc}s_{n}(p_{T})\\ c_{n}(p_{T})\end{array}\right)=\frac{\sum^{\rm ch.\ particles}_{\rm in\ p_{T}\ bin}\left(\begin{array}[]{cc}\sin(n\phi))\\ \cos(n\phi)\end{array}\right)}{\sum^{\rm ch.\ particles}_{\rm in\ p_{T}\ bin}}\,,\quad\phi\equiv\arctan\left(\frac{p_{y}}{p_{x}}\right)\,,$$

Full text

I Introduction

Experimental data on d+Au , p+Pb and 3He+Au collisions seems to indicate the presence of collective flow in an (at least partially) equilibrated system Abelev et al. (2013); Aad et al. (2013); Adare et al. (2013a); Chatrchyan et al. (2013); Huang (December 2014). While alternative explanation have been put forward Dumitru et al. (2014), the result from measurements of multi-particle cumulants are strongly suggestive of a hydrodynamic origin of this collectivity CMS (2014).

This study is motivated by the following, interrelated, questions:

Does hydrodynamics apply to systems created in light-heavy ion collisions? 2. 2.

At which system size or collision energy does hydrodynamics break down? 3. 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 event-by-event generalization of the SONIC model Habich et al. (2015), 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. Bozek (2012), Piotr Bozek studied p+Pb and d+Pb collisions at $\sqrt{s}=4.4$ TeV and $\sqrt{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+{\rm Pb}\,$ collisions at $\sqrt{s}=5.02$ TeV.

In Ref. Nagle et al. (2014), Nagle et al. studied p+Au , d+Au and 3He+Au collisions at $\sqrt{s}=200$ GeV and p+Pb collisions at $\sqrt{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. Bozek (2012), and first proposed 3He+Au collisions as an interesting handle on QCD transport properties. Elliptic and triangular flow in 3He+Au at $\sqrt{s}$ have since been measured experimentally Huang and Mitchell (2014).

In Refs. Schenke and Venugopalan (2014a, b), 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 ,3He+Au collisions (in agreement with experiment), but very little flow for p+Pb collisions (in disagreement with experiment).

In Ref. Kozlov et al. (2014), Kozlov et al. studied p+Pb collisions at $\sqrt{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 boost-invariant (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 Miller et al. (2007) rather than the IP-Glasma model Schenke et al. (2012), 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 , 3He+Au and p+Pb collisions at various energies in Sec.III, and the conclusions in Sec.IV.

II Methodology

II.1 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. Mrowczynski (1993); Arnold et al. (2005); Rebhan et al. (2005); Romatschke and Venugopalan (2006); Kurkela and Moore (2011); Epelbaum and Gelis (2013); Berges et al. (2014); Strickland (2013); Kurkela and Lu (2014)), 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. Nastase (2005); Kovchegov and Taliotis (2007); Amsel et al. (2008); Grumiller and Romatschke (2008); Chesler and Yaffe (2011); Wu and Romatschke (2011); van der Schee et al. (2013); Casalderrey-Solana et al. (2014); Bantilan and Romatschke (2014); Chesler and Yaffe (2015)).

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 Epelbaum and Gelis (2013); Kurkela and Lu (2014), results for hydrodynamic initial conditions are presently available only from the strong-coupling method van der Schee et al. (2013); Bantilan and Romatschke (2014); Chesler and Yaffe (2015), 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 ${\cal N}=4$ SYM it has been found that the radial fluid velocity is proportional to the gradient of the initial density distribution Romatschke and Hogg (2013); van der Schee et al. (2013); Chesler and Yaffe (2015). In this study, this finding is promoted to the fluid velocity in all directions, such that

[TABLE]

where $\tau\equiv\sqrt{t^{2}-z^{2}}$ and $R^{2}({\bf x})$ is the product of the particle densities at the time of collision ($\tau=0$). Note that this result is consistent with Ref. Vredevoogd and Pratt (2009), but the pre-factor is non-trivial (cf. Habich et al. (2015)). 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 $\tau_{0}$ in heavy-ion collisions van der Schee et al. (2013).

For the energy-density $\epsilon(\tau_{0},{\bf x})$, initial conditions for each event are constructed as follows. Using Woods-Saxon distribution functions for the heavy ions such as ${\rm Au,Pb}$ Jager et al. (1974); Vries et al. (1987), the Hulthen wavefunction for the deuteron (cf. Adare et al. (2014a)) and realistic calculations for the ${}^{3}{\rm He}$ wavefunction Carlson and Schiavilla (1998), probability distributions of the nucleons within the nuclei of interest (cf. Nagle et al. (2014)) are obtained. Using a Monte-Carlo Glauber Miller et al. (2007), these probability distributions are mapped to positions of individual nucleons in the transverse ($x,y$) plane on an event-by-event basis implementing a hard-core 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}({\bf 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

[TABLE]

with $E_{0}$ an overall constant (dependent on $\tau_{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. Luzum and Romatschke (2008); van der Schee et al. (2013). 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.

II.2 Thermal Stage: Viscous Hydrodynamics

In modern definitions, hydrodynamics is understood to be an effective theory of energy conservation at long wavelength Baier et al. (2008); Bhattacharyya et al. (2008). 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 $\tau_{\pi}$.

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 van der Schee et al. (2013)!

To test if second-order corrections are small compared to first-order gradients, the value of $\tau_{\pi}$ (more precisely the ratio $C_{\eta}=\frac{\tau_{\pi}(\epsilon+P)}{2\eta}$) is varied by 50 percent around the reference value $C_{\eta}=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_{\eta}$) by 50 percent. Conversely, if final results showed very little dependence on the value for $C_{\eta}$ (as is the case for heavy ion collisions, cf. Luzum and Romatschke (2008)), 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_{\eta}$, 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 $\eta/s$ and the bulk viscosity is set to identically zero ($\zeta=0$). For the equation of state, a parametrization of lattice QCD data, given in Ref. Borsanyi et al. (2014), 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) Luzum and Romatschke (2008), 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. Nagle et al. (2014)) 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}\ {\rm 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_{\rm 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.

II.3 Late Stage: Hadron Interactions in a Hadron Cascade

For the late-stage hadron interactions,the hadronic cascade code B3D Novak et al. (2014) 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_{\rm SW}$ with modifications of the momentum distribution to include deformations from viscous stress tensors (see Pratt and Torrieri (2010) 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 $\frac{dN_{\rm ch}}{dY}$, mean charged particle momentum $\langle p_{T}\rangle$ and flow coefficients $v_{n}(p_{T})$ for $n\geq 1$ from summing over individual particles with momenta ${\bf p}$. Specifically,

[TABLE]

where $\Delta_{T}=80$ MeV, $\Delta_{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 $\Delta_{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 $\frac{dN_{\rm ch}}{2\pi p_{T}dYdp_{T}}$ 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. 1. The first column in this table gives the system configuration, the second the collision energy, the third the inelastic cross-section (from Ref. Aad et al. (2011)) 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. 1 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 $\eta/s$ chosen in the hydrodynamic simulations, which were all started at an initial time $\tau_{0}=0.25$ fm/c with or without pre-equilibrium flow according to the sixth column of Tab. 1. For all runs, the energy scale factor $E_{0}$ that was chosen in order to match measured or expected values (Alver et al. (2011)) for the charged particle multiplicity, reported in column seven and eighth of Tab. 1. (Note that the calculated $\frac{dN_{\rm ch}}{dY}$ is converted to experimentally measured pseudo-rapidity distribution $\frac{dN_{\rm ch}}{d\eta}$ 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. 1 give the mean charged particle transverse momentum in the superSONIC simulation compared to experimental data where available.

As one can see from Tab. 1, simulated particle multiplicities range from $\frac{dN_{\rm ch}}{d\eta}\simeq 74$ down to $\frac{dN_{\rm ch}}{d\eta}\simeq 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 light-heavy 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_{\rm 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 3He+Pb collisions at $\sqrt{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_{\eta}=2-3$. From these plots, the first finding is that the hydrodynamic uncertainty range for light-heavy ion collisions at $\sqrt{s}=5.02$ TeV is rather small compared to the overall magnitude of the flow for $v_{2},v_{3},v_{4}$, thus providing solid theoretical support to the notion of true hydrodynamic behavior in light-heavy ion collisions. Furthermore, one finds that simulated flow coefficients are in overall agreement with experimental data, where available, when suitably adjusting the constant simulated shear viscosity over entropy density. Note that, somewhat surprisingly, it is hard to distinguish the cases of $\eta/s=0.16$ with pre-equilibrium flow and $\eta/s=0.08$ without pre-equilibrium flow using experimental data for $v_{2},v_{3},v_{4}$. For this reason alone, it would be highly desirable to have an experimental handle on the presence of pre-equilibrium flow. Finally, comparing simulations of 3He+Pb and p+Pb collisions at LHC energies one finds that one could expect $v_{2}$ to be 50 percent higher in 3He+Pb collisions than in p+Pb collisions, while the overall magnitude for $v_{3},v_{4}$ would be approximately the same. This presumably points to the fact that in 3He+Pb collision there is a sizable geometric component (other than event-by-event fluctuations) that drives $v_{2}$.

In figure 3, results for the flow coefficients $v_{n}$, n=2,3,4 minus $v_{5}$ are shown for p+Al , p+Au , d+Au and 3He+Au collisions at $\sqrt{s}=200$ GeV per nucleon pair (highest RHIC energies). As before, results with and without pre-equilibrium flow are shown and boxes indicate theoretical uncertainty. Again, one finds that simulated flow coefficients are in overall agreement with experimental data, where available. Interestingly, the simulated results, both with and without pre-equilibrium flow predict a $v_{3}$ component in d+Au collisions that is almost as large as in 3He+Au . While $v_{3}$ in central 3He+Au has been measured by the PHENIX experiment (Huang (December 2014)), a non-vanishing $v_{3}$ in d+Au central collisions has not been seen in any experiment yet. However, given the unambiguous presence of a sizable $v_{3}$ component in the superSONIC simulation package, as well as in other theory (hydrodynamic and non-hydrodynamic) simulations Schenke and Venugopalan (2014b); Koop et al. (2015), measuring or putting an upper bound on $v_{3}$ in d+Au collisions at $\sqrt{s}=200$ GeV could serve as a very useful experimental verification of theory ingredients to light-heavy ion collisions.

Similar to the finding for $\sqrt{s}=5.02$ TeV, for results at $\sqrt{s}=200$ shown in Fig. 3 one finds noticeably smaller $v_{2}$ component in p+Al and p+Au collisions compared to d+Au and 3He+Au . The comparison between different systems for $v_{3}$ is highly dependent on the presence or absence of pre-equilibrium flow. For instance, with pre-equilibrium flow, one finds that $v_{3}$ in p+Au collisions at $\sqrt{s}=200$ GeV is almost as large as for 3He+Au , while without the presence of pre-equilibrium flow $v_{3}$ in p+Au is about half as large as in 3He+Au . This suggests that $v_{3}$ in light-heavy ion collisions at $\sqrt{s}=200$ GeV could be a good experimental handle on pre-equilibrium QCD dynamics. Unlike the situation found for $\sqrt{s}=5.02$ TeV, simulation results predict $v_{4}$ to be consistent with zero except for p+Au ,d+Au and 3He+Au where some small, non-vanishing $v_{4}$ is found. Since $v_{4}$ is found to be so small at $\sqrt{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 $\sqrt{s}=200$ GeV. The system created in these collisions consists of an average multiplicity of only $\frac{dN_{\rm ch}}{d\eta}\sim 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 $\sqrt{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 $\sqrt{s}=200$ GeV.

In Fig. 4 results for the flow coefficients $v_{n}$, n=2,3,4 minus $v_{5}$ are shown for p+Au , d+Au and 3He+Au collisions at $\sqrt{s}=62.4$ GeV per nucleon pair. Comparing results in Fig. 4 to results at $\sqrt{s}=200$ GeV in Fig. 3, it is hard to find clear differences in any of the flow observables shown. All $v_{n}$ at $\sqrt{s}=62.4$ GeV in all systems are broadly consistent with results at $\sqrt{s}=200$ GeV, only slightly lower. Nevertheless, one should point out that this in particular implies a sizable $v_{2}$ component in p+Au at $\sqrt{s}=62.4$, and hence hydrodynamic behavior at these collision energies.

In order to answer the question of a breakdown of hydrodynamic applicability, results for flow coefficients in proton-nucleus collisions (p+Au and p+Pb ) at energies ranging from $\sqrt{s}=7.7$ GeV to $\sqrt{s}=5.02$ TeV are compared in figure 5. Most remarkable, the simulations predict a sizable $v_{2}(2.5\ {\rm p_{T}})\simeq 5-8$ percent even for p+Au collisions at $\sqrt{s}=7.7$ GeV, with a predicted multiplicity of only $\frac{dN_{\rm ch}}{d\eta}\simeq 4.5$. Again, according to the criterion of applicability of hydrodynamics from above, hydrodynamics still applies for systems created in p+Au collisions at $\sqrt{s}=7.7$ GeV. Failure of finding a break-down of hydrodynamics in p+Au collisions at $\sqrt{s}=7.7$ GeV begs the question of which value of $\sqrt{s}$, if any, one would have to study in order to truly see hydrodynamic no longer apply. Unfortunately, at collision energies below $\sqrt{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 $\sqrt{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 $\sqrt{s}=5.02$ TeV and lowering the collision energy down to $\sqrt{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 $\sqrt{s}$, even though at $\sqrt{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 $\sqrt{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 pre-equilibrium 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 $\sqrt{s}\leq 200$ GeV. Specifically, while $v_{3}$ in 3He+Au collisions has been measured at $\sqrt{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 $\sqrt{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 $\sqrt{s}=5.04$ TeV, both because 3He+Pb collisions at the LHC are not currently planned, and because 3He+Pb to p+Pb $v_{3}$ ratios are only sensitive to pre-equilibrium flow at $p_{T}\lesssim 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 $\sqrt{s}=7.7$ GeV to $\sqrt{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 $\eta/s$ and $\zeta/s$) over a large temperature range $T\lesssim 0.53$ GeV ($<T>=0.4$ GeV, while by contrast p+Au collisions at $\sqrt{s}=7.7$ GeV only probe $T<0.33$ GeV ($<T>=0.25$ GeV). Hence p+Au collisions at $\sqrt{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. $\eta/s$.

Finally, results on HBT Radii for selected collision systems are reported in Tab. 1 in Figures 7,8 and 9.

IV Summary and conclusions

In this work, central light on heavy ion collisions from $\sqrt{s}=7.7$ GeV to $\sqrt{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 $\sqrt{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 $\sqrt{s}$ is lowered from $5.02$ TeV to $\sqrt{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 $\sqrt{s}=7.7$ in order to potentially observe also $v_{2}$ collapse to zero in p+Au collisions.

For d+Au collisions at $\sqrt{s}=200$ GeV superSONIC simulations predicted $v_{3}$ results that were of the same order of magnitude as those measured in 3He+Au . Since similar results are observed in other theoretical models for d+Au collisions at $\sqrt{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 3He+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.

Acknowledgements.

I would like to thank S. Huang and J.L. Nagle for many helpful discussions on this topic. This work was supported by the Department of Energy, DOE award No. DE-SC0008027. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794) and the University of Colorado Boulder. The Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver and the National Center for Atmospheric Research. Janus is operated by the University of Colorado Boulder.

Appendix A Plot Collection

For better readability of the text, many of the plots for this article are collected in this appendix.