# The 3D structure of anisotropic flow in small collision systems
at the Relativistic Heavy Ion Collider

Wenbin Zhao  Department of Physics and Astronomy, Wayne State University,
Detroit, Michigan 48201, USA  Sangwook Ryu  Department of Physics and
Astronomy, Wayne State University, Detroit, Michigan 48201, USA  Chun Shen
Department of Physics and Astronomy, Wayne State University, Detroit, Michigan
48201, USA  RIKEN BNL Research Center, Brookhaven National Laboratory, Upton,
NY 11973, USA  Björn Schenke  Physics Department, Brookhaven National
Laboratory, Upton, NY 11973, USA

###### Abstract

We present (3+1)D dynamical simulations of asymmetric nuclear collisions at
the Relativistic Heavy Ion Collider (RHIC). Employing a dynamical initial
state model coupled to (3+1)D viscous relativistic hydrodynamics, we explore
the rapidity dependence of anisotropic flow in the RHIC small system scan at
200 GeV center of mass energy. We calibrate parameters to describe central
3He+Au collisions and make extrapolations to d+Au and p+Au collisions. Our
calculations demonstrate that approximately 50% of the $v_{3}(p_{T})$
difference between the measurements by the STAR and PHENIX Collaborations can
be explained by the use of reference flow vectors from different rapidity
regions. This emphasizes the importance of longitudinal flow decorrelation for
anisotropic flow measurements in asymmetric nuclear collisions, and the need
for (3+1)D simulations. We also present results for the beam energy dependence
of particle spectra and anisotropic flow in d+Au collisions.

##  I Introduction

Ultra-relativistic collisions of heavy ions are expected to create nucleus-
sized droplets of quark-gluon plasma (QGP), providing a unique opportunity to
study the properties of nuclear matter at extreme densities and temperatures.
Precise measurements of various flow observables performed at the Relativistic
Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) together with
the successful descriptions by hydrodynamic calculations have revealed that
the created QGP fireball behaves like a nearly perfect fluid with very small
specific shear viscosity Heinz and Snellings (2013); Gale  _et al._ (2013);
Shen and Heinz (2015); Shen and Yan (2020). Over the past ten years,
experiments and theorists have employed multiple techniques to assess whether
such QGP droplets are also formed in smaller collisions, such as p+A, p+p, and
even $\gamma^{*}$+A collisions Aad  _et al._ (2021); Zhao  _et al._ (2022);
Shen  _et al._ (2022) (see Dusling  _et al._ (2016); Loizides (2016);
Schlichting and Tribedy (2016); Nagle and Zajc (2018); Schenke  _et al._
(2020a); Schenke (2021) for reviews), and indeed, striking features of
collective expansion have been observed in high-multiplicity events of the
small collision systems studied at RHIC Aidala  _et al._ (2018, 2019) and the
LHC Li (2012); Dusling  _et al._ (2016); Nagle and Zajc (2018).

At RHIC, collisions of proton (p), deuteron (d), and helium-3 (3He)
projectiles on gold (Au) targets were proposed to discern whether “flow-like”
patterns are indeed geometry driven and possibly attributable to mini QGP
droplet formation. The PHENIX Collaboration Aidala  _et al._ (2019) has
reported a large set of elliptic ($v_{2}$) and triangular ($v_{3}$) azimuthal
anisotropy coefficients for p+Au, d+Au, and 3He+Au systems. PHENIX data
clearly shows the hierarchy of $v_{2}$ and $v_{3}$ expected based on the
differences in initial geometry between the three systems.

These collective features can be quantitatively reproduced by (2+1)D viscous
hydrodynamic calculations Habich  _et al._ (2015); Shen  _et al._ (2017a),
which translate the initial spatial anisotropies into final momentum
anisotropies of produced hadrons via the collective expansion of the bulk
matter. In contrast, calculations based solely on initial-state correlations
in the color glass condensate framework Mace  _et al._ (2018, 2019) are ruled
out by the PHENIX data as they predict the opposite hierarchy of $v_{2}$
between the systems. Recently, the STAR Collaboration has also reported
measurements of the differential elliptic flow and triangular flow
coefficients in these three small systems Lacey (2021); STA (2022). STAR
preliminary results show that the observed $v_{3}(p_{T})$ are system
independent, different from the PHENIX published data. In Lacey (2021); STA
(2022) sub-nucleon fluctuations are emphasized to play an important role in
the initial geometry and modify the expectation for the hierarchy between the
flow coefficients in the different systems.

An important difference between the PHENIX and STAR analyses is the use of
different pseudo-rapidity ranges in the two-particle correlation measurements
that yield the azimuthal anisotropy coefficients. PHENIX measures the two-
particle correlations between the backward rapidity ($-3.9<\eta<-3.1$) and
mid-rapidity regions ($-0.35<\eta<0.35$) with the event plane method Adare
_et al._ (2018); Aidala  _et al._ (2019). In contrast, STAR employs only the
mid-rapidity ($-0.9<\eta<0.9$) region, using the scalar-product method Lacey
(2021); STA (2022). It is already well known that in small asymmetric systems,
longitudinal decorrelations are significant Khachatryan  _et al._ (2015) and
boost-invariance is strongly broken Bozek  _et al._ (2015); Schenke and
Schlichting (2016); Shen and Alzhrani (2020); Jiang  _et al._ (2022); Wu and
Qin (2021); Zhao  _et al._ (2022). More recent calculations also indicate that
pre-equilibrium evolution, either in the weak Schenke  _et al._ (2020a) or
strong coupling Romatschke (2015) limit, has a significant effect on the small
and short-lived systems.

In this work, we include the effects of pre-equilibrium evolution and perform
full (3+1)D hydrodynamic simulations, which we argue are essential when
seeking meaningful quantitative comparisons to experimental data.

More specifically, we employ a (3+1)D framework combining a dynamic initial
state with hydrodynamics and hadronic transport Shen and Schenke (2022) to
explore the collective features in small asymmetric systems. The framework has
been shown to provide a unified and quantitative description of identified
hadron production in $\gamma^{*}$-nucleus, proton-proton, proton-nucleus, and
nucleus-nucleus collisions from center of mass energies of a few GeV to
several TeV Shen and Schenke (2022); Zhao  _et al._ (2022). Here we present
detailed results on $dN_{ch}/d\eta$, mean transverse momentum
($\left<p_{T}\right>$), and anisotropic flow in p+Au, d+Au, and 3He+Au
collisions in this (3+1)D hybrid hydrodynamic framework.

This paper is organized as follows: Section II briefly introduces the
3d-glauber+music+urqmd model for the description of the p/d/3He+Au systems.
Section III presents and discusses results from the 3d-glauber+music+urqmd
model on the $dN_{ch}/d\eta$, $\left<p_{T}\right>$, and collective flow.
Section IV concludes this paper with a summary.

##  II The theoretical framework

In this paper, we employ the 3d-glauber+music+urqmd hybrid model within the
open-source iebe-music framework iEB  to study paeticle spectra and collective
flow observables in p+Au, d+Au, and 3He+Au collisions at
$\sqrt{s_{\mathrm{NN}}}=200$ GeV, and for a range of lower beam energies in
d+Au collisions. The hybrid model uses a 3D Monte-Carlo (MC) Glauber initial
condition to dynamically deposit energy, momentum, and net baryon densities
into the evolving fluid system as the two colliding nuclei are crossing each
other Shen and Schenke (2018, 2022). The collective expansion of the QGP
fireball and the evolution of the conserved net-baryon current are described
by a (3+1)D viscous hydrodynamic model music Schenke  _et al._ (2011, 2010);
Paquet  _et al._ (2016); Denicol  _et al._ (2018); MUS . As the QGP expands
and becomes more dilute in the hadronic phase, the fluid dynamic description
is switched to a microscopic hadron cascade model, urqmd Bass  _et al._
(1998); Bleicher  _et al._ (1999); UrQ , to simulate the subsequent evolution
and dynamic decoupling of the hadronic matter.

Figure 1: (Color online) The charged hadron pseudo-rapidity distributions
$dN_{\rm ch}/d\eta$ in various multiplicity classes of p+Au, d+Au and 3He+Au
collisions at $\sqrt{s_{N}}$ = 200 GeV from the 3d-glauber+music+urqmd
simulations, compared to experimental data from the PHENIX Collaboration Adare
_et al._ (2018). Figure 2: (Color online) Identified particle mean transverse
momenta $\left<p_{T}\right>$ as functions of charged hadron multiplicity in
p+Au (solid lines), d+Au (dotted lines) and 3He+Au (dashed lines) collisions
from the 3d-glauber+music+urqmd framework. The d+Au results are compared to
experimental data from the STAR Collaboration Abelev  _et al._ (2009).

