###### Abstract

We present a comprehensive review of the theoretical and experimental progress
in the investigation of novel high-temperature quantum chromodynamics
phenomena in small systems at both the Relativistic Heavy Ion Collider and the
Large Hadron Collider. We highlight the challenges and opportunities
associated with studying small systems, by which we generally mean collision
systems that involve at least one light ion or even a photon projectile. We
discuss perspectives on possible future research directions to better
understand the underlying physics at work in the collisions of small systems.

##  Chapter 0 Progress and Challenges in Small Systems

###### Contents

  1. 0 Small Systems
    1. 1 Introduction
    2. 2 Multi-stage modeling for small system dynamics
      1. 1 Initial state models
      2. 2 The need for pre-equilibrium dynamics models
      3. 3 Relativistic hydrodynamics
    3. 3 Phenomenological studies in small systems
      1. 1 Collective flow
      2. 2 Multi-particle correlations
      3. 3 Partonic collectivity and strangeness enhancement in small systems
      4. 4 The $v_{n}$-$p_{T}$ correlations
      5. 5 Rapidity dependent anisotropic flow in small systems
      6. 6 Electromagnetic radiation in small systems
      7. 7 Hard probes and the $R_{pA}$ and $v_{2}(p_{T})$ puzzle in small systems
    4. 4 Further opportunities with small systems
      1. 1 Confronting ab initio nuclear structure with small system measurements
      2. 2 Accessing features of saturated QCD systems
      3. 3 Probing the boundaries of collective behavior
    5. 5 Instead of Conclusions
    6. 6 Acknowledgments

###  1 Introduction

For over two decades, ultrarelativistic heavy-ion collisions have produced and
systematically investigated the properties of the quark-gluon plasma (QGP), an
exotic phase of quantum chromodynamics (QCD) in which quarks and gluons are
not confined inside protons and neutrons and that existed microseconds after
the Big Bang. The current theoretical modeling [1, 2, 3] of these specks of
early universe matter produced in heavy-ion collisions is broadly consistent
with the statement that the quark-gluon plasma can flow [4] like a nearly
frictionless, strongly coupled quantum liquid over distance scales not much
larger than the proton radius. This makes the quark-gluon plasma formed in
colliders the hottest, smallest, densest, most perfect fluid known to
humanity.

Remarkable progress has been achieved over the years regarding the theoretical
description of heavy-ion collisions. A general framework consisting of initial
state modeling, followed by a pre-equilibrium phase and the subsequent
hydrodynamic evolution that is ultimately merged into hadronic transport,
forms the core of the so-called standard model of heavy-ion collisions [5, 6,
7]. This theoretical framework evolved from original qualitative success to
impressive predictive power [8, 9, 10, 11] and quantitative accuracy,
currently exemplified by the use of Bayesian inference techniques [12].

However, there are still many fundamental questions and potential pitfalls in
our understanding of heavy-ion collisions and the quark-gluon plasma, driven
by our general inability to perform ab initio real-time calculations in
quantum chromodynamics beyond weak coupling. In that context, the
investigation of the so-called small systems, such as proton-proton ($p$+$p$)
and proton-nucleus ($p$+A) collisions, has been particularly illuminating as
it requires pushing the boundaries of our understanding of the bulk collective
properties of the quark-gluon plasma.

Traditionally, relativistic hydrodynamics has been primarily employed in
studying large collision systems, such as those involving two heavy ions, e.g.
gold or lead. These systems exhibit an abundance of produced particles, a
long-lived quark-gluon plasma phase, and definite signatures of collective
behavior compatible with hydrodynamics. In the last decade, advances in
experimental capabilities have allowed us to explore experimental signatures
of collectivity in small collision systems. Those are characterized by more
significant pressure gradients and a shorter quark-gluon plasma phase compared
to that found in heavy-ion collisions. Nevertheless, despite their size, small
systems exhibit intriguing collective behavior, challenging the conventional
wisdom that collective effects are exclusive to large collision systems555We
note that some uncertainties appear already at the experimental level, as the
reduced number of particles produced in small systems introduces new
difficulties in experimental flow measurements [13, 14], which do not occur in
large systems.. However, one expects the uncertainties in theoretical modeling
to increase significantly as one moves down in system size (e.g., from Pb+Pb
to $p$+$p$ collisions).

For example, mirroring what is done in large systems, one may employ
relativistic hydrodynamics to investigate the collective properties of the
quark-gluon plasma encoded in the anisotropic flow of the final-state particle
distributions of small systems. Applying relativistic hydrodynamics to small
collision systems, however, presents several challenges. The reduced number of
particles produced in the events, and the shorter time for hydrodynamic
evolution, necessitate careful consideration of the initial conditions and the
freeze-out process. The large spatial gradients in the initial state can drive
the system significantly far from equilibrium, which can be challenging to
describe using fluid dynamics without violating fundamental physical
principles (such as relativistic causality). Due to their reduced system size,
in small systems, there is _a priori_ no clear separation between the scales
associated with the microscopic physics of QCD processes and the typical
scales related to gradients of conserved quantities. In this case, standard
arguments relying on the Knudsen number and inverse Reynolds number
expansions, widely used in the derivation of fluid dynamical models [15, 16],
are much less justifiable than in large heavy-ion collision systems (such as
central Pb+Pb, or Au+Au, collisions). Furthermore, stochastic thermal
fluctuations must become more relevant as the system decreases in size [17],
which introduces an inherent lack of determinism that has not yet been
included in state-of-the-art simulations of small systems.

On the other hand, applying relativistic hydrodynamics to small collision
systems represents a unique opportunity for progress in high-energy nuclear
physics [18, 19]. By studying flow observables and their hydrodynamic
interpretation in such extreme systems, one may uncover novel properties of
the quark-gluon plasma near and far from equilibrium, which can shed light on
the fundamental nature of strongly interacting matter. Continued theoretical
and experimental efforts in this field will deepen our understanding of the
quark-gluon plasma and its behavior across collision sizes.

In this chapter, we provide a concise overview of some recent developments and
challenges concerning the theoretical modeling of small systems. We intend to
be more illustrative than comprehensive, focusing on some theoretical results
and experimental measurements, which we believe give a reasonably good idea of
the progress achieved through the last decade and the remaining fundamental
questions. In this regard, we admit that the presentation is not entirely
balanced, with some results receiving more attention than others.

This paper is structured as follows. In Section 2, we introduce the models
employed to investigate the small quark-gluon plasma droplets formed in small
collision systems. This includes a discussion of the initial conditions, the
pre-equilibrium dynamics, and the hydrodynamic models. Section 3 presents the
model-to-data comparison, encompassing various observables. Section 4 outlines
additional theoretical advances that offer further insights into the physics
of small systems. Our final comments are presented in Section 5. _Notation_ :
We use a mostly minus metric signature $(+---)$ and natural units,
$\hbar=c=k_{B}=1$.

###  2 Multi-stage modeling for small system dynamics

####  1 Initial state models

Hydrodynamic simulations rely on models for the initial conditions. In
principle, these should provide the entire energy-momentum tensor and
conserved currents over the entire space at the initial time. However, many
models provide only the spatial energy density distribution, which for heavy
ion collisions with a long lifetime may be a good approximation but, for small
systems, is likely insufficient, as initial velocity and off-equilibrium
corrections can significantly affect the produced particle distributions.

Several studies have employed Monte Carlo Glauber type models [20, 21] with
nucleon degrees of freedom without substructure to initialize hydrodynamic
simulations of small systems [22, 23, 24, 25, 25, 26, 27, 28, 29, 30, 31, 32,
33, 34, 35, 36, 37]. These calculations generally reproduce azimuthal
anisotropies measured in experiments at Relativistic Heavy Ion Collider (RHIC)
and the Large Hadron Collider (LHC).

The success of these simulations depends on the choice of energy or entropy
deposition, which can be proportional to the sum of the thickness functions
666The thickness function is the integral of the nuclear density function over
the direction of the beamline. of the projectile and target or other
combinations, such as the product of thickness functions. Only the former
option yields good agreement with experimental data for nucleon-based initial
state models. For example, in $p$+A collisions, the proton projectile is
perfectly round when assuming no substructure, and fluctuations can only
originate from the fluctuating positions of nucleons in the target. In this
case, only the sum of thickness functions produces enough fluctuations to
explain the data [5, 38], while the product leads to overly smooth energy
density distributions [39].

Nevertheless, there are strong arguments for why the sum should be disfavored.
First, all Bayesian analyses using the T$\mathrel{\raisebox{-2.1pt}{{\rm
R}}}$ENTo model [38], which parametrizes the dependence of the initial entropy
or energy density distribution on the thickness functions, prefer an entropy
density proportional to the square root of the product of the two thickness
functions [40, 34, 41, 42, 7], rather than the sum, which is also included in
T$\mathrel{\raisebox{-2.1pt}{{\rm R}}}$ENTo’s parameter space. Further, the
effective theories and models that allow for calculating the energy deposition
without assuming the dependence on the thickness function all find
proportionality to (some power of) their product rather than the sum. This
includes the weak coupling calculations based on the color glass condensate
(CGC) effective theory [43, 44, 45, 46, 47] and strong coupling holographic
calculations [48, 2].

By adhering to an energy or entropy deposition description that utilizes the
product of thickness functions or a power thereof, an accurate representation
of anisotropic flow data in small systems can only be achieved by
incorporating subnucleon fluctuations. This has been explicitly demonstrated
in the CGC-based IP-Glasma [39] and T$\mathrel{\raisebox{-2.1pt}{{\rm
R}}}$ENTo [49] models.

The fluctuating spatial structure of a nucleon has not been computed from
first principles. To make progress, one can parametrize the substructure as a
collection of hot spots with certain sizes and then constrain the parameters
through Bayesian analyses of heavy-ion and small system collisions.[49, 7]
Although the precise number of hot spots could not be well constrained, a
finite number of hot spots was favored over a smooth nucleon. The width of the
subnucleonic hot spots was constrained to approximately $0.4\,\rm{fm}$,
similar to the typical width used for the _nucleon_ width in the IP-Glasma
model.

Another approach involves utilizing input from electron-proton ($e$+$p$)
scattering experiments. It has been demonstrated that the cross-sections for
incoherent diffractive vector meson production are sensitive to the proton
substructure and its fluctuations [50, 51, 52, 53]. By constraining the
parameters of a hot spot model using data from HERA [54, 55, 56, 57, 58], and
applying this model in proton-nucleus ($p$+A) collisions, successful
reproduction of flow observables has been achieved [59, 60, 61]. The most
sophisticated extraction of nucleon substructure parameters from HERA data
used Bayesian techniques [52]. While some parameters could be well
constrained, the number of hot spots, for example, was not well determined. A
global analysis using both $e$+$p$ and hadronic collisions could potentially
improve this situation.

As mentioned above, the initial state description should provide more than the
spatial distributions of energy density. For example, the spatial distribution
of the entire initial energy-momentum tensor can be computed within the IP-
Glasma model with the caveat that the system under consideration is purely
gluonic. In the classical approximation, the coupling is fixed, and the system
is conformally invariant. This means that the equation of state from this
calculation differs from the QCD equation of state [62, 63]. Furthermore,
conformal invariance implies that this system’s bulk viscosity transport
coefficient vanishes. Such differences introduce challenges when matching to
hydrodynamics[64, 65], for which one usually employs an equation of state
extracted from lattice QCD. For more on the early time dynamics that
transition the initial state to hydrodynamics, we refer the reader to the next
section.

The IP-Glasma initial energy-momentum tensor includes (shear) viscous
corrections and an initial flow profile. Interestingly, the initial transverse
flow velocity distribution is already anisotropic, a consequence of
correlations among the color fields in the colliding nuclei [66, 67]. Assuming
boost invariance of these correlations and performing a 2+1D (IP-Glasma-like)
Yang-Mills evolution, such correlations were shown to survive even after
hydrodynamic evolution for low multiplicity events [68]. However, performing
small $x$-evolution to obtain a rapidity-dependent initial state shows that
the correlations leading to such initial state momentum anisotropies decay
rapidly with the rapidity difference [69]. One may also expect that full 3+1D
Yang-Mills simulations, which include rapid growth of small fluctuations [70,
71, 72] and lead to at least partial pressure isotropization, will further
reduce the initial state anisotropies.

The extension to full 3D initial conditions has usually been done by
introducing a parametric longitudinal envelope folded with the transverse
energy density profile [73, 74, 25, 75, 76, 77, 78, 79]. In this family of
initial state models, one can further impose energy-momentum conservation when
mapping the initial-state density profiles to hydrodynamic fields [80, 81].
Different implementations result in different amounts of longitudinal
fluctuations in the initial energy density profiles, which can be constrained
by the longitudinal decorrelation measurements.

In Section 3 we will discuss in more detail the need for simulations that
explicitly consider all three spatial dimensions, i.e., that do not assume
boost invariance. That means that the initial state model needs to provide the
energy-momentum tensor as a function of both transverse space and space-time
rapidity. Several models that achieve this have been employed in the
description of small systems.[82, 83, 76, 84, 85].

####  2 The need for pre-equilibrium dynamics models

Immediately after the collisions between the nuclei, the matter is expected to
find itself in a gluon-dominated phase [44]. In principle, this idea can be
used to determine the spatial dependence of all the components of the system’s
energy-momentum tensor [39] and other conserved currents [86, 84, 87, 88] (in
this section, for simplicity, we focus on the energy-momentum tensor).
However, the system is usually still in a far-from-equilibrium state at such
early times, which implies that matching the initial state model directly to
hydrodynamics cannot be done without drastic, unjustified assumptions. This
can be better understood as follows.

Assume that some initial state model, say IP-Glasma, gives an initial
$T^{\mu\nu}$ at some early time $\tau_{0}$, which we intend to use to define
the initial conditions for hydrodynamic evolution directly. One can then use
Landau matching [17] to define a flow 4-velocity $u^{\mu}$ and energy density
$\varepsilon$ at $\tau_{0}$ as an eigenvector and its corresponding eigenvalue
of $T^{\mu\nu}$, i.e., $T^{\mu\nu}u_{\nu}=\varepsilon u^{\mu}$. We note that
one usually assumes that this procedure gives a flow velocity that obeys
$u_{\mu}u^{\mu}=1$ and $\varepsilon$ is non-negative. We note that these
properties are not guaranteed to hold in general[89, 90, 91]. Putting such
issues aside, one can unambiguously decompose the initial energy-momentum
tensor as

$T^{\mu\nu}=\varepsilon
u^{\mu}u^{\nu}-[p(\varepsilon)+\Pi]\Delta^{\mu\nu}+\pi^{\mu\nu}\,,$ (1)

where $p(\varepsilon)$ is the thermodynamic pressure related to the energy
density through the equation of state, $\Pi$ is the bulk pressure,
$\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$, and $\pi^{\mu\nu}$ is the shear
stress tensor defined by
$\pi^{\mu\nu}=\Delta^{\mu\nu}_{\alpha\beta}T^{\alpha\beta}$, where
$\Delta^{\mu\nu}_{\alpha\beta}=\frac{1}{2}\left(\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}+\Delta^{\mu}_{\beta}\Delta^{\nu}_{\alpha}\right)-\frac{1}{3}\Delta^{\mu\nu}\Delta_{\alpha\beta}$.

If hydrodynamics is a good approximation for the initial $T^{\mu\nu}$, then
one would expect that the deviations from equilibrium, quantified by the
inverse Reynolds numbers $\sim|\Pi|/p,\sqrt{|\pi^{\mu\nu}\pi_{\mu\nu}|}/p$,
are locally _small_ at the initial time $\tau_{0}$. It turns out that this is
generally not the case even in large heavy-ion collision systems [92, 93, 60,
94]. Therefore, one generally expects that, at sufficiently early times,
direct matching of the energy-momentum tensor from initial state models to
hydrodynamics is a very poor approximation. This should be even more pressing
in small systems, because of their reduced size compared to large heavy-ion
collisions.

