# A Decade of Collectivity in Small Systems

JAN FIETE GROSSE-OETRINGHAUS  CERN, 1211 Geneva, Switzerland
<EMAIL_ADDRESS>URS ACHIM WIEDEMANN  CERN, 1211 Geneva, Switzerland
<EMAIL_ADDRESS>
###### Abstract

Signatures of collectivity, including azimuthally anisotropic and radial flow
as well as characteristic hadrochemical dependencies, have been observed since
long in (ultra)relativistic nucleus–nucleus collisions. They underpin the
interpretation of these collision systems in terms of QGP formation and close-
to-perfect fluidity. Remarkably, however, essentially all these signatures of
collectivity have been identified within the last decade in collision systems
as small as $\mathrm{pp}$ and $\mathrm{p}$–$\mathrm{Pb}$, where collective
phenomena had been assumed to be absent traditionally. Precursor phenomena may
have been found even in ep and $\mathrm{e^{+}e^{-}}$ collisions. This article
provides a complete review of all data on small system collectivity. It
reviews model simulations of these data where available. However, in the
absence of a phenomenologically fully satisfactory description of collectivity
across all system sizes, we focus in particular on the theoretical basis of
all dynamical frameworks of collectivity invoked in heavy ion collisions, and
their expected scaling with system size. Our discussion clarifies to what
extent all dynamical explanations are challenged by the available data.

###### Contents

  1. 1 Introduction
  2. 2 Measures of collectivity and medium-modifications
    1. 2.1 Measures of anisotropic flow
    2. 2.2 Hadrochemical measures of collectivity
    3. 2.3 Hard processes in small systems
  3. 3 Experimental Overview
    1. 3.1 Long-range correlations
    2. 3.2 Bulk properties
    3. 3.3 Hard probes
    4. 3.4 Towards even smaller systems…
  4. 4 Theory overview
    1. 4.1 Dynamical frameworks of collectivity
    2. 4.2 Hydrodynamization and Thermalization - insights from simple models
    3. 4.3 Model studies and data comparisons
    4. 4.4 Model comparisons to hadrochemical measures of collectivity
  5. 5 Concluding discussion and outlook

> In high-energy physics, we have concentrated on experiments in which we
> distribute a higher and higher amount of energy into a region with smaller
> and smaller dimensions. But, in order to study the question of ‘vacuum’, we
> must turn to a different direction; we should investigate bulk phenomena by
> distributing high energy over a relatively large volume.
>
> T.D. Lee, 1975

##  1 Introduction

The LHC discovery of collectivity in small systems relates measurements in
proton–proton, proton–nucleus and nucleus–nucleus collisions in a remarkably
smooth way without any step-wise discontinuity. In addition to an intense
experimental effort at the LHC, it has motivated the study of dedicated small
collision systems at RHIC, the re-analysis and re-interpretation of data in
$\mathrm{e^{+}e^{-}}$ and ep collisions, and extensive theoretical work. As
this review will emphasize, this discovery impacts both: our understanding of
the dynamics of nucleus–nucleus collisions (in which signatures of
collectivity have been searched for traditionally) and our understanding of
the dynamics of proton–proton collisions (in which signatures of collectivity
had been assumed to be absent). To put these LHC measurements into a broader
scientific context, we recall the following:

With the discovery of QCD in 19731, 2 and with first insights into the QCD
phase structure of the strong interactions at finite temperature and density
3, 4, 5, confinement and the chiral condensate were understood to be
characteristics of the QCD vacuum. The strong interactions at finite
temperature became calculable, see e.g. Ref. 6. Their physical origin became
testable by creating the conditions for the QCD high-temperature phase — the
quark–gluon plasma (QGP) — in laboratory-based experiments. As captured
succinctly in the above quote of T.D. Lee 7, this led since the 1980s to
research in heavy-ion physics (HIP) that took a different direction from the
main efforts of high-energy physics (HEP).

HEP focused on “hard” high-momentum transfer processes that are calculable
perturbatively with controlled precision 8. Aimed initially at establishing
QCD as the theory of the strong interactions, this precision frontier is
advanced since the 1990s to control as precisely as possible the QCD
background on top of which new physics is searched for in collider
experiments. Most relevant for this program are processes in which the
systematically calculable “hard” short-distance physics factorizes from the
“soft” long-distance contributions which are process-independent and thus
independent of the partonic environment in which they may occur. Models of
soft multi-particle production in HEP multi-purpose event generators were
initially based on extrapolating this picture of collinearly factorized QCD
beyond its regime of guaranteed validity. Modern model implementations in
PYTHIA9, HERWIG10 and SHERPA11 go beyond a purely incoherent superposition of
multiple parton interactions (e.g. by implementing color reconnections in
hadronization 12) but HEP multi-purpose event generators are not partonic
transport codes: the partonic fragmentation processes and essential parts of
the hadronization process are modeled as in vacuum, i.e., independent of the
phase-space density within which they occur.

HIP, in contrast, focuses on testing QCD at finite phase-space density by
searching for collective phenomena in soft multi-particle production, and by
characterizing the modification of hard processes embedded in dense QCD
matter. Throughout the 1980s and 90s, fixed target experiments with ultra-
relativistic nuclear beams at the BNL AGS and at CERN SPS searched for signals
of the QGP phase transition in terms of step-like structures that would occur
in multiple classes of measurements at a common critical energy density 13,
14. These experiments stated circumstantial evidence for QGP formation based
on their discovery of strong sensitivities of soft particle production to the
nuclear environment in the collision region 15. However, the original simple
search strategy (that was akin of the methodology used to characterize phase
transitions in non-evolving solid state physics) had failed. No step-like
changes were discovered.

In the early 2000s, due to major theoretical advances 16, 17, 18 and in
comparison with RHIC data 19, 20, 21, 22, it became clear that dissipative
fluid dynamics provides a successful phenomenology of soft multi-particle
production in heavy-ion collisions 23, 24. This resulted in a radically
different approach towards testing hot and dense QCD. Since dissipative fluid
dynamics is formulated in terms of thermodynamic quantities that are
calculable from first principles in QCD, the rapid (hydro)dynamic evolution
from extreme initial conditions to a diluted final hadronic state was not
viewed any more as a confounding factor in testing equilibrium QCD but it
became the very tool via which QCD equilibrium properties would reveal
themselves. The resulting phenomenology supports the formation of a QCD fluid
of close to minimal dissipative properties in which the mean free path of
quasi-particles is so short that the very concept of a medium composed of
propagating quasi-particles becomes questionable 25. It cannot be
overemphasized that this HIP perfect fluid picture of nucleus–nucleus
collisions is the extreme opposite of the HEP picture for proton–proton
collisions. In the HIP case, no degree of freedom propagates through a sizable
fraction of the collision region without interaction with the surrounding
medium; in the HEP case, all degrees of freedom free-stream and fragment
through the entire collision region without interacting with any medium.

These maximally different default pictures for soft multi-particle production
in proton–proton (HEP) and nucleus–nucleus (HIP) collisions raise a set of
common sense questions, including: Can both pictures be correct? Is there a
critical size or density at which a QCD system produced in hadronic collisions
transits from free-streaming to fluid-dynamic behavior? Can precursors of
fluid dynamics be identified and could they inform us about the QCD non-
equilibrium dynamics that leads to fast hydrodynamization and equilibration?
All these questions point to the importance of understanding how collectivity
and medium-modifications arise as a function of systems size111Varying the
spatial extension of a system is not the only way of changing conditions of
collectivity. For instance, in simple conformal models, the degree to which a
collision system develops collective phenomena is known to depend on one
dimensionless product of energy density and transverse extension, the so-
called opacity26, 27, 28, 29, 30. In practice, however, a change of beam
particle from Pb to p changes the transverse spatial extension by a factor 10
– to change energy density by a comparable factor is more difficult because of
its approximately logarithmic $\sqrt{s}$-dependence. It is this point that is
informed by the LHC discovery of collectivity in small systems. These
measurements marked a paradigm shift as, historically, the experimental
programs with nuclear beams focused on central collisions of the heaviest
nuclei to distribute high energy over the largest achievable volume and to
create conditions most favorable for the formation of a QGP. In this context,
the study of small collision systems like $\mathrm{p}$–$\mathrm{A}$ or pp was
seen as a reference for large systems. The aim was to constrain the nuclear
dependence of parton distribution functions and to study so-called cold
nuclear matter effects 31, i.e. the particle production expected in the
absence of QGP formation. This logic prevailed throughout the CERN SPS and
RHIC experimental programs and it was the original motivation for
supplementing the LHC heavy-ion program with $\mathrm{p}$–$\mathrm{A}$
collisions 32. Heavy-ion experiments therefore always had an important small
systems program but without the expectation of gaining significant additional
understanding about QCD physics at extreme temperature and density. However,
with the LHC discovery of long-range correlation structures in $\mathrm{pp}$
collisions in 201033 and in $\mathrm{p}$–$\mathrm{Pb}$ collisions in 201234,
35, 36, for the first time, effects well-known from the study of the QGP, and
not accounted for by standard multi-purpose HEP event generators 37, occurred
in high-multiplicity $\mathrm{pp}$ and $\mathrm{p}$–$\mathrm{Pb}$ collisions.
That the near-side ridge in $\mathrm{pp}$ collisions was only the first
observation of a completely generic phenomenon, namely small system
collectivity, became gradually clear with the observation of similar effects
in $\mathrm{p}$–$\mathrm{Pb}$ collisions and with the equally unexpected
discovery of a second ridge structure on the away side. A new research
direction emerged. This manuscript refers to about 200 experimental
publications related to small systems, and it reviews the theoretical research
motivated by these data. For earlier reviews, see Refs. 38, 39, 40, 41.

This review separates the discussion of definitions of measurements in chapter
2, from the discussion of experimental data in chapter 3 and from the
discussion of theoretical concepts informed by these data in chapter 4. To the
extent possible, the chapters are self-contained. In particular, chapter 2
provides a detailed and necessarily technical discussion of how measures of
collectivity are defined. But one may want to skip chapter 2 in a first
reading if one is mainly interested in the state of art of experiment or
theory.

##  2 Measures of collectivity and medium-modifications

In the broader physics debate, it has been emphasized repeatedly that “more is
different”42 in the sense that macroscopic volumes of matter manifest
qualitatively new phenomena that are distinct from those that can be discerned
in interactions among few elementary constituents. In a broad sense,
collectivity refers to such phenomena that emerge with increasing size and
complexity of a physical system. Emergent collective phenomena pose novel
challenges for physics, as they often require distinct experimental techniques
for their characterization, distinct theoretical methods for their description
and as there is the fundamental challenge of understanding how collectivity in
many-body systems emerges from the fundamental degrees of freedom.

In the phenomenological practice of collider physics, a measure of
“collectivity” is more than a measurement. It is a measurement supplemented by
a baseline that is free of collective phenomena and on top of which
“collectivity” can be established. Baselines are chosen to separate
microscopic physics (that is thought to arise from the incoherent
superposition of elementary interactions) from those qualitatively novel
“collective” phenomena that can arise only in sufficiently macroscopic
systems. However, the simple idea of separating microscopic from macroscopic
physics faces obvious challenges if applied to small mesoscopic systems that
interpolate smoothly between the limiting cases of the small and the big.

###  2.1 Measures of anisotropic flow

Hadronic wavefunctions have a finite transverse extent. When they collide,
energy is deposited in a finite region of the plane transverse to the beam
direction. The shape of this initial spatial deposition depends on the kind of
hadron or nucleus that is collided, on the impact parameter of the collision
and on event-by-event fluctuations. In an azimuthal Fourier decomposition,
this spatial information is typically characterized in terms of spatial
eccentricities $\epsilon_{n}$. However, what is measured are certain
correlations of final azimuthal asymmetries $v_{n}$ in momentum space. If the
quanta produced in different positions in the transverse plane would not
interact with each other, the measured $v_{n}$ would not be sensitive to the
$\epsilon_{n}$’s. One important aspect discussing collectivity is to ask: What
is the non-trivial collective dynamics that turns spatial gradients and
asymmetries $\epsilon_{n}$ into momentum-anisotropies $v_{n}$,

$\underbrace{\\{\epsilon_{1},\epsilon_{2},\epsilon_{3},\dots\\}}_{\rm
spatial\,asymmetries}\qquad\xrightarrow[\rm dynamics?]{\rm
which\,collective}\qquad\underbrace{\\{v_{1},v_{2},v_{3},\dots\\}}_{\rm
momentum\,asymmetries}$ (1)

Experimentally, the initial conditions
$\\{\epsilon_{1},\epsilon_{2},\epsilon_{3},\dots\\}$ in (1) can be changed by
choosing different nuclear or hadron beams, by varying the centrality of the
collision, or by selecting particular event classes via event shape
engineering 43. As reviewed in the following section 3, a wealth of data maps
out the dynamical response (1) as a function of system size. Here, we start by
reviewing different experimental definitions of the anisotropic flow
coefficients $v_{n}$.

If $P(\vec{p}_{1},\vec{p}_{2})$ is the probability to produce a pair of
particles with momenta $\vec{p}_{1}$ and $\vec{p}_{2}$ in a collision, and if
$P(\vec{p}_{i})$ is the corresponding probability to produce a single particle
with momentum $\vec{p}_{i}$, then
$P(\vec{p}_{1},\vec{p}_{2})\not=P(\vec{p}_{1})P(\vec{p}_{2})$ signals
correlated particle production. However, not all correlated production is
“collective”. Measures of anisotropic flow aim at separating few-particle
correlations (that originate from microscopic sources like jets or resonance
decays, also referred to as non-flow) from collective phenomena (that involve
the correlation amongst essentially all particles in an event). Both sources
of correlations scale differently with charged event multiplicity222For
reasons of readability, we denote charged-particle multiplicity with $N$ and
without the usual subscript “ch”. $N$. If the correlation function contains a
collective phenomenon that affects with signal strength $S_{\rm collective}$
each particle pair, then it is multiplied by the number of pairs
$\tfrac{N(N-1)}{2}$. This collective contribution dominates parametrically by
$O(N)$ compared to the microscopic production, since

$P(1,2)=\frac{2}{N(N-1)}\left(\frac{N(N-1)}{2}S_{\rm
collective}(N)+\frac{N}{2}S_{\rm microscopic}\right)\,.$ (2)

Here, $S_{\rm collective}$ is a function of $N$ since it depends also on
system size. In contrast, the signal strength $S_{\rm microscopic}$ of
microscopic correlations is by definition independent of the environment it is
embedded in and it is thus independent of $N$.

As a next step we discuss what one would like to measure and what one would
like to remove from the measurement: Ideally, one would like to characterize
the azimuthal distribution of particles into Fourier coefficients

$v_{\rm n}\equiv\langle\langle
e^{in\varphi}\rangle\rangle=\langle\int\text{d}\varphi
e^{in\varphi}\tfrac{dN}{\text{d}\varphi}\rangle\,,$ (3)