More specifically, the dynamical initial condition is simulated by the 3D
Monte-Carlo Glauber model on an event-by-event basis Shen and Schenke (2018,
2022), and the space-time and momentum distributions of the initial energy-
momentum tensor and net baryon charge current are provided by the classical
string deceleration model Bialas  _et al._ (2018); Shen and Schenke (2018).
The transverse positions of valence quarks and soft partonic cloud in the
3d-glauber model are sampled from a 2D Gaussian,
$\exp\left[-\frac{x^{2}+y^{2}}{2}/B^{2}_{G}\right]$, with a transverse nucleon
width parameter $B_{G}=5$ GeV-2, corresponding to a nucleon width of $w=0.44$
fm. After the collision, the deposited energy density distribution has a
Gaussian profile in the transverse plane with a width of $w_{q}=0.2$ fm. We
parameterize the average rapidity loss function of the valence quarks and the
soft partonic cloud with their incoming rapidity $y_{\mathrm{init}}$ in the
collision pair rest frame as Shen and Schenke (2019, 2022),

$\langle
y_{\mathrm{loss}}\rangle(y_{\mathrm{init}})=Ay_{\mathrm{init}}^{\alpha_{2}}[\tanh(y_{\mathrm{init}})]^{\alpha_{1}-\alpha_{2}}.$
(1)

As shown in Fig. 1, when choosing the three parameters values
$A=1.32,\alpha_{1}=1.8$, and $\alpha_{2}=0.34$, we can reproduce the measured
pseudo-rapidity distributions of charged hadrons for p+Au, d+Au, and 3He+Au
collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV. The spatial topography of the
deuteron’s two-nucleon system is obtained from sampling the Hulthen wave
function Adare  _et al._ (2014), and fluctuating 3He configurations come from
results of Green’s function Monte Carlo calculations using the AV18+UIX model
interaction Carlson and Schiavilla (1998). The detailed implementation of this
initial condition model is discussed in Ref. Shen and Schenke (2022).

To incorporate the effects of pre-hydrodynamic flow, we include a finite
transverse initial velocity motivated by the blast-wave model Schnedermann
_et al._ (1993); Teaney  _et al._ (2001); Heinz (2004). We assume the initial
flow is radial and its strength is characterized by the transverse flow
rapidity,

$\eta_{\perp}=\alpha^{\rm pre-flow}r,$ (2)

where $r$ is the radial distance from the string center. The $\alpha^{\rm pre-
flow}$ parameter controls the strength of the pre-hydrodynamic flow. Here we
set $\alpha^{\rm pre-flow}=0.15$ to fit the $v_{n}(p_{T})$ in 3He+Au at 200
GeV. We will explore the effects of pre-hydrodynamic flow on hadronic
observables in p+Au collisions in Appendix B.

The produced strings from individual nucleon-nucleon collisions are treated as
dynamical source terms for the hydrodynamic evolution Okai  _et al._ (2017);
Shen  _et al._ (2017b); Shen and Schenke (2018); Du  _et al._ (2019); Shen and
Schenke (2022),

$\displaystyle\partial_{\mu}T^{\mu\nu}$ $\displaystyle=$ $\displaystyle
J^{\nu}$ (3) $\displaystyle\partial_{\mu}J_{B}^{\mu}$ $\displaystyle=$
$\displaystyle\rho_{B}.$ (4)

Here we use a crossover equation of state (NEOS-BQS) for the QCD matter at
finite chemical potentials, that is constructed using recent lattice data
Borsanyi  _et al._ (2012, 2014); Ding  _et al._ (2015); Bazavov  _et al._
(2017); Monnai  _et al._ (2019). We employ the strangeness neutrality
condition, $n_{s}=0$, and for simplicity set the net electric charge-to-baryon
density ratio to $n_{Q}/n_{B}=0.4$ Monnai  _et al._ (2019), a value most
appropriate for large nuclei. The specific shear and bulk viscosities are
parametrized as

$\frac{\eta T}{e+P}=(\eta/s)_{\rm norm}\left[1+(\eta/s)_{\rm
slope}\left(\frac{\mu_{B}}{\mu_{B,\rm scale}}\right)^{\alpha}\right]$ (5)
$\frac{\zeta T}{e+P}=(\zeta/s)_{\rm norm}{\rm exp}\left[-\left(\frac{T-T_{\rm
peak}}{T_{{\rm width},\lessgtr}}\right)^{2}\right],$ (6)

where $T_{\rm width,<}=0.015$ GeV for $T<T_{\rm peak}$, and $T_{\rm
width,>}=0.1$ GeV for $T>T_{\rm peak}$ with $T_{\rm peak}=0.17$ GeV and
$(\zeta/s)_{\rm norm}=0.08$. For the specific shear viscosity, we set
$(\eta/s)_{\rm norm}=0.13$, $(\eta/s)_{\rm slope}=1.0$, $\mu_{B,\rm
scale}=0.6$ GeV, and $\alpha=1.5$. The shear stress tensor and bulk viscous
pressure are evolved with the Denicol-Niemi-Molnár-Rischke (DNMR) theory with
spatial gradients up to the second order Denicol  _et al._ (2012, 2018). For
simplicity, the effect of charge diffusion is neglected in this work. These
transport coefficients are chosen to reproduce the mean transverse momentum
and the differential flow in central 3He+Au collisions at 200 GeV. In
addition, the $\mu_{B}$-dependence of $\eta/s$ is informed by calibrating flow
observables in Au+Au collisions from 7.7 GeV to 200 GeV San . Compared to Ref.
Shen and Schenke (2022), we use smaller sizes for nucleons and hotspots, which
are consistent with the sub-nucleonic structure constrained by a recent
Bayesian analysis Mäntysaari  _et al._ (2022). The small hotspot size
($w_{q}=0.2$ fm) and inclusion of pre-hydrodynamic flow require a non-zero QGP
bulk viscosity in the hydrodynamic evolution, which was neglected in Shen and
Schenke (2022).

In the iebe-music framework, the Cooper-Frye particlization of the fluid cells
is performed on a hyper-surface with a constant energy density of $e_{\rm
sw}=$0.50 GeV/fm3 using the Cornelius algorithm Huovinen and Petersen (2012)
and the open-source code package iSS Shen  _et al._ (2016); iSS . The non-
equilibrium shear and bulk viscosity corrections to the local equilibrium
distribution function are taken from Grad’s moments method with multiple
conserved charges San . The produced hadrons are then fed into the hadron
cascade model, urqmd, for further scatterings and decays until kinetic freeze-
out is achieved dynamically.

Figure 3: (Color online) The correlation between the flow vectors as a
function of pseudo-rapidity $\eta$ with the reference vector calculated within
$-3.9<\eta^{\rm Ref.}<-3.1$ in p+Au, d+Au and 3He+Au systems from the
3d-glauber+music+urqmd simulations. Figure 4: (Color online) Elliptic flow
$v_{2}$ as a function of pseudo-rapidity ($\eta$) in high-multiplicity 0-5%
central p+Au, d+Au and 3He+Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV
from the 3d-glauber+music+urqmd framework. The results are compared to the
experimental data from the PHENIX Collaboration Adare  _et al._ (2018). Figure
5: (Color online) The anisotropic flow $v_{n}(p_{T})$ as a function of $p_{T}$
in central p+Au, d+Au and 3He+Au collisions computed from the
3d-glauber+music+urqmd framework. The results are compared to the experimental
data from the PHENIX and STAR Collaborations Aidala  _et al._ (2019); Lacey
(2021). The PHENIX data are all 0-5% and STAR data are 0-2%, 0-10%, and 0-10%
for p+Au, d+Au, and 3He+Au collisions, respectively.

##  III RESULTS

In this section, we report results for observables, including charged hadron
pseudo-rapidity distributions, mean transverse momentum, and differential flow
harmonics for the p+Au, d+Au, and 3He+Au collision systems at RHIC energies.

###  III.1 Charged Hadron Multiplicity Rapidity Distribution and Average
Transverse Momentum