In recent years, pre-equilibrium models have been used to bridge the gap
between the initial state energy-momentum tensor and the energy-momentum
tensor at the start of the hydrodynamics phase. Such models are used to
describe the evolution of the system from $\tau_{0}$ up to
$\tau_{\text{hydro}}$, the time at which a hydrodynamic description is valid.

The simplest model one can consider is free streaming. In this case, one
assumes that on-shell non-interacting massless partons emerge isotropically
from an initial hard scattering at the initial time $\tau_{0}$ [95, 96].
Assuming longitudinal boost invariance, the energy-momentum tensor of such a
system can be obtained by integrating the initial number density of partons in
the transverse plane, effectively smoothing out the energy density of the
system. This free streaming dynamics is suddenly interrupted at
$\tau_{\text{fs}}>\tau_{0}$, where the system is suddenly assumed to attain
local thermal equilibrium. Clearly, this sudden transition from a zero to
finite (strong) coupling is unphysical. Therefore, free streaming pre-
equilibrium dynamics is necessarily a drastic assumption that, though
practical, is theoretically unjustified.

Progress toward a more realistic scenario was explored in [97, 98] using the
KøMPøST framework. There, the system’s evolution from $\tau_{0}$ to
$\tau_{\text{hydro}}$ is described by an effective kinetic theory model of
weakly coupled QCD [99]. In this approach, the system’s dynamics at early
times after the collision are assumed to be determined by classical Yang-Mills
equations. The system expands, its energy density drops, and at a time
$\tau_{\text{EKT}}$ the subsequent evolution of the system is modeled by
effective kinetic processes that drive the plasma towards a state where
inverse Reynolds numbers have sufficiently decreased such that relativistic
viscous fluid dynamics models may become a good approximation.

While this represents significant progress, current approximations in such
pre-equilibrium models can still be improved. For example, one can take the
trace of Eq. (1) to find that the total isotropic pressure is
$p(\varepsilon)+\Pi=\frac{\varepsilon-T^{\mu}_{\mu}}{3}$. One can see that for
conformal pre-equilibrium models, such as the original KøMPøST framework where
$T^{\mu}_{\mu}=0$, a simple calculation gives
$p(\varepsilon)+\Pi=\frac{\varepsilon}{3}$. This implies that an artificial
bulk viscous pressure component must be added at the beginning of the
hydrodynamic evolution because the conformal pressure $\epsilon/3$ differs
from the QCD pressure. Uncertainties associated with this issue have been
investigated [60, 64, 65, 100, 101]. Secondly, the dilute cold corona of the
collision requires a non-perturbative description of QCD confinement and
parton-hadron interactions [102, 103].

Pre-equilibrium models are needed in self-consistent hydrodynamic simulations
of large and small collision systems. Their inclusion in the standard modeling
of heavy-ion collisions emphasizes the need to investigate how the far-from-
equilibrium initial QCD fields decohere towards a state amenable to
hydrodynamic description. In fact, investigations into the physics of the pre-
equilibrium phase must consider the inherent non-Abelian physics of QCD in
out-of-equilibrium conditions, which is bound to lead to several new
theoretical developments. For example, the pre-equilibrium phase should be
crucial to understanding the onset of hydrodynamics and the presence of
hydrodynamic attractors [104, 105, 106], which have been proposed as a way to
extend the domain of applicability of hydrodynamics to the far-from-
equilibrium regime [107].

####  3 Relativistic hydrodynamics

Fluid dynamic behavior emerges in vastly different systems in nature.
Applications can be found across many scales, ranging from cosmology to the
fluid dynamic phenomena observed in everyday life. The ubiquitousness of
hydrodynamics stems from the presence of conservation laws and the assumption
of a significant separation between different characteristic length scales in
the system. In general, this hierarchy is quantified by the Knudsen number
$\mathrm{Kn}\sim\ell/L$, which broadly denotes the ratio between the relevant
microscopic scale $\ell$ and a macroscopic scale $L$ associated with the
variations of conserved quantities (such as energy, momentum, and baryon
number). Fluid dynamics is expected to be a good approximation to describe the
long-time, large wavelength physics of many-body systems in situations where
$\mathrm{Kn}\ll 1$ [17].

Relativistic hydrodynamics plays a crucial role in understanding the behavior
of matter and energy under extreme conditions [2]. The equations of motion of
relativistic hydrodynamics encode the conservation of energy, momentum, and
conserved charges:

$\partial_{\mu}T^{\mu\nu}=0,\qquad\partial_{\mu}J^{\mu}=0.$ (2)

Here, $T^{\mu\nu}$ is the system’s energy-momentum tensor, and the $J^{\mu}$
represents the 4-current of conserved charges (such as baryon number in the
case of heavy-ion collisions). In its simplest incarnation, known as ideal
fluid dynamics, dissipative processes are neglected, and these equations are
closed using the matter’s equation of state, e.g., $p(\varepsilon,n_{B})$,
where $n_{B}$ is the baryon density [108].

In many aspects, one may argue that the quark-gluon plasma produced in
ultrarelativistic collisions lies at the edge of the applicability of fluid
dynamics. In small systems, this is even more pressing since “macroscopic
scales” are, at the most, comparable to the proton radius (of the order of $1$
fm), and the relevant microscopic scales can be estimated in terms of the
inverse temperature $\ell\sim 1/T$. This exercise shows that for the
temperature range expected in such systems, $T\sim 150-450$ MeV, the Knudsen
number $\mathrm{Kn}$ is not necessarily small [92, 93, 32, 60, 3]. This has
important consequences for the hydrodynamic modeling of the quark-gluon
plasma, as discussed below.

In high-energy nuclear collisions, relativistic hydrodynamics has proven to be
an indispensable tool for investigating the intricate dynamics of the quark-
gluon plasma. Because the underlying QCD system is relativistic, and
dissipative processes may not be neglected (as $\mathrm{Kn}$ is not generally
small), a _causal_ formulation of relativistic _viscous_ hydrodynamics is
essential for its phenomenological applications in heavy-ion collisions. Given
that the standard relativistic generalization of Navier-Stokes equations,
derived by Eckart [109] and Landau and Lifshitz [17], lead to acausal
evolution [110] where small disturbances around the equilibrium state grow
exponentially [111], the development of other frameworks where dissipation is
not fundamentally at odds with relativistic causality became necessary. It is
now understood that the issues found in Eckart, and Landau and Lifshitz’s
theories follow from very general arguments, as proved recently[112, 113].
Basically, the main insight is that causality is necessary for the stability
of the equilibrium state in relativistic fluids. The previous formulations
from Eckart, and Landau and Lifshitz give rise to non-hyperbolic partial
differential equations (PDE)s [110] that necessarily violate relativistic
causality, creating unphysical instabilities of the equilibrium state.

The issues found in previous formulations can be fixed using the new approach
put forward by Bemfica, Disconzi, Noronha, and Kovtun (BDNK) [114, 115, 116,
112, 117]. Following standard effective theory reasoning [118], in BDNK non-
equilibrium corrections are taken into account in constitutive relations using
the most general expressions compatible with the symmetries involving
combinations of first-order spacetime derivatives of the standard hydrodynamic
variables (which vanish for the equilibrium state). This procedure creates
several new terms in the equations of motion, which were not considered before
and are crucial to restoring causality and stability. Though at this point,
this new formalism has not been applied to realistic quark-gluon plasma
simulations, nontrivial numerical solutions have already been developed in
other contexts[119, 120, 121, 122, 123].

Current hydrodynamic simulations of the quark-gluon plasma employ equations of
motion that can lead to causal and stable evolution within the framework
originally devised by Israel and Stewart [124]. This approach differs from the
BDNK reasoning mentioned above because in Israel-Stewart theories [124, 125,
15] the dissipative fluxes (such as the shear-stress tensor, $\pi_{\mu\nu}$,
and the bulk scalar, $\Pi$) are not defined via constitutive relations, i.e.,
they are not expressed only using the hydrodynamic variables and their
derivatives. Rather, the dissipative fluxes evolve according to additional
equations of motion of relaxation type (derived either from entropic
considerations [124], kinetic theory [15], effective theory arguments [125],
or anisotropic hydrodynamics [126]), which introduce a relaxation time that
characterizes how the dissipative fluxes relax towards their Navier-Stokes
values. The relaxation time effectively restores causality as, for example,
flow gradients are not automatically converted into acceleration as in the
relativistic Navier-Stokes equations. In-depth discussions of the different
theories of relativistic viscous fluid dynamics can be found in the
literature[16].

All state-of-the-art hydrodynamic simulations of the quark-gluon plasma formed
in heavy-ion collisions have been performed using Israel-Stewart-like
approaches. The equations of motion are now solved by many different groups,
using either Eulerian [127, 128, 129, 25, 32, 33, 59, 61] or Lagrangian
numerical algorithms [130, 131, 93], which pass standard tests of accuracy
[132]. However, despite this significant progress in solving these equations
of motion in the context of heavy-ion collisions, there are still
uncertainties in the modeling coming from the hydrodynamic description. This
is especially relevant in the case of simulations of small collision systems,
as discussed below.

##### Large Knudsen numbers and inverse Reynold’s numbers in small systems

Relativistic hydrodynamics describes the macroscopic dynamics of many-body
systems. The system’s entropy is packed in a small volume in space for high
multiplicity events of ultrarelativistic proton-nucleus and proton-proton
collisions. Despite small particle multiplicity, the local energy density can
be large enough such that, at least in principle, the system may behave as a
fluid.

In small collision systems, the limited size and short lifetime of the created
QGP pose unique challenges in the hydrodynamic description of the system’s
dynamics. One significant aspect relevant to the overall applicability of
hydrodynamics in small systems is that, in such systems, Knudsen numbers and
inverse Reynolds numbers can become large.

Let us first make a few general comments about the Knudsen number in small
systems. The initial transverse volume in small systems is about 50 times
smaller than in central Pb+Pb collisions. The reduction of the macroscopic
size significantly increases the Knudsen number at early times. This indicates
that the underlying collision rate is small, meaning that the collision time
is larger than the characteristic timescale of the hydrodynamic expansion. As
a result, the assumptions of local thermal equilibrium and continuous fluid
behavior, crucial to our standard understanding of the domain of applicability
of hydrodynamics, are progressively challenged.

The inverse Reynolds number ($\mathrm{Re}^{-1}$) serves as a complementary
measure, characterizing the importance of viscous effects in the system. It
describes the ratio of viscous forces to inertial forces, and it can be
estimated, for example, as the ratio of the viscous pressure contributions to
the fluid cell’s local equilibrium enthalpy. In practice, one may define the
inverse Reynolds number associated with shear stress as [15, 16]

$\mathrm{Re}^{-1}_{\pi}=\frac{\sqrt{|\pi_{\mu\nu}\pi^{\mu\nu}|}}{\varepsilon+p}\,.$
(3)

Near equilibrium, this quantity becomes small, as viscous stresses are then a
small contribution to the pressure/energy density.

We note that while $\mathrm{Re}^{-1}$ and $\mathrm{Kn}$ are certainly related,
in general they are not the same. To illustrate this point, consider the case
of shear-viscous effects. In that case, one may estimate the shear Knudsen
number as

$\mathrm{Kn}_{\sigma}\sim\frac{\sqrt{|\sigma_{\mu\nu}\sigma^{\mu\nu}|}}{T},$
(4)

where the shear tensor is
$\sigma^{\mu\nu}=\Delta^{\mu\nu\alpha\beta}\nabla_{\alpha}u_{\beta}$. The
$1/T$ factor denotes the microscopic scale777For a dilute gas, one may use the
mean free path as the microscopic scale $\ell\sim 1/n\sigma_{0}$, where $n$ is
the density and $\sigma_{0}$ the cross-section. For a strongly coupled system,
$\ell\sim 1/T$ is typically used in such estimates., while the other factor
describes the scale associated with gradients of the flow velocity. Now one
may compare Eq. (3) with Eq. (4). In the Navier-Stokes regime, where
$\pi_{\mu\nu}=2\eta\sigma_{\mu\nu}$, these quantities are completely
equivalent. However, in Israel-Stewart theories $\pi_{\mu\nu}$ evolves
according to a differential equation

$\tau_{\pi}\Delta^{\mu\nu}_{\alpha\beta}u^{\lambda}\partial_{\lambda}\pi^{\alpha\beta}+\pi^{\mu\nu}=2\eta\sigma^{\mu\nu}+\ldots\,,$
(5)

where $\tau_{\pi}$ is the shear relaxation time and $\ldots$ denote other
terms usually employed in simulations [16]. Thus, depending on the initial
conditions, the state of the system, and the values of the transport
coefficients, in general, $\pi^{\mu\nu}$ is not well approximated by
$2\eta\sigma^{\mu\nu}$. In that _transient_ case, $\mathrm{Re}^{-1}$ and
$\mathrm{Kn}$ should be treated as separate parameters that may be used in a
series expansion [15]. In small collision systems, the large pressure
gradients lead to a fast expansion rate, which drives the system significantly
out of equilibrium. Consequently, the viscous effects can dominate the
dynamics, leading to a large inverse Reynolds number, which is not reflected
by the corresponding $\mathrm{Kn}$.

The presence of large Knudsen and inverse Reynolds numbers in small systems
challenges the applicability of conventional relativistic viscous hydrodynamic
models that assume perturbative viscous corrections around local thermal
equilibrium. To address these challenges, efforts have been made to develop
hybrid approaches that combine relativistic hydrodynamics with effective
kinetic theory frameworks (which were discussed in Sec. 2). These hybrid
models aim to capture the non-equilibrium dynamics and viscous corrections
that become increasingly significant in small systems.

##### Possible cavitation under rapid expansion

The rapid expansion observed in small collision systems suggests that
cavitation may occur and play an important role in the hydrodynamic
description. Cavitation, commonly observed in fluid dynamics [133], refers to
the formation and subsequent collapse of vapor or gas-filled voids within a
liquid medium. It arises when the local pressure drops below the vapor
pressure, leading to the nucleation and growth of bubbles.

In heavy-ion collisions, the hydrodynamic evolution is characterized by an
intense and rapid expansion of the QGP, leading to a drop in pressure and
temperature. Under certain conditions, this rapid expansion may give rise to a
local total pressure that falls close to zero, potentially leading to
cavitation within the QGP [134, 135, 136, 137, 138, 60]. The idea behind it is
simple: the system’s total pressure has equilibrium and non-equilibrium
contributions, i.e., $p_{\mathrm{total}}=p+\Pi$. In Navier-Stokes theory
$\Pi=-\zeta\partial_{\mu}u^{\mu}$ [16], where $\zeta$ is the bulk viscosity.
For an expanding system $\partial_{\mu}u^{\mu}>0$, hence $\Pi<0$, and, if the
expansion is sufficiently large, there could be regions where the total
pressure vanishes. In practice, such arguments are qualitative at best, given
that in the transient hydrodynamic regime, $\Pi$ is not well approximated by
Navier-Stokes values, and consequently, the bulk inverse Reynolds number does
not equal the associated Knudsen number. In any case, if cavitation can occur
in heavy-ion collisions, it is reasonable to expect that it is more likely to
take place in small systems because of their large expansion rate/out-of-
equilibrium corrections. It is interesting to note that having zero local
pressure does not necessarily lead to fundamental problems. For example,
causality can still be preserved even if $p+\Pi<0$ as long as
$\varepsilon+p+\Pi>0$ [139].

Since the transition from the QGP phase to the hadron gas phase is a smooth
cross-over, cavitation in small collision systems may or may not have
significant consequences for interpreting experimental observables.
Nevertheless, these low-pressure voids can introduce non-equilibrium effects
and modify the system’s evolution. The expansion dynamics of the QGP, which is
typically described by viscous hydrodynamics, may be affected by the presence
of these low-pressure voids.