where the two brackets denote averages over all particles in an event and over
event samples. Then, one would like to remove from the measurement of $v_{\rm
n}$ all physics of microscopic origin, so that it becomes a measure of
collectivity. The problem here is two-fold. First, without reference to a
physically preferred azimuthal orientation (such as the orientation of the
reaction plane), the azimuthal integral will average to zero. Eq. 3 does not
correspond to an experimental measurement. To bypass this problem, one
constructs azimuthal multi-particle correlations that are independent of
overall azimuthal orientation. Second, any subtraction of microscopic physics
depends on a assumption of what this microscopic physics is and it is
therefore, strictly speaking, a choice. The notion of pseudorapidity gap plays
in this context a central role. Namely, the two main few-particle correlation
phenomena, jets and resonance decays, occur at small angular differences on
the near side, i.e. at small $\Delta\eta$ and $\Delta\varphi$, as we will see
later in the left panel of Fig. 2. Those can therefore be separated from long-
range contributions by considering only particle pairs with
$|\Delta\eta|>\eta_{\rm gap}$ ($\eta$-gap method). However, a jet has a back-
to-back component due to momentum conservation333In pp collisions, contrary to
elementary $\mathrm{e^{+}e^{-}}$ collisions, the parton–parton scattering
frame is not the lab frame and therefore the back-to-back side of the jet is
found at $\Delta\varphi=\pi$ but at different $\Delta\eta$. leading to a broad
away-side structure. This structure cannot be separated through an $\eta$-gap
and in consequence, a considerable amount of the jet correlation remains.

What one measures and what one removes from the measurement. In the following,
we summarize different procedures that are in use to quantify azimuthal
correlations and that are used to reduce few-particle correlations in these
measurements.

  1. 1.

Near-side ridge yields
Measurements of near-side ridge yields are arguably the simplest way of
bypassing the above-mentioned problem that jet-like correlations lead to a
broad away-side structure. As basic $\mathrm{e^{+}e^{-}}$ or pp inspired
calculations and Monte Carlo simulations do not have a mechanism which
produces near-side long-range correlation, the yield measurement carries a
strong physics message. Some uncertainties remain from the definition of the
baseline of the combinatorial background, and one may underestimate this yield
at low $p_{T}$ where the away-side peak is wide in azimuth. Also, from the
ridge yield alone one cannot extract $v_{\rm n}$ coefficients while the
converse is possible.

  2. 2.

Event-plane method
In this method 44, two different detectors are used. Typically, a forward
rapidity detector determines a preferred azimuthal orientation called event
plane $\Psi_{\rm EP}$ which in large collision systems is a proxy for the
orientation of the reaction plane. This angle is then used to extract Fourier
coefficients similar to Eq. 3 from particles measured in a mid-rapidity
detector:

$v_{\rm n}\equiv\langle\langle e^{in(\varphi-\Psi_{\rm
EP})}\rangle\rangle=\langle\int\text{d}\varphi e^{in(\varphi-\Psi_{\rm
EP})}\tfrac{dN}{\text{d}\varphi}\rangle\,.$ (4)

The rapidity separation of the two detectors allows to reduce non-flow
contributions. In small systems such as $\mathrm{p}$–$\mathrm{Pb}$ collisions,
this method has limitations, since there may not be a clear definition of the
reaction plane. For instance, if the system is not rapidity-invariant or if it
has small multiplicity, $\Psi_{\rm EP}$ can depend on rapidity. The precision
of the event plane depends on the number of considered particles and its
resolution is proportional to $\sqrt{N}$. If the event plane is misestimated,
the extracted signal (4) is then reduced as it is averaged out. In practice,
if the two used detector systems have a fixed rapidity gap between them by the
experiment’s design, an experimental study of the signal as function of the
rapidity separation is not possible and the associated uncertainties cannot be
experimentally assessed.

  3. 3.

Decomposing associated per-trigger yields into Fourier coefficients
Measuring the distance in azimuth ($\Delta\varphi$) and pseudorapidity
($\Delta\eta$) of a so-called trigger particle from a second particle
(commonly referred to as associated particle), one constructs the normalized
per-trigger yield $Y_{N}$ for events of multiplicity $N$

$\frac{1}{N_{\mathrm{trig}}}\frac{\text{d}N_{\mathrm{assoc}}}{\text{d}\Delta\varphi\text{d}\Delta\eta}\equiv
Y_{N}\equiv J_{N}(\Delta\varphi,\Delta\eta)+N\left[1+\sum_{n}2v_{\rm
n}^{2}\cos n\Delta\varphi\right]\,.$ (5)

This quantity is normalized as particle yield, i.e. an integral over a certain
phase space region expresses the number of associated particles to each
trigger particle. This method does not require to determine an event plane, as
pointed out e.g. in Ref. 45. To characterize a collective component in terms
of Fourier coefficients $v_{n}$ and to distinguish it from a jet-like
component $J_{N}(\Delta\varphi,\Delta\eta)$ due to microscopic interactions,
one then exploits that both contributions should scale differently with event
multiplicity. To this end, one chooses a factor $a$ such that the jet-like
near-side contribution cancels in the difference between two different
multiplicities $N$ and $M$

$\Delta Y(\Delta\varphi)=Y_{N}-a\,Y_{M}\,,$ (6)

and one decomposes this difference in terms of flow coefficients $v_{n}$. This
subtraction assumes that the away-side structures scale in the same way with
$a$ as the near-side structures. Two different subtraction techniques have
been used regularly:

    1. (a)

The low-multiplicity subtraction method35, 34
assumes that the collective component vanishes at low multiplicity $M$ and
extracts the $v_{n}$’s from the ansatz

$\Delta
Y(\Delta\varphi)=G^{\prime}+N\sum_{n}2v_{n}^{2}\,\cos(n\Delta\varphi)\,.$

    2. (b)

The template method46
assumes that the collective component is significant at low multiplicity and
extracts the $v_{n}$ from the ansatz444For pedagogical reasons we present here
the typical simultaneous fit of $a$, $G$ and $v_{\rm n}$ in several steps.

$\Delta
Y(\Delta\varphi)=G\left(1+\sum_{n}2v_{n}^{2}\,\cos(n\Delta\varphi)\right).$

To illustrate the subtle differences between both subtraction methods, one may
consider the case that $a=1$ and that the values of the flow coefficients
$v_{n,N}$ in (5) vary with multiplicity like
$v_{n,M}^{2}=\alpha_{N}v_{n,N}^{2}$. For the two methods, the measured
$v_{n,{\rm method}}^{2}$ deviates then from the true $v_{n,N}^{2}$ like

$\frac{v_{\rm
n,low\,mult\,subtraction}^{2}}{v_{n,N}^{2}}=\frac{N-\alpha_{N}M}{N}\qquad\hbox{versus}\qquad\frac{v_{\rm
n,template}^{2}}{v_{n,N}^{2}}=\frac{N-\alpha_{N}M}{N-M}.$

If one of these ratios equals unity then the “true” signal $v_{n,N}^{2}$ for
the $n$-th harmonic at multiplicity $N$ is returned by the chosen method. The
low-multiplicity subtraction is correct for $\alpha_{N}=0$, and the template
method is correct of $\alpha_{N}=1$. As illustrated in Fig. 1, the larger a
true collective component will be at low multiplicity, the more it will be
underestimated by the low-multiplicity subtraction method and the more it will
be overestimated by the template method. This also becomes apparent in
measurements if both methods are applied to the same dataset. For instance
this explains the differences seen at low event multiplicity in the $v_{2}$
measurements displayed later in Fig. 5. In the limit of high multiplicity, the
results of both methods converge. However, the discussion of what constitutes
collectivity at small multiplicity requires a particularly critical assessment
of the method employed since even qualitatively different conclusions may be
reached at face value. We therefore suggest to make the differences explicit
in the nomenclature, for instance by using $v_{\rm n,low\,mult\,subtraction}$
and $v_{\rm n,template}$.

Figure 1:  Ratio of the measured $v_{\rm 2}^{2}$ to the true $v_{\rm 2}^{2}$
as a function of relative multiplicity $N/M$ for different cases of
$\alpha_{N}$ ($v_{2,M}^{2}=\alpha_{N}v_{2,N}^{2}$). The solid lines in the
upper part show the template method where the $v_{\rm 2}$ is exact or
overestimated depending on $\alpha_{N}$, while the dashed lines in the lower
part show the low-multiplicity subtraction method where the $v_{\rm 2}$ is
exact or underestimated depending on $\alpha_{N}$. Each of the two methods is
exact only for a single $\alpha_{N}$ and the two values are maximally
different.

  4. 4.

Two-particle cumulant flow $v_{n}\\{2\\}$ and higher orders $v_{n}\\{2m\\}$
The azimuthal angles of all pairs $(1,2)$ in an event can be directly used to
compute a cumulant47, 48:

$\langle\langle
e^{in(\varphi_{1}-\varphi_{2})}\rangle\rangle=v_{n}\\{2\\}\,v_{n}\\{2\\}+\underbrace{\langle
e^{in(\varphi_{1}-\varphi_{2})}\rangle_{\rm micro}}_{\sim O(1/N_{s})}\,,$ (7)

where the brackets denote averages over all particles in an event and over
event samples. Non-flow effects are assumed to result only from correlations
between pairs of particles that both originate from the same microscopic
source. If there are $N_{s}$ such sources in an event, the second “non-flow”
contribution is parametrically $O(1/N_{s})$ smaller. The first collective
contribution dominates the measured two-particle correlation if the observed
signal satisfies $v_{n}\\{2\\}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}O(1/N_{s}^{1/2})$.

The cumulant method may be regarded as being theoretically preferred since it
has the important property that the sensitivity to “non-flow” few-particle
correlations can be systematically reduced by purely statistical arguments. To
this end, one studies connected $2m$-point correlation functions
(“cumulants”). For instance, the fourth order cumulant reads

$\displaystyle\langle\langle
e^{in(\varphi_{1}+\varphi_{2}-\varphi_{3}-\varphi_{4})}\rangle\rangle_{c}$
$\displaystyle\equiv$ $\displaystyle\langle\langle
e^{in(\varphi_{1}+\varphi_{2}-\varphi_{3}-\varphi_{4})}\rangle\rangle-\langle\langle
e^{in(\varphi_{1}-\varphi_{3})}\rangle\rangle\,\langle\langle
e^{in(\varphi_{2}-\varphi_{4})}\rangle\rangle$ (8)
$\displaystyle-\langle\langle
e^{in(\varphi_{1}-\varphi_{4})}\rangle\rangle\,\langle\langle
e^{in(\varphi_{2}-\varphi_{3})}\rangle\rangle$ $\displaystyle=$
$\displaystyle-(v_{n}\\{4\\})^{4}\,+O\left(\frac{1}{N_{s}^{3}}\right)\,.$

In general, cumulants are constructed such that a non-flow correlation that
manifests itself only on the level of $(2m-2)$-particle correlations does not
manifest itself in the $(2m)$-th cumulant. The signal strength needs to be
$v_{n}\\{2m\\}>O\left(1/N_{s}^{(2m-1)/m}\right)$ for flow contributions to
dominate the signal. In the phenomenological practice, if the signal
$v_{n}\\{2m\\}$ stabilizes with increasing $m$, this is taken as evidence for
collectivity. Mathematically, the limit $m\rightarrow\infty$ of these higher-
order cumulants corresponds to the so-called LYZ (Lee-Yang-Zero) method 49.
The parametric separation of flow- and non-flow effects via
$O\left(1/N_{s}^{(2m-1)/m}\right)$ power counting is powerful if $N_{s}$ is
large, i.e. if event multiplicities are large. Therefore, the capability of
the cumulant method to separate by purely statistical arguments flow from non-
flow is limited for small systems.

  5. 5.

Subevent technique
Higher-order cumulants are computationally intense and naively scale with
$O(N^{2m})$. Even for a 4-particle cumulant this is computationally
challenging if applied to large event samples. The complexity can be reduced
to $O({N})$ by the calculation of so-called $Q$-vectors whose expressions
instead of considering $2m$ particles sum over single particles47, 50. As
mentioned before, jet-like fragmentation patterns and resonance decays lead
mainly to correlations between particles that are close in
$\Delta\eta\times\Delta\varphi$. In order to exclude those when computing the
$Q$-vectors, where an $\eta$-gap between particles cannot be applied (as only
single particles are considered in the calculation), the event is split into
several subevents51, each of which is used for the computation of a
$Q$-vector. The extracted coefficient is then:

$\left(v_{n}\\{2,|\Delta\eta|>\gamma\\}\right)^{2}=\langle\langle
e^{in(\varphi_{i}-\varphi_{j})}\rangle\rangle\,,$ (9)

where all particles $i$ are in the range $\eta>\gamma/2$ while all particles
$j$ are in the range $\eta<-\gamma/2$ with $\gamma\geq 0$. Experimentally, one
can then test the dependence on $\gamma$ and if the signal stabilizes with
increasing $\gamma$, this is taken as evidence for collectivity. It should be
noted that regardless how large $\gamma$, the back-to-back jet structure is
not excluded from $v_{n}\\{2,|\Delta\eta|>\gamma\\}$, but can be suppressed by
extending this method into calculating higher-order cumulants from more than
two subevents51.

This enumeration demonstrates that there are numerous experimental measures
for $v_{\rm n}$ coefficients. Those use different definitions with different
underlying assumptions, in particular on what is the definition of the physics
in absence of collectivity. Only some of the different methods have separate
notation, e.g. $v_{n}\\{2m\\}$. Therefore, care has to be taken when comparing
different $v_{\rm n}$ between each other or with theoretical calculations.

###  2.2 Hadrochemical measures of collectivity

Hadrochemistry addresses the question with which relative abundance different
hadron species are produced in a collision. What is measured are identified
particle yields and their ratios, particle-identified transverse momentum
spectra and particle-identified flow coefficients. The baseline on top of
which collective effects are established is traditionally chosen to be set by
minimum bias $\mathrm{pp}$ collisions. There are some technical challenges
related to the question of how feed-down contributions from resonance decays
are included in or subtracted from particle-identified data52, and how the
experimental procedure is paralleled in theoretical descriptions. On the
modelling side, an additional challenge arises from the incomplete knowledge
of short-lived higher mass resonance states which feed down inevitably into
the measured yields, see e.g. Ref. 53. In general, however, what is measured
is easily stated and it seems unnecessary to describe the definition of
hadrochemical measures of collectivity in as much detail as we did for
measures of anisotropic flow. In the following, we therefore focus almost
exclusively on the physics motivations for which hadrochemical measures of
collectivity are studied:

  1. 1.

Hadrochemical equilibration
In QCD at sufficiently high energy density, the gluon fusion process
$gg\leftrightarrow s\bar{s}$ can drive the strangeness content to chemical
equilibrium. The resulting expectation of strangeness enhancement in heavy-ion
collisions counts amongst the very first proposed signatures of QGP formation
54, 55. In nucleus–nucleus collisions, strangeness enhancement (compared to a
baseline set by minimum bias $\mathrm{pp}$ collisions) is observed in the
phase-space integrated relative abundances of all identified strange and
multi-strange hadrons at all relativistic center of mass energies from BNL AGS
to the LHC 56. These observations can be accounted for by fitting to a
statistical hadronization model that involves two intensive quantities
(temperature and a baryo-chemical potential) and one extensive one (volume)
56, 53. Also the relative abundances of charm and beauty hadrons can be
included in the systematics of this statistical hadronization model 53, 57, if
their total production rate is set by initial hard (vacuum-like) processes.
That means that these rates are not thermally produced but statistically
distributed.

  2. 2.

