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

Over the last two decades, the properties of extremely hot and dense QCD matter, i.e., quark-gluon plasma (QGP), have been intensively studied by the relativistic heavy-ion programs at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC). It has been found that the QGP created in large collision systems, such as Au+Au and Pb+Pb collisions, behaves as an almost ''perfect liquid'' [1−15]. Its strong collective expansion and associated flow observables have been successfully described by hydrodynamic calculations with small specific shear viscosity, which is close to the lowest KSS bound [16−18]. Small collision systems, such as p+Pb and p+p collisions at the LHC, were originally intended to provide reference data for large collision systems. However, similar collective behaviour has been observed in high multiplicity events, including the ''double ridge'' structure in two particle correlations [19−29], multi-particle cumulants [22, 24, 29−32], and mass ordering of the anisotropic flow of identified particles [28, 31, 33]; see [14, 34−36] for recent reviews.

To understand these flow-like observables, a key question is whether QGP droplets are formed in small collision systems. The theoretical study of this question can be classified into two scenarios: the final state effect associated with the QGP fluid expansion and the initial state effect without QGP formation. In the final state scenario, the collective flow is related to the initial state geometry through non-linear evolution, where hydrodynamic or kinetic calculations can approximately reproduce the data with tuned parameters [37−59]. 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 experimental measurements such as two-particle and multi-particle correlations, long-range rapidity correlations, and mass ordering [60−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/³He + Au collisions with variation of the initial state geometry could provide useful information to identify the aforementioned two scenarios for flow-like signals in small systems [54]. More specifically, with the initial state geometry dominated by nucleon position fluctuations, models such as MC-Glauber produce the following eccentricity ordering: $ \varepsilon_2^{{p}+\text{Au}} < \varepsilon_2^{ {d}+\text{Au}} \simeq \varepsilon_2^{^3 \text{He}+\text{Au}} $, $\varepsilon_3^{ {p}+\text{Au}} \simeq \varepsilon_3^{{d}+\text{Au}} < \varepsilon_3^{^3 \text{He}+\text{Au}}$. Hydrodynamic evolution responses directly to such initial state eccentricities and generates the associated flow ordering as follows: $ v_2^{ {p}+\text{Au}} < v_2^{ {d}+\text{Au}} \approx v_2^{^3 \text{He}+\text{Au}} $, $ v_3^{{p}+\text{Au}} \approx v_3^{ {d}+\text{Au}} < v_3^{^3 \text{He}+\text{Au}} $ [37, 41, 71]. In contrast, CGC model calculations are not sensitive to such initial state geometry, predicting the flow anisotropies of much smaller magnitude for the three collision systems as follows: $ v_2^{{p}+\text{Au}} \gtrsim v_2^{{d}+\text{Au}} \gtrsim v_2^{^3 \text{He}+\text{Au}} $, $ v_3^{ {p}+\text{Au}} \gtrsim v_3^{ {d}+\text{Au}} \gtrsim v_3^{^3 \text{He}+\text{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/³He+Au collisions using the event-plane [54] and “3 $ \times $ 2PC” [73] methods, respectively. Both methods produced similar flow ordering, $ v_2^{p+\text{Au}} < v_2^{d+\text{Au}} \approx v_2^{^3\text{He}+\text{Au}} $ and $ v_3^{p+\text{Au}} \approx v_3^{d+\text{Au}} < v_3^{^3\text{He}+\text{Au}} $ [54], which agrees with the hydrodynamic predictions under MC-Glauber initial conditions [37, 41, 71, 74]. The PHENIX measurements also largely ruled out the CGC calculations with only initial state effects, producing 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 non-flow subtraction methods based on the template and Fourier expansion fits, respectively [76, 77]. Compared to the PHENIX results, the STAR measurements resulted in much larger values of $ v_3(p_T) $ in p+Au and d+Au systems with the TPC detector at mid-rapidity. Such system independent $ v_3 $ data might indicate the existence of sub-nucleon fluctuations, while the discrepancy between measurements may also originate from different non-flow subtractions or different detector pseudorapidity acceptances [73, 78, 79].

The flow signals of these small collision systems have been theoretically studied using various models and different initial conditions [37, 46, 56, 74, 80, 81]. The hydrodynamic calculations under MC-Glauber initial conditions can well describe the $ v_2(p_T) $ and $ v_3(p_T) $ data from PHENIX. However, they cannot explain the flow discrepancy between STAR and PHENIX, even with longitudinal decorrelations in 3+1-d simulations [74]. Using the initial conditions with sub-nucleon fluctuations, hydrodynamic simulations produce similar $ v_3 $ values in p/d/³He+Au collisions, which clearly differ 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 not well constrained yet.

In this study, we analyzed the flow observables in p/d/³He + Au collisions at $ \sqrt{s_{NN}} = 200 $ GeV using the 2+1-d $ {\mathrm{iEBE-VISHNU}}$ model under TRENTo initial conditions including both nucleon and sub-nucleon fluctuations. We tuned the model parameters to fit the elliptic and triangular flow data from PHENIX and STAR, respectively, and calculated the 4-particle cumulants $ c_2\{4\} $ in p+Au and d+Au systems. We also evaluated the validity of the hydrodynamics in terms of Knudsen numbers and depicted the predicted $ v_3/v_2(p_T) $ hydrodynamic band in p/d/³He + Au collisions with the relatively reliable region of 2+1-d hydrodynamics.

In this study, we implemented iEBE-VISHNU to study the flow observables in p/d/³He+Au collisions at $ \sqrt{s_{NN}} = 200 $ GeV. iEBE-VISHNU [82] is an event-by-event hybrid model that combines 2+1-d viscous hydrodynamics VISH2+1 [83−85] for QGP evolution, a particle sampler iSS [16, 86] for the particlization at a switching temperature, and the hadron cascade model UrQMD [87, 88] for subsequent hadronic evolution. According to [51, 52, 89−91], we used the HotQCD+HRG equation of state (EoS) [92−94] as input, and set temperature-dependent specific shear viscosity $ \eta/s $ and bulk viscosity $ \zeta/s $ [17].

We implemented TRENTo, a parameterized initial condition model, to generate the initial entropy density for hydrodynamic simulations starting at $ \tau_0 $ [94−96]. 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 ω:

By contrast, when considering sub-nucleon fluctuations, the nucleon is assumed to be composed of independent constituents, and the nucleon density is expressed as

where $ n_c $ is the constituent number, $ {\bf{x}}_i $ is the position of the i-th constituent, and the density $ \rho_{\text{constit}} $ is defined as $\rho_{\text{constit}} ({\bf{x}}) = {(2\pi v^2)^{-3/2}} \text{exp} \left(-{{\bf{x}}^2}/({2v^2})\right)$. The constituent Gaussian width v relates to the nucleon width ω with a standard deviation r as $ \omega = \sqrt{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 fluctuated thickness of the colliding nucleons is expressed as

Here, besides the nucleon/constituent position, the initial fluctuation is mainly controlled by the gamma random variable $ \gamma_i $, which is parameterized by the shape factor k. The resulting standard deviation of the initial fluctuation is denoted as

With the fluctuating thickness $ \tilde{T}_{A, B} $, the initial entropy density at mid-rapidity can be calculated by generalized means with a dimensionless parameter p:

Table 1 lists the model parameters used in the calculations for p+Au, d+Au and ³He+Au collisions at $ \sqrt{s_{NN}} $ = 200 GeV. Para-I is tuned to fit the $ v_2(p_T) $ and $ v_3(p_T) $ data from STAR, which include the initial nucleon sub-structure and a small constituent width to enlarge the fluctuation. Para-II and Para-III fit the published PHENIX data, tuned with/without sub-nucleon fluctuations, respectively. Para-III with sub-nucleon structure is similar to Bayesian analyses in p+Pb and Pb+Pb collisions [96], except that TRENTo paratemers k and ν were tuned to fit the charged particle multiplicity distribution in d+Au collisions at $ \sqrt{s_{NN}} = 200 $ GeV. Note that simultaneously reproducing the STAR $ v_2 $ and $ v_3 $ requires a large shear viscosity in Para-I, which lies outside the usual parameter range of hydrodynamic approaches. This will be further discussed in Sec. IV. Figure 1 represents the averaged eccentricities $ \varepsilon_{2, 3} $, calculated from the TRENTo initial conditions with the parameter sets listed in Table 1. Compared to the early results of the MC-Glauber model [41, 54], Fig. 1 shows a weaker ordering of $ \varepsilon_n $ for p+Au, d+Au, and ³He+Au collisions. Such weaker ordering in Para-I and Para-III can be explained by sub-nucleonic fluctuations. Owing to the imprinted multiplicity fluctuations, Para-II without subnucleonic fluctuations also shows a weaker ordering of $ \varepsilon_n $ for the TRENTo initial conditions [71].

	Parameter Set	p	k	$n_c$	ω	ν	$\tau_0$	$(\eta/s)_{\text{min}}$	$(\eta/s)_{\text{slope}}$	$(\zeta/s)_{\text{max}}$	${\rm{T}}_{\text{switch}}$
I	Sub-Nucl. Fluc.	0	0.6	5	0.5	0.2	0.6	0.28	1.6	0.022	154
II	Nucl. Fluc.	1	1.6	1	0.4	−	0.6	0.09	1.0	0.0	151
III	Sub-Nucl. Fluc.	0	0.28	6	0.92	0.55	0.37a	0.11	1.6	0.032	151
a with free-streaming

Table 1.  Parameter setups.

Figure 1.  (color online) Initial eccentricities $\varepsilon _2$ and $\varepsilon_3$calculated from TRENTo with impact parameter b < 2 fm.

In this section, we show the flow harmonic results calculated by iEBE-VISHNU with the parameters listed in Table 1. All parameter sets in Table 1 were well tuned to reproduce the multiplicity distribution in d+Au collisions. We also tuned 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 shows the differential flow harmonics $ v_2(p_T) $ and $ v_3(p_T) $ of all charged hadrons in 0−5% p+Au, d+Au, and ³He+Au collisions ¹. 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. According to [48, 97], we used the two-subevent cumulant method to calculate the two particle correlations with a kinematic cut $ 0.2<p_T<2.0 $ GeV/c and $ |\eta|<1.0 $ with a gap $ |\Delta \eta|>1.0 $.

Figure 2.  (color online) Flow harmonics $v_2(p_T)$ and $v_3(p_T)$ for all charged hadrons in 0−5% $p+$Au (a), $d+$Au (b), and ³He+Au (c) collisions at $\sqrt{s_{NN}} = 200$GeV. The theoretical curves were calculated by iEBE-VISHNU under TRENTo initial conditions with the three parameter sets listed in Table 1. The STAR and PHENIX data were taken from [76] and [54], respectively.

Concerning Para-I, iEBE-VISHNU properly fit the $ v_2(p_T) $ and $ v_3(p_T) $ data measured by STAR. We found that sub-nucleon fluctuations are essential to produce larger values of $ v_3 $, which are insensitive to the variation of the collision system. 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 have also been studied and discussed previously [56]. Regarding the PHENIX measurements, iEBE-VISHNU simulations with nucleon fluctuations in the initial state (Para-II) were able to approximately reproduce the $ v_2 $ and $ v_3 $ data within statistical error bars. Meanwhile, one can achieve similar results under sub-nucleonic initial conditions with a free streaming (Para-III). Here, Para-III was obtained from Bayesian analysis of Pb+Pb and p+Pb collisions at $ \sqrt{s_{NN}} = 5.02 $ TeV [96], except for the fluctuation parameter k and constituent width ν, which were tuned to fit the multiplicity fluctuation in top RHIC energy. We conclude that, for 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/³He+Au collisions.

Figure 3 depicts the 4-particle cumulant $ c_2\{4\} $ as a function of $ {\rm d}N_{\rm ch}/{\rm d}\eta $ for p+Au and d+Au collisions at $ \sqrt{s_{NN}} = 200 $ GeV. Panel (a) shows the iEBE-VISHNU predictions for Para-I, which was tuned to fit the STAR $ v_2(p_T) $ and $ v_3(p_T) $ data and generated positive $ c_2\{4\} $ for p+Au collisions and negative $ c_2\{4\} $ for d+Au collisions in high multiplicity events. This is qualitatively consistent with experimental measurements. In fact, large event-by-event fluctuations in Para-I lead to positive $ c_2\{4\} $ in p+Au collisions. Concerning d+Au collisions, the intrinsic geometry of the deutron produces a dominant contribution to the initial eccentricities, leading to a negative $ c_2\{4\} $ in high multiplicity events. Panels (b) and (c) show the iEBE-VISHNU results, calculated for Para-II and Para-III fitting the PHENIX $ v_2(p_T) $ and $ v_3(p_T) $ data. For p+Au collisions, $ c_2\{4\} $ is always close to zero for both parameter sets over the whole range of multiplicities owing to small flow fluctuations. Regarding d+Au collisions, $ c_2\{4\} $ values are always negative owing to the intrinsic geometry of the deutron.

Figure 3.  (color online) $c_2\{4\}$ as a function of ${\rm d}N_{\text{ch}}/{\rm d}\eta$ in p+Au and d+Au collisions at $\sqrt{s_{NN}} = 200$ GeV, calculated by iEBE-VISHNU under TRENTo initial conditions with three parameter sets. The PHENIX data were taken from [98].

Note that in the above calculations, the specific shear viscosity in some parameter sets tuned to fit the $ v_2(p_T) $ and $ v_3(p_T) $ data becomes notably large. To evaluate the validity of hydrodynamic simulations in small systems, we calculated the Knudsen number $ K_{n} $ defined as follows [99]:

where $ \tau_{\pi} $ is the relaxation time associated with the microscopic time scale and $ \theta = \partial_{\mu} u^\mu $ is the expansion rate associated with the macroscopic hydrodynamic time scale. $ K_n \rightarrow 0 $ is the perfect fluid limit at which the local equilibrium is maintained during hydrodynamic evolution. $ K_n \rightarrow \infty $ is the other limit, which corresponds to the case in which the fluid system breaks up into free-streaming particles. It is generally suggested that the hydrodynamics is relatively reliable for $ K_n < 1 $ [99, 100] ².

Figure 4 shows the time evolution of the averaged Knudsen number $ \langle K_n \rangle $ in the event-by-event hydrodynamic simulations for p/d/³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. Regarding Para-I with sub-nucleon fluctuations, we observed that the averaged Knudsen number $ \langle K_n \rangle $ is always larger than unity throughout the whole evolution for different collision systems. Such a large Knudsen number clearly indicates that the hydrodynamic simulations are beyond their applicable limit owing to the large specific shear viscosity $ \eta/s \sim 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 approximately equal or less than unity with a smaller specific shear viscosity $ \eta/s \sim 0.09 $. Regarding Para-III, the Knudsen number lies between those of Para-I and Para-II, which is large in early time instants owing to the free streaming evolution before thermalization. However, it drops below 1 after a certain interval of hydrodynamic evolution. In short, Fig. 4 suggests that hydrodynamic simulations with Para-I tuned to fit the STAR data are beyond the limit owing to a large Knudsen number.

Figure 4.  (color online) Evolution of the average Knudsen number $\langle K_n \rangle$ within the freeze-out hypersurface for 0–5% p/d/³He+Au collisions.

To further investigate whether iEBE-VISHNU could fit all the experimental data within its hydrodynamic limit, we explored the model parameter space as much as possible under the following constraint on the Knudsen number at the end of the evolution: $ \langle K_n \rangle<1 $. Our test parameter sets correspond to such initial conditions with/without the nucleon substructure and with/without the free-streaming effect. The range of free parameters is listed in Table 2. For $ n_c = 1 $, the initial conditions include only nucleon fluctuations, and for $ n_c = 2-9 $, the initial conditions include sub-nucleon fluctuations. In this study, the effective shear viscosity $ \eta/s $ and shape parameter k were fixed to reproduce the $ v_2(p_T) $ data in 0−5% ³He+Au collisions and the multiplicity distribution of d+Au collisions while neglecting the bulk viscosity.

Parameter	Description	Prior range
ω	Nucleon width	[0.4, 1.0] fm
ν	Constituent width ($< \omega$)	[0.3, 1.0] fm
$n_c$	Number of nucleon constituents	[1, 9]
$\tau_{\text{fs}}$	Free-streaming time	[0.2, 0.8] fm/c
$T_{\text{switching}}$	Switching temperature	[135, 165] MeV
Free-stream		on/off
Sub. fluct.		on/off

Table 2.  Free parameter range for the hybrid model.

Figure 5 represents $ p_T $ as a function of the $ v_3(p_T)/v_2(p_T) $ ratio in 0−5% p/d/³He+Au collisions; the theoretical band was calculated by iEBE-VISHNU according to the parameter range listed in Table 2 under the constraint $ \langle K_n \rangle<1 $. The experimental data were taken from STAR with the statistical uncertainty of $ v_2(p_T) $ and $ v_3(p_T) $ obtained from 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/³He+Au collisions can be reproduced by iEBE-VISHNU within the allowed parameter range $ \langle K_n \rangle<1 $. 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 three small collision systems with a reasonable parameter range within the hydrodynamic limit.

Figure 5.  (color online) Ratio of $v_3(p_T)$ and $v_2(p_T)$ for all charged hadrons in 0−5% p/d/³He+Au collisions calculated by iEBE-VISHNU using the parameter range listed in Table 2 under the constraint $\langle K_n \rangle<1$. The STAR data were taken from [76].

In this study, we implemented iEBE-VISHNU under TRENTo initial conditions to study the collective flow in p/d/³He+Au collisions at $ \sqrt{s_{NN}} = 200 $ GeV. For the PHENIX measurements, $ v_2(p_T) $ and $ v_3(p_T) $ data show evident hierarchies for different collision systems. These hierarchies can be reproduced by our hybrid model simulations under nucleon/sub-nucleon fluctuating initial conditions. The related simulations also reproduced a negative 4-particle cumulant $ c_2\{4\} $ for d+Au collisions, but an almost zero $ c_2\{4\} $ for p+Au collisions that does not describe the positive $ c_2\{4\} $ measured by PHENIX. Regarding the STAR measurements, the magnitude of $ v_3 $ is insensitive to the collision system. iEBE-VISHNU simulations with sub-nucleon fluctuating initial conditions can fit these $ v_2 $ and $ v_3 $ data, which can also approximately reproduce the positive and negative $ c_2\{4\} $ values measured in high multiplicity p+Au and d+Au collisions, respectively. However, owing to the large shear viscosity applied to fit experimental data, the hydrodynamic simulations were already beyond its limits with an average Knudsen number $ \langle K_n \rangle $notably larger than unity for these three collision systems. We also explored the model parameter space as much as possible under the following constraint on the Knudsen number: $ \langle K_n \rangle<1 $. We found that $ {\mathrm{iEBE-VISHNU}}$under TRENTo initial conditions with/without sub-nucleon fluctuations cannot describe the experimentally measured $ v_3 / v_2 $ ratio for p+Au collisions.

Our calculations demonstrate that for a meaningful evaluation of the collective flow in small systems, one should also evaluate the validity of hydrodynamics. As the collision system becomes smaller, the isotropization and thermalization conditions become progressively more difficult to reach. Besides applying full 3+1-d hydrodynamic simulations, other improved hydrodynamic theories such as anisotropic hydro [101−104] can be applied to small systems that may not reach equilibrium in an early stage. It was also found that the fragmentation and mini-jets effects become more important in small collision systems [53]. A comprehensive model that includes the core-corona effects [105, 106] is also required to further evaluate the flow signals in p/d/³He+Au collisions at $ \sqrt{s_{NN}} = 200 $ GeV.

We thank Wenbin Zhao for helpful discussions. We also acknowledge the extensive computing resources provided by the Supercomputing Center of Chinese Academy of Science (SCCAS), Tianhe-1A from the National Supercomputing Center in Tianjin, China and the High-performance Computing Platform of Peking University