The possibility of cavitation under rapid expansion in small systems presents
an intriguing avenue for research in relativistic hydrodynamics. Understanding
and quantifying the possible cavitation effects in small systems is
challenging. The dynamics of bubble formation and collapse in the QGP require
a detailed understanding of the equation of state, the thermodynamic
properties of the medium, and the transport coefficients. Furthermore, the
evolution of the cavitation bubbles must be considered in the context of the
overall expansion and freeze-out processes. Incorporating these complexities
into hydrodynamic models and analyzing their implications on experimental
observables represent active research areas.

##### Causality violations when simulating small systems

Causality is, in a broad sense, the statement that no information in physical
systems can be propagated faster than the speed of light. As mentioned before,
incorporating dissipative effects in a way that is consistent with
relativistic causality is still a matter of intense investigation [140, 141,
142]. Causality is a necessary requirement for thermodynamic stability in
relativistic many-body systems [113]. For small disturbances around
equilibrium, it is possible to obtain conditions involving the equation of
state and the transport coefficients of the system under which causality and
stability can be established. This is the case in all the theories employed to
investigate relativistic viscous fluids [16], such as BDNK [114, 115, 116,
112, 117] and Israel-Stewart theories [143, 144, 145, 146, 147, 148]. Once
nonlinear corrections are included, the equations of motion of hydrodynamics
become a nonlinear set of PDEs, and establishing causality (or hyperbolicity)
is a highly challenging task. In BDNK theories, where viscous corrections do
not require the introduction of additional fields besides those in
equilibrium, causality in the full nonlinear regime has been rigorously
established [114, 115, 116, 112, 117, 149]. For Israel-Stewart theory,
however, the analysis is much more complicated because the causality
conditions now explicitly include the dissipative fluxes, in addition to
transport coefficients and the equation of state. This implies that the
constraints needed to ensure causality in the linear regime around equilibrium
are not sufficiently powerful to guarantee causality in the nonlinear regime.
In that case, only a few general results, valid for general equations of
state, in 3+1 dimensions, etc, are known [139, 150].

When performing numerical simulations of heavy-ion collisions, it is crucial
to ensure that the evolution preserves causality. Violating causality
conditions can lead to significant consequences that impact the reliability
and interpretation of the results [151, 152]. This is particularly relevant
for small collision systems [153]. Causality violations can occur due to
several factors. One common issue is the choice of relaxation times for the
shear and bulk viscous tensors. Suppose they are too close to the causality
bounds derived in the linear perturbative region. In that case, there is
little room for the system to stay out of equilibrium, according to the
causality conditions recently derived for the Israel-Stewart hydrodynamic
theory in the full non-equilibrium regime [150]. In some cases, serious
violations of causality conditions can generate numerically unstable modes
that grow exponentially. However, in other cases, the signs of causality
violations are hard to trace, and the numerical simulations remain stable.

The violation of causality conditions in numerical simulations of small
systems demands careful consideration, because it can undermine the
reliability and accuracy of the obtained results. Unphysical artifacts
introduced by causality violations can distort the interpretation of the
system’s dynamics and compromise the comparison with experimental data [151,
152, 153]. This hampers the extraction of meaningful insights about the
underlying physics and, in particular, the transport properties of the QGP in
small collision systems. Moreover, violating causality conditions can hinder
developing and improving theoretical models and numerical algorithms.
Identifying and addressing the sources of unphysical behavior becomes
challenging when causality violations are present. Resolving these issues
requires careful investigations within the complex numerical simulations.

To mitigate the consequences of causality violations, ongoing research focuses
on developing more robust theories and implementing numerical algorithms that
explicitly account for the nonlinear causality conditions [139, 150] in small
system simulations [151, 152, 153].

##### Sensitivity to second-order transport coefficients

In large collision systems, the effects of viscosity on the collective
behavior of the QGP are well understood [1, 2]. They can be characterized by
the first-order transport coefficients, namely the shear viscosity ($\eta$)
and the bulk viscosity ($\zeta$). Recent studies have revealed that small
collision systems exhibit a heightened sensitivity to both the first-order and
second-order transport coefficients, such as the viscous relaxation times
$\tau_{\pi}$ and $\tau_{\Pi}$ [60, 154, 152].

In large collision systems, the collective behavior is dominated by the long-
lived QGP phase, reducing the sensitivity to the specific values of the
second-order transport coefficients. However, small collision systems have a
reduced QGP lifetime, leading to shorter hydrodynamic evolution timescales
close to the viscous relaxation times.

The sensitivity of small collision systems to the transport coefficients
manifests in various ways. First, the magnitude of the viscous effects becomes
more pronounced in small systems, leading to enhanced sensitivity to the
values of shear and bulk viscosity and their temperature dependence [60].
Second, the hydrodynamic evolution in small systems is more susceptible to
non-equilibrium and transient effects [155, 156]. The shorter duration of the
QGP phase reduces the time available for the system to relax towards local
thermal equilibrium, amplifying the influence of dissipative processes. In
this context, accurately determining the transport coefficients becomes
crucial, as their values directly affect the system’s ability to reach thermal
equilibrium and influence the subsequent particle production and flow
patterns.

The sensitivity of small collision systems to the second-order transport
coefficients poses challenges to theoretical modeling and interpretation of
experimental data. Precise knowledge of the transport coefficients is required
to accurately describe the dynamics and extract meaningful information from
experimental observables. This necessitates a deeper understanding of the
underlying microscopic processes and improved constraints on the transport
coefficients from experimental measurements.

To address these challenges, efforts are underway to systematically study the
second-order transport coefficients and their temperature dependence within
different theoretical frameworks [157, 158, 159, 160]. These studies involve
comparisons between theoretical calculations, hydrodynamic simulations, and
experimental data [7, 161, 152, 162].

Figure 1: Effective shear viscosity versus gradient strength for hydrodynamic
attractors in different theories (rBRSSS, kinetic theory and AdS/CFT). For
large gradient strength, $\frac{\eta_{B}}{s}\rightarrow 0$. This figure is
taken from published work. [163].

##### Resummed transport coefficients and attractors

Hydrodynamics is traditionally based on the assumption that the system is near
local equilibrium. A trivial rearrangement of this statement in terms of the
Knudsen number is that hydrodynamics is under control when $\mathrm{Kn}\ll 1$.
However, there is growing evidence supporting the capability of hydrodynamics
to describe small colliding systems that are not close to local equilibrium
and in situations where $\mathrm{Kn}\sim 1$. Recently, it has been argued that
part of this success may be attributed to the existence of hydrodynamic
attractors [104, 163, 105]. According to this view [107], the solution of the
dynamical equations converges quickly to this attractor through non-
hydrodynamic mode decay. Despite the fact that the nonlinear hydrodynamic
gradient series can diverge [164, 165], the series may be resummed to give
rise to the hydrodynamic attractor, which may be approximated by second-order
hydrodynamic solutions (at least, in very simplified cases).

Resummation of the gradient series may also provide a way to define transport
coefficients far from equilibrium [163, 166]. Figure 1 shows this effective
shear viscosity $\eta_{B}/s$ versus gradient strength [163]. For small
gradients (system near equilibrium), the effective viscosity reduces to the
viscosity computed using linear response theory. For large gradient strength
(system far from equilibrium), the effective viscosity $\eta_{B}/s$ tends to
zero. Thus, in systems out-of-equilibrium, it is possible that higher-order
viscous corrections effectively reduce the value of the viscosity entering the
viscous hydrodynamic equations [163, 167].

The idea that the far-from-equilibrium properties of the system may be taken
into account in new resummed transport coefficients can be potentially
relevant for phenomenology. For example, a unified hydrodynamic description of
the quark-gluon plasma generated in ultrarelativistic $p$+Pb and Pb+Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV has been proposed.[49, 61] The
authors employ a Bayesian parameter estimation model calibration framework to
conduct a global analysis of bulk observables in both $p$+Pb and Pb+Pb
collisions. Nevertheless, it is crucial to acknowledge that, as previously
discussed, smaller systems exhibit larger gradient strengths, leading to a
comparatively smaller effective viscosity than that observed in heavy-ion
collisions. A potential avenue for future research involves a more
comprehensive approach that integrates small and large systems. This framework
could involve determining the viscosity in small systems in the context of
gradient series resummation and resummed transport coefficients. The
anticipated outcome is a smaller effective viscosity and/or viscous effects in
small systems than it would be naively expected from standard bare
$\mathrm{Kn}$ arguments. This innovative approach holds promise for refining
our understanding of the nuanced dynamics within different collision systems.

##### Impacts from stochastic noise on the hydrodynamic description in small
systems

Dissipative systems near equilibrium must display thermal fluctuations – this
is the physics behind the celebrated fluctuation-dissipation theorem [168,
169]. Therefore, stochastic noise arises from the inherent fluctuations
present in the microscopic dissipative dynamics of the heavy-ion collision
system. In small systems, the reduced size and short lifetime are expected to
amplify the relative importance of these stochastic effects. The fluctuations
effectively contribute as a stochastic source term to the hydrodynamic
equations, converting the system to a nonlinear set of stochastic PDEs, which
is very challenging to directly evolve numerically [170, 171].

The impacts of stochastic noise on the hydrodynamic description in small
systems are expected to be twofold. First, the fluctuations can introduce an
additional source of event-by-event variations in the system’s evolution (for
the same initial state), resulting in broadened distributions for many
observables. In particular, stochastic noise can lead to additional event-by-
event fluctuations in the azimuthal anisotropy coefficients ($v_{n}$) and
other flow-related observables, which are not currently considered in state-
of-the-art simulations.

Second, the need for the inclusion of stochastic noise in small systems
provides another source of uncertainty in our current hydrodynamic models that
assume smooth and deterministic evolution. The fluctuations can induce
effective non-equilibrium behavior, which is particularly relevant in small
systems, where the QGP lifetime is shorter and the hydrodynamic evolution
timescale is comparable to or shorter than the characteristic relaxation
timescales.

Understanding and quantifying the impacts of stochastic noise on the
hydrodynamic description of relativistic systems is a very active area of
research. Several approaches and applications have been pursued in recent
years, which include work on hydro-kinetics [172, 173, 174, 175, 176],
fluctuations near a critical point [177, 178, 179, 180, 181, 182, 183, 184,
185, 186, 187], general studies of fluctuations in the context of heavy-ion
collisions [188, 189, 190, 171, 191, 192, 193], and more [194, 195, 196, 197,
118, 198, 170, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208]. Many of the
theoretical developments have focused on the Schwinger-Keldysh action on the
closed time path [209, 210, 211, 212, 213] where actions are constructed as an
effective field theory from the underlying quantum mechanical system using the
dynamical Kubo-Martin-Schwinger symmetry [169, 214, 215, 216, 217, 218, 219,
220, 221, 222, 223, 224]. New works have explicitly considered the recent
developments concerning causality and thermodynamic stability in relativistic
fluids [225, 226].

Experimental measurements as functions of the rapidity gap [227] can provide
valuable insights into the impacts of stochastic noise on small systems.
Measurements of higher-order cumulants and flow fluctuations, as well as the
study of event-by-event correlations in small systems with stochastic
fluctuations, are expected to be crucial for unraveling the underlying physics
and extracting more accurate information about the transport properties of the
QGP [228].

##### Particlization and freeze-out in small systems

To compute experimental observables, the fluid needs to be converted into
particles. This particlization stage is usually performed by the Cooper-Frye
procedure [229], which imposes continuity requirements for the energy-momentum
current across the conversion hyper-surface under the grand-canonical
ensemble. Each fluid cell on the conversion hyper-surface thermally emits
hadrons at its temperature and chemical potential, which are then boosted
according to the underlying flow velocity [5]. This prescription produces
particle thermal spectra, which also include perturbative out-of-equilibrium
corrections from the viscous part of the fluid’s energy-momentum tensor [230,
231, 232]. The details of such corrections to the particle spectra depend on
the underlying theory and, in practice, are an additional source of
theoretical uncertainty [233, 234, 162, 235]. Especially if the viscous
corrections are large at particlization, observables can be affected by the
underlying assumptions, especially for experimental observables with
$p_{T}\gtrsim 2$ GeV.

For small systems, the available energy in the system could be limited. As the
particle multiplicity becomes small, one would need to consider the effects
coming from the conservation of energy, momentum, and conserved charges [236,
237, 238].

Because the fluid cells in small systems particlize at a similar local energy
density compared to that in large collision systems [61], the system’s final
size is comparable with those in large collision systems at the same final-
state hadron multiplicity [239]. Due to the faster local expansion rate in
small systems, the collision rate in the hadronic scatterings was found to be
lower than that in the large systems at the same particle multiplicity [240,
32]. The hadronic rescatterings effects are not significant in small systems
[32].

##### The role of longitudinal structure in full (3+1)D simulations for small
systems

Many small systems are asymmetric light+heavy collisions, which explicitly
break the longitudinal boost-invariance at the collision geometry level.
Because of longitudinal momentum conservation, the initial-state energy
density profiles have non-trivial longitudinal structures, which result in
rapidity-dependent particle production and anisotropic flow coefficients in
the final state [241, 242, 243, 84, 85, 244, 245, 246]. We will discuss the
phenomenological studies concerning these effects in Sec. 5.

Measurements of rapidity decorrelation of flow observables study how the
anisotropic flow vector fluctuates and evolves as a function of the particle
rapidity [247, 248]

$r_{n}(\eta)=\frac{\Re\\{\langle
Q_{n}(-\eta)(Q_{n}^{\mathrm{ref}})^{*}\rangle\\}}{\Re\\{\langle
Q_{n}(\eta)(Q_{n}^{\mathrm{ref}})^{*}\rangle\\}},$ (6)

where the $n$-th order anisotropic flow vector $Q_{n}$ is defined as

$Q_{n}\equiv\sum_{j=1}^{M}e^{in\phi_{j}},$ (7)

and the sum runs over all $M$ final state particles and their transverse
momenta’s azimuthal angles $\\{\phi_{j}\\}$. The slope of $r_{n}$ as a
function of $\eta$ measures how fast the anisotropic flow vector
$Q_{n}(-\eta)$ decorrelates from $Q_{n}(\eta)$. The reference flow vector
$Q_{n}^{\mathrm{ref}}$ is usually chosen at forward or backward rapidity with
at least one unit of rapidity gap from both the $Q_{n}(\pm\eta)$ vectors to
suppress non-flow correlations. In light+heavy collisions, the asymmetric
collision geometry and event-by-event longitudinal fluctuations contribute to
the decorrelation of the flow vectors in $r_{n}(\eta)$. Full (3+1)D dynamic
evolution is essential to study the flow rapidity decorrelation measurements
and constrain the longitudinal dynamics of the collision systems [249, 250,
243].

The (3+1)D dynamical evolution of the relativistic heavy-ion collisions allows
us to interpret the final-state anisotropic flow vector $Q_{n}(\eta)$ as a
convolution of the initial state eccentricities $\varepsilon_{n}(\eta_{s})$
with the hydrodynamic response function $G_{n}(\eta-\eta_{s})$ [251, 252]

$Q_{n}(\eta)=\int d\eta_{s}G_{n}(\eta-\eta_{s})\varepsilon_{n}(\eta_{s}).$ (8)

This equation is a generalization of the simple linear response
$Q_{n}=k_{n}\varepsilon_{n}$ observed in boost-invariant simulations (for the
harmonic order $n=2,3$). Because the hydrodynamic response function
$G(\eta-\eta_{s})$ is a two-point function, thermal fluctuations in the
hydrodynamic phase also contribute to the hydrodynamic response function and
the rapidity decorrelation. A large flow rapidity decorrelation would reduce
the sensitivity of the final-state flow measurements to the initial geometry.

Performing (3+1)D dynamics is crucial when studying experimental observables
related to the evolution of longitudinal flow velocity, for example, the fluid
vorticity and $\Lambda$ polarization [253, 254]. In asymmetric light+heavy
collision systems, the light nucleus drills through the heavy nucleus,
resulting in a nontrivial longitudinal flow profile. The transverse coordinate
dependence of the longitudinal flow is directly related to the size of
initial-state flow vorticity in the collision system, which imprints itself to
the final-state polarization of the $\Lambda$ hyperons. Theoretical work
predicted the existence of “smoke ring” patterns for flow vorticity in
asymmetric collision systems [255]. Future measurements of such a pattern in
the $\Lambda$ polarization would verify and probe the fluid dynamic nature of
the produced collision system at the length scale of velocity gradients [256,
257].