Particle-identified flow
To the extent to which collective flow is built up in the partonic phase of a
collision, all hadron species emerge from the same common flow field
$u_{\mu}$. According to the Cooper-Frye freeze-out prescription, all hadrons
will decouple with statistical weights $\exp\left[-p^{\mu}u_{\mu}/T\right]$
from locally equilibrated patches of a globally expanding system 58. In
analytical blast-wave models59 and in fluid-dynamical simulations, this leads
to a characteristic mass hierarchy in the $p_{T}$-dependent particle-
identified flow coefficients $v_{n}(p_{T})$ and in the single inclusive
$p_{T}$ spectra (“radial flow”). These dependencies have been observed in
central and semi-peripheral nucleus–nucleus collisions at the CERN SPS, at
RHIC and at the LHC for sufficiently small transverse momenta (see Figure 8
and related discussion).
One central question is whether the onset of flow seen in the smallest pp and
$\mathrm{p}$–$\mathrm{Pb}$ collision systems has the same collective dynamical
origin as the flow in nucleus–nucleus collisions. Testing whether anisotropic
flow in large and small collision systems shows the same PID-dependence is of
obvious interest in this context.

  3. 3.

Hadronization mechanisms
The dynamics of hadronization in QCD is not understood. In high-momentum
transfer processes where QCD collinear factorization applies, high-$p_{T}$
partons hadronize the same irrespective of the hard partonic process in which
they are produced. Remarkably, also the Lund string60 or cluster hadronization
models10 implemented in multi-purpose $\mathrm{pp}$ event generators are
process-independent in this sense. On the other hand, the observation of
significant strangeness enhancement (as well as a suite of other measurements)
falsifies the assumption of a process-independent hadronization in
nucleus–nucleus collisions. Hadronization in a dense nuclear environment leads
to a characteristically different hadrochemical composition compared to that
in the vacuum. One important question is whether a process-independent
hadronization mechanism of bulk-particle production can also be falsified in
small pp and $\mathrm{p}$–$\mathrm{Pb}$ collision systems.
There are process-dependent hadronization models that are sensitive to the
environment in which hadronization occurs. In particular, in recombination
models hadrons result from combining valence-like partons that are produced
within the same collision in different partonic interactions61, 62. Signatures
of recombination are normally searched for at intermediate transverse momenta
($2\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$<$}}p_{T}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}5\,\mathrm{GeV}$). The reason is that
this mechanism cannot be dominant for soft bulk particle production as it
reduces the number of independent degrees of freedom and thus entropy. Also,
fragmentation from a high-$p_{T}$ power-law tail dominates over recombination
from a low-$p_{T}$ exponential slope61, 62. One qualitative prediction of
recombination is a characteristic baryon-meson splitting as a consequence of
the different number of valence-like quarks that need to be combined. It is of
obvious interest to study whether this baryon-meson separation persists across
system size.

###  2.3 Hard processes in small systems

In sufficiently central nucleus–nucleus collisions, essentially all
high-$p_{T}$ hadronic final states show numerically large medium-
modifications. Our understanding of these jet quenching 63 phenomena is based
on mechanisms of radiational and collisional parton energy loss that strongly
suppress high-$p_{T}$ colored probes in sufficiently large collision systems.
Jet quenching arises from interactions of jets with a medium that exhibits
collective phenomena. In theory, there is therefore a direct relation between
jet quenching and flow. If fluid dynamics or kinetic transport is the
dynamical origin of flow, then there must be medium interactions with colored
degrees of freedom at all momentum scales, and thus, there must be jet
quenching.

In this context, it is remarkable that in the smallest collision systems (pp,
$\mathrm{p}$–$\mathrm{Pb}$), LHC experiments have discovered sizable
anisotropic flow without finding evidence for jet quenching. Jet quenching in
small systems may have escaped detection so far since effects become smaller
with decreasing system size and since systematic uncertainties in measuring
jet quenching become more important. If this is so, it calls for studies based
on more sensitive observables and supported by more accurate calculations. We
therefore comment here on problems related to determining the system size
dependence of jet quenching in small systems.

In nucleus–nucleus collisions, one of the simplest and most widely used
measures of jet quenching is the nuclear modification factor

$\displaystyle\left.R^{h,j}_{\text{AA}}(p_{T},y)\right|_{\text{cent}}=\frac{1}{\langle
T_{AA}\rangle}\frac{\left.\frac{1}{N_{\text{ev}}}\frac{dN^{h,j}_{\text{AA}}}{dp_{T}dy}\right|_{\text{cent}}}{\frac{d\sigma_{pp}^{h,j}}{dp_{T}dy}}\,.$
(10)

This factor compares the number $dN^{h,j}_{\text{AA}}$ of charged hadrons
($h$) or calorimetrically defined jets ($j$) (produced in AA collisions in
some range of transverse momentum $p_{T}$ and rapidity $y$) to an equivalent
number of such hadrons or jets produced in pp collisions. Here, the equivalent
number of collisions $\langle N_{\rm coll}\rangle=\sigma_{\rm pp}^{\rm
inel}\,\langle T_{AA}\rangle$ is the product of the total inelastic pp cross
section and the nuclear overlap function $\langle T_{AA}\rangle$.
$R_{AA}(p_{T},y)$ is typically determined as a function of centrality.
Centrality is defined in terms of the fraction of highest multiplicity events
(or highest total transverse energy events) of the total inelastic
nucleus–nucleus cross section, and it is related to $\langle T_{AA}\rangle$ by
Glauber-type models 64.

For centrality to be a good proxy of system size, event-multiplicity should be
tightly correlated to the geometry of the collision and it should not be
correlated with the hard process. These assumptions are satisfied
approximately in sufficiently large systems but they become problematic in
small systems. Estimates of the uncertainties associated to $\langle
T_{AA}\rangle$ range from 3% in central to 15% in peripheral
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions 65. In addition, in peripheral
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions (centralities larger than $\sim$ 70%),
Eq. (10) is known to be affected by significant event selection and geometry
biases 66, 67. The centrality dependence of hard processes in
$\mathrm{p}$–$\mathrm{Pb}$ collisions reveals even larger biases. Some of
these biases arise from the larger influence of fluctuations in smaller
systems, other biases originate in correlations between the hard process and
the global event multiplicity68, 69. Finally, the nuclear modification factor
(10) measures jet quenching with respect to a baseline set by pp collisions,
and it is therefore not a good starting point for analyzing jet quenching in
pp collisions.

These findings indicate that the measurement of single inclusive jet or hadron
spectra as a function of event multiplicity (centrality) is of limited use for
the characterization of jet quenching in the smallest collision systems. For
small collision systems, it is clearly problematic to include geometrical
information such as $\langle T_{AA}\rangle$ in the definition of jet quenching
measures, as such information depends on soft physics models that are
difficult to constrain for small systems. One way out is the use of inclusive,
minimum bias spectra in light nucleus collisions for which a nuclear
modification factor can be defined without reference to soft physics modeling,
$\left.R^{h,j}_{\text{AA}}(p_{T},y)\right|_{\text{min
bias}}=\frac{1}{A^{2}}\frac{d\sigma^{h,j}_{\text{AA}}}{dp_{T}dy}\Big{/}\frac{d\sigma_{pp}^{h,j}}{dp_{T}dy}$.
It has been shown that such inclusive measurements have much smaller
uncertainties and that they therefore allow for the unambiguous detection of
smaller jet quenching signals in nucleus–nucleus collisions with $A$ as small
as oxygen 70, 71. A second general possibility of detecting signatures of jet
quenching in small collision systems is the use of self-normalizing processes
such as dijet asymmetries. We recall that the very first proposal of a jet
quenching effect in small p$\bar{\rm p}$ collisions made by Bjorken in 1982
was of that type 72.

##  3 Experimental Overview

###  3.1 Long-range correlations

The most prominent and from the point of view of analysis methodology rather
simple observation is the one of long-range effects in two-particle
correlations (5). QCD effects implemented in standard multi-purpose event
generators of proton–proton collisions lead to distinct patterns in these
correlations: At small angular distances ($\Delta\varphi\approx 0$ and
$\Delta\eta\approx 0$) a so-called near-side peak results from jet-like
fragmentation and resonance decays. A back-to-back structure from momentum
conservation is observed at $\Delta\varphi\approx\pi$, called away side. This
away-side structure is not concentrated around $\Delta\eta\approx 0$ since the
center-of-mass systems of any two partons in a hadronic collision is
distributed widely in rapidity. In addition to these expected patterns, Fig. 2
(left panel) shows also positive long-range rapidity correlations for
$|\Delta\eta|\gtrsim 2$ on the near side. This was surprising since the known
microscopic dynamics operational in proton–proton collisions does not lead to
long-range near-side correlations. Fig. 2 displays the first ever observed
“collective effect” in proton–proton collisions. At the time, it was based on
the first sizable LHC high-multiplicity sample although from today’s
perspective this sample is very small.

The second hallmark result was the observation of a double-ridge structure in
$\mathrm{p}$–$\mathrm{Pb}$ collisions34, 35, displayed on the right panel of
Fig. 2. For reasons which are detailed in section 2.1, the observation in
small systems (contrary to large systems) required to take the difference
between high-multiplicity and low-multiplicity collisions. Fig. 2 shows
clearly an almost rapidity-independent second harmonic $v_{2}$ imprinted on
the azimuthal two-particle correlation. This is the telltale sign of an
elliptic flow signal that had been observed in heavy-ion collisions at all
fixed-target and collider energies. Elliptic flow was thought to arise from
the collective dynamical response of a system to pressure gradients related to
initial spatial anisotropies $\varepsilon_{2}$. However, such a collective
effect had not been expected in the small $\mathrm{p}$–$\mathrm{Pb}$ collision
system. At this moment, the scientific community was convinced that something
extraordinary was going on and embarked on understanding if the observed
effects were truly of similar nature as in heavy-ion collisions. To this end,
studies focused first on characterizing higher order flow harmonics $v_{n}$
and their relation to finer details $\varepsilon_{n}$ of the initial
conditions (for motivation, see the discussion of (1)). Soon after, the
searches for high-density QCD phenomena in small systems extended to nearly
all areas of heavy-ion physics, see Table LABEL:table:smallsystems.

Figure 2:  Two-particle correlation showing the first observation of the long-
range correlation on the near side in pp collisions (left panel; Figure from
Ref. 33) and the first observation of the double-ridge structure in
$\mathrm{p}$–$\mathrm{Pb}$ collisions (right panel; Figure from Ref. 35). The
left panel shows the two-particle correlation without any subtraction. Hence
the near-side jet peak is clearly visible and the the away side is dominated
by the jet component. The right panel uses the low-multiplicity subtraction
described in section 2.1 which effectively reduces the jet component. This
results in an almost complete suppression of the near-side jet peak allowing
one to observe two ridges. The right panel is normalized as per-trigger yield,
while the left panel uses a relative normalization, see Ref. 33 for details.

Azimuthal anisotropies The observed long-range structures are quantified by
decomposing the per-trigger yield into Fourier coefficients $v_{n}$, see Eq.
5. The first measurement of non-zero $v_{\rm 2}$ and $v_{\rm 3}$ in
$\mathrm{p}$–$\mathrm{Pb}$ collisions, shown in Fig. 4, provided another
indication that the phenomena seen in small systems are heavy-ion like. Since
then, a large number of detailed studies have been done on different
coefficients of anisotropic flow. Today, measurements extend up to $v_{\rm 9}$
in $\mathrm{Pb}$–$\mathrm{Pb}$, $v_{\rm 5}$ in $\mathrm{p}$–$\mathrm{Pb}$,
$v_{\rm 3}$ in $\mathrm{d}$–$\mathrm{Au}$, and $v_{\rm 4}$ in $\mathrm{pp}$
collisions for charged particles. Fig. 4 shows the state-of-the art
measurement of $v_{\rm 2}$ in $\mathrm{pp}$, $\mathrm{p}$–$\mathrm{Pb}$ and
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions. These coefficients exhibit a weaker
multiplicity dependence in $\mathrm{pp}$ and $\mathrm{p}$–$\mathrm{Pb}$
collisions than in $\mathrm{Pb}$–$\mathrm{Pb}$ collisions where this is
closely related to the shape of the overlap region, while $v_{\rm 3}$ are
independent of the collision system as a function of multiplicity73, 74, 75,
34, 35, 36, 76, 77, 78, 46, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89. The
slight trend of increasing $v_{\rm 2}$ for decreasing multiplicity in pp
collisions is related to the extraction method used, see section 2.1. The
sensitivity to the method becomes clear when both methods are applied to the
same dataset, which is shown in Fig. 5. The $\eta$ dependence is overall weak
and similar in small and large systems 90, 91, 92, 93, 94, 82, 95, 96, 97, 98.

Figure 3:  $v_{\rm 2}$ and $v_{\rm 3}$ coefficients extracted from low-
multiplicity subtracted two-particle correlations in
$\mathrm{p}$–$\mathrm{Pb}$ collisions (Figure from Ref. 35). For details on
the methodology see Sec. 2.1.

Figure 4:  The $v_{\rm 2}$ coefficient extracted with the template method in
pp, $\mathrm{p}$–$\mathrm{Pb}$ and $\mathrm{Pb}$–$\mathrm{Pb}$ collisions
(Figure from Ref. 84). For details on the methodology see Sec. 2.1.

Figure 5:  The $v_{\rm 2}$ coefficient extracted with the template method
(left panel) and the low-multiplicity subtraction method (right panel). The
overestimation with the template method and underestimation with the low-
multiplicity subtraction method are clearly visible. Figure from Ref. 46. For
details on the methodology see Sec. 2.1, in particular Fig. 1

Figure 6:  Comparison of the $v_{\rm 3}$ measurement of PHENIX101 (red filled
symbols) and STAR102 (black open symbols) in three different collision systems
which favor round ($\mathrm{p}$–$\mathrm{Au}$), elliptical
($\mathrm{d}$–$\mathrm{Au}$) and triangular
($\mathrm{{}^{3}He}$–$\mathrm{Au}$) shapes in the initial state.

Figure 7:  Higher-order correlations of the $v_{\rm 2}$ coefficient measured
in pp collisions (Figure adapted from Ref. 83).

At RHIC, the influence of the initial geometry has been studied with special
collision configurations that favor round ($\mathrm{p}$–$\mathrm{Au}$),
elliptical ($\mathrm{d}$–$\mathrm{Au}$) and triangular
($\mathrm{{}^{3}He}$–$\mathrm{Au}$) shapes in the initial state. Fig. 7
compares the $v_{3}$ coefficients measured by the PHENIX101 and STAR102
collaborations. The larger signal strength in
$\mathrm{{}^{3}He}$–$\mathrm{Au}$ as compared to $\mathrm{p}$–$\mathrm{Au}$
and $\mathrm{d}$–$\mathrm{Au}$ observed by PHENIX101 supports the conclusion
that a larger triangular eccentricity leads to a larger triangular $v_{\rm 3}$
also in these small systems. However, this conclusion is not supported by the
STAR102 data which differ from those of PHENIX by up to a factor 3 in
$\mathrm{p}$–$\mathrm{Au}$ and $\mathrm{d}$–$\mathrm{Au}$, while they are
consistent for $\mathrm{{}^{3}He}$–$\mathrm{Au}$. In model studies, half of
this difference could be attributed to the different methods and rapidity
ranges555PHENIX uses the event-plane method with detectors at $-3.9<\eta<-3.1$
and $|\eta|<0.35$ while STAR uses the template method with both particles
within $|\eta|<0.9$ and an $\eta$ gap of 1 unit. used in the extraction
methods285 but this does not reconcile the different data. The discrepancy is
of a qualitative, not only a quantitative nature. To draw lasting physics
conclusions from the unique RHIC opportunity of studying these different small
collision systems, a resolution of the apparent inconsistencies between STAR
and PHENIX remains important.

