Collective flow and the fluid behavior in p/d/$^3$He+Au collisions at $\sqrt{s_{NN}} = 200$ GeV

By varying the intrinsic initial geometry, the p/d/$^3$He+Au collisions at the Relativistic Heavy Ion Collider (RHIC) provide a unique opportunity to understand the collective behavior and probe the possible sub-nucleon fluctuations in small systems. In this paper, we employ the hybrid model iEBE-VISHNU with Trento initial conditions to study the collective flow and the fluid behavior in p/d/$^3$He+Au collisions. With fine-tuned parameters, iEBE-VISHNU can describe the $v_2(p_T)$ and $v_3(p_T)$ data from the PHENIX and STAR collaborations. However, for some certain parameter sets with initial sub-nucleon fluctuation, the hydrodynamic simulations have already beyond their limits with the average Knudsen number $\langle K_n \rangle$ obviously larger than one. Our calculations demonstrate that, for a meaningful evaluation of the fluid behavior in the small systems, model simulations should also pay attention to the validity range of hydrodynamics.

To understand these flow-like observables, a key question is whether the QGP droplets were formed in the small collision systems.The theoretical study can be classified into two scenarios, the final state effect associated with the QGP fluid expansion and the initial state effect without the QGP formation.In the final state scenario, the collective flow is related to the initial state geometry through the non-linear evolution, where hydrodynamic or kinetic calculations can roughly reproduce the data with tuned parameters .In contrast, the initial state effect describes the observed anisotropies by momentum correlation of the initially produced particles in the color field domain.One typical model is the Color Glass Condensate(CGC), which can qualitatively describe the experimental measurements such as two-particle and multi-particle correlations, long-range rapidity correlations and mass ordering [60][61][62][63][64][65][66][67][68][69][70].Compared to the final state scenario, a major difference is that the correlations are expected to be weaker for larger collision systems with more uncorrelated domains involved.
It was believed that comparative runs of p/d/ 3 He + Au collisions with variation of the initial state geometry could provide useful information to identify the above two scenarios for flow-like signals in the small systems [54].More specifically, with the initial state geometry dominated by the nucleon position fluctuations, models such as MC-Glauber show an eccentricity ordering as: ε p+Au [69,72].
With this motivation, the PHENIX Collaboration has measured the flow anisotropies of v 2 (p T ) and v 3 (p T ) in p/d/ 3 He+Au collisions using the event-plane method [54] and the "3 × 2PC" method [73], respectively.Both methods observed similar flow ordering: 3 He+Au 3 [54], which consists with the hydrodynamic predictions with MC-Glauber initial conditions [37,41,71,74].The PHENIX measurements also largely ruled out the CGC calculations with only initial state effects, which show different v 2 and v 3 orderings of smaller magnitude [72,73,75].Recently, the STAR collaboration has also measured the flow coefficients of these three small systems using nonflow subtraction methods based on the template fit and the Fourier expansion fit, respectively [76,77].Compared On the theoretical side, the flow signals of these small collision systems have been studied by various models with different initial conditions [37,46,56,[80][81][82].The hydrodynamic calculations with MC-Glauber initial conditions can well describe the v 2 (p T ) and v 3 (p T ) data from PHENIX, but cannot explain the flow discrepancy between STAR and PHENIX, even with the longitudinal decorrelations in the 3+1-d simulations [80].Using the initial conditions with sub-nucleon fluctuations, hydrodynamic simulations produce similar v 3 values in p/d/ 3 He+Au collisions, which is obviously different from the PHENIX results and the calculations with nucleon fluctuations.However, these hydrodynamic calculations cannot simultaneously describe the STAR v 2 and v 3 data and the initial sub-nucleon structure is still not well constrained.
In this paper, we study the flow observables in p/d/ 3 He + Au collisions at √ s N N = 200 GeV, using 2+1-d iEBE-VISHNU model with TRENTo initial conditions including both nucleon and sub-nucleon flucuations.We tune the model parameters to fit the elliptic and triangular flow data from PHENIX and STAR, respectively, and calculate the 4-particle cumulants c 2 {4} in p+Au and d+Au systems.We also evaluate the validity of the hydrodynamics by Knudsn numbers and draw the hydrodynamic predicted v 3 /v 2 (p T ) band in p/d/ 3 He + Au collisions with the relatively reliable region of 2+1-d hydrodynamics.
We implement TRENTo, a parameterized initial condition model, to generate the initial entropy density for the hydrodynamic simulations starting at τ 0 [95][96][97].In the case without sub-nucleon structure, the fluctuations come from the distribution of nucleon center position.For each nucleon, its density distribution is parameterized as a Gaussian function with nucleon Gaussian width ω: On the other hand, when considering the sub-nucleon fluctuation, the nucleon is assumed to be composed of independent constituents, and the nucleon density is written as: where n c is the constituent number, x i is the position of the i-th constituent, density ρ constit defined as ).The constituent Gaussian width v relates to the nucleon width ω with a standard deviation r: ω = √ r 2 + v 2 and in this case ω is defined as the root mean square radius of a nucleon.
After obtaining the nucleon density distribution, the  I.The STAR and PHENIX data are taken from [76] and [54],respectively.fluctuated thickness of the colliding nucleons is then written as: Here, besides the nucleon/constituent position, the initial fluctuation is controlled mainly by the gamma random variable γ i , which is parameterized by the shape factor k. The resulting standard devaition of the initial fluctuation is denoted as With the fluctuating thickness TA,B , the initial entropy density at mid-rapidity can be calculated by the generalized means with a dimensionless parameter p: and v 3 at the same time requires a large shear viscosity in Para-I, which lies outside the usual parameter range of the hydrodynamic approaches.This will be further discussed in section IV.Fig. 1 plots the averaged eccentricities ε 2,3 , calculated from the TRENTo initial conditions with the parameter sets in Tab.I. Compared to the early results of the MC-Glauber model [41,54], Fig. 1 shows a weaker ordering of ε n for p+Au, d+Au and 3 He+Au collisions.Such weaker ordering in Para-I and Para-III can be explained by sub-nucleonic fluctuations.Due to the imprinted multiplicity fluctuations, Para-II without subnucleonic fluctuations also shows a weaker ordering of ε n for the TRENTo initial conditions [71].