Figure 1 shows the charged hadron pseudo-rapidity distributions in different
centrality bins for p+Au, d+Au, and 3He+Au collisions at
$\sqrt{s_{\mathrm{NN}}}=200$ GeV. We note that all parameters were tuned to
reproduce observables in 3He+Au collisions. Keeping all the model parameters
fixed, our model’s extrapolations provide a reasonable description of the
experimental data from the PHENIX Collaboration Adare  _et al._ (2018) for the
p+Au and d+Au systems. In particular, the increase of the asymmetry around
midrapidity from peripheral to central collisions is well captured as a
consequence of the increasing number of participant nucleons from the Au
target towards central collisions. Here we have used the charged hadron
multiplicity in the pseudo-rapidity range $-3.9<\eta<-3.1$ (in the Au-going
direction) to determine the collision centrality bins, which is consistent
with the method used by PHENIX Adare  _et al._ (2018). As pointed out in Ref.
Shen and Schenke (2022), the correlation of particle production between
forward and mid-rapidity is crucial to reproduce the centrality dependence of
charged particles in these asymmetric collision systems. This is highlighted
in Fig. 10 in Appendix A, where we show that using different centrality
definitions generates different shapes for $dN_{ch}/d\eta$.

Within the hydrodynamic paradigm, the mean transverse momentum of identified
hadrons strongly constrains the radial expansion of the fireball. This radial
expansion is very sensitive to the bulk viscosity, initial hotspot size, and
the strength of pre-hydrodynamic flow. Consequently, the mean transverse
momentum provides a strong constraint on $(\zeta/s)(T)$ Ryu  _et al._ (2015,
2018). In Fig. 2 we show the identified particle mean transverse momentum as a
function of charged hadron multiplicity in p+Au, d+Au, and 3He+Au collisions
at 200 GeV and compare to the STAR measurements for d+Au collisions Abelev
_et al._ (2009). Using the temperature-dependent bulk viscosity $(\zeta/s)(T)$
in Eq. (6) with $(\zeta/s)_{\rm norm}$ = 0.08, the calculation shows good
agreement with the STAR data in d+Au collisions Abelev  _et al._ (2009). This
indicates that our hybrid model simulations generate proper amounts of radial
flow and the transverse expansion has been well constrained.

Figure 6: (Color online) The differential elliptic flow $v_{2}(p_{T})$ in
3He/d/p+Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV determined by the
3$\times$2PC method using pseudo-rapidity ranges corresponding to the FVTXS-
CNT-FVTXN (BF) and BBCS-FVTXS-CNT (BB) detector combinations. The model
results are compared to experimental data from the PHENIX Collaboration
Acharya  _et al._ (2022a, b)

The flow harmonics describe the momentum anisotropy of the final produced
hadrons. They are commonly measured using the two-particle correlation method
with a large pseudo-rapidity gap between the particles to suppress non-flow
correlations Bilandzic  _et al._ (2011, 2014). The presence and size of the
gap will affect the measured value of $v_{n}\\{2\\}$ if the correlation
between the flow vectors depends on the pseudo-rapidity separation of two
particles. Therefore, we calculate the $\eta$-dependent correlations between
the flow vectors $Q_{n}$ defined by

${\rm Corr}\left(Q_{n}(\eta),~{}Q^{\rm Ref.}_{n}\right)=\frac{\left\langle
Q^{*}_{n}(\eta)Q_{n}^{\rm
Ref.}\right\rangle}{\sqrt{\left\langle|Q_{n}(\eta)|^{2}\right\rangle}\sqrt{\left\langle|Q_{n}^{\rm
Ref.}|^{2}\right\rangle}},$ (7)

to characterize the longitudinal decorrelation of flow vectors. Here,
$Q_{n}(\eta)$ is the $n$-th order flow vector defined in the $\eta$ bin of
interest, and $Q_{n}^{\rm Ref.}$ is the $n$-th order reference flow vector
within the assigned reference $\eta$ range. Figure 3 shows the correlation
${\rm Corr}\left(Q_{n}(\eta),~{}Q^{\rm Ref.}_{n}\right)$ of the second and
third order flow vectors with their reference flow vectors calculated within
$-3.9<\eta<-3.1$ in central 0-10% p+Au, d+Au, and 3He+Au systems, following
PHENIX’s detector coverage Aidala  _et al._ (2019).

The values of ${\rm Corr}\left(Q_{n}(\eta),~{}Q^{\rm Ref.}_{n}\right)$ for
$|\eta|<0.35$ show the strength of flow vector correlations relevant for the
PHENIX measurement of $v_{n}\\{2\\}$ Aidala  _et al._ (2019). For elliptic
flow, we find that the correlation drops slowly with pseudo-rapidity and is
approximately 0.9 for $|\eta|<0.35$ in d+Au and 3He+Au collisions. However,
the flow vector correlation drops quite rapidly in p+Au collisions, reaching
$\sim 0.6$ in the $|\eta|<0.35$ bin, which indicates that the $v_{2}\\{2\\}$
measurements could depend strongly on the rapidity regions used in the
analysis. For the triangular flow coefficients, the flow correlations are
significantly below unity for all three systems and drop more rapidly with
$\eta$ compared to the case of $n=2$. This result demonstrates that it is non-
trivial to extract the fireball’s initial triangularity hierarchy from
$v_{3}\\{2\\}$ measurements. Because the boost-invariant (2+1)D simulations
can not capture these longitudinal flow correlations, their results for
triangular flow in (p,d,3He)+Au collisions are highly idealized and contain
substantial theoretical uncertainties.

The three panels of Figure 4 show the elliptic flow $v_{2}(\eta)$ as a
function of pseudo-rapidity in the 0–5% centrality bin of
$\sqrt{s_{\mathrm{NN}}}=200$ GeV p+Au, d+Au, and 3He+Au collisions,
respectively. Following the same method as used by the PHENIX Collaboration
Adare  _et al._ (2018), we calculate the second flow harmonic, $v_{2}(\eta)$,
using the event plane (EP) method with the event plane determined in the range
$-3.9<\eta<-3.1$, corresponding to using the Beam Beam Counter (BBC) in the
Au-going direction, and with centrality defined by the charged hadron
multiplicity measured in the same range. The flow harmonics for the event
plane method are calculated using