Multi-particle nature of the observed correlation In heavy-ion collisions,
collective phenomena are a correlation between all particles (through common
symmetry planes) and therefore correlations are in principle present at all
orders. To test this assumption, higher-order cumulants (see Eq. (8) and
following text) have been measured using the Lee-Yang Zeros (LYZ) method and
multi-particle correlation techniques with up to 8 particles in both
$\mathrm{Pb}$–$\mathrm{Pb}$ and $\mathrm{p}$–$\mathrm{Pb}$ collisions and up
to 6 particles in $\mathrm{pp}$ collisions (shown in Fig. 7) for $v_{\rm 2}$
111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 83, 84, 121, 96, 122, 88,
123. The measurements of the cumulants at different orders ($n\geq 4$) are
similar within 10% within each collision system. This is a tell-tale sign of
collectivity, but it does not identify unambiguously the dynamical mechanism
underlying collectivity. From fluid dynamics (see section 4.1.1) to one-hit
kinetic theory (see Section 4.1.2), very different dynamical frameworks can
lead to correlations amongst essentially all particles in an event.

Figure 8:  Compilation of identified-particle $v_{\rm 2}$ in
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions (left panel) compared to
$\mathrm{p}$–$\mathrm{Pb}$ collisions (right panel). Figures from Ref. 171).

Particle-mass dependence Studying the $v_{\rm n}$ coefficients as a function
of $p_{\scriptstyle\mathrm{T}}$ for different particle species shows a
characteristic mass dependent pattern up to $v_{\rm 5}$ in
$\mathrm{Pb}$–$\mathrm{Pb}$ (left panel of Fig. 8), $v_{\rm 3}$ in
$\mathrm{p}$–$\mathrm{Pb}$ (right panel of Fig. 8) and other small systems, as
well as $v_{\rm 2}$ in $\mathrm{pp}$ collisions where heavier particles are
depleted at low $p_{\scriptstyle\mathrm{T}}$ 103, 104, 105, 106, 78, 81, 83,
107, 87, 108, 110. If all hadron species were emitted with thermal weight
$\exp\left[-p_{\mu}u^{\mu}/T\right]$ from the same common flow field
$u^{\mu}$, then this would lead to a characteristic
$m_{T}=\sqrt{m^{2}+p_{T}^{2}}$-scaling of radial and anisotropic flow (see
section 2.2). If, on the other hand, hadron species would result from the
recombination of valence-like quarks, then a baryon-meson grouping of
observables may be expected. We emphasize that $m_{T}$-scaling and baryon-
meson grouping are distinct dependencies that cannot be realized
simultaneously. Existing data show approximate $m_{T}$-scaling at low and
approximate meson-baryon grouping at intermediate $p_{T}$. Remarkably,
essentially the same hadrochemical and $p_{T}$-dependencies are observed
across all systems from $\mathrm{Pb}$–$\mathrm{Pb}$ to
$\mathrm{p}$–$\mathrm{Pb}$ to $\mathrm{pp}$ although both mechanisms may be
expected to depend on system size or phase space density.

The extent to which signatures of collectivity show close commonalities
between small and large collision systems had not been expected. Beyond the
early findings described above, commonalities have been found by now in a
large number of refined measurements listed in Table LABEL:table:smallsystems.
In the following, we summarize the main experimental conclusions as follows:

Non-linear response Correlations between different anisotropic flow harmonics
characterize the non-linear coupling between different collective excitations.
In the modelling of heavy-ion collisions, non-linear mode-mode coupling is a
generic feature of fluid dynamics286, 287, 288, 289 where it is found to be
sensitive to transport properties. However, non-linear mode-mode coupling is
also a generic feature of kinetic transport theory where it persists even in
the dilute one-hit approximation 290, 291, see section 4.1.2 for details.
Mixed harmonics or so-called symmetric cumulants have been measured up to
$\rm{SC}(5,3)$ in $\mathrm{Pb}$–$\mathrm{Pb}$ and $\rm{SC}(4,2)$ in
$\mathrm{Au}$–$\mathrm{Au}$, $\mathrm{Xe}$–$\mathrm{Xe}$,
$\mathrm{p}$–$\mathrm{Pb}$ and $\mathrm{pp}$ collisions 124, 125, 126, 127,
86, 128, 88, 129, 130. A further class of observables, linear ($v_{\rm
n}^{\rm{L}}$) and non-linear ($v_{\rm n}^{\rm{NL}}$) flow modes, have been
investigated in $\mathrm{Pb}$–$\mathrm{Pb}$ collisions up to the seventh order
flow harmonic 131, 132, 89 but not yet in $\mathrm{pp}$ and
$\mathrm{p}$–$\mathrm{Pb}$ collisions. In general, the measured non-linear
response varies only mildly with system size.

More detailed $v_{n}$-characterizations Numerous other characterizations
constrain the phenomenology of $v_{n}$. The factorization of long-range
azimuthal two-particle correlations into a product of single particle
anisotropies and the breaking of this factorization due to event-plane angle
decorrelations in $p_{\scriptstyle\mathrm{T}}$ and $\eta$ has been measured in
both $\mathrm{Pb}$–$\mathrm{Pb}$ and $\mathrm{p}$–$\mathrm{Pb}$ collisions
133, 85, 135, 136. Recent measurements show this behavior using four-particle
correlators 137. With the existing data such measurements are not yet possible
in $\mathrm{pp}$ collisions. It would also be interesting to study event-by-
event distributions of $v_{\rm n}$ coefficients in both
$\mathrm{p}$–$\mathrm{Pb}$ and $\mathrm{pp}$ collisions, which have so far
been done in $\mathrm{Pb}$–$\mathrm{Pb}$ collisions 138, 139, 140. The
interplay of flow harmonics and radial expansion is studied with correlations
of $v_{\rm n}$ coefficients and the mean $p_{\scriptstyle\mathrm{T}}$ event-
per-event. The correlation has been measured through modified Pearson
correlation coefficient as a function of event activity up to fourth order in
$\mathrm{Pb}$–$\mathrm{Pb}$, third order in $\mathrm{Xe}$–$\mathrm{Xe}$
collisions and second order in $\mathrm{p}$–$\mathrm{Pb}$ collisions 141, 142.
First measurements of higher-order correlations, e.g. $v_{\rm 2}$–$v_{\rm
3}$–$p_{\scriptstyle\mathrm{T}}$ have been made 142.

Directed flow In $\mathrm{Pb}$–$\mathrm{Pb}$ collisions, directed flow of
charged particles at mid-rapidity was measured relative to the collision
symmetry plane defined by the spectator nucleons, and evidence for dipole-like
initial-state density fluctuations in the overlap region was found143. In
small systems, the concept of directed flow is less clear, especially in
$\mathrm{pp}$ collisions. If there is collectivity in $\mathrm{pp}$
collisions, one could also expect a non-zero directed flow measurement. This
is technically challenging since the measurement of the spectator plane is not
feasible in small systems and, hence, $v_{\rm 1}$ could only be measured using
higher-order ($n\geq 4$) cumulants which has not been achieved to date.

Charge-dependent correlations The balance function probes the charge creation
time and the development of collectivity in the produced system. Its width,
$\langle\Delta\eta\rangle$ and $\langle\Delta\varphi\rangle$, has been
measured for charged particles in $\mathrm{pp}$, $\mathrm{p}$–$\mathrm{Pb}$
and $\mathrm{Pb}$–$\mathrm{Pb}$ collisions144, 145. A picture emerges where
the system exhibits larger radial flow with increasing multiplicity but also
whose charges are created at the later stages of the collision. Charge-
dependent azimuthal correlations are measured in both
$\mathrm{Pb}$–$\mathrm{Pb}$ and $\mathrm{p}$–$\mathrm{Pb}$ collisions172, 173,
174, 177. These studies assess the chiral magnetic effect (CME) and the chiral
magnetic wave (CMW) on the produced particles. Their interpretation is still
today influenced by strong background contributions, for example from local
charge conservation and possibly radial and anisotropic flow.

###  3.2 Bulk properties

Particle spectra and production Particle production yields have been measured
for numerous particle species in small and large collision systems167, 153,
154, 157, 158, 168, 52, 161, 169, 170. Figure 9 shows the production yields of
different identified non-strange, strange and multiple-strange hadrons
normalized to the pion yield. As a function of multiplicity, these ratios
increase more for strange than for non-strange hadrons, and they increase more
for multi-strange hadrons than for strange hadrons 168, 52, 169. That
strangeness in nucleus–nucleus collisions is enhanced compared to minimum-bias
proton–proton collisions has been observed at all relativistic collision
energies at the BNL AGS, the CERN SPS, the RHIC56 and now the LHC. Up until
today, the factor 25 enhancement of the $\Omega/\pi$ ratio observed by the
NA57 Collaboration at the CERN SPS in comparing yields in
$\mathrm{Pb}$–$\mathrm{Pb}$ and $\mathrm{p}$–$\mathrm{Pb}$ collisions is
arguably one of the numerically largest medium-modification ever observed in
heavy-ion collisions 292. The qualitatively novel information of Figure 9 is
that this strangeness enhancement sets in smoothly with increasing
multiplicity so that even $\mathrm{pp}$ collisions show an enhanced
strangeness content for higher multiplicity, and $\mathrm{p}$–$\mathrm{Pb}$
collisions must not be considered as a baseline in which strangeness
enhancement is absent. Moreover, strangeness enhancement shows a universal
dependence on event multiplicity irrespective of the collision system. This is
in marked contrast to kinetic signatures of collectivity such as elliptic flow
(see e.g. Fig. 4) that depend not only on event multiplicity as they contain
also spatial information about the collision systems and their centrality
dependence.

The $p_{\scriptstyle\mathrm{T}}$ spectra of identified particles harden with
increasing multiplicity for all studied collision systems. The hardening is
stronger for heavier particles153, 162, 163, 156, 164, 165, 166, 160. In
central nucleus–nucleus collisions, the slopes of soft charged pion, kaon and
(anti-)proton $p_{\scriptstyle\mathrm{T}}$ spectra are consistent with models
in which all hadron spectra are blue-shifted by the same transverse radial
flow 151. In this framework, a larger radial flow is observed52 in
$\mathrm{pp}$ and $\mathrm{p}$–$\mathrm{Pb}$ collisions as compared to
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions at the same multiplicity.

Figure 9:  Particle production yields of strange and non-strange particles
normalized to the pion yield as a function of multiplicity in pp,
$\mathrm{p}$–$\mathrm{Pb}$, $\mathrm{Pb}$–$\mathrm{Pb}$ and
$\mathrm{Xe}$–$\mathrm{Xe}$ collisions (Figure from Ref. 171).

Source size The measurement of quantum-statistic correlations between pairs of
same-charge particles in particular pions at low-momentum transfer allows to
assess the size of the emitting source 293. These freeze-out radii are
measured in three orthogonal directions (“out”, “side”, “long”). While clearly
the expectation is that these are vastly different in small and large
collision systems, their dependence on event and particle pair variables can
be compared. The radii scale with $\sqrt[3]{N}$ in all collision systems
indicating a constant density at freeze-out and decrease with increasing pair
momentum $k_{\scriptstyle\mathrm{T}}$. This is consistent with the expectation
that velocity gradients in the collision region should lead to smaller
“homogeneity regions” with increasing $k_{\scriptstyle\mathrm{T}}$ 293. The
size along the emission direction is similar to the geometric size of the
system ($R_{\rm out}/R_{\rm side}\approx 1$) in $\mathrm{Pb}$–$\mathrm{Pb}$
collisions173, 175, 176, 178, 178, 179, 180 while this ratio is smaller than
unity for both $\mathrm{p}$–$\mathrm{Pb}$ and $\mathrm{pp}$ collisions172,
173, 174, 177.

Direct photons Direct-photon measurements have been performed in
$\mathrm{Pb}$–$\mathrm{Pb}$ and $\mathrm{pp}$ collisions in the
$p_{\scriptstyle\mathrm{T}}$ region sensitive to thermal photon production.
For central $\mathrm{Pb}$–$\mathrm{Pb}$ collisions, an effective temperature
(averaged over the evolution) of about 300 MeV has been extracted from the
slope of the photon transverse momentum spectrum 181, 183. Fluid models that
reproduce these spectra start from an initial QGP temperature of more than 400
MeV181 in central $\mathrm{Pb}$–$\mathrm{Pb}$ collisions. In $\mathrm{pp}$
collisions, no significant direct-photon signal has been identified at low
momentum182. The reported upper bounds lie well above expectations of next-to-
leading order QCD. A significant further increase in experimental accuracy
would be needed to gain access to a possible thermal component.

###  3.3 Hard probes

Parton energy loss In nucleus–nucleus collisions, jet quenching and its
centrality dependence has been quantified through measurements of the nuclear
modification factor $R_{\rm AA}$ (10) 189, 190, 191, 192, 194, 195, 196, 65,
67, 200, 201, 202, 203, 204 as well as correlation measurements213, 214, 215,
217, 218, 219, 220, 222, 204. Fig. 10 (left panel) shows a strong suppression
($R_{\rm AA}$ as small as 0.13) in central $\mathrm{Pb}$–$\mathrm{Pb}$
collisions for charged particles and various identified particles while the
$R_{\rm AA}$ for $\gamma$ is consistent with unity. Also, consistent with the
qualitative expectations of parton energy loss, a large asymmetry in back-to-
back jet $p_{\scriptstyle\mathrm{T}}$ accompanied by slightly modified jet
fragmentation functions inside small jet cone sizes ($\rm{R}=0.4$) has been
observed. The radiated energy appears mostly at large angles
($\rm{R}>0.8$)205, 206, 208, 209, 210, 211, 212. Furthermore, in semi-
peripheral collisions, reduced but non-zero $v_{\rm n}$ coefficients have been
found to persist up to the highest transverse momenta 223, 224 and di-jet
events have been found correlated with elliptic flow 225. These data are
interpreted in terms of path-length dependent parton energy loss which in
azimuthally asymmetric collision region leads to an aziumthal dependence of
high-$p_{T}$ hadron yields. They illustrate the connection between jet
quenching and flow emphasized in section 2.3.