###  3 Phenomenological studies in small systems

In this section, we will review the diverse signals associated with quark-
gluon plasma formation in the soft sector of small systems. We will conduct a
thorough comparison between the model calculations and experimental data.

####  1 Collective flow

Collective flow is one of the most important observables in relativistic
nuclear collisions. Its study provides valuable information on the initial
state and several QGP properties. Collective flow observables determine the
anisotropy in particle transverse momentum distributions correlated with the
flow symmetry plane $\Psi_{n}$ [258]. The various characteristic patterns of
the anisotropic flow can be obtained from a Fourier expansion of the azimuthal
particle distribution in a given event [259]:

$\frac{{\rm d}N}{{\rm d}\varphi}\propto
1+2\sum_{n=1}^{\infty}v_{n}\,e^{in(\varphi-\Psi_{n})}$ (9)

where $v_{n}=\langle\cos\,[n(\varphi-\Psi_{n})]\rangle$ is the anisotropic
flow coefficient and $\Psi_{n}$ the corresponding flow symmetry plane, both
fluctuating event by event. Here, the angular average $\langle O\rangle$ is
defined as

$\langle O\rangle\equiv\frac{\int d\phi\,O\,dN/d\phi}{\int d\phi\,dN/d\phi}.$
(10)

In general, the first coefficient, $v_{1}$, is called _directed flow_ , the
second coefficient, $v_{2}$, is called _elliptic flow_ and the third
coefficient $v_{3}$, is called _triangular flow_. For $n\geq 3$, we refer to
the $v_{n}$ as higher-order flow harmonics. Because the global rotation of the
collision orientation introduces an arbitrary random phase event by event, the
flow symmetry plane cannot be directly measured. Therefore, the anisotropic
flow coefficient $v_{n}$ cannot be measured using single-particle information.
A popular approach is the event-plane method [260], where the azimuthal
correlation of emitted particles is measured with respect to an event-plane.
However, it was found that the results from the event-plane method strongly
depend on the resolution of the event-plane (limited mainly by the finite
number of measured particles), which introduces an uncontrolled bias in the
measurement [261]. As an alternative approach, multi-particle azimuthal
correlation measurements [262, 263] have been employed as they allow for a
more robust measurement of the underlying anisotropic flow.

The $Q$-cumulant method measures the flow harmonics $v_{n}$ from multi-
particle correlations without knowing the event plane. With the flow
$Q_{n}$-vector defined in Eq. (7), the 2-and 4-particle azimuthal correlations
in a single event can be calculated as [262]:

$\displaystyle\langle 2\rangle_{n,-n}$
$\displaystyle=\frac{|Q_{n}|^{2}-M}{M(M-1)},$ (11) $\displaystyle\langle
4\rangle_{n,n,-n,-n}$
$\displaystyle=\frac{|Q_{n}|^{4}+|Q_{2n}|^{2}-2\cdot{\rm{Re}}[Q_{2n}Q_{n}^{*}Q_{n}^{*}]}{M(M-1)(M-2)(M-3)}$
$\displaystyle\ \
~{}~{}-2\frac{2(M-2)\cdot|Q_{n}|^{2}-M(M-3)}{M(M-1)(M-2)(M-3)},$

where $M$ is the particle multiplicity for the flow vector. Here, we have used
the general notation of the single-event $k$-particle correlators $\langle
k\rangle_{n_{1},n_{2},...,n_{k}}\equiv\langle\cos(n_{1}\varphi_{1}\\!+\\!n_{2}\varphi_{2}\\!+\\!\cdots\\!+\\!n_{k}\varphi_{k})\rangle\,(n_{1}\geq
n_{2}\geq\cdots\geq n_{k})$ and $\langle...\rangle$ means an average over all
the particles in a single event. After averaging over all events within the
selected centrality bin, the obtained 2- and 4-particle cumulants are:

$\displaystyle c_{n}\\{2\\}$ $\displaystyle=\langle\langle
2\rangle\rangle_{n,-n},$ (12) $\displaystyle c_{n}\\{4\\}$
$\displaystyle=\langle\langle 4\rangle\rangle_{n,n,-n,-n}-2\cdot\langle\langle
2\rangle\rangle^{2}_{n,-n},$

Then, the 2- and 4-particle integrated flow harmonics can be calculated as
[262]:

$\displaystyle v_{n}\\{2\\}$ $\displaystyle=\sqrt{c_{n}\\{2\\}},\quad
v_{n}\\{4\\}$ $\displaystyle=\sqrt[4]{-c_{n}\\{4\\}}.$ (13)

In general, the 4-particle correlations used to determine the flow harmonics
$v_{n}\\{4\\}$ can largely suppress the non-flow effects from jets, resonance
decays, etc. Non-flow effects significantly influence $v_{n}\\{2\\}$ obtained
from the 2-particle correlations. To suppress such non-flow effects, one
divides the whole event into sub-events with a certain pseudorapidity gap
$|\Delta\eta|$ and then calculates the 2-particle azimuthal correlations to
further suppress the non-flow contamination.

The high energy proton–lead ($p$+Pb) and proton–proton ($p$+$p$) collisions at
the LHC were originally aimed to provide the reference data (the “no-QGP”
reference) for high energy nucleus-nucleus collisions. However, various
unexpected phenomena have been observed in these small systems, especially in
the high multiplicity region. One surprising discovery was the long-range
“ridge” structure in two-particle azimuthal correlations with a large pseudo-
rapidity separation in high multiplicity $p$+Pb and $p$+$p$ collisions. Such
long-range correlation structures were first discovered in Au+Au and Pb+Pb
collisions and interpreted as a signature of the collective expansion. This
has also triggered the study of collectivity in other small systems, such as
$p$+Au, $d$+Au, and 3He+Au at RHIC [264, 265, 60, 85], as well as the system
size scan of various collision systems at the LHC [266]. Keeping in mind the
caveats discussed in Section 2, the observed flow-like signals in the small
systems can be quantitatively or semi-quantitatively described by hydrodynamic
calculations, which translate initial spatial anisotropies into final momentum
anisotropies of produced hadrons with the collective expansion of the bulk
matter.

Figure 2: Anisotropy coefficients $v_{2}\\{2\\}$ (a) and $v_{3}\\{2\\}$ (b)
for charged hadrons vs. charged hadron multiplicity in various collision
systems at LHC, compared to experimental data from the ALICE Collaboration
[267]. The figures are taken from published work.[61].

Fig. 2 shows hydrodynamic calculations with IP-Glasma initial conditions [61]
of the charged hadron $v_{2}\\{2\\}$ (left) and $v_{3}\\{2\\}$ (right) as a
function of charged hadron multiplicity from small to large collision systems
at LHC and compares to experimental data [267]. Overall, the hydrodynamic
model reproduces the multiplicity dependence of the integrated $v_{2}\\{2\\}$
and $v_{3}\\{2\\}$ from $p$+Pb to Pb+Pb collisions well. In particular,
hydrodynamics seems to reproduce well the $v_{2}\\{2\\}$ in Xe+Xe and Pb+Pb
collisions from peripheral to central collisions. For $p$+Pb collisions, the
agreement becomes worse, and for $p$+$p$ collisions, this hydrodynamic
simulation misses the experimentally observed magnitude and multiplicity
dependence. The situation is similar for $v_{3}\\{2\\}$. The agreement is best
for more central larger systems. The $v_{3}\\{2\\}$ is rather insensitive to
the system’s average geometry, as it is driven solely by fluctuations. Similar
to $v_{2}\\{2\\}$, hydrodynamic calculations underestimate the experimental
data for $v_{3}\\{2\\}$ both in $p$+$p$ and $p$+Pb collisions, which requires
further study in the future.

Figure 3: Four-particle correlations, $c_{2}\\{4\\}$, as a function of
$N_{ch}$ for charged hadrons in pp collisions 13 TeV, calculated by
hydrodynamic model. The CMS data and the ATLAS data are taken from [268] and
[269], respectively. This figure is taken from published work. [270].

####  2 Multi-particle correlations

Compared with the 2-particle correlations, multi-particle cumulants are less
influenced by the non-flow effects and thus are key observables to evaluate
the anisotropic collectivity in small systems. Additionally, the multi-
particle and 2-particle cumulants show different sensitivities to event-by-
event flow fluctuations [271, 272, 273]. An extensive measurement of these
different cumulants could provide tight constraints on the initial state
fluctuations. To extract real values of the flow coefficients, the 2-, 4-, 6-,
and 8-particle cumulants are expected to carry positive, negative, positive,
and negative signs, respectively. Such “changing sign pattern” has been
observed in the measured 2- and multi-particle cumulants in Pb+Pb collisions
at the LHC [274, 275]. Similarly, 2- and multi-particle cumulants have been
measured in high multiplicity $p$+$p$ and $p$+Pb collisions at the LHC.
However, it was found that the standard multi-particle cumulants in the small
systems are still largely affected by the residual short-range non-flow and
its fluctuations, which even presents fake flow signals with the “right sign”.
To further remove the residual non-flow from jets, 2- and 3-subevent methods
for the multi-particle cumulants were developed.[276, 277] With this
technique, several experimental collaborations have confirmed the observations
of positive 2-particle cumulants and negative 4-particle cumulants in high
multiplicity $p$+$p$ and $p$+Pb collisions at the LHC [278].

With properly tuned parameters, hydrodynamic simulations can quantitatively or
at least semi-quantitatively describe the two-particle correlations in $p$+Pb
and $p$+$p$ and multi-particle correlations in $p$+Pb collisions at the LHC
and RHIC (see reviews [279, 280, 281, 18, 60, 3, 19]). However, the measured
negative $c_{2}\\{4\\}$ in $p$+$p$ collisions, which could naively be
interpreted as evidence of hydrodynamic flow, could not be reproduced by any
of the hydrodynamic calculations on the market [270, 61, 282]. As shown in
Fig. 3, hydrodynamic model calculations with different initial conditions
yield positive values for $c_{2}\\{4\\}$ in $p$+$p$ collisions. It is still
unknown if the wrong sign of $c_{2}\\{4\\}$ is because of the incorrect
fluctuation spectrum in initial conditions or due to issues with the
applicability of the hydrodynamic model to $p$+$p$ collisions.

####  3 Partonic collectivity and strangeness enhancement in small systems

The origins of the observed collective behavior in small systems are still
under debate. For the soft hadron data measured in high-multiplicity $p$+Pb
collisions, hydrodynamics or transport models based on final-state effects can
describe many of the flow-like signals. On the other hand, color glass
condensate models [283, 284, 285, 67] including the IP-Glasma model [66, 286],
which by themselves produce momentum anisotropies as initial-state effects,
have also been used to explain some features of the observed collectivity.
However, generally, the initial state models produce the wrong multiplicity
dependence of the anisotropic flow. As discussed in Section 1, the usual
assumption of boost-invariance of the signal has been challenged, as has the
probability of its survival if 3+1 dimensional early time evolution is
included.

In large systems, hard probes have been used as signatures for QGP formation.
In small systems, due to their limited sizes and lifetime, the energy lost by
energetic partons no longer leads to obvious signatures at high $p_{T}$ to
discern if the QGP is formed. The relatively small nuclear modification
effects for large $p_{T}$ light and heavy flavor hadrons and jets measured in
small systems are consistent with the expectations of cold nuclear matter
effects [287, 288, 289, 290, 291, 292]. For more information about the hard
probes in small systems, we refer the reader to subsection 7.

Besides collective flow and the quenching of energetic jets, the number of
constituent quark (NCQ) scaling of elliptic flow and strangeness enhancement
was proposed long ago as signals of QGP formation in high-energy nuclear
collisions. The ATLAS, CMS, and ALICE Collaborations have measured the
$v_{2}(p_{T})$ of charged and identified hadrons with high precision in high
multiplicity events of small systems [293, 294, 295]. The resulting data show
a similar approximate NCQ scaling of $v_{2}$ at intermediate $p_{T}$ as
observed in heavy ion collisions. Also, strangeness enhancement in small
colliding systems was reported by the ALICE Collaboration [296]. Yield ratios
of (multi)strange hadrons to charged pions exhibit a monotonic and continuous
increase as functions of multiplicity across collision systems. In high-
multiplicity events of small systems, the strange hadron yield ratios almost
reach similar values to those in heavy-ion collisions.

Figure 4: (Color online) Differential elliptic flow $v_{2}(p_{T})$ of pions,
kaons and protons in 0-20% p-Pb at $\sqrt{s_{NN}}=5.02$ TeV measured in
experiments (left window) and calculated with the Hydro-Coal-Frag model
(middle window) and the Hydro-Frag model (right window). The ALICE, CMS, and
ATLAS data are taken from [297], [293] and [294], respectively. This figure is
taken from published work. [298]

The role of partonic degrees of freedom in high multiplicity $p$+Pb collisions
was investigated in the quark coalescence model [298]. The partonic
collectivity will manifest itself at intermediate $p_{T}$ of hadrons after the
quark coalescence process. In this model, the mesons and baryons are coalesced
from the quark and antiquark employing the overlaps of the quarks’ phase-space
distributions and the hadrons’ Wigner functions [299, 300]. In the calculation
shown in Fig. 4 [298], the coalescence process includes thermal-thermal,
thermal-hard and hard-hard parton recombinations, where the thermal quarks are
taken from the hydrodynamic medium and hard partons are taken from the linear
Boltzmann transport energy loss model [301, 302, 303]. The remnant hard
partons are grouped into strings and fragmented to hard hadrons using the
“hadron standalone mode” of PYTHIA8 [304].

As shown in Fig. 4, the newly developed hybrid hadronization model nicely
describes the measured proton-to-pion ratio and the $v_{2}(p_{T})$ of pions,
kaons and protons from 0 to 6 GeV in high multiplicity $p$+Pb collisions at
5.02 TeV. Specifically, the low $p_{T}$ mass ordering of $v_{2}(p_{T})$ of
identified hadrons is reproduced by the hydrodynamic part. For $p_{T}>2.5$
GeV, the $v_{2}(p_{T})$ of protons becomes larger than that of pions and
kaons, which is attributed to the quark coalescence process. In contrast,
without the quark coalescence process, the model fails to describe the
$v_{2}(p_{T})$ of identified hadrons at intermediate $p_{T}$. This
demonstrates that including the quark coalescence contribution to the
production of hadrons is essential in reproducing the measured $v_{2}(p_{T})$
of identified hadrons at intermediate $p_{T}$. It thus provides a strong
indication for the existence of partonic degrees of freedom and the possible
formation of the QGP in high multiplicity $p$+Pb collisions at 5.02 TeV.

Figure 5: (Color online) Hadron yield ratios of (a) cascades ($\Xi^{-}$ and
$\bar{\Xi}^{+}$) and (b) lambdas ($\Lambda$ and $\bar{\Lambda}$) to charged
pions ($\pi^{-}$ and $\pi^{+}$) as functions of multiplicity at midrapidity,
$\left<dN_{\mathrm{ch}}/d\eta\right>$, from $p$+$p$ to Pb+Pb collisions at LHC
energies. The ALICE data is taken from [305, 306, 307, 308, 309, 310, 296].
This figure is taken from published work. [102].