$v_{n}(\eta,p_{T})\\{\mathrm{EP}\\}=\frac{\langle\cos(n(\phi(\eta,p_{T})-\Psi^{\mathrm{Ref.}}_{n})\rangle}{R(\Psi^{\mathrm{Ref.}}_{n})},$
(8)

where $\phi(\eta,p_{T})$ is the azimuthal angle of particles of interest for
specific $\eta$ and/or $p_{T}$ bins, $\Psi^{\rm Ref.}_{n}$ is the $n$-th order
azimuthal reference event plane, and ${\rm R(\Psi^{Ref.}_{n})}$ is the
resolution of $\Psi^{\rm Ref.}_{n}$. Following the PHENIX measurements Aidala
_et al._ (2018, 2017); Adare  _et al._ (2018), the event plane resolution is
estimated by:

${\rm
R(\Psi^{Ref.}_{n})}=\sqrt{\frac{\langle\frac{Q_{nA}}{|Q_{nA}|}\frac{Q^{*}_{nB}}{|Q_{nB}|}\rangle\langle\frac{Q_{nA}}{|Q_{nA}|}\frac{Q^{*}_{nC}}{|Q_{nC}|}\rangle}{\langle\frac{Q_{nB}}{|Q_{nB}|}\frac{Q^{*}_{nC}}{|Q_{nC}|}\rangle}}$
(9)

with the $n$th-order flow vector,

$Q_{n}=\sum_{k}e^{in\phi_{k}},$ (10)

with the $\sum$ summing over particles in one event. The $Q_{nA,B,C}$ vectors
are defined in the $-3.9<\eta<-3.1$ (south BBC, Au-going side),
$-3.0<\eta<-1.0$ (south Forward Silicon Vertex (FVTX)) and $|\eta|<0.35$
(central - CNT), respectively. The 3d-glauber+music+urqmd model gives a
reasonable description of the $\eta$-dependent elliptic flow in d+Au and
3He+Au systems. Both the experimental data and our model simulations have an
increasing flow coefficient $v_{2}$ towards backward rapidity up to
$\eta=-3.0$ in the Au-going direction. At $\eta<-3.0$, our model predicts a
decrease of $v_{2}(\eta)$. The $v_{2}(\eta)$ of 3He+Au is slightly larger than
that of d+Au collisions because of the longer fireball lifetime in the
hydrodynamic phase in 3He+Au collisions. For the p+Au system, the PHENIX data
shows a pronounced increase of the $v_{2}(\eta)$ towards negative rapidities,
especially for $\eta<-2.0$. Our model fails to describe the $v_{2}(\eta)$ for
$\eta<$ 1.0 in p+Au collisions. We expect $v_{2}(\eta)$ in p+Au collisions to
be smaller than in the other systems because of the stronger decorrelation and
in p+Au collisions, compared to d+Au and 3He+Au collisions, as shown in Fig.
3. Also the shorter lifetime in p+Au collisions can lead to smaller
$v_{2}(\eta)$.

Figure 5 shows the computed $p_{T}$-differential flow $v_{n}(p_{T})$ $(n=2,3)$
for charged hadrons compared to experimental data from the PHENIX and STAR
Collaborations. Here, we adopt the same methods and kinematic cuts as employed
in the PHENIX or STAR analyses, respectively. Specifically, when comparing to
the PHENIX data Aidala  _et al._ (2019), we select collision events for 0-5%
centrality using the number of charged hadrons in $\eta\in[-3.9,-3.1]$. The
$v_{n}(p_{T})\\{\rm EP\\}$ (PHENIX) for charged hadrons in the midrapidity
region covering $|\eta|<0.35$ (central CNT) are calculated with the event
plane method by Eq. (8), where the second-order event plane is determined in
the Au-going direction, in the pseudo-rapidity range $-3.0<\eta<-1.0$ (south
FVTX) in p/d+Au and $-3.9<\eta<-3.1$ (south BBC) in 3He+Au collisions, and the
third order event planes are determined in the range of the south BBC for all
systems Aidala  _et al._ (2019). The event plane resolutions are calculated
with the three-subevent method, which correlates measurements in the south
BBC, south FVTX, and central CNT arms Aidala  _et al._ (2019); Adare  _et al._
(2018) as given explicitly in Eq. (9). The centrality classes are defined
using the multiplicity of charged hadrons in the pseudo-rapidity range
$-3.9<\eta<-3.1$.

To compare with the experimental data from the STAR Collaboration Lacey
(2021); STA (2022), we select the 0-2% most central p+Au collisions and 0-10%
most central d+Au and 3He+Au collisions, based on the charged hadron
multiplicity in $\eta\in[-0.9,0.9]$. The STAR flow measurements imposed a
rapidity gap $|\Delta\eta|>1$ for charged particles in $\eta\in[-0.9,0.9]$.
For convenience we compute the two-particle flow coefficients using the
scalar-product (SP) method with particles of interest (POIs) and reference
particles (RPs) selected from two sub-events within $-0.9<\eta<-0.5$ and
$0.5<\eta<0.9$, respectively. We confirmed that this approximates STAR’s
method very well.

As mentioned above, we tune parameters to fit observables, including
$v_{2}(p_{T})$ and $v_{3}(p_{T})$, in 3He+Au collisions, and then predict
observables in p+Au and d+Au collisions. Fig. 5 shows that the
3d-glauber+music+urqmd model gives an overall good description of the PHENIX
$v_{n}(p_{T})$ data for d+Au and 3He+Au collisions. We underestimate the
$v_{2}(p_{T})$ in p+Au collisions by about 20-30%, possibly because of a too
large longitudinal flow decorrelation in the model. Compared to the STAR
elliptic flow, our model provides a reasonable description for $p_{T}<1$ GeV
for (d, 3He)+Au collisions, but underestimates the data for p+Au collisions.
Our calculation overestimates the STAR $v_{2}(p_{T})$ at higher $p_{T}$. The
STAR $v_{3}(p_{T})$ measurements show a much weaker system size dependence
than the model, which gives a good description of $v_{3}(p_{T})$ in 3He+Au but
significantly underestimates $v_{3}(p_{T})$ in (p, d)+Au collisions.

Both in the model calculations and experimental data, the $v_{3}(p_{T})$ with
the STAR definition (stars) are systemically larger than those determined
using the PHENIX definition (squares). In our calculations, this difference is
mainly caused by the different magnitudes of the longitudinal decorrelation of
flow vectors of $v_{3}$ between the different pseudo-rapidity bins. That is,
the decorrelation between the flow vectors at $-3.9<\eta<-3.1$ and
$-0.35<\eta<0.35$ is greater than that between $-0.9<\eta<-0.5$ and
$0.5<\eta<0.9$.

However, our model still underestimates the STAR data of $v_{3}(p_{T})$ in
p+Au and d+Au collisions, indicating that the decorrelations present in the
model do not fully account for the differences between the STAR and PHENIX
results. We leave a more systematic calibration of fitting both STAR and
PHENIX data for future studies using a Bayesian framework. This study should
shed more light on the role of any residual non-flow in the experimental data.

Recently, the PHENIX Collaboration measured anisotropic flow coefficients
using the $3\times$2PC method, whose two-particle azimuthal correlations are
constructed with three different sets of pairs Acharya  _et al._ (2022a, b).
Taking the $p_{T}$-integrated flow vector $Q_{nA,B,C}$, and the
$p_{T}$-differential flow vector $q_{nA,B,C}(p_{T})$ of particles of interest
(using finite $p_{T}$ bins) from different sub-events $A$, $B$, and $C$ that
indicate different pseudo-rapidity bins, the differential $v_{n}(p_{T})$ of
the particles of interest in the single sub-event $C$ can be determined as

$v_{n}^{C}(p_{T})=\sqrt{\frac{c_{n}^{AC}(p_{T})c_{n}^{BC}(p_{T})}{C_{n}^{AB}}}.$
(11)

where the $p_{T}$-differential $c_{n}(p_{T})$ and $p_{T}$-integrated $C_{n}$
coefficients are defined as

$\displaystyle C_{n}^{AB}$ $\displaystyle=\left\langle
Q_{nA}Q^{*}_{nB}\right\rangle,$ (12) $\displaystyle c_{n}^{AC}(p_{T})$
$\displaystyle=\left\langle Q_{nA}q^{*}_{nC}(p_{T})\right\rangle,$ (13)
$\displaystyle c_{n}^{BC}(p_{T})$ $\displaystyle=\left\langle
Q_{nB}q^{*}_{nC}(p_{T})\right\rangle.$ (14)

Here $\left\langle...\right\rangle$ represents event averaging and taking the
real part of the correlator.

Figure 6 shows the differential elliptic flow $v_{2}(p_{T})\\{3\times 2{\rm
PC}\\}$ at mid-rapidity $|\eta|<0.35$ in 3He+Au, d+Au, and p+Au systems,
respectively. Following the PHENIX measurements Acharya  _et al._ (2022a, b),
we calculate the $v_{2}(p_{T})$ from the 3×2PC method with two combinations,
“$-3.0<\eta<-1.0,~{}-0.35<\eta<0.35~{}{\rm and}~{}1.0<\eta<3.0$”, FVTXS-CNT-
FVTXN, detector combination (BF) and “$-3.9<\eta<-3.1,~{}-3.0<\eta<-1.0~{}{\rm
and}~{}-0.35<\eta<0.35$”, BBCS-FVTXS-CNT, detector combination (BB).

The 3d-glauber+music+urqmd model reproduces the $v_{2}(p_{T})\\{3\times 2{\rm
PC}\\}$ in d+Au and 3He+Au collisions quite well up to 60% in centrality. The
model calculations produce a slowly decreasing trend in $v_{2}(p_{T})$’s
magnitude from central to peripheral collisions for both the BF and BB
combinations. This trend is caused by the shorter hydrodynamic lifetimes in
more peripheral collisions. More importantly, our model quantitatively
reproduces the difference between $v_{2}(p_{T})({\rm BF})$ and
$v_{2}(p_{T})({\rm BB})$, which suggests that our model predicts the right
amount of longitudinal flow decorrelation in d+Au and 3He+Au collisions from
most central to middle central collisions. Our model results underestimate the
measurements in the peripheral collisions, for example in 40-60% bin, where we
expect sizable non-flow correlations in the measurements.

Similar to the previous $v_{n}\\{2\\}$ comparison, our calculations
underestimate the p+Au $v_{2}(p_{T})$ by 40% across all centrality bins.
Interestingly, our model can reasonably reproduce the difference between
$v_{2}(p_{T})({\rm BF})$ and $v_{2}(p_{T})({\rm BB})$ in 0-5% and 5-10%
central p+Au collisions. Furthermore, this $v_{2}(p_{T})$ difference between
BF and BB is larger in p+Au collisions than in (d, 3He)+Au collisions,
suggesting a large longitudinal flow decorrelation in p+Au collisions, as
present in our model. The PHENIX $v_{2}(p_{T})({\rm BF})$ for 20% central or
more peripheral p+Au collisions show rapid increases at high $p_{T}$,
suggesting sizable residual non-flow correlations.

Figure 7: (Color online) The charged hadron pseudo-rapidity distributions
$dN_{\rm ch}/d\eta$ of central d+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ =
200, 62.4, 39 and 19.6 GeV from the 3d-glauber+music+urqmd simulations,
compared to experimental data from the PHENIX Collaboration Aidala  _et al._
(2017). Figure 8: (Color online) Charged hadron elliptic flow $v_{2}\\{2\\}$,
$v_{2}\\{4\\}$ and $v_{2}\\{2,|\Delta\eta|>2.0\\}$ as a function of charged
hadron multiplicity $N_{\rm Track}$ of d+Au collisions at
$\sqrt{s_{\mathrm{NN}}}=200$, 62.4, 39 and 19.6 GeV from the
3d-glauber+music+urqmd framework. The results are compared to experimental
data from the PHENIX Collaboration Aidala  _et al._ (2018). Figure 9: (Color
online) Elliptic flow $v_{2}$ as a function of pseudo-rapidity in high-
multiplicity central d+Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$, 62.4, 39
and 19.6 GeV from the 3d-glauber+music+urqmd framework. The results are
compared to the experimental data from the PHENIX Collaboration Aidala  _et
al._ (2017).

###  III.2 Results for the d+Au Beam Energy Scan

In this sub-section, we extrapolate our model to d+Au collisions at lower
collision energies. To do so we only reduce the center-of-mass collision
energy and keep all model parameters fixed.

Figure 7 shows our model predictions for the collision energy dependence of
the charged hadron pseudorapidity distribution in central d+Au collisions. Our
model gives an overall good description of the asymmetric charged hadron
pseudorapidity distributions compared with the PHENIX data from
$\sqrt{s_{\mathrm{NN}}}=$ 200 GeV to 19.6 GeV. The peaks of $dN_{ch}/d\eta$ in
the theoretical calculations are shifted more towards midrapidity compared
with the PHENIX measurements at low collision energies, suggesting a small
overestimation of the rapidity loss in the initial-state model for the lower
collision energies. Compared to the results shown in Ref. Shen and Schenke
(2022), the current $dN_{ch}/d\eta$ at 19.6 GeV is slightly smaller because a
higher switching energy density $e_{\rm sw}=0.5$ GeV/fm3 is used here for all
collision energies.

Figure 8 shows the model-to-data comparisons for $v_{2}\\{2\\}$,
$v_{2}\\{4\\}$ and $v_{2}\\{2,|\Delta\eta|>2.0\\}$ in d+Au collisions at 200,
62.4, 39, and 19.6 GeV, using the PHENIX FVTX ($1.0<|\eta|<3.0$) acceptance
Aidala  _et al._ (2018). Here, following the PHENIX measurements Aidala  _et
al._ (2018); Adare  _et al._ (2019), the $N_{\rm Track}$ is estimated by
$N_{\rm Track}=p_{0}\times N_{ch}$, with $p_{0}=0.67$ being the tracking
efficiency and acceptance in the FVTX detector. The $N_{ch}$ is the charged
hadron multiplicity within the FVTX acceptance $1<|\eta|<3$ computed in our
model. With a rapidity gap of $|\Delta\eta|>2.0$, non-flow correlations can be
effectively suppressed.

We find our calculations agree nicely with the PHENIX
$v_{2}\\{2,|\Delta\eta|>2.0\\}$ data for all four collision energies. The
difference between $v_{2}\\{2\\}$ and $v_{2}\\{2,|\Delta\eta|>2.0\\}$ data
indicate the contributions from event plane decorrelation and non-flow
correlations. In our calculations, the difference between these two
correlation observables only comes from the event plane decorrelations.
Therefore, the difference between the model and PHENIX $v_{2}\\{2\\}$ can
serve as an estimation of non-flow correlations in the measurements.

The elliptic flow measured using the four-particle cumulant method,
$v_{2}\\{4\\}$, in general, contains fewer non-flow correlations. Our model
results overestimate $v_{2}\\{4\\}$ for central d+Au collisions at 200 GeV.
For the lower collision energies, our $v_{2}\\{4\\}$ underestimates the PHENIX
data. At the current stage, it is not clear whether the discrepancy originates
from possible residual non-flow correlation in the measurements or different
variances in the elliptic flow distributions.

Finally, Figure 9 shows the rapidity-dependent $v_{2}(\eta)\\{\rm EP\\}$ in
central d+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200, 62.4, 39 GeV, and
19.6 GeV, calculated using the same event plane method as that in Fig. 4. Our
model produces results close to the PHENIX measurements at the top three
collision energies. Our $v_{2}(\eta)\\{\rm EP\\}$ at 19.6 GeV serves as a
model prediction. The $v_{2}(\eta)\\{\rm EP\\}$ at all four energies have
similar shapes with respect to $dN_{\rm ch}/d\eta$ that decreases with
increasing $\eta$ between $0<\eta<3$ in the proton-going side, and peaks at
the Au-going side. This correlation between $v_{2}(\eta)$ and $dN_{\rm
ch}/d\eta$ reflects that the size of the elliptic flow is proportional to the
fireball lifetime in the hydrodynamic phase Aidala  _et al._ (2017).

##  IV Summary

We explored the collectivity in asymmetric (p, d, 3He)+Au collisions at 200
GeV and d+Au collisions at 19.6 - 200 GeV within a (3+1)D hybrid theoretical
framework. We find that the elliptic flow vector correlations in the PHENIX
kinematic regions remain $\sim 0.9$ for central (d, 3He)+Au collisions at 200
GeV, ensuring a good signal of the geometric response from initial 2D
ellipticity to final-state charged hadron elliptic flow. However, the elliptic
flow correlation drops to 0.6 in central p+Au collisions. For triangular flow,
the flow correlations are significantly below unity for all three collision
systems. The presence of substantial flow decorrelations leads to the
conclusion that full (3+1)D simulations are required for any quantitative
comparisons to experimental data from these asymmetric collision systems.

With parameters calibrated using the central 3He+Au measurements, the
3d-glauber+music+urqmd model reproduces the PHENIX two-particle $v_{n}(p_{T})$
in central d+Au collisions, but underestimates the p+Au $v_{n}(p_{T})$ by 30%.
By computing $v_{n}(p_{T})$ using the STAR definitions, we find a 30% larger
$v_{3}(p_{T})$ with the STAR definition compared to the triangular flow from
the PHENIX definition. Compared to the STAR measurements, our calculations
suggest that approximately 50% of the difference between the PHENIX and STAR
$v_{3}(p_{T})$ originates from the different flow correlations between
different rapidity regions. The remaining discrepancy between our calculations
and the experimental data could come from different origins, including pre-
hydrodynamic flow and non-flow in the experimental data that is not present in
the model.

We further obtain a reasonable agreement with PHENIX measurements of
differential elliptic flow using the $3\times$2PC method in d+Au and 3He+Au
collisions up to 60% in centrality. We underestimate the $v_{2}(p_{T})$ in
p+Au by approximately 40%. However, our model quantitatively reproduces the
difference between the $v_{2}(p_{T})$ measured using the BF and BB
combinations of detectors for all three systems. The larger observed
difference between $v_{2}(p_{T})({\rm BF})$ and $v_{2}(p_{T})({\rm BB})$ in
p+Au collisions compared to (d, 3He)+Au collisions suggests a larger
longitudinal flow decorrelation in p+Au collisions, as present in our model.

Lastly, we extrapolate our model calculations to d+Au collisions at 19.6, 39,
and 62.4 GeV. With all the model parameters fixed by the 3He+Au collisions at
200 GeV, our model predictions give good descriptions of PHENIX measurements
of $v_{2}\\{2\\}$ with a rapidity gap and the rapidity dependent elliptic flow
in d+Au at the lower collision energies.

In summary, our study explicitly demonstrates that longitudinal flow
decorrelations play a central role in anisotropic flow measurements in
asymmetric nuclear collisions. We conclude that full (3+1)D hybrid simulations
are essential for making quantitative comparisons with observables from small
system collisions at RHIC.

## Acknowledgements

We thank Ron Belmont for providing the PHENIX data and Mashhood Munir for
estimates of the pre-hydrodynamic flow. We thank Shengli Huang, Ulrich W.
Heinz, Jiangyong Jia and Zhengbu Xu for helpful discussions and suggestions.
W.B.Z. is supported by the National Science Foundation (NSF) under grant
numbers ACI-2004571 within the framework of the XSCAPE project of the JETSCAPE
collaboration. S.R. is supported by the U.S. Department of Energy(DOE) under
grant number DE-SC0013460. B.P.S. and C.S. are supported by the U.S.
Department of Energy, Office of Science, Office of Nuclear Physics, under DOE
Contract No. DE-SC0012704 and Award No. DE-SC0021969, respectively. C.S.
acknowledges support from a DOE Office of Science Early Career Award. This
work is in part supported within the framework of the Beam Energy Scan Theory
(BEST) Topical Collaboration and under contract number DE-SC0013460. This
research was done using resources provided by the Open Science Grid (OSG)
Pordes  _et al._ (2007); Sfiligoi  _et al._ (2009), which is supported by the
National Science Foundation award #2030508.

##  Appendix A Centrality Selection Dependence of $dN_{ch}/d\eta$ in
Asymmetric Collisions

We study the effects of different centrality definitions on the charged hadron
pseudo-rapidity distributions in 0-5% central (p, d, 3He)+Au collisions at 200
GeV. Figure 10 shows the comparison of the $dN_{ch}/d\eta$ between the PHENIX
centrality definition (determined by charged hadron multiplicity within
$-3.9<\eta<-3.1$) and the STAR centrality definition (determined by charged
hadron multiplicity within $-0.9<\eta<0.9$). The STAR centrality definition
gives a larger yield at mid-rapidity than that using PHENIX’s definition,
while the central events according to the PHENIX definition have larger yields
in the backward (negative) rapidity region. Figure 10 illustrates that it is
crucial to adopt the same centrality definition as the experiment for model-
to-data comparisons in these asymmetric collisions, emphasizing again the need
for full (3+1)d hybrid simulations.

##  Appendix B Importance of Pre-Hydrodynamic Evolution for Anisotropic Flow
in p+Au Collisions

We investigate the importance of the pre-hydrodynamic flow for hadronic
observables in p+Au collisions. Figs. 11 and 12 show the identified particle
mean transverse momenta and the $v_{n}(p_{T})$ ($n=2,3$) with different values
of the $\alpha^{\rm pre-flow}$ defined in Eq. (2) when simulating central p+Au
collisions at 200 GeV. Here the $\alpha^{\rm pre-flow}$ parameter controls the
strength of the pre-hydrodynamic flow. Figure 11 shows that a larger pre-
hydrodynamic flow leads to a faster transverse expansion, which results in
larger $\langle p_{T}\rangle$ for identified hadrons. The proton’s mean
$p_{T}$ shows a stronger sensitivity to $\alpha^{\rm pre-flow}$ than that of
pions and kaons. It’s also noteworthy that the faster fireball expansion with
a larger $\alpha^{\rm pre-flow}$ shortens the overall fireball lifetime and
its longitudinal expansion. Therefore, there is slightly more energy left at
mid-rapidity to produce more charged hadrons.

Figure 12 shows that the pre-hydrodynamic flow has significant effects on the
charged hadron anisotropic flow in p+Au collisions. At $p_{T}\sim 1$ GeV, the
elliptic flow increases by 50% and triangular flow increases by 300% when we
change from $\alpha^{\rm pre-flow}=0.15$ to $\alpha^{\rm pre-flow}=1.5$. The
stronger sensitivity of higher order anisotropic flow to pre-hydrodynamic flow
is in qualitative agreement with early results from boost-invariant
simulations Romatschke (2015). With $\alpha^{\rm pre-flow}=1.5$, our
simulations can reproduce the STAR $v_{3}(p_{T})$ data but they overestimate
the $v_{2}(p_{T})$ for $p_{T}>0.6$ GeV at the same time. This result is in
qualitative agreement with results shown in Ref. Schenke  _et al._ (2020b).

Figure 10: (Color online) The charged hadron pseudo-rapidity distributions
$dN_{\rm ch}/d\eta$ with different centrality selections of central p+Au, d+Au
and 3He+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV from the
3d-glauber+music+urqmd simulations. Figure 11: (Color online) Identified
particle mean transverse momenta $\left<p_{T}\right>$ as functions of charged
hadron multiplicity in p+Au collisions with different values of $\alpha^{\rm
pre-flow}$ from the 3d-glauber+music+urqmd framework. The results are compared
to experimental d+Au data from the STAR Collaboration Abelev  _et al._ (2009).
Figure 12: (Color online) The charged hadron anisotropic flow $v_{n}(p_{T})$
as a function of $p_{T}$ in central p+Au collisions with different values of
$\alpha^{\rm pre-flow}$ computed from the 3d-glauber+music+urqmd framework.
The results are compared to experimental data from the STAR Collaboration
Lacey (2021).

## References

  * Heinz and Snellings (2013) Ulrich Heinz and Raimond Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions,” Ann. Rev. Nucl. Part. Sci. 63, 123–151 (2013), arXiv:1301.2826 [nucl-th] .
  * Gale  _et al._ (2013) Charles Gale, Sangyong Jeon, and Bjoern Schenke, “Hydrodynamic Modeling of Heavy-Ion Collisions,” Int. J. Mod. Phys. A 28, 1340011 (2013), arXiv:1301.5893 [nucl-th] .
  * Shen and Heinz (2015) Chun Shen and Ulrich Heinz, “The road to precision: Extraction of the specific shear viscosity of the quark-gluon plasma,” Nucl. Phys. News 25, 6–11 (2015), arXiv:1507.01558 [nucl-th] .
  * Shen and Yan (2020) Chun Shen and Li Yan, “Recent development of hydrodynamic modeling in heavy-ion collisions,” Nucl. Sci. Tech. 31, 122 (2020), arXiv:2010.12377 [nucl-th] .
  * Aad  _et al._ (2021) Georges Aad _et al._ (ATLAS), “Two-particle azimuthal correlations in photonuclear ultraperipheral Pb+Pb collisions at 5.02 TeV with ATLAS,” Phys. Rev. C 104, 014903 (2021), arXiv:2101.10771 [nucl-ex] .
  * Zhao  _et al._ (2022) Wenbin Zhao, Chun Shen, and Björn Schenke, “Collectivity in Ultra-Peripheral Pb+Pb Collisions at the Large Hadron Collider,” (2022), arXiv:2203.06094 [nucl-th] .
  * Shen  _et al._ (2022) Chun Shen, Wenbin Zhao, and Björn Schenke, “Collectivity in ultra-peripheral heavy-ion collisions,” in  _20th International Conference on Strangeness in Quark Matter 2022_ (2022) arXiv:2209.15065 [nucl-th] .
  * Dusling  _et al._ (2016) Kevin Dusling, Wei Li, and Björn Schenke, “Novel collective phenomena in high-energy proton–proton and proton–nucleus collisions,” Int. J. Mod. Phys. E 25, 1630002 (2016), arXiv:1509.07939 [nucl-ex] .
  * Loizides (2016) Constantin Loizides, “Experimental overview on small collision systems at the LHC,” Nucl. Phys. A 956, 200–207 (2016), arXiv:1602.09138 [nucl-ex] .
  * Schlichting and Tribedy (2016) Sören Schlichting and Prithwish Tribedy, “Collectivity in Small Collision Systems: An Initial-State Perspective,” Adv. High Energy Phys. 2016, 8460349 (2016), arXiv:1611.00329 [hep-ph] .
  * Nagle and Zajc (2018) James L. Nagle and William A. Zajc, “Small System Collectivity in Relativistic Hadronic and Nuclear Collisions,” Ann. Rev. Nucl. Part. Sci. 68, 211–235 (2018), arXiv:1801.03477 [nucl-ex] .
  * Schenke  _et al._ (2020a) Bjoern Schenke, Chun Shen, and Prithwish Tribedy, “Hybrid Color Glass Condensate and hydrodynamic description of the Relativistic Heavy Ion Collider small system scan,” Phys. Lett. B 803, 135322 (2020a), arXiv:1908.06212 [nucl-th] .
  * Schenke (2021) Björn Schenke, “The smallest fluid on Earth,” Rept. Prog. Phys. 84, 082301 (2021), arXiv:2102.11189 [nucl-th] .
  * Aidala  _et al._ (2018) C. Aidala _et al._ (PHENIX), “Measurements of Multiparticle Correlations in $d+\mathrm{Au}$ Collisions at 200, 62.4, 39, and 19.6 GeV and $p+\mathrm{Au}$ Collisions at 200 GeV and Implications for Collective Behavior,” Phys. Rev. Lett. 120, 062302 (2018), arXiv:1707.06108 [nucl-ex] .
  * Aidala  _et al._ (2019) C. Aidala _et al._ (PHENIX), “Creation of quark–gluon plasma droplets with three distinct geometries,” Nature Phys. 15, 214–220 (2019), arXiv:1805.02973 [nucl-ex] .
  * Li (2012) Wei Li, “Observation of a ’Ridge’ correlation structure in high multiplicity proton-proton collisions: A brief review,” Mod. Phys. Lett. A 27, 1230018 (2012), arXiv:1206.0148 [nucl-ex] .
  * Habich  _et al._ (2015) M. Habich, J. L. Nagle, and P. Romatschke, “Particle spectra and HBT radii for simulated central nuclear collisions of C + C, Al + Al, Cu + Cu, Au \+ Au, and Pb + Pb from $\sqrt{s}=62.4$ \- $2760$ GeV,” Eur. Phys. J. C 75, 15 (2015), arXiv:1409.0040 [nucl-th] .
  * Shen  _et al._ (2017a) Chun Shen, Jean-François Paquet, Gabriel S. Denicol, Sangyong Jeon, and Charles Gale, “Collectivity and electromagnetic radiation in small systems,” Phys. Rev. C 95, 014906 (2017a), arXiv:1609.02590 [nucl-th] .
  * Mace  _et al._ (2018) Mark Mace, Vladimir V. Skokov, Prithwish Tribedy, and Raju Venugopalan, “Hierarchy of Azimuthal Anisotropy Harmonics in Collisions of Small Systems from the Color Glass Condensate,” Phys. Rev. Lett. 121, 052301 (2018), [Erratum: Phys.Rev.Lett. 123, 039901 (2019)], arXiv:1805.09342 [hep-ph] .
  * Mace  _et al._ (2019) Mark Mace, Vladimir V. Skokov, Prithwish Tribedy, and Raju Venugopalan, “Systematics of azimuthal anisotropy harmonics in proton–nucleus collisions at the LHC from the Color Glass Condensate,” Phys. Lett. B 788, 161–165 (2019), [Erratum: Phys.Lett.B 799, 135006 (2019)], arXiv:1807.00825 [hep-ph] .
  * Lacey (2021) Roy A. Lacey (STAR), “Long-range collectivity in small collision-systems with two- and four-particle correlations @ STAR,” Nucl. Phys. A 1005, 122041 (2021), arXiv:2002.11889 [nucl-ex] .
  * STA (2022) “Measurements of the elliptic and triangular azimuthal anisotropies in central 3He+Au, $d$+Au and $p$+Au collisions at $\sqrt{s_{\rm NN}}$ = 200 GeV,” (2022), arXiv:2210.11352 [nucl-ex] .
  * Adare  _et al._ (2018) A. Adare _et al._ (PHENIX), “Pseudorapidity Dependence of Particle Production and Elliptic Flow in Asymmetric Nuclear Collisions of $p+$Al, $p+$Au, $d+$Au, and 3He$+$Au at $\sqrt{s_{{}_{NN}}}=200$ GeV,” Phys. Rev. Lett. 121, 222301 (2018), arXiv:1807.11928 [nucl-ex] .
  * Khachatryan  _et al._ (2015) Vardan Khachatryan _et al._ (CMS), “Evidence for transverse momentum and pseudorapidity dependent event plane fluctuations in PbPb and pPb collisions,” Phys. Rev. C 92, 034911 (2015), arXiv:1503.01692 [nucl-ex] .
  * Bozek  _et al._ (2015) Piotr Bozek, Adam Bzdak, and Guo-Liang Ma, “Rapidity dependence of elliptic and triangular flow in proton–nucleus collisions from collective dynamics,” Phys. Lett. B 748, 301–305 (2015), arXiv:1503.03655 [hep-ph] .
  * Schenke and Schlichting (2016) Bjoern Schenke and Soeren Schlichting, “3D glasma initial state for relativistic heavy ion collisions,” Phys. Rev. C 94, 044907 (2016), arXiv:1605.07158 [hep-ph] .
  * Shen and Alzhrani (2020) Chun Shen and Sahr Alzhrani, “Collision-geometry-based 3D initial condition for relativistic heavy-ion collisions,” Phys. Rev. C 102, 014909 (2020), arXiv:2003.05852 [nucl-th] .
  * Jiang  _et al._ (2022) Ze-Fang Jiang, Shanshan Cao, Xiang-Yu Wu, C. B. Yang, and Ben-Wei Zhang, “Longitudinal distribution of initial energy density and directed flow of charged particles in relativistic heavy-ion collisions,” Phys. Rev. C 105, 034901 (2022), arXiv:2112.01916 [hep-ph] .
  * Wu and Qin (2021) Xiang-Yu Wu and Guang-You Qin, “Asymmetric longitudinal flow decorrelations in proton-nucleus collisions,” (2021), arXiv:2109.03512 [hep-ph] .
  * Romatschke (2015) Paul Romatschke, “Light-Heavy Ion Collisions: A window into pre-equilibrium QCD dynamics?” Eur. Phys. J. C 75, 305 (2015), arXiv:1502.04745 [nucl-th] .
  * Shen and Schenke (2022) Chun Shen and Björn Schenke, “Longitudinal dynamics and particle production in relativistic nuclear collisions,” Phys. Rev. C 105, 064905 (2022), arXiv:2203.04685 [nucl-th] .
  * (32) The iEBE-MUSIC is a general-purpose numerical framework to simulate dynamical evolution of relativistic heavy-ion collisions event-by-event. This work uses v0.5 of this framework, which can be downloaded from https://github.com/chunshen1987/iEBE-MUSIC.
  * Shen and Schenke (2018) Chun Shen and Björn Schenke, “Dynamical initial state model for relativistic heavy-ion collisions,” Phys. Rev. C 97, 024907 (2018), arXiv:1710.00881 [nucl-th] .
  * Schenke  _et al._ (2011) Bjorn Schenke, Sangyong Jeon, and Charles Gale, “Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics,” Phys. Rev. Lett. 106, 042301 (2011), arXiv:1009.3244 [hep-ph] .
  * Schenke  _et al._ (2010) Bjoern Schenke, Sangyong Jeon, and Charles Gale, “(3+1)D hydrodynamic simulation of relativistic heavy-ion collisions,” Phys. Rev. C 82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Paquet  _et al._ (2016) Jean-François Paquet, Chun Shen, Gabriel S. Denicol, Matthew Luzum, Björn Schenke, Sangyong Jeon, and Charles Gale, “Production of photons in relativistic heavy-ion collisions,” Phys. Rev. C 93, 044906 (2016), arXiv:1509.06738 [hep-ph] .
  * Denicol  _et al._ (2018) Gabriel S. Denicol, Charles Gale, Sangyong Jeon, Akihiko Monnai, Björn Schenke, and Chun Shen, “Net baryon diffusion in fluid dynamic simulations of relativistic heavy-ion collisions,” Phys. Rev. C 98, 034916 (2018), arXiv:1804.10557 [nucl-th] .
  * (38) MUSIC is the numerical implementation of (3+1)D relativistic viscous hydrodynamic simulations for high energy heavy-ion collisions. Its official website is http://www.physics.mcgill.ca/music. This work uses v2.5 of this framework, which can be downloaded from https://github.com/MUSIC-fluid/MUSIC.
  * Bass  _et al._ (1998) S. A. Bass _et al._ , “Microscopic models for ultrarelativistic heavy ion collisions,” Prog. Part. Nucl. Phys. 41, 255–369 (1998), arXiv:nucl-th/9803035 .
  * Bleicher  _et al._ (1999) M. Bleicher _et al._ , “Relativistic hadron hadron collisions in the ultrarelativistic quantum molecular dynamics model,” J. Phys. G 25, 1859–1896 (1999), arXiv:hep-ph/9909407 .
  * (41) We use the official UrQMD v3.4 and set it up to run as the afterburner mode, https://bitbucket.org/Chunshen1987/urqmd_afterburner/src/master/.
  * Abelev  _et al._ (2009) B. I. Abelev _et al._ (STAR), “Systematic Measurements of Identified Particle Spectra in $pp,d^{+}$ Au and Au+Au Collisions from STAR,” Phys. Rev. C 79, 034909 (2009), arXiv:0808.2041 [nucl-ex] .
  * Bialas  _et al._ (2018) Andrzej Bialas, Adam Bzdak, and Volker Koch, “Stopped nucleons in configuration space,” Acta Phys. Polon. B 49, 103 (2018), arXiv:1608.07041 [hep-ph] .
  * Shen and Schenke (2019) Chun Shen and Björn Schenke, “Dynamical initialization and hydrodynamic modeling of relativistic heavy-ion collisions,” Nucl. Phys. A 982, 411–414 (2019), arXiv:1807.05141 [nucl-th] .
  * Adare  _et al._ (2014) A. Adare _et al._ (PHENIX), “Centrality categorization for $R_{p(d)+A}$ in high-energy collisions,” Phys. Rev. C 90, 034902 (2014), arXiv:1310.4793 [nucl-ex] .
  * Carlson and Schiavilla (1998) J. Carlson and R. Schiavilla, “Structure and Dynamics of Few Nucleon Systems,” Rev. Mod. Phys. 70, 743–842 (1998).
  * Schnedermann  _et al._ (1993) Ekkard Schnedermann, Josef Sollfrank, and Ulrich W. Heinz, “Thermal phenomenology of hadrons from 200-A/GeV S+S collisions,” Phys. Rev. C 48, 2462–2475 (1993), arXiv:nucl-th/9307020 .
  * Teaney  _et al._ (2001) D. Teaney, J. Lauret, and Edward V. Shuryak, “Flow at the SPS and RHIC as a quark gluon plasma signature,” Phys. Rev. Lett. 86, 4783–4786 (2001), arXiv:nucl-th/0011058 .
  * Heinz (2004) Ulrich W. Heinz, “Concepts of heavy ion physics,” in  _2nd CERN-CLAF School of High Energy Physics_ (2004) pp. 165–238, arXiv:hep-ph/0407360 .
  * Okai  _et al._ (2017) Michito Okai, Koji Kawaguchi, Yasuki Tachibana, and Tetsufumi Hirano, “New approach to initializing hydrodynamic fields and mini-jet propagation in quark-gluon fluids,” Phys. Rev. C 95, 054914 (2017), arXiv:1702.07541 [nucl-th] .
  * Shen  _et al._ (2017b) Chun Shen, Gabriel Denicol, Charles Gale, Sangyong Jeon, Akihiko Monnai, and Bjoern Schenke, “A hybrid approach to relativistic heavy-ion collisions at the RHIC BES energies,” Nucl. Phys. A 967, 796–799 (2017b), arXiv:1704.04109 [nucl-th] .
  * Du  _et al._ (2019) Lipei Du, Ulrich Heinz, and Gojko Vujanovic, “Hybrid model with dynamical sources for heavy-ion collisions at BES energies,” Nucl. Phys. A 982, 407–410 (2019), arXiv:1807.04721 [nucl-th] .
  * Borsanyi  _et al._ (2012) Szabolcs Borsanyi, Zoltan Fodor, Sandor D. Katz, Stefan Krieg, Claudia Ratti, and Kalman Szabo, “Fluctuations of conserved charges at finite temperature from lattice QCD,” JHEP 01, 138 (2012), arXiv:1112.4416 [hep-lat] .
  * Borsanyi  _et al._ (2014) Szabocls Borsanyi, Zoltan Fodor, Christian Hoelbling, Sandor D. Katz, Stefan Krieg, and Kalman K. Szabo, “Full result for the QCD equation of state with 2+1 flavors,” Phys. Lett. B 730, 99–104 (2014), arXiv:1309.5258 [hep-lat] .
  * Ding  _et al._ (2015) H. T. Ding, Swagato Mukherjee, H. Ohno, P. Petreczky, and H. P. Schadler, “Diagonal and off-diagonal quark number susceptibilities at high temperatures,” Phys. Rev. D 92, 074043 (2015), arXiv:1507.06637 [hep-lat] .
  * Bazavov  _et al._ (2017) A. Bazavov _et al._ , “The QCD Equation of State to $\mathcal{O}(\mu_{B}^{6})$ from Lattice QCD,” Phys. Rev. D 95, 054504 (2017), arXiv:1701.04325 [hep-lat] .
  * Monnai  _et al._ (2019) Akihiko Monnai, Björn Schenke, and Chun Shen, “Equation of state at finite densities for QCD matter in nuclear collisions,” Phys. Rev. C 100, 024907 (2019), arXiv:1902.05095 [nucl-th] .
  * Denicol  _et al._ (2012) G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, “Derivation of transient relativistic fluid dynamics from the Boltzmann equation,” Phys. Rev. D 85, 114047 (2012), [Erratum: Phys.Rev.D 91, 039902 (2015)], arXiv:1202.4551 [nucl-th] .
  * (59) S. Ryu, C. Shen, B. Schenke, in preparation.
  * Mäntysaari  _et al._ (2022) Heikki Mäntysaari, Björn Schenke, Chun Shen, and Wenbin Zhao, “Bayesian inference of the fluctuating proton shape,” Phys. Lett. B 833, 137348 (2022), arXiv:2202.01998 [hep-ph] .
  * Huovinen and Petersen (2012) Pasi Huovinen and Hannah Petersen, “Particlization in hybrid models,” Eur. Phys. J. A 48, 171 (2012), arXiv:1206.3371 [nucl-th] .
  * Shen  _et al._ (2016) Chun Shen, Zhi Qiu, Huichao Song, Jonah Bernhard, Steffen Bass, and Ulrich Heinz, “The iEBE-VISHNU code package for relativistic heavy-ion collisions,” Comput. Phys. Commun. 199, 61–85 (2016), arXiv:1409.8164 [nucl-th] .
  * (63) The iSS code package is an open-source particle sampler based on the Cooper-Frye freeze-out prescription. It converts fluid cells to particle samples. This work uses v1.0 of the iSS, which can be downloaded from https://github.com/chunshen1987/iSS/releases.
  * Ryu  _et al._ (2015) S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, “Importance of the Bulk Viscosity of QCD in Ultrarelativistic Heavy-Ion Collisions,” Phys. Rev. Lett. 115, 132301 (2015), arXiv:1502.01675 [nucl-th] .
  * Ryu  _et al._ (2018) Sangwook Ryu, Jean-François Paquet, Chun Shen, Gabriel Denicol, Björn Schenke, Sangyong Jeon, and Charles Gale, “Effects of bulk viscosity and hadronic rescattering in heavy ion collisions at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider,” Phys. Rev. C 97, 034910 (2018), arXiv:1704.04216 [nucl-th] .
  * Acharya  _et al._ (2022a) U. A. Acharya _et al._ (PHENIX), “Kinematic dependence of azimuthal anisotropies in $p+Au,$ $d+Au,$ and ${}^{3}He+Au$ at $\sqrt{s}_{NN}=200$ GeV,” Phys. Rev. C 105, 024901 (2022a), arXiv:2107.06634 [hep-ex] .
  * Acharya  _et al._ (2022b) U. A. Acharya _et al._ (PHENIX), “Measurements of second-harmonic Fourier coefficients from azimuthal anisotropies in $p+p,p$+Au $d$+Au, and 3He + Au collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV,” (2022b), arXiv:2203.09894 [nucl-ex] .
  * Bilandzic  _et al._ (2011) Ante Bilandzic, Raimond Snellings, and Sergei Voloshin, “Flow analysis with cumulants: Direct calculations,” Phys. Rev. C 83, 044913 (2011), arXiv:1010.0233 [nucl-ex] .
  * Bilandzic  _et al._ (2014) Ante Bilandzic, Christian Holm Christensen, Kristjan Gulbrandsen, Alexander Hansen, and You Zhou, “Generic framework for anisotropic flow analyses with multiparticle azimuthal correlations,” Phys. Rev. C 89, 064904 (2014), arXiv:1312.3572 [nucl-ex] .
  * Aidala  _et al._ (2017) C. Aidala _et al._ (PHENIX), “Measurements of azimuthal anisotropy and charged-particle multiplicity in $d+$Au collisions at $\sqrt{s_{{}_{NN}}}=$200, 62.4, 39, and 19.6 GeV,” Phys. Rev. C 96, 064905 (2017), arXiv:1708.06983 [nucl-ex] .
  * Adare  _et al._ (2019) A. Adare _et al._ (PHENIX), “Multiparticle azimuthal correlations for extracting event-by-event elliptic and triangular flow in Au$+$Au collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV,” Phys. Rev. C 99, 024903 (2019), arXiv:1804.10024 [nucl-ex] .
  * Pordes  _et al._ (2007) Ruth Pordes _et al._ , “The Open Science Grid,” J. Phys. Conf. Ser. 78, 012057 (2007).
  * Sfiligoi  _et al._ (2009) Igor Sfiligoi, Daniel C. Bradley, Burt Holzman, Parag Mhashilkar, Sanjay Padhi, and Frank Wurthwrin, “The pilot way to Grid resources using glideinWMS,” WRI World Congress 2, 428–432 (2009).
  * Schenke  _et al._ (2020b) Bjoern Schenke, Chun Shen, and Prithwish Tribedy, “Running the gamut of high energy nuclear collisions,” Phys. Rev. C 102, 044905 (2020b), arXiv:2005.14682 [nucl-th] .