Early data on $R_{\rm AA}$ in $\mathrm{Pb}$–$\mathrm{Pb}$ collisions show that
a significant suppression $R_{\rm AA}<1$ persists even in the most peripheral
event classes. By now, this phenomenon is understood as an intrinsic selection
bias (events with less activity are biased towards larger nucleon–nucleon
impact parameter which in turn have less high-$Q^{2}$ scatterings) 67. There
is only a tight upper bound but no direct evidence for parton energy loss in
the most peripheral $\mathrm{Pb}$–$\mathrm{Pb}$ collisions 67, 294.

Figure 10:  Left panel: $R_{\rm AA}$ in $\mathrm{Pb}$–$\mathrm{Pb}$ collisions
for charged particles as well as $\pi$, K, p and $\phi$ compared to $\gamma$
(Figure from Ref. 171). Right panel: $R_{\rm AA}$ in
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions for charged particles up to high
$p_{\scriptstyle\mathrm{T}}$ compared to the result in
$\mathrm{p}$–$\mathrm{Pb}$ collisions (Figure from Ref. 65).

In $\mathrm{p}$–$\mathrm{Pb}$ collisions, measurements of inclusive
high-$p_{\scriptstyle\mathrm{T}}$ hadron and inclusive jet yields in minimum-
bias $\mathrm{p}$–$\mathrm{Pb}$ collisions at the LHC are within the current
accuracy of approximately 20% consistent with $R_{\rm pPb}=1$193, 197, 198,
65, 199, 201, 203, see Fig. 10 (right panel). As explained in section 2.3, in
small systems, the binning of $R_{\rm AA}$ and other measures of parton energy
loss as a function of event activity (centrality) does not necessarily inform
us about the system size dependence of jet quenching. In particular, inclusive
jet yields $R_{\rm pPb}$ are strongly suppressed relative to unity in
“central” $\mathrm{p}$–$\mathrm{Pb}$ collisions, and strongly enhanced in
“peripheral” $\mathrm{p}$–$\mathrm{Pb}$ collisions 193, but this is attributed
to selection biases 216 and it does not inform us about the onset of jet
quenching in “sufficiently central” $\mathrm{p}$–$\mathrm{Pb}$ collisions.
Also, signatures of jet quenching are neither observed in inclusive hadron
production nor in the dijet transverse momentum imbalance207, 197, 201. No
significant jet quenching effects have been found in
$\mathrm{p}$–$\mathrm{Pb}$ collisions using correlations by either studying
hadron production relative to the jet direction221 or by measuring the semi-
inclusive yield of jets recoiling from a high-$p_{\scriptstyle\mathrm{T}}$
trigger hadron216. Ref. 216 sets an upper limit of 400 MeV (at 90% CL) on
medium-induced energy transport outside a jet cone with $R=0.4$. This is about
a factor 20 smaller than the magnitude of energy transport from within the jet
to outside the jet cone in $\mathrm{Pb}$–$\mathrm{Pb}$ collisions 213. On the
other hand, as mentioned, the presence of energy loss leads to sizable $v_{\rm
n}$ coefficients at high-$p_{\scriptstyle\mathrm{T}}$ 223 in
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions. Surprisingly, non-zero $v_{\rm 2}$ are
measured above 9 $~{}\text{GeV}\\!/\mspace{-2.0mu}c$ also in
$\mathrm{p}$–$\mathrm{Pb}$ collisions despite the lack of other energy-loss
signals184. This indicates that our current understanding of the interplay
between jet quenching and flow phenomena is incomplete.

$v_{\rm n}$ in presence of a hard probe. The $v_{\rm 2}$ measured in events
with a $Z$ candidate in $\mathrm{pp}$ collisions are identical to the
measurement in inclusive events185. Similarly, selecting events with or
without jets in pp collisions does not influence the measured $v_{\rm 2}$186.
Measurements in $\mathrm{p}$–$\mathrm{Pb}$ collisions in presence of a
high-$p_{\scriptstyle\mathrm{T}}$ jet are similar to the ones in inclusive
events for $v_{\rm 2}$ and $v_{\rm 3}$ 184.

Jet constituent $v_{\rm n}$ Azimuthal correlations have also been measured
between jet constituents and the bulk. No significant correlation is found
between constituents of high-$p_{\scriptstyle\mathrm{T}}$
($p_{\scriptstyle\mathrm{T}}>15$ GeV) jets and the remaining hadrons in the
event in pp collisions186. On the contrary, attempts to isolate jet
constituent $v_{\rm 2}$ in jet-like three-particle correlation in
$\mathrm{p}$–$\mathrm{Pb}$ show significant non-zero values up to
$p_{\scriptstyle\mathrm{T}}\sim
8\,\mathrm{~{}\text{GeV}\\!/\mspace{-2.0mu}c}$187.

$v_{\rm n}$ in jets. In an attempt to identify collective phenomena in the
smallest systems, there has been a study of high-multiplicity high-energy
($p_{\scriptstyle\mathrm{T}}>550\,\mathrm{~{}\text{GeV}\\!/\mspace{-2.0mu}c}$)
jets with up to 100 charged constituents in the final state188. Two-particle
correlations and $v_{\rm 2}$ are extracted in a coordinate system defined
relative to the jet direction. The measured $v_{\rm 1}$, $v_{\rm 2}$ and
$v_{\rm 3}$ show a distinct trend as a function of multiplicity which is
reproduced by Monte Carlo simulations without collective phenomena except at
multiplicities above 80 where the trend changes and departs from the Monte
Carlo expectation. The authors interpret this result as possible collective
phenomena occurring in the dense environment of the jet.

Open heavy flavor and quarkonia Heavy flavor quarks (charm and beauty) are
produced in the initial hard scattering and then interact and rescatter with
the medium. Thus, also heavy-flavor particles exhibit finite anisotropies as
shown with non-zero $v_{\rm 2}$ measurements for $D$ mesons, J/$\psi$ and HF
decay electrons and muons in both $\mathrm{Pb}$–$\mathrm{Pb}$ and
$\mathrm{p}$–$\mathrm{Pb}$ collisions226, 227, 228, 229, 230, 231, 232, 233,
234, 235, 236, 237, 108, 238, 239, 240, 241, 242, 243, 244. An example is
shown in the left panel of Fig. 11. Non-zero $v_{\rm 2}$ of particles from
decays involving bottom quarks are measured in $\mathrm{Pb}$–$\mathrm{Pb}$ and
$\mathrm{p}$–$\mathrm{Pb}$ collisions but with large uncertainties240 in
$\mathrm{p}$–$\mathrm{Pb}$ collisions, see also the left panel of Fig. 11.
$v_{\rm 2}$ measurements for $\Upsilon$ are consistent with zero in
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions246, 248 and $\mathrm{p}$–$\mathrm{Pb}$
collisions245. In $\mathrm{pp}$ collisions, significant $v_{\rm 2}$ for charm
are measured at high multiplicity247, 240 while the $v_{\rm 2}$ for muons from
b decays is consistent with zero247, see the right panel of Fig. 11. Amongst
measurements of collectivity seen in $\mathrm{Pb}$–$\mathrm{Pb}$ collisions
but not yet in smaller systems, there are data on non-zero $v_{\rm 2}$ for
$\psi(2S)$ as well as non-zero $v_{\rm 3}$ for $D$ mesons and J/$\psi$.238,
243, 244

Figure 11:  Left panel: $v_{\rm 2}$ coefficient of prompt $D^{0}$ mesons and
J/$\psi$ as well as non-prompt $D^{0}$ compared to strange particles as a
function of $p_{\scriptstyle\mathrm{T}}$ in high-multiplicity
$\mathrm{p}$–$\mathrm{Pb}$ collisions (Figure adapted from Ref. 240). Right
panel: $v_{\rm 2}$ coefficient of muons from charm and beauty decays as a
function of multiplicity in pp collisions (Figure from Ref. 247).

The nuclear modification factor $R_{AA}$ for $J/\psi$ in large systems is
enhanced at LHC with respect to RHIC energies249, 250, 257, 259, 261, 232,
267, 273, 276, 282, 283. This is qualitatively different from the
$\sqrt{s}$-dependence of $R_{AA}^{J/\psi}$ from CERN SPS to RHIC energies.
J/$\psi$ suppression has been expected due to the effect of the medium on the
J/$\psi$ binding energy, while the enhancement at LHC energies can be
incorporated as a J/$\psi$ regeneration component from deconfined charm
quarks282. $\Upsilon(1S)$, $\Upsilon(2S)$ and $\Upsilon(3S)$ production is
equally found suppressed in $\mathrm{Pb}$–$\mathrm{Pb}$ collisions with
respect to the pp reference250, 251, 272, 274, 277, 280, 284.

In $\mathrm{p}$–$\mathrm{Pb}$ collisions, J/$\psi$ and $\Upsilon(1S)$ are
suppressed relative to $\mathrm{pp}$ collisions252, 254, 256, 265, 266, 269,
270, 271, 275, 279, 281 which can be attributed to nuclear modification of the
gluon PDF171. The production of the excited charmonium state, $\Psi(2S)$ as
well as excited bottomonium states $\Upsilon(nS)$, $n\geq 2$ have been
measured in both $\mathrm{Pb}$–$\mathrm{Pb}$ and $\mathrm{p}$–$\mathrm{Pb}$
collisions253, 254, 255, 257, 258, 260, 262, 263, 264, 268, 273, 278, 281, 284
which shows a suppression w.r.t. the ground state.

###  3.4 Towards even smaller systems…

Table 1: Summary of observables measured in $\gamma$–p (realized in $\mathrm{p}$–$\mathrm{Pb}$ and ep colliders), $\gamma$–Pb, low multiplicity pp, $\mathrm{e^{+}e^{-}}$ and ep collision systems. Low $N$ refers to values close to the minimum-bias multiplicity. | Observable
---|---
|  Near-side |  Azimuthal |  Higher-order
System |  ridge yields |  anisotropy ($v_{\rm n}$) |  cumulants
$\gamma$–Pb |  observed (template) |  up to $v_{\rm 3}$ (template) |  not studied
(in $\mathrm{Pb}$–$\mathrm{Pb}$ UPC) |  |  |
$\gamma$–p |  not observed |  up to $v_{\rm 3}$ $\approx$ MC |  not studied
(in $\mathrm{p}$–$\mathrm{Pb}$ UPC) |  |  (template) |
$\gamma$–p |  not studied |  up to $v_{\rm 2}$ $\approx$ MC |  n = 4 $\approx$ MC
(in ep) |  |  |  wrong sign
ep (DIS) |  not studied |  up to $v_{\rm 2}$ $\approx$ MC |  not studied
pp (low $N$) |  observed |  up to $v_{\rm 2}$ (template) |  not studied
ee |  limits + hint at $1.02\sigma$ |  up to $v_{\rm 3}$ $\approx$ MC |  not studied
Refs. |  295, 296, 297, 298, 100, 299 |  46, 83, 300, 301 302, 298, 299 |  301

Studies of the mentioned phenomena have been extended into even smaller
systems, by studying low multiplicity pp collisions and by using ultra-
peripheral (UPC) heavy-ion collisions where $\gamma$–p and $\gamma$–Pb
processes occur. Furthermore, data of former $\mathrm{e^{+}e^{-}}$ and ep
experiments have been re-analyzed where the latter gives access to deep
inelastic scattering (DIS) as well as resolved photoproduction resembling
$\gamma$–p processes. While these studies only cover a fraction of the
observables discussed in this section, they provide a number of interesting
findings. Table 1 summarizes the current status and a brief overview is given
in the following.

Figure 12:  $v_{\rm 2}$ coefficient measured in ultra-peripheral $\gamma$–Pb
collisions compared to the measurement in $\mathrm{p}$–$\mathrm{Pb}$ and pp
collisions (Figure adapted from Ref. 302).

Figure 13:  Near-side ridge yield measured in pp collisions compared to upper
limits set in $\mathrm{e^{+}e^{-}}$ collisions (Figure from Ref. 100).

Ultraperipheral collisions The electromagnetic cloud traveling with Pb ions in
a particle accelerator leads to collisions of photons with beam particles.
This allows one to study $\gamma$–p ($\gamma$–Pb) collisions in a
$\mathrm{p}$–$\mathrm{Pb}$ ($\mathrm{Pb}$–$\mathrm{Pb}$) collider
configuration in so-called ultraperipheral collisions. Although the study of
$\gamma$–p processes does not show any near-side ridge structure298, azimuthal
Fourier coefficients up to third order have been measured. The positive
$v_{\rm 2}$ is reasonably described by simulations without involving
collective phenomena298. The study of $\gamma$–Pb collisions shows a different
picture: significant non-zero $v_{\rm 2}$ and $v_{\rm 3}$ are measured, albeit
with the template subtraction method. The $v_{\rm n}$ are lower than the
equivalent measurement at the same multiplicity $\mathrm{p}$–$\mathrm{Pb}$
collisions but similar to the ones in pp collisions, see Fig. 13. Contrary to
the $\gamma$–p result, simulations without a collective component do not
reproduce the result, even qualitatively302.

ep collisions Studies in ep collisions in the DIS region find no evidence of a
$v_{\rm 2}$ signal beyond the one expected from momentum conservation and are
reproduced by MC generators without involving collective phenomena300, 301. In
the photoproduction region $v_{\rm 2}$ is extracted with 2 and 4-particle
correlations. No evidence of effects beyond momentum conservation are
observed301. The 4-particle signal at second order has the opposite sign301 to
the one observed in high-multiplicity pp to $\mathrm{Pb}$–$\mathrm{Pb}$
collisions. The authors of Ref. 301 report their $v_{\rm 2}$ signal to be
consistent with the one in $\gamma$–Pb collisions302 and attribute it to jet
production unrelated to hydrodynamical behavior.

Low-multiplicity pp collisions While high-multiplicity ridge yields in pp
collisions have been the first observation of collective phenomena in pp
collisions, see section 3.1, studies at low multiplicity are intrinsically
difficult as the signal scales with the square of the particle multiplicity
and is dominated by jet fragmentation and resonance decays at low
multiplicity. Recently, near-side ridge yields have been measured with great
precision to multiplicities lower than the minimum-bias multiplicity100.
Interestingly, significant long-range correlations emerge already at the
minimum-bias multiplicity, see Fig. 13. Studies beyond ridge yields, for
instance of $v_{\rm n}$ coefficients are scarce at low multiplicity. Only a
$v_{\rm 2}$ measurement exists in pp collisions just above the minimum-bias
multiplicity46, 83 applying a template subtraction. Higher-order cumulants are
only extracted at a multiplicity of 3–4 times the minimum-bias multiplicity at
present.

$\mathrm{e^{+}e^{-}}$ collisions These collisions are free of hadronic effects
in the initial state. As such, several mechanisms for collectivity invoked in
other collision systems may be difficult to realize in $\mathrm{e^{+}e^{-}}$.
This makes it particularly interesting to check whether or not precursors of
collective effects seen in other collision systems extend to
$\mathrm{e^{+}e^{-}}$. These collisions have a randomly oriented event
topology in the detector, therefore the mentioned studies are performed in the
so-called thrust axis which is a proxy for the direction of the hard process
occurring in e.g. $\mathrm{e^{+}e^{-}}\rightarrow q\bar{q}$. Archived ALEPH
data as well as data from Belle have been analyzed and limits on near-side
ridge yields have been set at collisions energies of about 10.5 GeV296, 297,
91 GeV295 and at 183–209 GeV299. In the highest multiplicity studied at
183–209 GeV a non-zero near-side ridge yield has been found at
$1.02\sigma$299, see rightmost ALEPH point in Fig. 13. A special case is the
study of $\Upsilon(4S)$ on-resonance data. Here a non-zero near-side yield is
found which is however attributed to the special topology in the decay of the
$\Upsilon\rightarrow B\bar{B}$ system and reproduced by Monte Carlo
simulations without collective phenomena. Furthermore, $v_{\rm n}$
coefficients are extracted and compared to MC simulations without collective
phenomena which agree except at the largest studied multiplicities299. The
statistical precision does not yet allow the conclusion that such phenomena
are also found in $\mathrm{e^{+}e^{-}}$ collisions.