III. RESULTS AND DISCUSSIONS
In this section, we show the flow harmonic results calculated by iEBE-VISHNU with the parameters listed in Table I.All parameter sets in Table I are well tuned to reproduce the multiplicity distribution in d+Au collisions.We also tune Para-I to fit v 2 (p T ) and v 3 (p T ) from STAR, and Para-II and Para-III to fit v 2 (p T ) and v 3 (p T ) from PHENIX.
Figure 2  TRENTo initial condition with three parameter sets.The PHENIX data are taken from [98].
and v 3 (p T ) of all charged hadrons in 0-5% p+Au, d+Au and 3 He+Au collisions [99].For the triangular flow v 3 (p T ), the STAR and PHENIX measurements are not consistent, particularly for p+Au and d+Au collisions, where the STAR data are larger than the PHENIX data by a factor of 3.Such apparent discrepancies may be due to the different rapidity region and non-flow substraction methods used by these two collaborations [73,76].Therefore, we fit the STAR and PHENIX data, respectively.Following [48,100], we use the two-subevent cumulant method to calculate the two particle correlations with a kinematic cut 0.2 < p T < 2.0 GeV/c and |η| < 1.0 with a gap |∆η| > 1.0.
With the Para-I, iEBE-VISHNU nicely fits the v 2 (p T ) and v 3 (p T ) data measured by STAR.We find that subnucleon fluctuations are essential to produce larger v 3 , which are insensitive to the collision systems variation.Meanwhile a large shear viscosity is also required to simultaneously fit the v 2 and v 3 data from STAR.The validity of hydrodynamic simulations with such large shear viscosities will be discussed in the next section.Note that the effects of sub-nucleon fluctuations on the flow of small systems also have been studied and discussed in the earlier paper [56].For the PHENIX measurements, iEBE-VISHNU simulations with nucleon fluctuations in the initial state (Para-II) were able to roughly reproduce the v 2 and v 3 data within the statistical error bars.Meanwhile, one can achieve similar results with sub-nucleonic initial conditions with a free streaming (Para-III).Here, Para-III is the one obtained from the Bayesian analysis for Pb+Pb and p+Pb collisions at √ s N N = 5.02 TeV [97] except for the fluctuation parameter k and the constituent width ν, which are tuned to fit the multiplicity fluctuation in top RHIC energy.We conclude that, for the PHENIX measurements, iEBE-VISHNU simulations with both nucleonic and sub-nucleonic initial state fluctuations can fit the v 2 and v 3 hierarchies in p/d/ 3 He+Au collisions.