Concerning strangeness enhancement, since the colliding nuclei do not contain
a strange valence quark, yields of strange hadrons should be sensitive to
details of the reaction dynamics. A dynamical core–corona framework to study
the event activity dependence of the hadron yield ratios in various collision
systems and over a wide range of collision energies was developed.[102] In
this approach, the system generated in high-energy nuclear collisions is
described with two components: equilibrated matter (core) and non-equilibrated
matter (corona). Non-equilibrated partons can act as sources of QGP fluid as
they lose energy and momentum. The core/fluid is described by hydrodynamics.
The low-density corona, where partons suffer only a few collisions, is treated
microscopically, and hadrons are produced via string fragmentation. Fig. 5
shows that, because of the interplay between core and corona components, the
model reasonably describes the strangeness of hadron yield ratios to charged
pions as functions of multiplicity ranging from $p$+$p$ to $p$+Pb to Pb+Pb
collisions. The ratios increase up to $\left<dN_{ch}/d\eta\right>\sim 100$ and
saturate in high multiplicity events. This tendency implies that the
contribution of the fluid becomes large and dominant in high multiplicity
events. Importantly, at a given multiplicity value, the high multiplicity
small systems have similar values of the strange hadron ratios to large
systems, where the system reaches the chemically equilibrated fluid limit.
This result provides a strong indication of partial QGP generation in high
multiplicity small colliding systems.

####  4 The $v_{n}$-$p_{T}$ correlations

An interesting observable involving the correlation between the elliptic
momentum anisotropy $v_{2}$ and the average transverse momentum $[p_{T}]$ has
been proposed to help constrain initial state models.[311] This observable has
been used to extract information on the nucleon size and on the origin of the
observed momentum anisotropy [312, 313]. The correlation between $v_{n}^{2}$
and $[p_{t}]$ expressed by the Pearson correlation coefficient is [311]

$\hat{\rho}_{n}\equiv\hat{\rho}(v_{n}^{2},[p_{t}])=\frac{\langle\delta
v_{n}^{2}\delta[p_{t}]\rangle}{\sqrt{\langle(\delta
v_{n}^{2})^{2}\rangle\langle(\delta[p_{t}])^{2}\rangle}},$ (14)

where $\delta O=O-\langle O\rangle$ for any observable $O$ and the flow
anisotropic coefficient $v_{n}$ is defined in Eq. (9).

Figure 6:  Top: The correlator $\hat{\rho}(v_{2}^{2},[p_{t}])$ (circles, solid
lines) together with estimators based on the initial geometry
($\hat{\rho}_{\rm est}(\varepsilon_{2}^{2},[s])$, stars) and the initial
momentum anisotropy ($\hat{\rho}_{\rm est}(\varepsilon_{p}^{2},[s])$, squares)
in d+Au collisions at $\sqrt{s}=200\,{\rm GeV}$, and the Pearson coefficients
between $v_{2}$ and the initial ellipticity (stars) and the initial momentum
anisotropy (squares), respectively. Bottom: Results of the IP-
Glasma+MUSIC+UrQMD framework for $\rho(v_{2}^{2},[p_{t}])$ with different
types of shaded bands correspond to different values of the nucleon width in
Pb+Pb collisions. These figures are taken from published works. [312, 313].

Two predictors for the $\hat{\rho}$-correlator as a function of multiplicities
are presented in the top panel of Fig. 6 for $200\,{\rm GeV}$ $d$+Au
collisions. The $\hat{\rho}_{\rm est}(\epsilon_{2}^{2},[s])$, is based
entirely on the initial geometry using the initial spatial eccentricity
$\epsilon_{2}$, which is a good estimator for $v_{2}$. The $\hat{\rho}_{\rm
est}(\epsilon_{p}^{2},[s])$ uses the initial momentum anisotropy
$\epsilon_{p}$ as estimator for $v_{2}$. It is computed from the initial
energy-momentum tensor of the IP-Glasma model as [60]

$\mathcal{E}_{p}\equiv\epsilon_{p}e^{i2\psi_{2}^{p}}\equiv\frac{\langle
T^{xx}-T^{yy}\rangle+i\langle 2T^{xy}\rangle}{\langle
T^{xx}+T^{yy}\rangle}\,,$ (15)

evaluated at $\tau=\,0.1\,{\rm fm}$, where here $\langle\cdot\rangle$ is
defined as average over the transverse plane [60]. The $[p_{T}]$ is estimated
from the average initial entropy density $[s]$, where $[s]=[e^{3/4}]$ for
ideal parton gas at $\tau=\,0.1\,{\rm fm}$. The top panel in Fig. 6 also
displays the Pearson coefficients of $\mathcal{E}_{2}$ with $V_{2}$ and
$\mathcal{E}_{p}$ with $V_{2}$, which is defined as [314, 315, 316]

$Q(\mathcal{E},V_{2})=\frac{{\rm
Re}\langle\mathcal{E}V_{2}^{*}\rangle}{\sqrt{\langle|\mathcal{E}|^{2}\rangle\langle|V_{2}|^{2}\rangle}}\,,$
(16)

where $V_{2}$ is the complex valued $2^{\rm nd}$ order flow harmonic.

One can see that for higher multiplicities $\hat{\rho}(v_{2}^{2},[p_{T}])$
approaches the geometric estimator, while at lower multiplicities, the initial
momentum anisotropy predicts $\hat{\rho}(v_{2}^{2},[p_{T}])$ better. Based on
the color domain interpretation of the initial state momentum anisotropy [317,
318, 319, 320, 67], the $\hat{\rho}_{\rm est}(\epsilon_{p}^{2},[s])$ is
expected to be positive, because at fixed multiplicity, a larger $[p_{T}]$
selects events with smaller transverse size. This reduces the number of color
domains with an average size of $1/Q_{s}$, which enhances the magnitude of
initial momentum anisotropy in the CGC description [67].

The Pearson coefficients $Q(\mathcal{E},V_{2})$ show that the behavior of
$\hat{\rho}(v_{2}^{2},[p_{T}])$ is a result of the geometry dominating the
elliptic flow in high multiplicity events, and the initial momentum anisotropy
driving the final $v_{2}$ at low multiplicity. The Pearson coefficients were
studied already before [60], but they are not experimentally observable. In
contrast, the $\hat{\rho}(v_{2}^{2},[p_{T}])$ is an observable whose sign
change as a function of multiplicity is an indicator of the origin of the
elliptic flow in small systems, and the presence of initial state momentum
anisotropies as predicted from the color glass condensate.

We note, however, that as discussed earlier, the above study assumed boost
invariance and that recent studies going beyond that assumption indicate that
the initial state momentum anisotropies are short-range in rapidity
($|\Delta\eta|<1$) [69]. Given the rapidity gaps used in the experimental
determination of $\hat{\rho}(v_{2}^{2},[p_{T}])$, this means that the initial
state momentum contribution should be much smaller than shown in Fig. 6.
Furthermore, a complete CGC calculation beyond the classical limit, including
the rapid growth of small fluctuations [70], is also expected to wash out
initial momentum anisotropies.

The bottom panel of Fig. 6 shows the $v_{n}^{2}$-$[p_{t}]$ correlations,
modeled using different nucleon sizes in Pb+Pb collisions at 5.02 TeV. Note
that $\hat{\rho}_{2}$ for a smooth nucleon with $w=0.4$ fm differs above 50%
centrality from the result with sub-nucleon constituents, showing that the
details of the sub-nucleon structure are important. The figure reports as well
experimental data from the ATLAS [321] and the ALICE [322] collaborations at
LHC. Data is in qualitative disagreement with calculations implementing
$w=0.8$ fm or larger and suggests that it is indeed possible to constrain the
size of nucleons (or their constituents) from such observations in heavy-ion
collisions. It also shows the results from a low-viscosity run, which
demonstrates that $\hat{\rho}_{n}$ are driven by initial-state properties and
are largely insensitive to medium effects. The sensitivity of
$\hat{\rho}(v_{3}^{2}$-$[p_{t}])$ to variations in $w$ is also noteworthy.
Nucleon structure properties significantly influence the observable of final-
state two-particle correlations, particularly in peripheral collisions.
Notably, the sign of $\hat{\rho}(v_{3}^{2}$-$[p_{t}])$ can flip depending on
the value of $w$ in peripheral collisions.

####  5 Rapidity dependent anisotropic flow in small systems

Figure 7:  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 3+1D hydrodynamic
model. The PHENIX and STAR data are from [265, 323]. This figure is taken from
published work. [85].  Figure 8:  The pseudorapidity dependence of $p_{T}$
integrated $v_{2}$ in $p$+Pb at $\sqrt{s_{NN}}=5.02$ TeV. Comparison of the
ALICE data [324] (black circles) with a calculation within a 3+1D
hydrodynamical model (blue band) [243] for the 0–5% (top left), 5–10% (top
right), 10–20% (bottom left), and 20–40% (bottom right) centrality classes.
This figure is taken from published work [324].

The (3+1)D dynamics is also crucial to explore the collectivity in central
p+Au, d+Au and $\rm{{}^{3}He}$+Au collisions at RHIC and high-multiplicity
$\gamma^{*}$-nucleus collisions in ultra-peripheral collisions (UPC) at the
LHC [243, 85, 325]. Figure 7 shows the $v_{n}(p_{T})$ ($n=2,3$) for charged
hadrons compared between (3+1)D hydrodynamic model simulations and the
experimental data from the PHENIX and STAR Collaborations. The (3+1)D model
gives an overall good description of the PHENIX data for d+Au and
$\rm{{}^{3}He}$+Au collisions. The model has some tension with the STAR data
for $v_{3}(p_{T})$ in p+Au and d+Au collisions. [85, 325]

It is crucial to note that both in the model calculations and experimental
data, the $v_{3}(p_{T})$ with the STAR definition are systemically larger than
those determined using the PHENIX definition. In (3+1)D hydrodynamic
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 used by the two collaborations. Figure 8 shows
the pseudorapidity dependence of $p_{T}$ integrated $v_{2}(\eta)$ in $p$+Pb at
5.02 TeV. In the (3+1)D model calculations, the longitudinal decorrelation of
flow vectors plays a significant role in shaping the pseudorapidity dependence
of $v_{2}$. For UPCs, the incoming energies between the virtual photons and
the nucleus are highly unbalanced and fluctuate event by event. The large
global rapidity shift between the center-of-mass frame and the lab frame
introduces the non-trivial longitudinal decorrelations of flow vectors. This
underscores the significance of employing (3+1)D model simulations in small
systems.

####  6 Electromagnetic radiation in small systems

Black-body radiation is a characteristic phenomenon for any finite-temperature
system in equilibrium. Electromagnetic (EM) radiation from relativistic heavy-
ion collisions is clean because photons and dileptons suffer negligible final-
state interactions after production [326, 327, 328]. An observation of thermal
EM radiation is a smoking gun for producing thermalized hot matter in
collisions. Furthermore, their transverse momentum spectra encode information
on local temperature and blue shift from flow velocity. Therefore, searching
for thermal EM radiation in small systems can provide direct evidence of
producing a nearly equilibrated QCD medium in these collision systems [240,
32, 329, 330].

Detecting thermal photon production from relativistic nuclear collisions has
been challenging experimentally because of overwhelming background radiation
from hadronic decays and prompt photon production. Nevertheless, the PHENIX
Collaboration measured the direct (including thermal) photon radiation in
small systems, such as ($p$, $d$, 3He)+Au collisions at 200 GeV. A hint of
thermal photon enhancement above the prompt photon background was observed in
0-5% most central collision events [331]. While the measurements were
consistent with the theoretical predictions based on hydrodynamic evolution
and thermal photon emission rates, they still contain considerable
uncertainty, preventing any decisive conclusion.

Future systematic studies of thermal photon radiation as a function of charged
hadron multiplicity across multiple collision systems [332] would shed
additional light on characterizing the nature of the QCD matter produced in
small systems.

####  7 Hard probes and the $R_{pA}$ and $v_{2}(p_{T})$ puzzle in small
systems

Motivated by the significant azimuthal anisotropies measured in small systems
for low-$p_{T}$ hadrons with patterns similar to those observed in heavy-ion
collisions, measurements aimed at detecting signs of jet quenching in these
small collision systems were performed. However, they have found no such
effect. Analyses of high $p_{T}$ hadron and jet $p_{T}$ spectra indicate
production yields consistent with those in $p$+$p$ collisions scaled by the
expected nuclear thickness in $p$+Pb [333, 287, 334] and $d$+Au collisions
[335]. On the other hand, the experimental observations of the hadron
azimuthal anisotropy up to $p_{T}\approx 12$ GeV indicate a non-zero
anisotropy extending into the region beyond the usual hydrodynamic
interpretation and into the regime of jet quenching [333, 334]. In heavy ion
collisions, this effect was interpreted as resulting from directionally
dependent jet quenching, driven by the anisotropic shape of the underlying
medium. However, it is unlikely that there can be differential jet quenching
as a function of orientation relative to the QGP geometry if there is no jet
quenching in $p$+Pb collisions in the first place.[336]

Thus, there are two related outstanding puzzles, one being the lack of jet
quenching observed in the spectra, if indeed small droplets of QGP are formed,
and the other being the mechanism that leads to high-$p_{T}$ hadron
anisotropies other than differential jet quenching. In recent developments,
some progress has been achieved in addressing this issue [337, 338].

###  4 Further opportunities with small systems

Figure 9:  Charged hadron anisotropic flow coefficients $v_{2}\\{2\\}$ and
$v_{3}\\{2\\}$ as functions of charged hadron multiplicity $N_{\rm ch}$ in
$p$+Pb (dashed lines) and $\gamma^{*}$+Pb (solid lines) collisions at LHC
energies from 3+1D hydrodynamics simulations. The ATLAS data is from [294,
339]. This figure is taken from published work. [243].

####  1 Confronting ab initio nuclear structure with small system measurements

The collision geometry of small systems is dominated by the shape of light
nuclei, which can be computed using ab initio methods in nuclear structure
theory [340]. Therefore, small collision systems, especially light+heavy
collisions, can serve as excellent cameras for the light colliding nuclei
[341]. Unlike low-energy nuclear experiments, relativistic nuclear collisions
can probe the multi-nucleon wavefunction at the collision impact [342, 343,
53]. Therefore, precision measurements on anisotropic flow coefficients can
set complementary constraints on the multi-particle correlations inside the
nuclei. On the other hand, the precision nuclear structure inputs to
relativistic nuclear collisions can significantly reduce the theoretical
uncertainty in modeling the collisions of small systems.

####  2 Accessing features of saturated QCD systems

An alternative approach to unraveling the nature of long-range azimuthal
anisotropies in small systems emerges within the framework of the color glass
condensate effective theory. In this framework, the effective degrees of
freedom are color sources at large $x$ and gauge fields at small $x$. At high
energies, the former are frozen color configurations, randomly distributed
event-by-event due to time dilation. The latter are dynamic fields coupled to
the static color sources. In the CGC, small $x$ gluons saturate with large
occupation numbers, and their typical momenta peak at the saturation scale
$Q_{s}$. These gluon fields exhibit correlations within a characteristic
length of $1/Q_{s}$, which can be interpreted as color field domains of that
size. Within this picture, given an initial non-perturbative distribution of
sources at an initial scale, one can compute systematically $n$-point gluon
correlation functions and their evolution with $x$ order by order in
perturbation theory [344].

For the anisotropic flow in small systems, it was demonstrated that multi-
gluon production from the CGC gives rise to notable correlations in azimuthal
anisotropies across long-range rapidities in momentum space without the need
for final state interactions.[345, 346, 347] Recent calculations also have
demonstrated that CGC effective theory can capture certain collective features
in small systems, such as multiparticle cumulants [348] and the mass ordering
of anisotropy coefficients [286, 349]. However, as alluded to earlier, the CGC
calculations struggle with reproducing the magnitudes or the correct
systematics of the experimental data with multiplicity or collision system
(e.g., $p$+Au, $d$+Au, 3He+Au at RHIC [265]). This suggests that final state
effects are necessary to describe the data, at least for high multiplicities
regions in small systems.

####  3 Probing the boundaries of collective behavior

Intriguing experimental results on two-particle azimuthal correlations in
ultra-peripheral Pb+Pb collisions (UPCs) at the LHC have recently appeared
[339]. These UPCs involve appreciable rates of photo-nuclear interactions
[350, 351], and the ATLAS measurements of such photo-nuclear ($\gamma^{*}$+Pb)
interactions in Pb+Pb UPCs indicate the persistence of collective phenomena
with the strength of correlations comparable to that observed in proton-proton
and proton-lead collisions in similar multiplicity ranges [339].