This data allows for a clear comparison to previously mentioned low-
multiplicity ridge yields measured in pp collisions. The latter are at least
$5\sigma$ larger100 than the ones in $\mathrm{e^{+}e^{-}}$ collisions 666To be
specific: The multiplicity range from 8 to 24 charged particles contributes
most to this global assessment since ALEPH data allow one to establish the
tightest upper bounds in this region. At high multiplicity, the upper bounds
are much weaker. While this leaves more room for a potential signal, it must
not be interpreted as a signal per se. It needs to be mentioned that the
comparison as a function of multiplicity of these different collision systems
at very different center-of-mass energies is intrinsically difficult and
requires modeling assumptions resulting in a 5–6.3$\sigma$ difference100., see
Fig. 13. This demonstrates that additional phenomena are observed in low-
multiplicity pp collisions compared to $\mathrm{e^{+}e^{-}}$ collisions.

##  4 Theory overview

In the phenomenological discussion of ultra-relativistic nucleus–nucleus
collisions, it has been emphasized repeatedly that fluid dynamic models with
close-to-minimal dissipative properties can account for a broad range of
collective phenomena303. Data comparisons with iEBE-VISHNU304, superSONIC305,
MUSIC306 and other fluid models307, 308 support this conclusion. However, the
physics encoded in such phenomenologically successful “fluid” models goes
significantly beyond fluid dynamics777Here and in the following, we put the
notion “fluid” in quotation marks to highlight - where required in content -
that “fluid” models invoke more than fluid dynamics.. It includes suitably
chosen initial conditions, hadronization prescriptions, phenomenological
choices for kinetic transport in the pre-equilibrium dynamics, and hadronic
rescattering. Also the modeling of the fluid dynamic phase itself includes
more than fluid dynamics, as explained in section 4.1.1. Early efforts to
extend such fluid modeling to small collision systems are reviewed in section
4.3. The current state of the art in model comparisons is set by Bayesian
inference techniques 309.

The earliest Bayesian analyses of fluid models did not include any data from
small systems and they were based on relatively few $p_{T}$-integrated
measurements 310, 311, 312, 313, 314. Qualitatively different initial
conditions were also studied315. Early analyses that included some
$\mathrm{p}$–$\mathrm{Pb}$ data remained restricted to $p_{T}$-integrated
quantities 316. The most complete analysis as of today is given by
Trajectum317, 318. It includes more than 500 data points on
$\mathrm{p}$–$\mathrm{Pb}$ and $\mathrm{Pb}$–$\mathrm{Pb}$ to overconstrain
$O(20)$ model parameters, it extends to $p_{T}$-differential information and
it has been used to illustrate the physics opportunities of future light ion
beams by generating mock data for oxygen–oxygen collisions at the LHC 319.

Remarkably, the inclusion of data on small collision systems has led so far to
no significant additional constraints on the dynamical mechanisms and material
properties invoked to describe successfully central nucleus–nucleus
collisions. The reason for this is thought to lie in increased uncertainties
associated to the initial conditions of small systems and in an increased
sensitivity to final state effects, see section 4.3.2. Within these complex
multi-parameter models, our currently limited understanding of such
uncertainties seems to offset the discriminatory power which a wider variation
of system size could provide 320, 317, 318. A similar statement applies to
other complex simulation tools of ultra-relativistic heavy ion collisions such
as A Multi Phase Transport (AMPT) model 321, 322, 323, BAMPS 324, 325 or EPOS
326, 327. We summarize in section 4.3 to what extent these codes have been
compared to measures of collectivity in small systems, but we mention already
here that more systematic studies are missing and that no significant further
constraints on the dynamical mechanisms underlying collectivity have been
reported.

The situation is thus rather peculiar: On the one hand, experiments at the LHC
have established the persistence of almost all soft
$p_{\scriptstyle\mathrm{T}}$ signatures of collectivity across all system
sizes from central $\mathrm{Pb}$–$\mathrm{Pb}$ down to the low multiplicity
$\mathrm{pp}$ collisions in which collectivity had not been expected. On the
other hand, this wealth of qualitatively novel data does not seem to challenge
any of the complex model frameworks used traditionally to simulate soft
physics in heavy-ion collisions. The conclusion here is not that fluid-dynamic
or transport models of ultra-relativistic heavy-ion collisions naturally
explain collectivity in small systems. The conclusion is rather that the
predictivity of these models for small systems is sufficiently weak to
accommodate existing data.

In this situation, a discussion of what we learn from small system
collectivity profits from going beyond mere data-model comparisons. One should
ask whether the experimentally established phenomena of collectivity are more
generic than the concepts invoked for their explanation. To be specific:

  1. I

The phenomenology of nucleus–nucleus collisions invokes the concept of rapid
equilibration – but can this equilibration be rapid enough to explain
collectivity down to $\mathrm{pp}$?

  2. II

The phenomenology of nucleus–nucleus collisions invokes fluid dynamics – but
what is the smallest droplet to which fluid dynamics applies?

  3. III

The phenomenology of nucleus–nucleus collisions invokes close-to-minimal
dissipative properties as the sine qua non for numerically large collective
effects – but can this be the main cause for collectivity in
$\mathrm{p}$–$\mathrm{Pb}$, in $\mathrm{pp}$ or even in $\mathrm{e^{+}e^{-}}$?

The field of ultra-relativistic heavy-ion collisions knows of precedents where
substantial insights were achieved by going beyond a mere model-data
comparison. For instance, two decades ago, a calculation of the viscosity over
entropy ratio in a maximally supersymmetric field theory had contributed to
revolutionizing heavy-ion phenomenology not because it led itself to a direct
data-model comparison, but because it indicated that a hitherto unimagined
dynamics could be realized in quantum field theory. In the same spirit, the
LHC discovery of collectivity in small systems has triggered since 2010 many
studies in theoretically clean simplified set-ups to understand what is
imaginable in terms of fast hydrodynamization, fast kinetic thermalization and
fast chemical equilibration, and to what extent fluid dynamic or non-fluid
dynamic degrees of freedom are at work in such fast processes. These lines of
research are an intrinsic part of the LHC legacy of small system collectivity
even if they address often “only” the conceptual and not the phenomenological
part of the problem. The present chapter summarizes the most pertinent
findings of this broad theoretical research effort. Finally, the last
subsection will discuss data-model comparisons.

###  4.1 Dynamical frameworks of collectivity

Fluid dynamics, kinetic theory and the evolution of initial-state correlations
(as well as various non-perturbative mechanisms) are the dynamical frameworks
considered for the explanation of collective phenomena in ultra-relativistic
nucleus–nucleus collisions. Here, we recall the assumptions on which they are
based. This will enable a discussion about their applicability to collective
effects in the smallest collision systems.

####  4.1.1 Fluid dynamic descriptions of collectivity

In nuclear physics, first fluid dynamic formulations of collectivity can be
traced back to the liquid-drop model in the 1930s – a very informative account
of these earliest developments is given in Ref. 40. In ultra-relativistic
heavy ion collisions, it was Bjorken’s work 328 on longitudinally boost-
invariant systems that initiated fluid dynamic formulations of the collective
expansion.

Fluid dynamics is formulated in terms of thermodynamic quantities, and it
emerges in the long wavelength limit of any more complete dynamics. Different
derivations of relativistic fluid dynamics exist.

  1. 1.

The thermodynamic (or entropy-wise) derivation
Fluid dynamics can be derived directly from conservation laws in a gradient
expansion. In the simplest case (expansion to zeroth order in gradient), the
energy-momentum tensor
$T_{\mu\nu}=(\varepsilon+p)u_{\mu}u_{\nu}-p\,g_{\mu\nu}$ depends only on
energy density $\varepsilon$, pressure $p$ and three independent components of
the local flow field $u_{\mu}$. The collective dynamics of these five unknowns
is fully determined in terms of four constraints from energy-momentum
conservation $\nabla_{\mu}T^{\mu\nu}=0$ supplemented by the equation of state
$\varepsilon=\varepsilon(p)$.
To first order in gradients, the tensor decomposition of $T_{\mu\nu}$ with
respect to the flow field $u_{\mu}$ contains a shear viscous term

$\Pi_{\mu\nu}^{\rm
constitutive}=\eta\left[\Delta^{\mu\alpha}\nabla_{\alpha}u^{\nu}+\Delta^{\nu\alpha}\nabla_{\alpha}u^{\mu}-\tfrac{2}{3}\Delta^{\mu\nu}\nabla_{\alpha}u^{\alpha}\right]\,,$
(11)

as well as a bulk viscous and a heat conducting one. Here, $\eta$ denotes the
shear viscosity and $\Delta^{\mu\nu}$ is the projector on the subspace
orthogonal to $u^{\mu}$. Remarkably, requiring that entropy increases locally
($\nabla_{\mu}S^{\mu}\geq 0$) closes the dissipative fluid equations to first
order in gradients 329. In this sense, dissipative fluid dynamics results
solely from the fundamental laws of thermodynamics supplemented by material
properties ($\varepsilon=\varepsilon(p)$, $\eta$, …) that are calculable in an
equilibrium field theory.
However, implementing the spatial gradients $\nabla_{\mu}u^{\alpha}$ in (11)
instantaneously would violate causality. Therefore, causal viscous fluid
dynamics888Other causal formulations of fluid dynamics are possible 330, 331
but these have not been employed in phenomenological studies so far. specifies
how the physical shear viscous tensor $\Pi_{\mu\nu}$ relaxes to its first-
order constitutive form $\Pi_{\mu\nu}^{\rm constitutive}$. This is done by ad
hoc assumptions that lie outside the realm of fluid dynamics 999 As this
relaxation dynamics must preserve transversality and orthogonality of
$\Pi_{\mu\nu}$, its precise formulation is a bit technical and it will not be
given here. For details, see e.g. Eq. (3.12) of Ref. 18.. Israel-
Stewart(IS)-type formulations of causal viscous fluid dynamics invoke a simple
relaxation time $\tau_{\Pi}$; the temporal decay of a shear viscous excitation
takes then the form

$A_{\rm excitation\,decay}(t)\simeq c_{\rm
hyd}\exp\left[-\frac{\eta}{\varepsilon+p}k^{2}\,t\right]+c_{\rm non-
hyd}\exp\left[-\frac{t}{\tau_{\Pi}}\right]\,.$ (12)

This time dependence reflects the precise statement that in IS fluid dynamics,
the shear viscous channel of the retarded propagator of $T_{\mu\nu}$ has poles
at $-i\tfrac{\eta}{\varepsilon+p}k^{2}$ and $-i\tfrac{t}{\tau_{\Pi}}$, see
e.g. Ref.27 for details. Since hydrodynamics is a gradient expansion, the
lifetime of fluid dynamic excitations must grow with the inverse of the
squared wave number $k^{2}$. The $k$-independent decay with lifetime
$\tau_{\Pi}$ is thus explicitly non-fluid dynamic. Israel-Stewart dynamics is
fluid dynamics only for sufficiently large wavelengths $\lambda$

$\frac{1}{\tau_{\Pi}}>\frac{\eta}{\varepsilon+p}k^{2}=\frac{\eta}{s}\frac{(2\pi)^{2}}{T\,\lambda^{2}}\,.$
(13)

  2. 2.

Derivation of relativistic fluid dynamics from transport theory
In weakly coupled theories, there is a scale separation between the typical
size of the wave packet $1/T$ and the mean free path between individual
scatterings. Therefore, for time separations much larger than $1/T$, when
interference effects between subsequent scatterings can be neglected, the
evolution is determined by Boltzmann transport in which collision kernels are
given by in-medium scattering processes in the field theory. In QCD, this is
the basis of AMY effective kinetic theory (EKT).332, 333
In general, Boltzmann transport theory evolves one-point probability densities
$f(x,p,t)$ in response to collision kernels. 101010The evolution of $f(x,p,t)$
is coupled to that of correlated $n$-point density functions in the so-called
BBGKY-hierarchy, and Boltzmann transport results from the truncation of that
hierarchy. Fluid dynamics can be derived by restricting Boltzmann transport to
the lowest momentum moments of $f(x,p,t)$, namely the local conserved currents
$N_{\mu}=\int p_{\mu}f(x,p,t)$ and the energy-momentum tensor $T_{\mu\nu}=\int
p_{\mu}p_{\nu}f(x,p,t)$. For relativistic dissipative fluid dynamics, this
derivation was pioneered by Israel and Stewart 334 (see also Refs. 335, 336,
337) and it was adapted to heavy-ion collisions subsequently 338, 339, 340,
341.
For fluid-dynamic excitations, the kinetic and entropy-wise derivations yield
the same equations of motion. However, the spectrum of non-fluid dynamic
excitations depends on how the collision kernel is modeled 342, 343, 344, 345
and it differs from that of Israel-Stewart dynamics. In the simplest
relaxation time approximation (RTA) of kinetic theory, the non-propagating
non-fluid dynamic excitation of IS fluid dynamics in (12) is replaced by
propagating quasi-particle excitations with lifetime
$\tau_{\Pi}=a\,\tfrac{\eta}{s\,T}$, $a=5$. In close analogy to IS dynamics,
this kinetic evolution will not be governed by fluid dynamic degrees of
freedom on length scales $\lambda$ that are too short to satisfy (13).

  3. 3.

Derivation of relativistic fluid dynamics from quantum field theory
In general, dissipative transport coefficients can be calculated via the
Green-Kubo formula as certain long-wavelength limits of correlation functions
in the thermal field theory 346, 347. For QCD, such transport coefficients
have been calculated in finite temperature perturbation theory to (almost)
next-to-leading order 348, and there are some exploratory studies in lattice
QCD 349, 350. However, nothing is known rigorously about the spectrum of non-
fluid excitations in QCD, except that they must exist, and that they may be
more complicated than a few well-isolated poles 351. In this sense, the
microscopic structure of the QGP is unknown to date.
In the strong coupling limit of non-abelian quantum field theories with
gravity dual, both fluid dynamics and the microscopic dynamics beyond fluid
dynamics is known. Dissipative fluid dynamics can be derived from the gradient
expansion of gravitational perturbations in the dual classical supergravity
theory 18, 352. String-theory inspired derivations have helped to identify the
complete set of allowed second-order transport terms 18, including transport
coefficients that can arise only in anomalous fluid dynamics 353 and that can
lead to qualitatively novel phenomena such as the chiral magnetic effect 354.
These anomalous transport coefficients can also be identified in an entropy-
wise derivation 355.
For ${\cal N}=4$ SYM theory in the limit of strong coupling
$\lambda=g^{2}N_{c}\to\infty$ and large number of colors $N_{c}\to\infty$, the
relaxation time in (12) is $\tau_{\Pi}=a\,\tfrac{\eta}{s\,T}$ with
$a=4-\log(4)\approx 2.61$ (see Ref. 25 for a more complete review of
relaxation times). The excitation spectrum of non-fluid (so-called quasi-
normal) modes is fully known in the strong-coupling limit 356, 357. The
question of how the non-fluid dynamic properties of these non-abelian plasmas
change with coupling strength is subject of recent research358, 359, 360, 361.