IV. APPLICABILITY OF HYDRODYNAMIC SIMULATION
We have noticed that in the above calculation, the specific shear viscosity in some parameter sets tuned to fit the v 2 (p T ) and v 3 (p T ) data become quite large.In order to evaluate the validity of the hydrodynamic simulations in small systems, we calculate the Knudsen number K n defined as [101]: where τ π is the relaxation time associated with the microscopic time scale and θ = ∂ µ u µ is the expansion rate associated with the macroscopic hydrodynamic time scale.K n → 0 is the perfect fluid limit where the local equilibrium is maintained during the hydrodynamic evolution.
K n → ∞ is the other limit, which corresponds to the case that the fluid system breaks up into free-streaming particles.It is generally suggested that the hydrodynamics is relatively reliable with Fig. 4 shows the time evolution of the averaged Knudsen number ⟨K n ⟩ in the event-by-event hydrodynamic simulations for p/d/ 3 He+Au collisions at 0-5% centrality.The average is taken within the freeze-out hypersurface with the local energy density as the weight for each time step.For the Para-I with sub-nucleon fluctuations we observed that the averaged Knudsen number ⟨K n ⟩ is always larger than 1 throughout the whole evolution for different collision systems.Obviously, such a large Knudsen number indicates that the hydrodynamic simulations are beyond their applicable limit, which is mainly due to the large specific shear viscosity η/s ∼ 0.28 and the large initial gradients introduced by fluctuations to fit the v 3 data.In contrast, the average Knudsen number for Para-II is about or less than 1 with a smaller specific shear viscosity η/s ∼ 0.09.For Para-III, the Knudsen number lies between those of Para-I and Para-II, which is large in the early time due to the free streaming evolution before thermalization, but drops below 1 after certain time of hydrodynamic evolution.In short, Fig. 4 suggests that hydrodynamic simulations with Para-I that tuned to fit the STAR data are beyond the limit due to the large Knudsen number.To further investigate whether iEBE-VISHNU could fit all the experimental data within its hydrodynamic limit, we explore the model parameter space as far as possible but with the constraint of the Knudsen number ⟨K n ⟩ < 1 at the end of the evolution.Our test parameter sets correspond such initial conditions with/without nucleon substructure and with/without the free-streaming effect.The range of the free parameters is listed in Tab.II.With n c = 1, the initial conditions include only nucleon fluctuations and with n c = 2 − 9, the initial conditions include sub-nucleon fluctuations.In our investigation, the effective shear viscosity η/s and the shape parameter k are fixed to reproduce the v 2 (p T ) data in 0-5% 3 He+Au collisions and the multiplicity distribution of d+Au collisions with neglecting the bulk viscosity.
Fig. 5 shows the p T dependent v 3 (p T )/v 2 (p T ) ratio in 0-5% p/d/ 3 He+Au collisions, where the theoretical band is calculated by iEBE-VISHNU with the parameter range listed in Tab.II, together with the constraint ⟨K n ⟩ < 1.The experimental data are taken from STAR with the statistical uncertainty of v 2 (p T ) and v 3 (p T ) using the error propagation formula.As shown in panels (b) and (c), the flow harmonic ratio v 3 (p T )/v 2 (p T ) in d/ 3 He+Au collisions can be reproduced by iEBE-VISHNU within the allowed parameter range ⟨K n ⟩ < 1.While panel (a) shows that the upper limit of the v 3 (p T )/v 2 (p T ) ratio in p+Au collisions, calculated from iEBE-VISHNU simulations, is clearly below the experimental data.These results indicate that the current hybrid model calculations are not able to simultaneously describe the STAR flow data in the three small collision systems with reasonable parameter range within the hydrodynamic limit.