On the theoretical front, hydrodynamic simulations with the appropriate
initial state model have proven effective in describing the systematic
behavior of measured azimuthal anisotropies. These simulations span a wide
range of systems, from high-energy heavy-ion collisions to small hadronic
collision systems producing only a limited number of charged hadrons per unit
of rapidity. Collective QGP signatures in $\gamma^{*}$+Pb collisions at the
LHC were explored using a full (3+1)D dynamical framework[84] with
hydrodynamics and hadronic transport.[243] Fig. 9 shows the multiplicity
dependence of the $p_{T}$-integrated anisotropic flow coefficients
$v_{2}\\{2\\}$ and $v_{3}\\{2\\}$ computed with two subevents with the
$|\Delta\eta|>2.0$ for $\gamma^{*}$+Pb and $p$+Pb collisions. Extrapolating
from $p$+Pb to $\gamma^{*}+$Pb collisions, hydrodynamic calculations, assuming
strong final-state interactions, reproduce the hierarchy observed for the
elliptic flow coefficient $v_{2}$ in $N_{\mathrm{ch}}\in[20,60]$ in the ATLAS
data. Within this framework, the elliptic flow hierarchy between $p$+Pb and
$\gamma^{*}$+Pb collisions is dominated by the difference in longitudinal flow
decorrelations. The hydrodynamic model predicts triangular flow in
$\gamma^{*}+$Pb collisions smaller than that in $p$+Pb collisions at the same
charged hadron multiplicity, again because of the larger longitudinal
decorrelation, which does not agree with the ATLAS data. The magnitude of
$v_{3}\\{2\\}$ in $\gamma^{*}+$Pb collisions may be sensitive to the vector
meson’s detailed substructure fluctuations.

The color glass condensate effective theory also generates significant
correlations in the initial state of $\gamma^{*}$+Pb collision systems. It has
also been shown to mimic collective behavior to a certain degree [352].
Notably, the hydrodynamic model and the CGC effective theory predict the
opposite dependence of the collectivity on the $\gamma^{*}$’s virtuality
$Q^{2}$. In $\gamma^{*}$+Pb collisions, the average size of the projectile
vector meson decreases with increasing $Q^{2}$. Thus, hydrodynamics predicts
that vector mesons with smaller virtuality lead to greater elliptic flow
coefficients. This is attributed to the increased transverse space for
geometric fluctuations, resulting in larger average ellipticities. In
contrast, the CGC model predicts that vector mesons with smaller virtuality
result in weaker elliptic flow coefficients. This is because the increased
number of independent color domains in the transverse space reduces two-
particle correlations on average. Future experiments at an Electron-Ion
Collider (EIC) will provide direct access to the photon’s virtuality, allowing
for systematic tests of predictions from both the hydrodynamic and CGC
frameworks.

More recently, the CMS Collaboration has reported a new measurement of two-
particle correlations of particles within a jet cone in $p$+$p$ collisions at
13 TeV [353]. For low multiplicity events, the extracted $v_{2}\\{2\\}$ was
found to be consistent with the PYTHIA8 or SHERPA model results, where the
observed correlations arise from short-range correlations associated with the
presence of jets or mini-jets. Surprisingly, for very high multiplicity jets,
the CMS data reveal an increasing trend in the value of $v_{2}$. The PYTHIA8
and SHERPA models, which do not include the long-range collective effects,
fail to describe this data. This discrepancy raises intriguing questions about
QGP droplets within a jet and calls for further investigations in the future.

###  5 Instead of Conclusions

Research on the physics of small collision systems is ongoing, with several
open theoretical and phenomenological questions. Therefore, instead of
conclusions, in this section, we briefly summarize the progress made to date
and list some important open questions that need to be studied in the future.

In small systems, hydrodynamic + hadronic transport models, together with the
proper initial conditions and early time evolution, can describe most of the
features of the experimental data at a largely quantitative level.
Nonetheless, challenges persist in hydrodynamic simulations of small systems.
Large Knudsen and inverse Reynolds numbers pose significant challenges to
applying hydrodynamics in small systems. Additionally, driven by far-from-
equilibrium effects, numerical simulations may violate causality, potentially
introducing unphysical artifacts.

Furthermore, a multitude of open questions remains. For example, the sign of
the four-particle correlation from hydrodynamic simulations ($c_{2}\\{4\\}>0$)
conflicts with experimental data in proton-proton collisions. Jet quenching
does not seem to be observed in small systems, and the $R_{pA}$ and
$v_{2}(p_{T})$ cannot simultaneously be described within any model framework
in small systems. Clear evidence of enhanced production of electromagnetic
probes, indicating thermal production, is also still outstanding.

Hydrodynamics is being pushed to the edge of validity in these small colliding
systems. This has led to probably the most significant advances in
relativistic fluid dynamic theory in recent decades. Yet, more theoretical
work is still needed to provide a complete understanding of how collectivity
emerges in small systems and if and how it can be consistently described
within a hydrodynamic theory. Experimentally, we are moving to even more
extreme cases, such as ultraperipheral collisions and even the interior of
jets, and the future EIC is likely also to provide new insights into the
origins of collective behavior in multi-particle systems governed by quantum
chromodynamics.

###  6 Acknowledgments

This work is supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, under DOE Contract No. DE-SC0012704 (B.P.S.) and
Award No. DE-SC0021969 (C.S.). J.N. is partially supported by the U.S.
Department of Energy, Office of Science, Office for Nuclear Physics under
Award No. DE-SC0023861. C.S. acknowledges a DOE Office of Science Early Career
Award. 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 and US DOE under Contract No. DE-AC02-05CH11231, and
within the framework of the Saturated Glue (SURGE) Topical Theory
Collaboration.