In summary, these findings support the following qualitative statements:

  * •

The fluid dynamic evolution equations for relativistic collective dynamics are
universal: they take the same form, irrespective of whether they are derived
from conservation laws (entropy-wise derivation)362, from specific kinetic
theories via Grad’s 14 moment expansion 363 or with string-theoretical
techniques in strongly coupled QFTs 18, 352.

  * •

Any dynamics limited to fluid dynamic excitations is acausal. Any causal
relativistic collective dynamics exhibits also non-fluid dynamic exitations.
These are non-universal and they depend on details of the microscopic
dynamics. As non-fluid dynamic exitations must exist, the relevant question is
not the qualitative one (to be a fluid or not to be a fluid?), but the
quantitative one: to what extent is collectivity dominated by fluid-dynamic or
non-fluid dynamic excitations?

  * •

Fluid-dynamic excitations dominate if they are more long-lived than the non-
fluid ones. For this to be the case, the wavelengths $\lambda$ of excitations
should satisfy 111111Conditions for the applicability of fluid dynamics have
been formulated alternatively in terms of Reynolds number and Knudsen number,
as reviewed succinctly in Ref. 41. Reynolds number is an intrinsically fluid
dynamic concept used in the wider rheological practice to discriminate laminar
from turbulent flow. Knudsen number measures system size in units of mean free
path and thus requires knowledge about which degrees of freedom propagate with
which cross section. Both concepts can be generalized41. Instead, we prefer
here to formulate the applicability of fluid dynamics directly in terms of the
decay time of non-fluid dynamic modes in (14), as this brings to the forefront
what is really experimentally at stage: to learn about the nature of the non-
fluid degrees of freedom which are characteristic features of the QGP, to
learn about their decay which governs hydrodynamization, and to be open to the
logical possibility that over a significant period in the evolution, non-fluid
degrees of freedom may be non-negligible and coexistent with the fluid dynamic
ones.

$\lambda>\sqrt{\frac{\eta}{s}\frac{(2\pi)^{2}}{T}\tau_{\Pi}}=(2\pi)\frac{\eta}{s}\frac{\sqrt{a}}{T}\,.$
(14)

For an almost perfect plasma with $\tfrac{\eta}{s}=\tfrac{1}{4\pi}$ and
varying $a$ from $a\approx 2.61$ (in ${\cal N}=4$ SYM) to $a=5$ (in kinetic
theory), one finds $\lambda>(0.81-1.12)\tfrac{1}{T}$. In particular,
$\lambda>1\,{\rm fm}$ for $T=200$ MeV. For a less perfect fluid, this bound is
more stringent.

According to these numerical estimates, the smallest $\mathrm{pp}$ and
$\mathrm{p}$–$\mathrm{Pb}$ collision systems have an initial transverse
extension into which fluid dynamic excitations of wavelength (14) barely fit.
For a predominantly fluid dynamic explanation of $\mathrm{pp}$ and
$\mathrm{p}$–$\mathrm{Pb}$ collisions, however, one would require that scales
significantly smaller than the transverse system size can be resolved by
excitations that propagate fluid dynamically. According to the estimate (14),
this must not be taken for granted. Apparently supportive of this parametric
statement, there is some evidence that if “fluid” modeling is applied to
smaller and smaller systems, then the fluid dynamic aspects of that modeling
become less and less relevant (see section 4.3.2 below).

####  4.1.2 Kinetic theory and transport models

The application of kinetic theory to ultra-relativistic heavy-ion collisions
is as old as that of fluid dynamics. Around the same time at which Bjorken
formulated boost-invariant fluid dynamics, Baym formulated boost-invariant
kinetic theory to understand “the approach to thermodynamic equilibrium” of
“the excitations present” 364. As indicated by this wording, kinetic theory
had been viewed from the very beginning as a tool for understanding how and on
what time scale equilibrium is achieved (“thermalization”) and how fluid
dynamic behavior emerges (“hydrodynamization”). The potential relevance and
the potential limitations of formulating small system collectivity in terms of
kinetic transport theory may be summarized as follows:

  * Pro

Unlike fluid dynamics, kinetic theory applies to arbitrarily small systems. It
is a candidate dynamics that interpolates smoothly between free-streaming in
sufficiently small and dilute systems (assumed to be realized in multi-purpose
event generators for $\mathrm{pp}$ collisions) and dissipative fluid dynamics
in sufficiently large and dense systems (assumed to be realized in fluid
simulations for $\mathrm{Pb}$–$\mathrm{Pb}$ collisions). These qualitative
features have motivated many recent studies of the onset of collective
phenomena in small collision systems (see section 4.2.1).

  * Con

The formulation of kinetic transport theory relies on a scale separation
between the size $\sim 1/T$ of quantum mechanical wavepackets and the mean
free path between individual scatterings. In the strong coupling limit of
quantum field theories with gravity duals, this scale separation is not
realized. In high temperature QCD, the scale separation is perturbatively
realized and it allows one to derive an effective kinetic theory, called AMY
EKT 332, 333. For phenomenologically relevant values of the running QCD
coupling (say $\alpha_{s}(Q^{2}=2{\rm GeV}^{2})\simeq 0.3$), a parametrically
clean scale separation does not exist. In this case, the use of transport
models may still promise qualitative insights but it pushes kinetic theory
beyond its region of guaranteed applicability.

In kinetic theory, collective dynamics arises from the microscopic dynamics
implemented in the collision kernels of the transport equations. In the
phenomenological practice, any formulation of this microscopic dynamics is
necessarily incomplete. Without entering what would be an inevitably technical
and model-dependent discussion of how different choices of collision kernels
are motivated, we highlight here two qualitative insights that have emerged
within the last decade:

  1. 1.

Rapid kinetic hydrodynamization and thermalization
QCD effective kinetic theory, AMY EKT 332, 333 propagates one-point
probability densities in response to in-medium $2\to 2$ elastic and $1\to 2$
LPM collision kernels. It implements the so-called BMSS bottom-up
thermalization 365 according to which a longitudinally Bjorken-expanding
initially oversaturated system dilutes and becomes undersaturated while $1\to
2$ splittings build up a reservoir of soft scattering centers. Subsequent
interactions with these soft scattering centers drive the system then to
equilibrium. Numerical studies with phenomenologically realized values of
$\alpha_{s}$ have demonstrated 366, 367 that the resulting evolution towards
equilibrium can be sufficiently fast to hydrodynamize initial conditions in
ultra-relativistic heavy ion collisions efficiently.
Beyond conceptual insights, this line of research has resulted by now in a
simulation tool. A linearized version of AMY EKT underlies the KoMPoST code
368, 369 that sets the current state of the art of modeling the pre-
equilibrium stage of “fluid” models. The duration over which pre-equilibrium
evolution deviates significantly from viscous fluid dynamics in AMY EKT is
comparable to the transverse extent of a $\mathrm{pp}$ collision (see
discussion of Figure 14 below). This is another indication that if “fluid”
models are applied to small systems, the fluid part of their evolution plays
an increasingly smaller role for smaller system size. It is in line with the
conclusions summarized at the end of section 4.1.1.
Beyond the linearized treatement of KoMPoST, there is a first proof-of-
principle study of a Lorentz invariant parton cascade ALPACA that solves the
AMY EKT Boltzmann equations 370.

  2. 2.

Onset of collectivity in “one-hit” kinetic theory
Evolving with $2\to 2$ Boltzmann transport an expanding mixture of several
relativistic massive particle species, it was found early on that a single
collision per particle in average can already lead to sizable elliptic flow,
with mass ordering between the species 371. Similar findings were reported for
the phenomenologically more complete AMPT code372, emphasizing that the
contribution from hydrodynamic-type mechanisms is small, see section 4.3.
In a sufficiently small and dilute system, the first correction to free-
streaming is obtained by expanding the collective dynamics to first order in
the collision term. This is an expansion to first order in opacity, where
opacity measures a specific combination of transverse geometric extension and
local transverse energy density 291, 373. In simple kinetic models
(isotropization time approximation – ITA), the one-hit approximation leads
naturally and for generic reasons to a hierarchy of the linear response
coefficients $v_{n}/\varepsilon_{n}$ in (1) that is in line with experimental
observations. Remarkably, also non-linear response coefficients result in one-
hit kinematics as a consequence of simple scaling behaviors290. In ITA kinetic
theory, non-linear response coefficients like $v_{4}/\varepsilon_{2}^{2}$ or
$v_{5}/(\varepsilon_{2}\varepsilon_{3})$ are found to have the
phenomenologically required order of magnitude in one-hit dynamics291, 373.

In contrast to fluid dynamics, the region of validity of kinetic theory is not
limited to quanta of low momentum. Kinetic theory is thus a framework in which
bulk dynamics and the jet quenching of hard processes could emerge as
different manifestations of the same underlying dynamics. This is particularly
apparent in AMY EKT which is known to give rise to fast hydrodynamization
while its $1\to 2$ Landau-Pomeranchuk–Migdal collision kernel is known 374 to
implement the BDMPS-Z formulation375 of radiative parton energy loss that
underlies many jet quenching formulations.

####  4.1.3 Initial state quantum interference as a source of collectivity

Both fluid dynamic and kinetic models attribute the observed phenomena of
collectivity to final state partonic and hadronic scattering. In contrast, it
has been suggested in the context of saturation physics that initial state
correlations could be responsible for the flow-like phenomena observed in
small systems. With the so-called IP-Glasma model376, 377, saturation physics
has motivated an ab initio formulation of subnucleon fluctuations. Models
within this framework were shown to reproduce di-hadron correlation data in
$\mathrm{pp}$ 378, 379 and $\mathrm{p}$–$\mathrm{Pb}$ 379 collisions, and in
various small collision systems studied at RHIC 380. However, a reassessment
of this saturation physics explanation in an independent model implementation
381 called into question some of the essential elements reported in Ref. 380.
In this context, the results of Ref. 380 were referred to as
“counterintuitive”.

A short succinct review of the initial state perspective on collectivity in
small collision systems can be found in Ref. 382. What is referred to as
“glasma graphs” in saturation physics can be viewed as the calculation of
quantum interference effects in the superposition of emission amplitudes from
different sources. It follows directly from Heisenberg’s uncertainty principle
that quantum interference can map spatially anisotropic distributions of these
initial sources into momentum anisotropies in the final state. In this sense,
non-vanishing elliptic flow from sources correlated in the initial state may
be intuitively expected even though the size of their contribution may be
difficult to calculate and even though current model implementations may be
unsatisfactory.

The phenomenology of ultra-relativistic heavy ion collisions knows of
collective effects that are of a pure quantum nature. In particular, HBT
correlations are collective in the sense that they are correlations amongst
all identical particles in an event and that they persist in the limit of a
large number of sources; they also translate spatio-temporal anisotropies into
momentum anisotropies293. Explicit calculations in saturation physics383, 379
have identified such HBT-like correlations. Different calculational set-ups
indicate that momentum anisotropies $v_{n}$ resulting from quantum
interference can persist in the limit of a large number of sources384, 385,
386 irrespective of whether the conditions for saturation physics are
satisfied384, 387. However, the size of possible effects depends on
assumptions, e.g., about the orientation of color fields 388. Also, saturation
physics alone cannot be expected to explain all important manifestations of
collectivity. In particular, the formalism tends to lead to relatively short-
range correlations in rapidity389 and in transverse momenta388, it decreases
with increasing system size390, 380 and it can be washed out by final-state
interactions320, 391. In addition, it does not address the observed mass
ordering of collective flow and it does not explain naturally the system size
dependence of flow from small to large systems. In summary, the observed
collective phenomena are clearly more generic than their explanation within
the framework of saturation physics although, from a theoretical viewpoint,
quantum interference can be expected to contribute to momentum anisotropies at
some (as yet not fully understood) scale.

####  4.1.4 Strings, ropes, shoving, pomerons and all that

The collective phenomena observed in small systems relate to physics at small
momentum transfers which is nominally non-perturbative. Fluid dynamic
descriptions aim at capturing this non-perturbative physics in terms of a few
transport coefficients; saturation physics aims at capturing it by doing
perturbation theory on top of a non-perturbative high-density background;
kinetic theory views non-perturbative information as included in the collision
kernels. All three frameworks have important limitations if applied to small
systems. Also, they do not exhaust the range of non-perturbative concepts that
can be motivated within QCD. This subsection reviews other non-perturbative
approaches that have been explored in the context of small system
collectivity.

In QCD, partonic scatterings separate color charges. Separated color charges
are connected via flux tubes that fragment eventually. The description of this
non-perturbative dynamics is model-dependent, but it is an essential aspect of
QCD. In the so-called Lund string, it is at the basis of a successful
hadronization model 60 that is employed in modern multi-purpose event
generators 9. It has been argued that this string fragmentation must be
supplemented by thermal effects to account for the observation of heavy-ion
like behavior in the hadrochemical composition of high-multiplicity pp
collisions at the LHC 37. In inelastic hadronic or nuclear collisions, flux
tubes stretch between partons at projectile and target rapidities. With any
interaction that shifts the distribution of such flux tubes in the transverse
plane, these strings become candidates for long-range rapidity correlations
392, 393. Such a string shoving mechanism was shown to implement hadrochemical
and kinematic signatures of collectivity in small systems 394. A related
mechanism has been suggested but not worked out in Ref. 395. Even the founding
father of fluid dynamic applications to ultra-relativistic heavy ion
collisions has argued that while fluid dynamics may account for the observed
ridge in nucleus–nucleus collisions, the apparently similar ridge effect in
high-multiplicity pp collisions could be due to “a quite different physical
mechanism”, namely the collision of aligned flux tubes 396.

In general, elevating strings and more complex, spatially extended colored
structures (such as ropes, baryon junctions etc.) to dynamical degrees of
freedom has a long history in the phenomenological modeling of nucleus–nucleus
collisions. Early implementations include the dual parton model 397, RQMD 398,
URQMD 399 and extensions of HiJING 400. Also AMPT enhances its initial energy
density signficantly with a string melting mechanism 401 and it is thus more
than a parton cascade. Another non-perturbative concept that emerges naturally
in QCD is that of pomeron exchanges 402. Such mechanisms have been argued to
give rise to fluid dynamic behavior 403.

###  4.2 Hydrodynamization and Thermalization - insights from simple models

In general, collectivity is thought to arise from far out-of-equilibrium
initial conditions in a dynamics which thermalizes locally and which thus
lends itself eventually to a fluid dynamic description. But in a smaller
collision system, any collective dynamics is at work for a shorter time. In
the following, we summarize efforts to understand in simplified model set-ups
the central question I whether thermalization and hydrodynamization can occur
sufficiently quickly for explaining collectivity in the smallest collision
systems.