V. SUMMARY
In this paper, we implemented iEBE-VISHNU with TRENTo initial condition to study the collective flow in p/d/ The ratio of v3(pT ) and v2(pT ) for all charged hadrons in 0-5% p/d/ 3 He+Au collisions, calculated by iEBE-VISHNU using the the parameter range listed in Tab.II, together with the constraint ⟨Kn⟩ < 1.The STAR data are taken from [76].
PHENIX measurements, v 2 (p T ) and v 3 (p T ) data show obvious hierarchies for different collision systems, which can be reproduced by our hybrid model simulations with nucleon/sub-nucleon fluctuating initial conditions.The related simulations also reproduce a negative 4-particle cumulant c 2 {4} for d+Au collisions, but give an almost zero c 2 {4} for p+Au collisions, which can not describe the positive c 2 {4} measured by PHENIX.For the STAR measurements, the magnitude of v 3 are insensitive to the collision systems.iEBE-VISHNU simulations with sub-nucleon fluctuating initial conditions can fit these v 2 and v 3 data, which can also roughly reproduce the positive and negative c 2 {4} measured in the high multiplicity p+Au and d+Au collisions, respectively.However, due to the large shear viscosity applied to fit experimental data, the hydrodynamic simulations has already beyond its limits with the average Knudsen number ⟨K n ⟩ obviously larger than one for these three collision systems.We also explore the model parameter space as far as possible with the constraint of the Knudsen number ⟨K n ⟩ < 1, and found that iEBE-VISHNU with the TRENTo initial condition with/without sub-nucleon fluctuations cannot describe the experimental measured v 3 /v 2 ratio for p+Au collisions.
Our calculations demonstrate that for a meaningful evaluation of the collective flow in the small systems, one should also evaluate the validity of hydrodynamics.As the collision system become smaller, the isotropization and thermalization conditions are harder and harder to reach.Besides applying a full 3+1-d hydrodynamic simulations, other improved hydrodynamic theories like anisotropic hydro [104][105][106][107] should be implemented to the small systems that may not reach equilibrium in the early stage.It was also found that the fragmentation and mini-jets effects become more important in small collision systems [53,108], a comprehensive model includes the core-corona effects [109,110] is also required to further evaluate the flow signals in p/d/ 3 He+Au collisions at √ s N N = 200 GeV.
list the model parameters used in the calculations for p+Au, d+Au and 3 He+Au collisions at √ s N N = 200 GeV.The Para-I is tuned to fit the v 2 (p T ) and v 3 (p T ) data from STAR, which includes the initial nucleon sub-structure and a small constituent width to enlarge the fluctuation.The Para-II and Para-III fit the published PHENIX data, which are tuned with/without sub-nucleon fluctuation, respectively.Para-III with subnucleon structure is similar to the Bayesian analyses in the p+Pb and Pb+Pb collisions [97], except that TRENTo paratemr k and ν are tuned to fit the charged particle multiplicity distribution in d+Au collisions at √ s N N = 200 GeV.Note that, reproducing the STAR v 2

Fig. 3
Fig. 3 plots the 4-particle cumulant c 2 {4} as a function of dN ch /dη for p+Au and d+Au collisions at √ s N N =

TABLE I .
Parameters setups

TABLE II .
Free parameter range for the hybrid model.

3
He+Au collisions at √ s N N = 200 GeV.For the