## References

  * 1. C. Gale, S. Jeon, and B. Schenke, Hydrodynamic Modeling of Heavy-Ion Collisions, _Int. J. Mod. Phys. A_. 28, 1340011 (2013).  10.1142/S0217751X13400113.
  * 2. P. Romatschke and U. Romatschke, _Relativistic Fluid Dynamics In and Out of Equilibrium_. Cambridge Monographs on Mathematical Physics, Cambridge University Press (5, 2019).  ISBN 978-1-108-48368-1, 978-1-108-75002-8.  10.1017/9781108651998.
  * 3. C. Shen and L. Yan, Recent development of hydrodynamic modeling in heavy-ion collisions, _Nucl. Sci. Tech._ 31(12), 122 (10, 2020).  10.1007/s41365-020-00829-z.
  * 4. U. Heinz and R. Snellings, Collective flow and viscosity in relativistic heavy-ion collisions, _Ann. Rev. Nucl. Part. Sci._ 63, 123–151 (2013).  10.1146/annurev-nucl-102212-170540.
  * 5. C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, The iebe-vishnu code package for relativistic heavy-ion collisions, _Comput. Phys. Commun._ 199, 61–85 (2016).  10.1016/j.cpc.2015.08.039.
  * 6. J. H. Putschke et al., The JETSCAPE framework (3, 2019).
  * 7. G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Bayesian analysis of heavy ion collisions with the heavy ion computational framework Trajectum, _Phys. Rev. C_. 103(5), 054909 (2021).  10.1103/PhysRevC.103.054909.
  * 8. C. Shen, U. Heinz, P. Huovinen, and H. Song, Radial and elliptic flow in Pb+Pb collisions at the Large Hadron Collider from viscous hydrodynamic, _Phys. Rev. C_. 84, 044903 (2011).  10.1103/PhysRevC.84.044903.
  * 9. U. Heinz, Z. Qiu, and C. Shen, Fluctuating flow angles and anisotropic flow measurements, _Phys. Rev. C_. 87(3), 034913 (2013).  10.1103/PhysRevC.87.034913.
  * 10. H. Niemi, K. J. Eskola, R. Paatelainen, and K. Tuominen, Predictions for 5.023 TeV Pb + Pb collisions at the CERN Large Hadron Collider, _Phys. Rev. C_. 93(1), 014912 (2016).  10.1103/PhysRevC.93.014912.
  * 11. J. Noronha-Hostler, M. Luzum, and J.-Y. Ollitrault, Hydrodynamic predictions for 5.02 TeV Pb-Pb collisions, _Phys. Rev. C_. 93(3), 034912 (2016).  10.1103/PhysRevC.93.034912.
  * 12. J.-F. Paquet, Applications of emulation and Bayesian methods in heavy-ion physics (10, 2023).
  * 13. M. Aaboud et al., Measurement of long-range multiparticle azimuthal correlations with the subevent cumulant method in $pp$ and $p+Pb$ collisions with the ATLAS detector at the CERN Large Hadron Collider, _Phys. Rev. C_. 97(2), 024904 (2018).  10.1103/PhysRevC.97.024904.
  * 14. K. Křížková Gajdošová, Probing QGP with flow: An experimental overview, _Nucl. Phys. A_. 1005, 121802 (2021).  10.1016/j.nuclphysa.2020.121802.
  * 15. G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, Derivation of transient relativistic fluid dynamics from the Boltzmann equation, _Phys. Rev._ D85, 114047 (2012).  10.1103/PhysRevD.85.114047, 10.1103/PhysRevD.91.039902.  [Erratum: Phys. Rev.D91,no.3,039902(2015)].
  * 16. G. S. Rocha, D. Wagner, G. S. Denicol, J. Noronha, and D. H. Rischke, Theories of Relativistic Dissipative Fluid Dynamics (11, 2023).
  * 17. L. D. Landau and E. M. Lifshitz, _Fluid Mechanics_. Pergamon (1987).  https://doi.org/10.1016/B978-0-08-033933-7.50001-5.
  * 18. J. L. Nagle and W. A. Zajc, Small System Collectivity in Relativistic Hadronic and Nuclear Collisions, _Ann. Rev. Nucl. Part. Sci._ 68, 211–235 (2018).  10.1146/annurev-nucl-101916-123209.
  * 19. B. Schenke, The smallest fluid on Earth, _Rept. Prog. Phys._ 84(8), 082301 (2021).  10.1088/1361-6633/ac14c9.
  * 20. M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Glauber modeling in high energy nuclear collisions, _Ann. Rev. Nucl. Part. Sci._ 57, 205–243 (2007).  10.1146/annurev.nucl.57.090506.123020.
  * 21. D. d’Enterria and C. Loizides, Progress in the Glauber Model at Collider Energies, _Ann. Rev. Nucl. Part. Sci._ 71, 315–344 (2021).  10.1146/annurev-nucl-102419-060007.
  * 22. P. Bozek, Collective flow in p-Pb and d-Pd collisions at TeV energies, _Phys. Rev. C_. 85, 014911 (2012).  10.1103/PhysRevC.85.014911.
  * 23. P. Bozek and W. Broniowski, Correlations from hydrodynamic flow in p-Pb collisions, _Phys. Lett. B_. 718, 1557–1561 (2013).  10.1016/j.physletb.2012.12.051.
  * 24. P. Bozek and W. Broniowski, Size of the emission source and collectivity in ultra-relativistic p-Pb collisions, _Phys. Lett. B_. 720, 250–253 (2013).  10.1016/j.physletb.2013.02.014.
  * 25. P. Bozek and W. Broniowski, Collective dynamics in high-energy proton-nucleus collisions, _Phys. Rev. C_. 88(1), 014903 (2013).  10.1103/PhysRevC.88.014903.
  * 26. P. Bozek, W. Broniowski, and G. Torrieri, Mass hierarchy in identified particle distributions in proton-lead collisions, _Phys. Rev. Lett._ 111, 172303 (2013).  10.1103/PhysRevLett.111.172303.
  * 27. A. Bzdak, B. Schenke, P. Tribedy, and R. Venugopalan, Initial state geometry and the role of hydrodynamics in proton-proton, proton-nucleus and deuteron-nucleus collisions, _Phys. Rev. C_. 87(6), 064906 (2013).  10.1103/PhysRevC.87.064906.
  * 28. G.-Y. Qin and B. Müller, Elliptic and triangular flow anisotropy in deuteron-gold collisions at $\sqrt{s_{NN}}=200$ GeV at RHIC and in proton-lead collisions at $\sqrt{s_{NN}}=5.02$ TeV at the LHC, _Phys. Rev. C_. 89(4), 044902 (2014).  10.1103/PhysRevC.89.044902.
  * 29. K. Werner, M. Bleicher, B. Guiot, I. Karpenko, and T. Pierog, Evidence for Flow from Hydrodynamic Simulations of $p$-Pb Collisions at 5.02 TeV from $\nu_{2}$ Mass Splitting, _Phys. Rev. Lett._ 112(23), 232301 (2014).  10.1103/PhysRevLett.112.232301.
  * 30. I. Kozlov, M. Luzum, G. Denicol, S. Jeon, and C. Gale, Transverse momentum structure of pair correlations as a signature of collective behavior in small collision systems (5, 2014).
  * 31. P. Romatschke, Light-Heavy Ion Collisions: A window into pre-equilibrium QCD dynamics?, _Eur. Phys. J. C_. 75(7), 305 (2015).  10.1140/epjc/s10052-015-3509-3.
  * 32. C. Shen, J.-F. Paquet, G. S. Denicol, S. Jeon, and C. Gale, Collectivity and electromagnetic radiation in small systems, _Phys. Rev. C_. 95(1), 014906 (2017).  10.1103/PhysRevC.95.014906.
  * 33. R. D. Weller and P. Romatschke, One fluid to rule them all: viscous hydrodynamic description of event-by-event central p+p, p+Pb and Pb+Pb collisions at $\sqrt{s}=5.02$ TeV, _Phys. Lett. B_. 774, 351–356 (2017).  10.1016/j.physletb.2017.09.077.
  * 34. G. Giacalone, J. Noronha-Hostler, and J.-Y. Ollitrault, Relative flow fluctuations as a probe of initial state fluctuations, _Phys. Rev. C_. 95(5), 054910 (2017).  10.1103/PhysRevC.95.054910.
  * 35. M. D. Sievert and J. Noronha-Hostler, CERN Large Hadron Collider system size scan predictions for PbPb, XeXe, ArAr, and OO with relativistic hydrodynamics, _Phys. Rev. C_. 100(2), 024904 (2019).  10.1103/PhysRevC.100.024904.
  * 36. N. Summerfield, B.-N. Lu, C. Plumberg, D. Lee, J. Noronha-Hostler, and A. Timmins, 16O 16O collisions at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider comparing $\alpha$ clustering versus substructure, _Phys. Rev. C_. 104(4), L041901 (2021).  10.1103/PhysRevC.104.L041901.
  * 37. R. Katz, C. A. G. Prado, J. Noronha-Hostler, and A. A. P. Suaide, System-size scan of $D$ meson $R_{AA}$ and $v_{n}$ using PbPb, XeXe, ArAr, and OO collisions at energies available at the CERN Large Hadron Collider, _Phys. Rev. C_. 102(4), 041901 (2020).  10.1103/PhysRevC.102.041901.
  * 38. J. S. Moreland, J. E. Bernhard, and S. A. Bass, Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions, _Phys. Rev. C_. 92(1), 011901 (2015).  10.1103/PhysRevC.92.011901.
  * 39. B. Schenke and R. Venugopalan, Eccentric protons? Sensitivity of flow to system size and shape in p+p, p+Pb and Pb+Pb collisions, _Phys. Rev. Lett._ 113, 102301 (2014).  10.1103/PhysRevLett.113.102301.
  * 40. J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Applying Bayesian parameter estimation to relativistic heavy-ion collisions: simultaneous characterization of the initial state and quark-gluon plasma medium, _Phys. Rev. C_. 94(2), 024907 (2016).  10.1103/PhysRevC.94.024907.
  * 41. D. Everett et al., Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging, _Phys. Rev. Lett._ 126(24), 242301 (2021).  10.1103/PhysRevLett.126.242301.
  * 42. D. Everett et al., Multisystem Bayesian constraints on the transport coefficients of QCD matter, _Phys. Rev. C_. 103(5), 054904 (2021).  10.1103/PhysRevC.103.054904.
  * 43. L. D. McLerran and R. Venugopalan, Computing quark and gluon distribution functions for very large nuclei, _Phys. Rev. D_. 49, 2233–2241 (1994).  10.1103/PhysRevD.49.2233.
  * 44. L. D. McLerran and R. Venugopalan, Gluon distribution functions for very large nuclei at small transverse momentum, _Phys. Rev. D_. 49, 3352–3355 (1994).  10.1103/PhysRevD.49.3352.
  * 45. A. Kovner, L. D. McLerran, and H. Weigert, Gluon production from nonAbelian Weizsacker-Williams fields in nucleus-nucleus collisions, _Phys. Rev. D_. 52, 6231–6237 (1995).  10.1103/PhysRevD.52.6231.
  * 46. E. Iancu and R. Venugopalan, _The Color glass condensate and high-energy scattering in QCD_ , In eds. R. C. Hwa and X.-N. Wang, _Quark-gluon plasma 4_ , pp. 249–3363 (3.  , 2003).  10.1142/9789812795533_0005.
  * 47. T. Lappi, Energy density of the glasma, _Phys. Lett. B_. 643, 11–16 (2006).  10.1016/j.physletb.2006.10.017.
  * 48. P. Romatschke and J. D. Hogg, Pre-Equilibrium Radial Flow from Central Shock-Wave Collisions in AdS5, _JHEP_. 04, 048 (2013).  10.1007/JHEP04(2013)048.
  * 49. J. S. Moreland, J. E. Bernhard, and S. A. Bass, Bayesian calibration of a hybrid nuclear collision model using p-Pb and Pb-Pb data at energies available at the CERN Large Hadron Collider, _Phys. Rev. C_. 101(2), 024911 (2020).  10.1103/PhysRevC.101.024911.
  * 50. H. Mäntysaari and B. Schenke, Evidence of strong proton shape fluctuations from incoherent diffraction, _Phys. Rev. Lett._ 117(5), 052301 (2016).  10.1103/PhysRevLett.117.052301.
  * 51. H. Mäntysaari and B. Schenke, Revealing proton shape fluctuations with incoherent diffraction at high energy, _Phys. Rev. D_. 94(3), 034042 (2016).  10.1103/PhysRevD.94.034042.
  * 52. H. Mäntysaari, B. Schenke, C. Shen, and W. Zhao, Bayesian inference of the fluctuating proton shape, _Phys. Lett. B_. 833, 137348 (2022).  10.1016/j.physletb.2022.137348.
  * 53. H. Mäntysaari, B. Schenke, C. Shen, and W. Zhao, Multiscale Imaging of Nuclear Deformation at the Electron-Ion Collider, _Phys. Rev. Lett._ 131(6), 062301 (2023).  10.1103/PhysRevLett.131.062301.
  * 54. S. Chekanov et al., Measurement of proton dissociative diffractive photoproduction of vector mesons at large momentum transfer at HERA, _Eur. Phys. J. C_. 26, 389–409 (2003).  10.1140/epjc/s2002-01079-0.
  * 55. A. Aktas et al., Diffractive photoproduction of $J/\psi$ mesons with large momentum transfer at HERA, _Phys. Lett. B_. 568, 205–218 (2003).  10.1016/j.physletb.2003.06.056.
  * 56. A. Aktas et al., Elastic J/psi production at HERA, _Eur. Phys. J. C_. 46, 585–603 (2006).  10.1140/epjc/s2006-02519-5.
  * 57. S. Chekanov et al., Exclusive photoproduction of J / psi mesons at HERA, _Eur. Phys. J. C_. 24, 345–360 (2002).  10.1007/s10052-002-0953-7.
  * 58. C. Alexa et al., Elastic and Proton-Dissociative Photoproduction of J/psi Mesons at HERA, _Eur. Phys. J. C_. 73(6), 2466 (2013).  10.1140/epjc/s10052-013-2466-y.
  * 59. H. Mäntysaari, B. Schenke, C. Shen, and P. Tribedy, Imprints of fluctuating proton shapes on flow in proton-lead collisions at the LHC, _Phys. Lett. B_. 772, 681–686 (2017).  10.1016/j.physletb.2017.07.038.
  * 60. B. Schenke, C. Shen, and P. Tribedy, Hybrid Color Glass Condensate and hydrodynamic description of the Relativistic Heavy Ion Collider small system scan, _Phys. Lett. B_. 803, 135322 (2020).  10.1016/j.physletb.2020.135322.
  * 61. B. Schenke, C. Shen, and P. Tribedy, Running the gamut of high energy nuclear collisions, _Phys. Rev. C_. 102(4), 044905 (2020).  10.1103/PhysRevC.102.044905.
  * 62. S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, Full result for the QCD equation of state with 2+1 flavors, _Phys. Lett. B_. 730, 99–104 (2014).  10.1016/j.physletb.2014.01.007.
  * 63. A. Bazavov et al., Equation of state in ( 2+1 )-flavor QCD, _Phys. Rev. D_. 90, 094503 (2014).  10.1103/PhysRevD.90.094503.
  * 64. T. Nunes da Silva, D. Chinellato, M. Hippert, W. Serenone, J. Takahashi, G. S. Denicol, M. Luzum, and J. Noronha, Pre-hydrodynamic evolution and its signatures in final-state heavy-ion observables, _Phys. Rev. C_. 103, 054906 (2021).  10.1103/PhysRevC.103.054906.
  * 65. T. N. da Silva, D. D. Chinellato, A. V. Giannini, M. N. Ferreira, G. S. Denicol, M. Hippert, M. Luzum, J. Noronha, and J. Takahashi, Prehydrodynamic evolution in large and small systems, _Phys. Rev. C_. 107(4), 044901 (2023).  10.1103/PhysRevC.107.044901.
  * 66. B. Schenke, S. Schlichting, and R. Venugopalan, Azimuthal anisotropies in p$+$Pb collisions from classical Yang–Mills dynamics, _Phys. Lett. B_. 747, 76–82 (2015).  10.1016/j.physletb.2015.05.051.
  * 67. T. Lappi, B. Schenke, S. Schlichting, and R. Venugopalan, Tracing the origin of azimuthal gluon correlations in the color glass condensate, _JHEP_. 01, 061 (2016).  10.1007/JHEP01(2016)061.
  * 68. G. Giacalone, B. Schenke, and C. Shen, Observable signatures of initial state momentum anisotropies in nuclear collisions, _Phys. Rev. Lett._ 125(19), 192301 (2020).  10.1103/PhysRevLett.125.192301.
  * 69. B. Schenke, S. Schlichting, and P. Singh, Rapidity dependence of initial state geometry and momentum correlations in p+Pb collisions, _Phys. Rev. D_. 105(9), 094023 (2022).  10.1103/PhysRevD.105.094023.
  * 70. J. Berges, B. Schenke, S. Schlichting, and R. Venugopalan, Turbulent thermalization process in high-energy heavy-ion collisions, _Nucl. Phys. A_. 931, 348–353 (2014).  10.1016/j.nuclphysa.2014.08.103.
  * 71. A. Ipp and D. I. Müller, Progress on 3+1D Glasma simulations, _Eur. Phys. J. A_. 56(9), 243 (2020).  10.1140/epja/s10050-020-00241-6.
  * 72. S. McDonald, S. Jeon, and C. Gale, 3+1D initialization and evolution of the glasma, _Phys. Rev. C_. 108(6), 064910 (2023).  10.1103/PhysRevC.108.064910.
  * 73. T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey, and Y. Nara, Hadronic dissipative effects on elliptic flow in ultrarelativistic heavy-ion collisions, _Phys. Lett. B_. 636, 299–304 (2006).  10.1016/j.physletb.2006.03.060.
  * 74. P. Bozek, W. Broniowski, and J. Moreira, Torqued fireballs in relativistic heavy-ion collisions, _Phys. Rev. C_. 83, 034911 (2011).  10.1103/PhysRevC.83.034911.
  * 75. W. Ke, J. S. Moreland, J. E. Bernhard, and S. A. Bass, Constraints on rapidity-dependent initial conditions from charged particle pseudorapidity densities and two-particle correlations, _Phys. Rev. C_. 96(4), 044912 (2017).  10.1103/PhysRevC.96.044912.
  * 76. M. Barej, A. Bzdak, and P. Gutowski, Wounded nucleon, quark, and quark-diquark emission functions versus experimental results from the BNL Relativistic Heavy Ion Collider at $\sqrt{s_{NN}}$ =200 GeV, _Phys. Rev. C_. 100(6), 064902 (2019).  10.1103/PhysRevC.100.064902.
  * 77. X.-Y. Wu, G.-Y. Qin, L.-G. Pang, and X.-N. Wang, (3+1)-D viscous hydrodynamics at finite net baryon density: Identified particle spectra, anisotropic flows, and flow fluctuations across energies relevant to the beam-energy scan at RHIC, _Phys. Rev. C_. 105(3), 034909 (2022).  10.1103/PhysRevC.105.034909.
  * 78. L. Du, C. Shen, S. Jeon, and C. Gale, Probing initial baryon stopping and equation of state with rapidity-dependent directed flow of identified particles, _Phys. Rev. C_. 108(4), L041901 (2023).  10.1103/PhysRevC.108.L041901.
  * 79. Z.-F. Jiang, X.-Y. Wu, S. Cao, and B.-W. Zhang, Hyperon polarization and its relation with directed flow in high-energy nuclear collisions, _Phys. Rev. C_. 108(6), 064904 (2023).  10.1103/PhysRevC.108.064904.
  * 80. C. Shen and S. Alzhrani, Collision-geometry-based 3D initial condition for relativistic heavy-ion collisions, _Phys. Rev. C_. 102(1), 014909 (2020).  10.1103/PhysRevC.102.014909.
  * 81. D. Soeder, W. Ke, J. F. Paquet, and S. A. Bass, Bayesian parameter estimation with a new three-dimensional initial-conditions model for ultrarelativistic heavy-ion collisions (6, 2023).
  * 82. Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, A Multi-phase transport model for relativistic heavy ion collisions, _Phys. Rev. C_. 72, 064901 (2005).  10.1103/PhysRevC.72.064901.
  * 83. P. Bozek and W. Broniowski, Collective flow in ultrarelativistic 3He-Au collisions, _Phys. Lett. B_. 739, 308–312 (2014).  10.1016/j.physletb.2014.11.006.
  * 84. C. Shen and B. Schenke, Longitudinal dynamics and particle production in relativistic nuclear collisions, _Phys. Rev. C_. 105(6), 064905 (2022).  10.1103/PhysRevC.105.064905.
  * 85. W. Zhao, S. Ryu, C. Shen, and B. Schenke, 3D structure of anisotropic flow in small collision systems at energies available at the BNL Relativistic Heavy Ion Collider, _Phys. Rev. C_. 107(1), 014904 (2023).  10.1103/PhysRevC.107.014904.
  * 86. P. Carzon, M. Martinez, M. D. Sievert, D. E. Wertepny, and J. Noronha-Hostler, Monte Carlo event generator for initial conditions of conserved charges in nuclear geometry, _Phys. Rev. C_. 105(3), 034908 (2022).  10.1103/PhysRevC.105.034908.
  * 87. L. Du, C. Shen, S. Jeon, and C. Gale.  Constraints on initial baryon stopping and equation of state from directed flow.  In _30th International Conference on Ultrarelativstic Nucleus-Nucleus Collisions_ (12, 2023).
  * 88. G. Pihan, A. Monnai, B. Schenke, and C. Shen.  Tracing baryon and electric charge transport in isobar collisions (12.  , 2023).
  * 89. P. Arnold, P. Romatschke, and W. van der Schee, Absence of a local rest frame in far from equilibrium quantum matter, _JHEP_. 10, 110 (2014).  10.1007/JHEP10(2014)110.
  * 90. R. Rougemont, J. Noronha, W. Barreto, G. S. Denicol, and T. Dore, Violation of energy conditions and entropy production in holographic Bjorken flow, _Phys. Rev. D_. 104(12), 126012 (2021).  10.1103/PhysRevD.104.126012.
  * 91. R. Rougemont, W. Barreto, and J. Noronha, Hydrodynamization times of a holographic fluid far from equilibrium, _Phys. Rev. D_. 105(4), 046009 (2022).  10.1103/PhysRevD.105.046009.
  * 92. H. Niemi and G. S. Denicol, How large is the Knudsen number reached in fluid dynamical simulations of ultrarelativistic heavy ion collisions? (4, 2014).
  * 93. J. Noronha-Hostler, J. Noronha, and M. Gyulassy, Sensitivity of flow harmonics to subnucleon scale fluctuations in heavy ion collisions, _Phys. Rev. C_. 93(2), 024909 (2016).  10.1103/PhysRevC.93.024909.
  * 94. G. Inghirami and H. Elfner, The applicability of hydrodynamics in heavy ion collisions at $\sqrt{s_{\mathrm{NN}}}$ = 2.4–7.7 GeV, _Eur. Phys. J. C_. 82(9), 796 (2022).  10.1140/epjc/s10052-022-10718-x.
  * 95. W. Broniowski, W. Florkowski, M. Chojnacki, and A. Kisiel, Free-streaming approximation in early dynamics of relativistic heavy-ion collisions, _Phys. Rev. C_. 80, 034902 (2009).  10.1103/PhysRevC.80.034902.
  * 96. J. Liu, C. Shen, and U. Heinz, Pre-equilibrium evolution effects on heavy-ion collision observables, _Phys. Rev. C_. 91(6), 064906 (2015).  10.1103/PhysRevC.91.064906.  [Erratum: Phys.Rev.C 92, 049904 (2015)].
  * 97. A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, and D. Teaney, Effective kinetic description of event-by-event pre-equilibrium dynamics in high-energy heavy-ion collisions, _Phys. Rev. C_. 99(3), 034910 (2019).  10.1103/PhysRevC.99.034910.
  * 98. A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, and D. Teaney, Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions to Hydrodynamics with QCD Kinetic Theory, _Phys. Rev. Lett._ 122(12), 122302 (2019).  10.1103/PhysRevLett.122.122302.
  * 99. P. B. Arnold, G. D. Moore, and L. G. Yaffe, Effective kinetic theory for high temperature gauge theories, _JHEP_. 01, 030 (2003).  10.1088/1126-6708/2003/01/030.
  * 100. D. Liyanage, D. Everett, C. Chattopadhyay, and U. Heinz, Prehydrodynamic evolution and its impact on quark-gluon plasma signatures, _Phys. Rev. C_. 105(6), 064908 (2022).  10.1103/PhysRevC.105.064908.
  * 101. S. Jaiswal, C. Chattopadhyay, L. Du, U. Heinz, and S. Pal, Nonconformal kinetic theory and hydrodynamics for Bjorken flow, _Phys. Rev. C_. 105(2), 024911 (2022).  10.1103/PhysRevC.105.024911.
  * 102. Y. Kanakubo, Y. Tachibana, and T. Hirano, Unified description of hadron yield ratios from dynamical core-corona initialization, _Phys. Rev. C_. 101(2), 024912 (2020).  10.1103/PhysRevC.101.024912.
  * 103. K. Werner, Core-corona procedure and microcanonical hadronization to understand strangeness enhancement in proton-proton and heavy ion collisions in the EPOS4 framework (6, 2023).
  * 104. M. P. Heller and M. Spalinski, Hydrodynamics Beyond the Gradient Expansion: Resurgence and Resummation, _Phys. Rev. Lett._ 115(7), 072501 (2015).  10.1103/PhysRevLett.115.072501.
  * 105. M. Strickland, J. Noronha, and G. Denicol, Anisotropic nonequilibrium hydrodynamic attractor, _Phys. Rev. D_. 97(3), 036020 (2018).  10.1103/PhysRevD.97.036020.
  * 106. G. Giacalone, A. Mazeliauskas, and S. Schlichting, Hydrodynamic attractors, initial state energy and particle production in relativistic nuclear collisions, _Phys. Rev. Lett._ 123(26), 262301 (2019).  10.1103/PhysRevLett.123.262301.
  * 107. J. Jankowski and M. Spaliński, Hydrodynamic attractors in ultrarelativistic nuclear collisions, _Prog. Part. Nucl. Phys._ 132, 104048 (2023).  10.1016/j.ppnp.2023.104048.
  * 108. A. Monnai, B. Schenke, and C. Shen, QCD Equation of State at Finite Chemical Potentials for Relativistic Nuclear Collisions, _Int. J. Mod. Phys. A_. 36(07), 2130007 (2021).  10.1142/S0217751X21300076.
  * 109. C. Eckart, The thermodynamics of irreversible processes III. Relativistic theory of the simple fluid, _Physical Review_. 58, 919–924 (1940).
  * 110. G. Pichon, Étude relativiste de fluides visqueux et chargés, _Annales de l’I.H.P. Physique théorique_. 2(1), 21–85 (1965).  URL http://www.numdam.org/item/AIHPA_1965__2_1_21_0.
  * 111. W. A. Hiscock and L. Lindblom, Generic instabilities in first-order dissipative fluid theories, _Phys. Rev. D_. 31(4), 725–733 (1985).
  * 112. F. S. Bemfica, M. M. Disconzi, and J. Noronha, First-Order General-Relativistic Viscous Fluid Dynamics, _Phys. Rev. X_. 12(2), 021044 (2022).  10.1103/PhysRevX.12.021044.
  * 113. L. Gavassino, Can We Make Sense of Dissipation without Causality?, _Phys. Rev. X_. 12(4), 041001 (2022).  10.1103/PhysRevX.12.041001.
  * 114. F. S. Bemfica, M. M. Disconzi, and J. Noronha, Causality and existence of solutions of relativistic viscous fluid dynamics with gravity, _Phys. Rev._ D98(10), 104064 (26 pages) (2018).  10.1103/PhysRevD.98.104064.
  * 115. P. Kovtun, First-order relativistic hydrodynamics is stable (2019).
  * 116. F. S. Bemfica, F. S. Bemfica, M. M. Disconzi, M. M. Disconzi, J. Noronha, and J. Noronha, Nonlinear Causality of General First-Order Relativistic Viscous Hydrodynamics, _Phys. Rev. D_. 100(10), 104020 (2019).  10.1103/PhysRevD.100.104020.  [Erratum: Phys.Rev.D 105, 069902 (2022)].
  * 117. R. E. Hoult and P. Kovtun, Stable and causal relativistic Navier-Stokes equations, _JHEP_. 06, 067 (2020).  10.1007/JHEP06(2020)067.
  * 118. K. Kovtun, Lectures on hydrodynamic fluctuations in relativistic theories, _J. Phys._ A45, 473001 (2012).  10.1088/1751-8113/45/47/473001.
  * 119. A. Pandya and F. Pretorius, Numerical exploration of first-order relativistic hydrodynamics, _Phys. Rev. D_. 104(2), 023015 (2021).  10.1103/PhysRevD.104.023015.
  * 120. A. Pandya, E. R. Most, and F. Pretorius, Conservative finite volume scheme for first-order viscous relativistic hydrodynamics, _Phys. Rev. D_. 105(12), 123001 (2022).  10.1103/PhysRevD.105.123001.
  * 121. A. Pandya, E. R. Most, and F. Pretorius, Causal, stable first-order viscous relativistic hydrodynamics with ideal gas microphysics, _Phys. Rev. D_. 106(12), 123036 (2022).  10.1103/PhysRevD.106.123036.
  * 122. H. Bantilan, Y. Bea, and P. Figueras, Evolutions in first-order viscous hydrodynamics, _JHEP_. 08, 298 (2022).  10.1007/JHEP08(2022)298.
  * 123. G. S. Rocha, G. S. Denicol, and J. Noronha, Perturbative approaches in relativistic kinetic theory and the emergence of first-order hydrodynamics, _Phys. Rev. D_. 106(3), 036010 (2022).  10.1103/PhysRevD.106.036010.
  * 124. W. Israel and J. M. Stewart, Transient relativistic thermodynamics and kinetic theory, _Ann. Phys._ 118, 341–372 (1979).
  * 125. R. Baier, P. Romatschke, D. T. Son, A. O. Starinets, and M. A. Stephanov, Relativistic viscous hydrodynamics, conformal invariance, and holography, _JHEP_. 04, 100 (2008).  10.1088/1126-6708/2008/04/100.
  * 126. M. Alqahtani, M. Nopoush, and M. Strickland, Relativistic anisotropic hydrodynamics, _Prog. Part. Nucl. Phys._ 101, 204–248 (2018).  10.1016/j.ppnp.2018.05.004.
  * 127. M. Luzum and P. Romatschke, Conformal relativistic viscous hydrodynamics: Applications to rhic results at s(nn)**(1/2) = 200-gev, _Phys. Rev._ C78, 034915 (2008).  10.1103/PhysRevC.78.034915, 10.1103/PhysRevC.79.039903.  [Erratum: Phys. Rev.C79,039903(2009)].
  * 128. B. Schenke, S. Jeon, and C. Gale, (3+1)D hydrodynamic simulation of relativistic heavy-ion collisions, _Phys. Rev. C_. 82, 014903 (2010).  10.1103/PhysRevC.82.014903.
  * 129. B. Schenke, S. Jeon, and C. Gale, Elliptic and triangular flow in event-by-event (3+1)d viscous hydrodynamics, _Phys. Rev. Lett._ 106, 042301 (2011).  10.1103/PhysRevLett.106.042301.
  * 130. J. Noronha-Hostler, G. S. Denicol, J. Noronha, R. P. G. Andrade, and F. Grassi, Bulk viscosity effects in event-by-event relativistic hydrodynamics, _Phys. Rev._ C88(4), 044916 (2013).  10.1103/PhysRevC.88.044916.
  * 131. J. Noronha-Hostler, J. Noronha, and F. Grassi, Bulk viscosity-driven suppression of shear viscosity effects on the flow harmonics at energies available at the bnl relativistic heavy ion collider, _Phys. Rev._ C90(3), 034907 (2014).  10.1103/PhysRevC.90.034907.
  * 132. H. Marrochio, J. Noronha, G. S. Denicol, M. Luzum, S. Jeon, and C. Gale, Solutions of conformal israel-stewart relativistic viscous fluid dynamics, _Phys. Rev._ C91(1), 014903 (2015).  10.1103/PhysRevC.91.014903.
  * 133. C. E. Brennen, _Cavitation and Bubble Dynamics_. Cambridge University Press (2013).  10.1017/CBO9781107338760.
  * 134. G. Torrieri, B. Tomasik, and I. Mishustin, Bulk viscosity driven clusterization of quark-gluon plasma and early freeze-out in relativistic heavy-ion collisions, _Phys. Rev._ C77, 034903 (2008).  10.1103/PhysRevC.77.034903.
  * 135. K. Rajagopal and N. Tripuraneni, Bulk viscosity and cavitation in boost-invariant hydrodynamic expansion, _JHEP_. 03, 018 (2010).  10.1007/JHEP03(2010)018.
  * 136. M. Habich and P. Romatschke, Onset of cavitation in the quark-gluon plasma, _JHEP_. 12, 054 (2014).  10.1007/JHEP12(2014)054.
  * 137. G. S. Denicol, C. Gale, and S. Jeon, The domain of validity of fluid dynamics and the onset of cavitation in ultrarelativistic heavy ion collisions, _PoS_. CPOD2014, 033 (2015).
  * 138. M. Byres, S. H. Lim, C. McGinn, J. Ouellette, and J. L. Nagle, Bulk viscosity and cavitation in heavy ion collisions, _Phys. Rev. C_. 101(4), 044902 (2020).  10.1103/PhysRevC.101.044902.
  * 139. F. S. Bemfica, M. M. Disconzi, and J. Noronha, Causality of the Einstein-Israel-Stewart Theory with Bulk Viscosity, _Phys. Rev. Lett._ 122(22), 221602 (2019).  10.1103/PhysRevLett.122.221602.
  * 140. L. Gavassino, M. M. Disconzi, and J. Noronha, Dispersion relations alone cannot guarantee causality (7, 2023).
  * 141. D.-L. Wang and S. Pu, Stability and causality criteria in linear mode analysis: stability means causality (9, 2023).
  * 142. R. E. Hoult and P. Kovtun, Causality and classical dispersion relations (9, 2023).
  * 143. W. A. Hiscock and L. Lindblom, Stability and causality in dissipative relativistic fluids, _Annals of Physics_. 151(2), 466–496 (1983).
  * 144. T. S. Olson, Stability and causality in the Israel-Stewart energy frame theory, _Annals Phys._ 199, 18 (1990).  10.1016/0003-4916(90)90366-V.
  * 145. G. S. Denicol, T. Kodama, T. Koide, and P. Mota, Stability and causality in relativistic dissipative hydrodynamics, _J. Phys._ G35, 115102 (2008).  10.1088/0954-3899/35/11/115102.
  * 146. S. Pu, T. Koide, and D. H. Rischke, Does stability of relativistic dissipative fluid dynamics imply causality?, _Phys. Rev._ D81, 114039 (2010).  10.1103/PhysRevD.81.114039.
  * 147. L. Gavassino, Applying the Gibbs stability criterion to relativistic hydrodynamics, _Class. Quant. Grav._ 38(21), 21LT02 (2021).  10.1088/1361-6382/ac2b0e.
  * 148. L. Gavassino, M. Antonelli, and B. Haskell, Thermodynamic Stability Implies Causality, _Phys. Rev. Lett._ 128(1), 010606 (2022).  10.1103/PhysRevLett.128.010606.
  * 149. N. Abboud, E. Speranza, and J. Noronha, Causal and stable first-order chiral hydrodynamics (8, 2023).
  * 150. F. S. Bemfica, M. M. Disconzi, V. Hoang, J. Noronha, and M. Radosz, Nonlinear Constraints on Relativistic Fluids Far from Equilibrium, _Phys. Rev. Lett._ 126(22), 222301 (2021).  10.1103/PhysRevLett.126.222301.
  * 151. C. Plumberg, D. Almaalol, T. Dore, J. Noronha, and J. Noronha-Hostler, Causality violations in realistic simulations of heavy-ion collisions, _Phys. Rev. C_. 105(6), L061901 (2022).  10.1103/PhysRevC.105.L061901.
  * 152. C. Chiu and C. Shen, Exploring theoretical uncertainties in the hydrodynamic description of relativistic heavy-ion collisions, _Phys. Rev. C_. 103(6), 064901 (2021).  10.1103/PhysRevC.103.064901.
  * 153. R. Krupczak et al., Causality violations in simulations of large and small heavy-ion collisions (11, 2023).
  * 154. C. Shen, Studying QGP with flow: A theory overview, _Nucl. Phys. A_. 1005, 121788 (2021).  10.1016/j.nuclphysa.2020.121788.
  * 155. M. Strickland, Anisotropic hydrodynamics: Three lectures, _Acta Phys. Polon._ B45(12), 2355–2394 (2014).  10.5506/APhysPolB.45.2355.
  * 156. M. Strickland, Small system studies: A theory overview, _Nucl. Phys. A_. 982, 92–98 (2019).  10.1016/j.nuclphysa.2018.09.071.
  * 157. G. S. Denicol, S. Jeon, and C. Gale, Transport Coefficients of Bulk Viscous Pressure in the 14-moment approximation, _Phys. Rev. C_. 90(2), 024912 (2014).  10.1103/PhysRevC.90.024912.
  * 158. G. S. Rocha, G. S. Denicol, and J. Noronha, Novel Relaxation Time Approximation to the Relativistic Boltzmann Equation, _Phys. Rev. Lett._ 127(4), 042301 (2021).  10.1103/PhysRevLett.127.042301.
  * 159. G. S. Rocha, M. N. Ferreira, G. S. Denicol, and J. Noronha, Transport coefficients of quasiparticle models within a new relaxation time approximation of the Boltzmann equation, _Phys. Rev. D_. 106(3), 036022 (2022).  10.1103/PhysRevD.106.036022.
  * 160. C. V. P. de Brito and G. S. Denicol, Third-order relativistic dissipative fluid dynamics from the method of moments, _Phys. Rev. D_. 108(9), 096020 (2023).  10.1103/PhysRevD.108.096020.
  * 161. G. Nijs and W. van der Schee, Predictions and postdictions for relativistic lead and oxygen collisions with the computational simulation code Trajectum, _Phys. Rev. C_. 106(4), 044903 (2022).  10.1103/PhysRevC.106.044903.
  * 162. D. Everett et al., Multisystem Bayesian constraints on the transport coefficients of QCD matter, _Phys. Rev. C_. 103(5), 054904 (2021).  10.1103/PhysRevC.103.054904.
  * 163. P. Romatschke, Relativistic Fluid Dynamics Far From Local Equilibrium, _Phys. Rev. Lett._ 120(1), 012301 (2018).  10.1103/PhysRevLett.120.012301.
  * 164. M. P. Heller, R. A. Janik, and P. Witaszczyk, Hydrodynamic gradient expansion in gauge theory plasmas, _Phys. Rev. Lett._ 110(21), 211602 (2013).  10.1103/PhysRevLett.110.211602.
  * 165. A. Buchel, M. P. Heller, and J. Noronha, Entropy production, hydrodynamics, and resurgence in the primordial quark-gluon plasma from holography, _Phys. Rev._ D94(10), 106011 (2016).  10.1103/PhysRevD.94.106011.
  * 166. G. S. Denicol and J. Noronha, Connecting far-from-equilibrium hydrodynamics to resummed transport coefficients and attractors, _Nucl. Phys. A_. 1005, 121748 (2021).  10.1016/j.nuclphysa.2020.121748.
  * 167. J.-P. Blaizot and L. Yan, Fluid dynamics of out of equilibrium boost invariant plasmas, _Phys. Lett. B_. 780, 283–286 (2018).  10.1016/j.physletb.2018.02.058.
  * 168. H. B. Callen and T. A. Welton, Irreversibility and generalized noise, _Phys. Rev._ 83, 34–40 (1951).  10.1103/PhysRev.83.34.
  * 169. R. Kubo, Statistical mechanical theory of irreversible processes. 1. General theory and simple applications in magnetic and conduction problems, _J. Phys. Soc. Jap._ 12, 570–586 (1957).  10.1143/JPSJ.12.570.
  * 170. C. Young, J. I. Kapusta, C. Gale, S. Jeon, and B. Schenke, Thermally Fluctuating Second-Order Viscous Hydrodynamics and Heavy-Ion Collisions, _Phys. Rev. C_. 91(4), 044901 (2015).  10.1103/PhysRevC.91.044901.
  * 171. M. Singh, C. Shen, S. McDonald, S. Jeon, and C. Gale, Hydrodynamic Fluctuations in Relativistic Heavy-Ion Collisions, _Nucl. Phys. A_. 982, 319–322 (2019).  10.1016/j.nuclphysa.2018.10.061.
  * 172. Y. Akamatsu, A. Mazeliauskas, and D. Teaney, A kinetic regime of hydrodynamic fluctuations and long time tails for a Bjorken expansion, _Phys. Rev. C_. 95(1), 014909 (2017).  10.1103/PhysRevC.95.014909.
  * 173. Y. Akamatsu, A. Mazeliauskas, and D. Teaney, Bulk viscosity from hydrodynamic fluctuations with relativistic hydrokinetic theory, _Phys. Rev. C_. 97(2), 024902 (2018).  10.1103/PhysRevC.97.024902.
  * 174. X. An, G. Basar, M. Stephanov, and H.-U. Yee, Relativistic Hydrodynamic Fluctuations, _Phys. Rev. C_. 100(2), 024910 (2019).  10.1103/PhysRevC.100.024910.
  * 175. M. Martinez and T. Schäfer, Stochastic hydrodynamics and long time tails of an expanding conformal charged fluid, _Phys. Rev. C_. 99(5), 054902 (2019).  10.1103/PhysRevC.99.054902.