####  4.2.1 Studies in boost-invariant 1+1 dimensional systems

Figure 14:  Evolution of the ratio of longitudinal pressure over energy
density in different 1+1D models. Left panel: All microscopic theories (dashed
curves) have the same late-time attractor which is relativistic viscous fluid
dynamics. All microscopic theories hydrodynamize (i.e. approach this
attractor) on the same time scale $\sim(4\pi)\tfrac{\eta}{sT}$, including
variations of initial conditions within the same model (grey lines). Figure
taken from Ref. 404. Right panel: Not all microscopic theories have an early
time attractor. For the strongly coupled QFT shown here, initial conditions
set at five different times and two different choices all hydrodynamize at the
time scale $\sim(4\pi)\tfrac{\eta}{sT}$ without displaying universal behavior
at earlier times. Figure taken from Refs. 405. The hydrodynamization scale
identified in these studies is consistent with the simple parametric estimate
in (14).

A dynamical system with boost-invariant initial conditions evolves in a boost-
invariant way. Also, if initial conditions are translationally invariant in
the transverse plane, this symmetry is preserved throughout the evolution.
Invoking both symmetries, the resulting dynamics is essentially
1+1-dimensional. As first realized by Bjorken 328 for hydrodynamics and by
Baym 364 for kinetic theory, this offers the possibility of formulating
technically simple, analytically tractable toy models that share important
commonalities with physically realized almost boost-invariant collision
systems.

Timescale of hydrodynamization: Due to the symmetries in 1+1D systems, a very
limited set of out-of-equilibrium conditions can be explored: Transverse
gradients and density fluctuations are forbidden, conditions with initial
eccentricity cannot be initialized and the dynamical mapping (1) from
$\epsilon_{n}$ to $v_{n}$ cannot be explored. However, initial conditions that
deviate locally from equilibrium can be imposed. In fluid dynamics 406, 407,
408, 405 (in the absence of conserved currents), the only observables are the
longitudinal pressure $p_{L}$ and the energy density $\varepsilon$; both
depend only on time, and their ratio can be initialized arbitrarily far away
from the equilibrium value $p_{L}/\varepsilon=\tfrac{1}{3}$. In Boltzmann
kinetic theory 409, 407, 410, 411, 405, 412, also higher moments of the one-
particle distribution function $f$ can be studied. However, the timescales for
their relaxation are rather trivially related and numerically comparable to
those governing $p_{L}$ and $\varepsilon$ 410, 413, 405. In numerically
tractable, strongly coupled field theories with gravity dual 407, 414, 415,
405, the evolution is not limited to a one-particle probability distribution
and a larger variety of initial conditions may be initiated 405.

Technically, the question of hydrodynamization amounts to asking on which
timescale arbitrarily initialized components of the energy-momentum tensor
relax to their constitutive expressions, (such as (11)) in a fluid dynamic
gradient expansion. In the simplest 1+1D models (which have only $\varepsilon$
and $p_{L}$ as observables), this reduces to the question on which time scale
$p_{L}/\varepsilon$ reaches the viscous hydrodynamic correction (see black
line in the left panel of Fig. 14) to the equilibrium value
$p_{L}/\varepsilon=\tfrac{1}{3}$. The answer is: in a Bjorken expanding
system, hydrodynamization occurs on timescales
$O\left((4\pi)\tfrac{\eta}{sT}\right)$. This is so far all microscopic
dynamics explored in Fig. 14, irrespective of whether it is strongly coupled
or weakly coupled. This answer fully supports the parametric estimate (14) for
the limited range of validity of fluid dynamics. In this sense, the problem of
hydrodynamization is fully solved in 1+1D.

Early-time and late-time attractors: In the recent literature on
hydrodynamization, the more general concept of attractor solutions plays an
important role, see Ref. 416 for a broad review with relations to many
neighboring fields. Attractors are the particular solutions of a dynamical
system to which arbitrary initial conditions within the bassin of attraction
relax at sufficiently late times. The late-time hydrodynamic attractor is
clearly visible in Fig. 14 and it is reached for arbitrary initial conditions.
Remarkably, in 1+1D boost-invariant fluid dynamics and kinetic transport,
attractor solutions are found not only at late times, but also at arbitrarily
early times 409, 407, 410, 411, 405, 412. However, the existence of an early-
time attractor is not the tell-tale sign of a particularly efficient
hydrodynamization. Rather, it is a consequence of the interplay between local
interaction and the rapid longitudinal expansion in boost-invariant systems
405. In kinetic theory, this is a non-interacting free-streaming dynamics that
drives the system to negligible $p_{L}$, in fluid dynamics, the early time
attractor risks being unphysical since longitundinal gradients become
arbitrarily large at arbitrarily early time and they thus lie outside a fluid
dynamic description. In strongly coupled quantum field theories, no early time
attractor exists and the system hydrodynamizes on the natural time scale
$O(1/T)$ on which interactions drive the system to local equilibrium, see Fig.
14. Related findings have been corroborated in QCD effective kinetic theory
without 367, 368 and with 417, 418 quark degrees of freedom. In summary, only
the late-time attractor is independent of the microscopic dynamics, i.e.
universal. It reflects the onset of fluid dynamics. In contrast, the early-
time attractor is non-universal and it is not related to the onset of fluid
dynamics. Fluid dynamics is not unreasonably effective. Rather, it sets in
when standard textbook arguments expect it to set in, namely on scales on
which fluid dynamic excitations live longer as non-fluid dynamic ones (see
discussion leading to (14)).

The approach of particular initial conditions to known hydrodynamic late-time
attractor solutions has also been explored to exhibit the asymptotic nature of
the hydrodynamic gradient expansion, to develop resummation techniques 419,
420, 421, 343, 422, 411, 405 and it has been used to illustrate the role of
non-hydrodynamic poles in the relaxation to the attractor 415.

####  4.2.2 Studies in boost-invariant higher-dimensional systems

The fact that arbitrary initial conditions relax to a universal late-time
fluid-dynamic attractor solution is equivalent to saying that initial
conditions are forgotten. In the boost-invariant 1+1-dimensional systems
discussed so far, essentially all measurable quantities forget initial
conditions on a hydrodynamization time scale which is significantly shorter
than the thermalization time scale at which the system approaches local
equilibrium (here $p_{L}/\varepsilon=1/3$), see Figure 14. Having identified
this universal hydrodynamization time scale in multiple models based on
qualitatively different physics assumptions is arguably one of the main
insights of these 1+1-dimensional studies. However, while thermalization is
about forgetting, hydrodynamization is not. A fluid dynamic evolution that
implements the mapping (1) from initial spatial eccentricities $\epsilon_{n}$
to final momentum anisotropies $v_{n}$ is a dynamics that remembers the
initial conditions $\epsilon_{n}$. The fact that the evolution of 1+1D models
forgets any initial condition at late times is not a general physics feature
but it points to an intrinsic limitation of 1+1D studies: these models simply
do not allow for the description of measurable quantitites that are not
forgotten in the evolution. Higher dimensional systems need to be studied to
this end.

While most studies of hydrodynamization and thermalization are carried out in
1+1D models because of their technical simplicity, there have been studies in
higher dimensions 411, 423, 27, 405, 424, 425, 426, 29, 427. Early numerical
studies of radially symmetric colliding show waves mimicking small systems in
strongly coupled field theory corroborated the importance of the time and
length scale $O(1/T)$ for driving small systems toward fluid dynamic
evolution428, 429, 430. More generally, once radial symmetry is not imposed,
one can study how efficiently initial conditions $\epsilon_{n}$ are
remembered, and to what extent the response depends on the assumed dynamics.
For instance, evolving a system with initial spatial eccentricity
$\epsilon_{2}$, first with kinetic theory up to a switching time $\tau_{s}$
and then with the corresponding Israel-Stewart dynamics to later times informs
us about the sensitivity of collective response to the nature of non-fluid
dynamic degrees of freedom. In such model studies, one third to one half of
the elliptic flow signal seen in fully hydrodynamized systems can be built up
in small collision systems that extend over only one mean free path and that
do not hydrodynamize26. Along a different line of investigation, there has
been a significant computational effort to quantify in simple conformal
kinetic theories the full opacity dependence of all relevant linear and non-
linear response coefficients in the mapping
$\\{\varepsilon_{n}\\}\longrightarrow\\{v_{n}\\}$ 27, 424, 425, 29, 427. This
opacity dependence interpolates smoothly between free-streaming and almost
perfect fluid dynamics as a function of opacity. Despite their simplicity,
these models have led to statements about the (in)applicability of
hydrodynamics in $\mathrm{pp}$, $\mathrm{p}$–$\mathrm{Pb}$ and light nucleus
collisions28.

###  4.3 Model studies and data comparisons

Up to this point, our discussion has focused mainly on the description of
basic physics concepts that may underlie the formulation of collectivity in
small systems (section 4.1), and on the illustration of these concepts in
highly simplified model studies whose conceptual transparency offers unique
opportunities for understanding how and how efficiently collectivity may be
built up in small systems. However, such simple models do not reflect the full
complexity of ultra-relativistic hadronic and nuclear collisions. Here, we
finally summarize data comparisons with phenomenologically more complete model
frameworks.

####  4.3.1 Theory expectations prior to LHC data

Prior to the experimental observations of collectivity in small systems,
different works had contemplated the possibility that azimuthal asymmetries
may be seen in $\mathrm{pp}$ collisions at LHC as a consequence of various
nuclear-like effects 431, 432, as a consequence of fluid dynamic behavior 433,
434, 435, 436 or more generally as a consequence of a collective dynamical
response to initial spatial eccentricities in the $\mathrm{pp}$ overlap 437,
438. The very idea that $\mathrm{pp}$ collisions may develop fluid dynamic
behavior is not new and had been considered early on, see e.g. Ref. 439.
However, predictions of the size of $v_{2}$ in $\mathrm{pp}$ varied from zero
433 to values of $O(0.1)$, depending on model assumptions, and mechanisms
ranged from collective explosion 395 to initial state effects in saturation
physics 440. In short, these early works formulated interesting expectations,
but none of them constitutes a well-motivated prediction that can be regarded
as having received detailed experimental support.

####  4.3.2 Small system collectivity in fluid models

Following the phenomenological success in ultra-relativistic nucleus–nucleus
collisions, several works aimed at including small-system collectivity within
the same “fluid” modeling in small and large systems. This working hypothesis
“one fluid to rule them all” has been supported by a viscous fluid dynamic
description of high-multiplicity $\mathrm{pp}$, $\mathrm{p}$–$\mathrm{Pb}$ and
$\mathrm{Pb}$–$\mathrm{Pb}$ collisions at $\sqrt{s}=5.02$ TeV 441, 442.
However, it also became clear that “one fluid might not rule them all” 443.
While the wordings cited here sound like contradicting each other, the
scientific facts reported in these studies are consistent with each other. The
main message is simply that “fluid” models are more than fluid models: they
include additional physics. Additional physics found numerically relevant for
small collision systems includes the physics of hadronization (recombination,
see section 2.2 for details) 444, the physics of hadronic445, 446 non-fluid
contributions to flow, and the physics of initial state fluctuations447, 320,
448, 449. These additional physical mechanisms appear to become numerically
more relevant for smaller collision systems. In this sense, fluid dynamics
does not become abruptly invalid below a certain system size, but it becomes
gradually less and less relevant for the description of the system. The
studies of “fluid” models cited above compare to a small subset of the
measures of collectivity in Table LABEL:table:smallsystems only (mainly
measurements of $v_{2}(p_{T})$ and $v_{3}(p_{T})$). Where more detailed sets
of data were included in the analysis, issues of fine-tuning and of a possibly
limited range of validity of hydrodynamics became more important 443, 450.

The EPOS event generator employs a core-corona model. The initial conditions
are generated in a model of pomeron-type multiple scattering326. The corona of
hadronic and nuclear collisions is modeled in a microscopic-dynamic picture
that is assumed to be free of collective dynamics. The dynamics in the core
(which becomes more important in large systems) is described by viscous fluid
dynamics in recent versions of EPOS327. In this sense, EPOS is a particular
variant of a viscous “fluid” model. Hadronic rescattering is simulated with
UrQMD, but the fluid dynamic phase seems essential for the dynamical
translation (1) from spatial to momentum anisotropies. Different versions of
EPOS reproduce the ridge in pp collisions451, anisotropic flow phenomena in
$\mathrm{p}$–$\mathrm{Pb}$ collisions452, and radial flow phenomena in
$\mathrm{pp}$ and $\mathrm{p}$–$\mathrm{Pb}$ collisions including their mass
hierarchy453. But EPOS also predicts a characteristic nuclear modification of
heavy-flavor hadrons454 and a characteristic enhancement of thermal photon
production455 in $\mathrm{pp}$ collisions which have not been observed, yet.

####  4.3.3 Small system collectivity in transport models

In ultra-relativistic heavy ion collisions, the modeling of the bulk of the
collective evolution with partonic transport codes has a long history 456,
457, 458, 459, 460, 461, 462, 463, 321, 324. For all these simulation
packages, the code is the model in the sense that the code is more than a tool
for solving an easily stated set of transport equations of motion. The
interpretational question of why a transport code (that implement much more
physics than transport physics) does or does not agree with data is as old as
the comparison of transport codes to data and it can be difficult to address.
Amongst the above-mentioned partonic transport codes, mainly A Multi Phase
Transport (AMPT) model 321, 322, 323 and BAMPS 324, 325 have been compared
recently to data from the LHC.

AMPT successfully reproduces elliptic and triangular flow in
$\mathrm{p}$–$\mathrm{Pb}$, peripheral $\mathrm{Pb}$–$\mathrm{Pb}$ 464, 465,
466 and $\mathrm{{}^{3}He}$–$\mathrm{Au}$ collisions467. It has also been
shown to satisfactorily reproduce mass ordering of $v_{2}$ and $v_{3}$ in
settings in which coalescence was enabled 468. The question for why AMPT is
successful has focused in particular on the problem that the AMPT code
performs very differently from Molnar’s Parton Cascade MPC462 although it uses
nominally similarly small partonic cross sections. This links the discussion
of AMPT to the time-honored problem of understanding the so-called “opacity
puzzle” 462 at the beginning of the RHIC heavy-ion program. As Molnar had
realized at the time, to account for the size of the elliptic flow at RHIC,
MPC had to be initialized either with unphysically high parton densities per
unit rapidity, or it had to be evolved with unphysically large partonic cross
sections 469. There is a detailed analysis in Ref. 401 of which model-
dependent choices in AMPT may underlie its phenomenological success where
other transport codes fail.

We mention in passing the role of hadronic transport in heavy-ion
phenomenology. In simulations of nucleus–nucleus collisions that involve fluid
dynamic modeling, the system is hadronized along Cooper-Frye 58 freeze-out
hypersurfaces. The resulting hadron gas is sufficiently dilute to be evolved
with Boltzmann transport codes. Interactions in this hadronic phase affect
kinetic and hadrochemical distributions. RQMD 398, URQMD 399 and most recently
SMASH 470 are the main hadronic transport codes used to this end.

###  4.4 Model comparisons to hadrochemical measures of collectivity