The RHIC measurements featured a similar suppression level to that observed at SPS energies and a significantly larger suppression in the forward ($1.2 < |y| < 2.2$) than in the central rapidity region ($|y| < 0.35$). Detailed studies were performed to interpret the RHIC results, with a possible explanation assuming that the direct suppression of the tightly bound $J/\psi$ is partly compensated by a significant (re)generation effect [691,706–709]. However, the rather complex energy dependent interplay between hot and cold matter effects [710–712] does not allow for a significant conclusion on the magnitude of either the $J/\psi$ suppression or (re)generation.

In this situation, the studies at LHC energies are of paramount importance to settle these ambiguities in the interpretation of the results. The very large charm-quark multiplicity ($\sim 100$ c$\bar{c}$pairs in a central Pb–Pb collision at $\sqrt{s_{NN}} = 5.02$ TeV, a > 10 factor larger with respect to central Au–Au collisions at top RHIC energy) provides a test for the (re)generation mechanism which, if unambiguously demonstrated, would implicitly provide a strong evidence for the existence of a deconfined phase of nuclear matter, implying that coloured partons can roam freely over distances much larger than the hadronic scale.

Another fundamental area where LHC is expected to provide a decisive step forward is the study of the bottomonium sector, since the increase of the production cross sections with collision energy make such a measurement feasible with good enough statistical precision. The binding energies of the bottomonium vector states range from ~1.1 down to about 0.2 GeV and consequently $Υ(1S, 2S, 3S)$ represent an excellent test for the sequential suppression mechanism [713]. The (re)generation component, contrary to charmonia, is expected to be small, due to the much lower multiplicity of b$\bar{b}$ pairs ($5-10$ in central Pb–Pb collisions) compared to c$\bar{c}$. In addition, bottomonia are very promising from a theoretical point of view, thanks to a much more pronounced separation of scales between the heavy-quark mass, the typical bottomonium size and the radial or orbital angular-momentum excitation energy ($M \gg Mv \gg Mv^2$, where M is the heavy-quark mass and v the relative velocity of the pair in the bound state). That makes them good candidates for a direct application of effective field theory approaches [681].

Finally, at lower energies, competing non-QGP effects like quarkonium break-up by the nucleons of the colliding nuclei played a significant role in the description of the results [714]. At the LHC, due to the extremely short crossing time of the nuclei, the only sizeable effect not related to the medium is represented by nuclear modification of gluon densities (shading [46,47,715] or by CGC-related effects [716–718]), or coherent energy loss [719]. These effects can be investigated in the study of p–Pb collisions.

In the following sections we review the main results obtained by ALICE for various charmonium and bottomonium states, comparing them to experimental results at lower energies and to various theory calculations.

### 2.5.1 Study of the charmonium ground state: evidence for the (re)generation and demonstration of deconfinement at LHC energies

As outlined in the previous section, results at LHC energies were expected to clarify, in the charmonium sector, the presence of a (re)generation mechanism related to (re)combination of deconfined charm quarks. In ALICE, charmonia are measured at both central ($|y| < 0.9$) and forward rapidity ($2.5 < y < 4$), and down to zero $p_T$. This phase space coverage is ideal for such studies because the amount of (re)generation is expected to depend on rapidity and $p_T$, with a stronger effect around $y = 0$ and at small $p_T$, due to the higher charm-quark multiplicity in these kinematic regions. It should be remarked that the results shown in these sections refer to inclusive production, i.e., the contribution of decays of hadrons containing a b quark was not subtracted. Such a contribution has a small effect on the nuclear modification factor at moderate $p_T$ ($< 10\%$ for $p_T < 5$ GeV/c) [720] and increases at higher $p_T$ [721].

One of the main results on $J/\psi$ production is shown in the top panel of Fig. 55, where the nuclear modification factor for inclusive $J/\psi$ measured by ALICE in Pb–Pb collisions, at $\sqrt{s_{NN}} = 5.02$ TeV at central rapidity in the $e^+e^-$ decay channel, is shown as a function of the charged hadron pseudorapidity density $\langle dN_{ch}/d\eta \rangle$ at midrapidity [438]. The latter quantity is directly related, for a certain collision system at a given energy, to the centrality of the interaction and, for different collision systems, is roughly proportional to the initial energy density [201]. The ALICE result is compared with the $R_{AA}$ values measured at RHIC by STAR (Au–Au at $\sqrt{s_{NN}} = 0.2$ TeV) in $|y| < 0.5$ [705] and those obtained from NA50 results at the SPS in Pb–Pb at $\sqrt{s_{NN}} = 0.017$ TeV for $0 < y < 1$ [694], using a pp reference cross section extrapolated from NA60 pA measurements ranging from p–Be to p–U [695]. The selection $p_T > 0.15$ GeV/c is meant to remove a non-negligible contribution due to $J/\psi$ photoproduction in peripheral hadronic Pb–Pb collisions [722]. The prominent feature of this set of results is the strong decrease of the $J/\psi$ suppression moving from low- to high-energy experiments and, at the LHC, the disappearance of suppression effects when going towards central collisions. Both effects provide a strong indication of the presence of (re)generation effects on the $J/\psi$. Similar observations were first carried out by ALICE at forward rapidity [721]. As discussed later in Sect. 2.5.3, the contribution of non-QGP effects, in particular of nuclear shadowing, does not alter this conclusion.

A more quantitative assessment of the physics mechanisms at play requires a comparison of ALICE results with models. In the bottom panel of Fig. 55, calcula-
<figure><img src="image_67.png" /><figcaption>**Fig. 55** The nuclear modification factor as a function of the charged hadron multiplicity. For ALICE and STAR results, the selection $p_T > 0.15$ GeV/c minimises the contribution of photoproduced $J/\psi$ [722]. (Top) Comparison between SPS (NA50) [694], RHIC (STAR) [705] and LHC (ALICE) [438] results; (bottom) the ALICE results compared with theoretical calculations (SHM [723], Comovers [724], TM-TAMU [725], TM-Tsinghua [690])</figcaption></figure>

tions of $R_{AA}$ carried out in the frame of transport models (TM), of a comover model and of a statistical hadronisation model (SHM) are presented. The transport models (TM-TAMU [725], TM-Tsinghua [690]) rely on the solutions of macroscopic transport equations. Dissociation/(re)generation thermal rates for quarkonium states in the QGP are calculated taking into account a lattice-QCD inspired evaluation of the dependence of their spectral properties on the evolving thermodynamical properties of the medium. The two models discussed here mainly differ in the choice of the rate equation and of the open charm cross section. In the comover model [724], the scattering of the nascent charmonium states with “comovers” (partonic or hadronic) produced in the same kinematic region is at the origin of the suppression, with the comover density being tuned on the measured hadron yields. The dissociation cross

section, extracted from results at lower energy, is assumed to be energy independent and the (re)generation effects are implemented through a gain and loss differential transport equation. Finally, in the SHM [723], charmonium yields are assumed to be determined at chemical freeze-out according to their statistical weights, introducing a charm fugacity factor related to charm conservation and obtained from the charm production cross section. The comparison between data and models shows a fair agreement with most of the approaches, with TM-TAMU and SHM giving the best reproduction of the $\langle dN_{ch}/d\eta \rangle$ dependence. The uncertainties of the models are clearly large when compared to the experimental data, mainly as a direct consequence of the uncertainty on the total charm cross section for Pb–Pb collisions, which is a fundamental input of the calculations. Getting a precise evaluation of this quantity is not easy, due to the problem of measuring open charm mesons and baryons down to zero $p_T$ in nuclear collisions. Recent measurements by ALICE [726] have shown strong modifications of the fragmentation fractions of the c quark to the various final states from $e^+e^-$ to pp collisions and allowed a fairly precise (~10% uncertainty) estimate of the total charm cross section for the latter collision system (see also Sects. 2.3.1 and 6.2.2 for more details on the measurements). However, the extrapolation from pp collisions is not trivial as it involves an estimate of the non-negligible shadowing effects on the initial state [727].

Other valuable insights on the mechanisms at play can be obtained by considering the $p_T$ dependence of the $R_{AA}$ and $v_2$ of the $J/\psi$. In the top panel of Fig. 56 results on the nuclear modification factor are presented for the two rapidity ranges accessed by ALICE, $|y| < 0.9$ ($e^+e^-$ decay) and $2.5 < y < 4$ ($\mu^+\mu^-$ decay), for central Pb–Pb events [438]. The significant rise of $R_{AA}$ at low-$p_T$, where the bulk of charm production occurs, with this effect being stronger at midrapidity, strongly hints at the presence of (re)generation as a dominant mechanism in this transverse momentum range. Remarkably, at higher $p_T$ ($> 8–10$ GeV/c) there is no rapidity dependence for the strong suppression, and the $R_{AA}$ values become compatible with those measured for charged hadrons (see Fig. 45 in Sect. 2.4.2). This observation may suggest parton energy loss, as in the case of $J/\psi$ production from the splitting of a hard gluon, as a significant effect in this kinematic range [730,731], even if a quantitative description inspired to NRQCD underestimate the measured $R_{AA}$ [732]. In the bottom panel of Fig. 56 the $p_T$ dependence of $v_2$ for Pb–Pb collisions in the centrality interval 20–40% is shown [728]. The measurement is performed using the scalar product method with a pseudorapidity gap $|\Delta\eta| > 1.1$ between the $J/\psi$ and the event flow vector, to suppress non-flow effects. A large $v_2$ reaching ~0.1 at $p_T \sim 5$ GeV/c is observed, which could be explained if a large fraction of the detected $J/\psi$ originates from the (re)combination of charm quarks, which acquire their anisotropy by taking part in the
<figure><img src="image_68.png" /></figure>

<figure><img src="image_69.png" /><figcaption>**Fig. 56 (Top)** The nuclear modification factor as a function of the transverse momentum, for the central ($|y|<0.9$) and forward ($2.5 < y < 4$) rapidity regions [438], compared to TM-TAMU and SHM model calculations. (Bottom) The $p_T$ dependence of the inclusive $J/\psi v_2$, compared to TM-TAMU model calculations [728,729]</figcaption></figure>

collective expansion of the system. At the same time, the large $J/\psi v_2$ values, according to Ref. [725], favour a late $J/\psi$ formation time, since the charm-quark anisotropy requires time to build up.

The results of Fig. 56 are compared with theoretical calculations from TM-TAMU and SHM (only $R_{AA}$ for the latter). For $R_{AA}$, the models reproduce the different size of the rise at low-$p_T$. At large $p_T$ only TM-TAMU is in agreement with data, due to the absence in the SHM of additional production mechanisms such as $J/\psi$ production from gluon fragmentation in jet, which are expected to contribute in that range. In the high-$p_T$ range the slight rise of $R_{AA}$ in TM-TAMU is mainly due to the longer resonance formation time, that is Lorentz-dilated, coupled to the expected lower suppression rate for the pre-resonant state [706,708]. For $v_2$, TM-TAMU is able to reproduce the data over the explored $p_T$ range. In the low- and intermediate-$p_T$ region, where the $J/\psi$ is pro-

duced mostly via (re)combination, space-momentum correlations of the diffusing charm and anti-charm quarks in the hydrodynamically expanding fireball represent an important ingredient of the model [729]. At high-$p_T$, where the $v_2$ is likely determined by path-length dependent effects [612], the calculations describe the data within the large statistical uncertainties. It is worth noting that no energy loss effects are implemented in the TM-TAMU calculations, so the positive $v_2$ in the high-$p_T$ region can only be originated in this model from path-length dependent dissociation.

### 2.5.2 Study of the bottomonium ground state: strong suppression and small (re)generation effects

A remarkable observation at the LHC was that the $\Upsilon$ states are strongly suppressed, in a hierarchy apparently dictated by their binding energies [682]. First observed by the CMS Collaboration at midrapidity [733], this hierarchy of suppression was confirmed by the ALICE Collaboration at forward rapidity [734,735]. In Fig. 57 the $\langle N_{part} \rangle$ dependence of the $\Upsilon(1S)$ $R_{AA}$ measurements performed by ALICE at $\sqrt{s_{NN}} = 5.02$ TeV is shown and indicates a quick onset of the suppression of the $\Upsilon(1S)$ already in semi-peripheral collisions. The results are compared with various models. For each of them two lines are shown, corresponding to upper and lower uncertainty limits. All the models describe the data quite well, although the theory uncertainties are rather large. First, the same comover model [736] introduced in Sect. 2.5.1 is applied to the bottomonia, but, in contrast to the charmonium case, no (re)generation component is introduced. The model uncertainties are related to the used nPDF parameterisation and to the comover-$\Upsilon$ dissociation cross section. Second, the transport model introduced above for the $J/\psi$ (TM-TAMU [737]) is also used for the $\Upsilon(1S)$. The model results are shown with and without a (small) (re)generation component. The former is favoured by the measured $R_{AA}$ in central collisions. Third, a calculation of another transport model based on the framework of coupled Boltzmann equations [738] is also shown. In this approach the (re)generation is dominated by real-time (re)combination of correlated heavy-quark pairs, and the uncertainties correspond to those on the EPPS16 nPDF parametrisation [47]. Fourth, the model labelled “aHydro” in the figures implements a thermal modification of a complex heavy-quark potential inside an anisotropic plasma [739]. The survival probability of bottomonia is evaluated based on the local energy density, integrating a rate equation over the proper time of each state. The uncertainties represent the envelope of calculations performed with values of the viscosity-to-entropy density ratio $\eta/s$ ranging from $1/(4\pi)$ to $3/(4\pi)$. No modifications of nuclear PDFs or any (re)generation phenomenon are included. Finally, the “Open Quantum Systems” model results are obtained using poten-
<figure><img src="image_70.png" /><figcaption>**Fig. 57** The nuclear modification factor for inclusive Υ(1S) as a function of $N_{\text{part}}$ [735], compared to model calculations [736–739,741]</figcaption></figure>

tial Non-Relativistic QCD (pNRQCD) and the formalism of open quantum systems to treat the interaction of the Υ(1S) state with the QGP [740,741]. The uncertainties shown here correspond to variations of the non-perturbative transport coefficient $\kappa$ (heavy-quark momentum diffusion coefficient), keeping the value of its dispersive counterpart $\gamma$ fixed at its central value.

In pp collisions, up to 30–50% of the measured Υ(1S) yield results from the feed-down from other states [68,696, 742]. Consequently, a significant amount of Υ(1S) suppression may arise from the very strong suppression of the excited states (see Sects. 2.5.3 and 2.5.4). The models mentioned above follow similar approaches for the evaluation of the feed-down contribution, based on PDG values [37] for the branching ratios of the excited states, and information on the production cross sections of those states from pp data taken at LHC energy.

The $p_T$ and rapidity dependence of the $R_{AA}$ for Υ(1S) are shown in Fig. 58 [735]. The data exhibit a very weak (if any) $p_T$ dependence and a weak rapidity dependence, apparent only for the most forward rapidities $3 < y < 4$ covered by ALICE. The $p_T$ dependence is reproduced by both TM-TAMU and aHydro models [737,744]. It is worth noting that in the TM-TAMU model [737], the only one to implement (re)generation, its (small) effect becomes visible at $p_T$ of the order of the Υ(1S) mass, where the uncertainties are still rather large. More interesting is the rapidity dependence of the Υ(1S) $R_{AA}$. Complemented by the CMS results [743], the measurement spans four units of rapidity. While both the hydrodynamic [744] and the transport model [738] describe the Υ(1S) $R_{AA}$ at midrapidity rather well, they predict an opposite trend at forward rapidity than suggested by the data. In the models, the rapidity dependence of the suppression is determined by the rapidity dependence

of the energy density in the QGP phase, which is expected to decrease slightly towards forward $y$. The opposite trend suggested by the data may find a natural explanation in case a significant (re)generation component is present, even if the size of this effect, as visible in the left panel of Fig. 58, is predicted to be relatively small in TM-TAMU calculations. Another effect that would, in principle, lead to the observed behaviour for $R_{AA}$ is the presence of significant energy loss effects that may lead to a shrinking of the rapidity distribution in Pb–Pb collisions. Models that implement coherent energy loss effects in cold nuclear matter [745] suggest a decrease of $R_{AA}$ at forward rapidity, but the effect is much smaller than observed.

The bottomonium data were recently used for a first explicit extraction of the heavy-quark potential in the hot QCD medium [746], found to be characterised by significant remnants of the long-range confining force in the QGP. The bottomonium family is of utmost relevance for the next advances in the theoretical treatment, based on a quantum evolution in the deconfined medium, where first quantitative results are just becoming available [740,741].

### 2.5.3 Study of non-QGP effects on ground state suppression

The $J/\psi$ and Υ(1S) production yields are, to various extents, modified in Pb–Pb collisions. As detailed in Sects. 2.5.1 and 2.5.2, the $J/\psi$ nuclear modification factor presents a suppression of the charmonium ground state significantly increasing with $p_T$ (see Fig. 56 top), while for the Υ(1S) the $p_T$ dependence is weaker (see Fig. 58 left).

Even if QGP effects are significant at LHC energies, other non-QGP mechanisms have influence on the production of quarkonium. Commonly defined as “cold nuclear matter” (CNM) effects, in contrast to the hot-matter effects related to the presence of the QGP, they can be relevant in collisions involving light and heavy nuclei. Therefore, the quantification of their size is important to correctly define the net impact of QGP effects on the production of quarkonia. Experimentally, CNM effects are studied in minimum-bias pA collisions, where hot-matter mechanisms are not expected to dominate.

ALICE has studied quarkonium production in p–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV and 8.16 TeV [540, 727, 748, 749, 752–755]. The $p_T$ dependence of the nuclear modification factor of the $J/\psi$ and Υ(1S), at the highest energy, is shown in Fig. 59 at forward (p-going) and backward (Pb-going) rapidity. For both $J/\psi$ and Υ(1S) a significant $p_T$ dependence, with a stronger suppression in the low-$p_T$ region, can be seen. Furthermore, for the $J/\psi$ the effect is more important in the p-going ($2.03 < y_{\text{cms}} < 3.53$) than in the Pb-going direction ($-4.46 < y_{\text{cms}} < 3.53$), while for the Υ(1S) no significantly different behaviour, within the statistical and systematic uncertainties, is observed in the two accessible rapidity
<figure><img src="image_71.png" /><figcaption>**Fig. 58** The transverse momentum (left, data from ALICE [735]) and rapidity (right, data from ALICE [735] and CMS [743]) dependence of the inclusive Υ(1S) nuclear modification factor, compared to model calculations [737–739]</figcaption></figure>

domains. The comparison of the $J/\psi$ and Υ(1S) $R_{pA}$ with theoretical models that implement only nuclear modification of the gluon PDFs [750] confirms the negligible role of QGP-induced mechanisms in pA collisions, and it underlines the importance of CNM effects related to gluon saturation and/or nuclear shadowing.

Significantly different nuclear modification factors are measured in p–Pb and Pb–Pb collisions, with stronger suppressions observed in Pb–Pb, except in the low-$p_T$ region. At high-$p_T$, the $J/\psi$ suppression measured in Pb–Pb collisions can certainly not be due to the CNM effects alone and, even in the limits of this qualitative comparison, this is a clear confirmation of the role played by mechanisms related to the hot matter formation. At low-$p_T$, $R_{AA}$ and $R_{pPb}$ become quite similar. This does not imply that $R_{AA}$ is dominated by CNM effects, but rather that (re)combination, which was shown to be strong in this kinematic region (see left panel of Fig. 56) completely balances suppression effects, leading to an approximate recovery of the binary collision scaling once CNM effects are taken into account. We also remark that the comparison of p–Pb and Pb–Pb results at two different energies is meaningful when considering the shadowing effects, since the variation in the Bjorken-$x$ ranges due to the different $\sqrt{s_{NN}}$ is compensated by the different $y_{cms}$ coverage (about half a unit rapidity shift).

For the Υ(1S), the difference in the $p_T$ dependence of $R_{AA}$ and $R_{pA}$, in particular at high-$p_T$, gives an indication that hot-matter final state effects are significant. At low-$p_T$, the $R_{AA}$ and $R_{pPb}$ values become more similar. Since (re)combination effects are negligible for bottomonium, this observation implies that a significant fraction of the observed Υ(1S) suppression in Pb–Pb collisions might be ascribed to CNM effects.

A more quantitative estimate of the magnitude of CNM effects in nucleus–nucleus collisions can be performed assuming that such effects are dominated by nuclear shad-

owing. Under this hypothesis, the magnitude of shadowing effects in Pb–Pb collisions can be obtained as the product $R_{pPb}^{\text{forw}} \times R_{pPb}^{\text{backw}}$ [727,752]. By normalising the measured $R_{AA}$ to this product, the quantity $S = R_{AA}/(R_{pPb}^{\text{forw}} \times R_{pPb}^{\text{backw}})$ can be computed. This quantity, shown in Fig. 60 for both $J/\psi$ and Υ(1S), can be considered as the Pb–Pb nuclear modification factor, corrected for shadowing effects. For the $J/\psi$, values larger than unity are obtained at low-$p_T$. This indicates a net $J/\psi$ enhancement, clearly related to (re)generation effects [727]. At high-$p_T$, values as low as $S_{J/\psi} \sim 0.3$ are measured, corresponding to a strong suppression.

Also $S_{\Upsilon(1S)}$ shows a significant $p_T$ dependence. At high-$p_T$ its value is similar to $S_{J/\psi}$, while at low-$p_T$ it becomes compatible with unity. Taking into account that for Pb–Pb collisions the contribution from the (re)combination of b quarks should be negligible, this result could imply that once shadowing effects are divided out, only a modest suppression effect survives for the strongly bound Υ(1S). This tentative conclusion can be further quantified by taking into account the feed-down contribution from the excited bottomonium states. A thorough investigation of these effects can be found in Sec. 2.1 of Ref. [696]. Indeed, both Υ(2S) and Υ(3S), as well as various $\chi_b$ states have a non-negligible branching ratio towards the $1^{--}$ ground state. The feed-down fractions were measured in pp collisions by ATLAS, CMS and LHCb, with the resulting contribution from all these states being $\sim 30\%$ at relatively low-$p_T$, increasing to $\sim 50\%$ at high-$p_T$ [696]. A precise evaluation of the direct Υ(1S) $R_{AA}$ would require the knowledge of the corresponding quantity for the excited states, which is available for the Υ(2S) and Υ(3S) [735,743], but not for the various $\chi_b$. Assuming $R_{AA}(\chi_b(1P)) \sim R_{AA}(\Upsilon(2S))$ and $R_{AA}(\chi_b(2P)) \sim R_{AA}(\Upsilon(3S))$, as it may be suggested by the similar binding energies, one may evaluate the direct Υ(1S) nuclear modification factor by subtracting out the contribution of
<figure><img src="image_72.png" /><figcaption>**Fig. 59** J/ψ (left) and Υ(1S) (right) $R_{AA}$ as a function of $p_T$ in centrality-integrated Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV in $2.5 < y_{cms} < 4$ (blue circles) [735,747]. Results are compared to the corresponding $R_{pA}$ values obtained in p–Pb collisions at $\sqrt{s_{NN}} = 8.16$ TeV in $2.03 < y_{cms} < 3.53$ (light gray circles) and $-4.46 < y_{cms} < -2.96$ (dark gray circles) [748,749]. The gray bands correspond to model calculations based on a pure nuclear shadowing scenario using EPS09LO as set of nuclear parton distribution functions (filled bands) [750]. For the J/ψ case, a theory calculation based on a CGC approach coupled with a NRQCD production mechanism is also shown (dashed band) [751]</figcaption></figure>

the excited states as $R_{AA}(\Upsilon(1S))^{\text{dir}} = (R_{AA}(\Upsilon(1S))^{\text{incl}} - R_{AA}(\Upsilon(2S)) F_{\Upsilon(2S)\to\Upsilon(1S)} - \dots)/F_{\Upsilon(1S)^{\text{dir}}}$. A numerical evaluation for the lowest $p_T$ interval ($p_T < 2$ GeV/c) in Fig. 59 (bottom) gives $R_{AA}(\Upsilon(1S))^{\text{dir}} \sim 0.45 \pm 0.08$, to be compared with the product of the forward and backward $R_{pPb}$ which yields $0.42 \pm 0.11$. This approximate estimate suggests that shadowing and feed-down effects could be responsible for most of the suppression observed at low-$p_T$ for inclusive $\Upsilon(1S)$ production. At larger $p_T$ CNM effects become weaker but the feed-down fractions increase. The corresponding results, affected by larger uncertainties, are $R_{AA}(\Upsilon(1S))^{\text{dir}} \sim 0.68 \pm 0.22$ while $R_{pPb}^{\text{forw}} \times R_{pPb}^{\text{backw}} \sim 0.86 \pm 0.26$.

<figure><img src="image_73.png" /><figcaption>**Fig. 60** $S_{J/\psi}$ (green symbols) and $S_{\Upsilon(1S)}$ (orange symbols) as a function of $p_T$. See text for further details on the definition of this quantity</figcaption></figure>

### 2.5.4 Excited quarkonium states

Both c̄c and b̄b states have complex spectroscopies. The various states have very different binding energies, ranging from about 50 MeV for the $\psi(2S)$ to about 650 MeV for the J/ψ in the charmonium sector and from about 200 MeV for the loosely bound Υ(3S) up to about 1 GeV for the tightly-bound Υ(1S) in the bottomonium sector. Hence the comparison of their behaviours in a hot and dense medium can shed light on various properties of the created system. In particular, by linking the sequential disappearance of these resonances to the melting temperature, as evaluated by lattice QCD studies [686], one might also have a tool to access the temperature reached in the collision. In addition, also the details of the production process itself in nuclear collisions, with the contribution of direct and regenerated states, can be different for the various quarkonium resonances. Thus, it is expected that a precise measurement of both ground and excited states will allow distinguishing between statistical hadronisation [756], which assumes that all quarkonium states are produced according to thermal weights determined at the chemical freezeout, and microscopic transport models [725]. In the latter, the production of quarkonia is continuous through the lifetime of the fireball and can thus occur out of chemical equilibrium.

ALICE has measured both the J/ψ and the $\psi(2S)$ at forward rapidity ($2.5 < y < 4$) and down to zero transverse momentum [757]. These results extend the measurements of the CMS experiment, which has studied the J/ψ and the $\psi(2S)$ in the kinematic ranges $|y| < 1.6, 6.5 < p_T < 30$ GeV/c and $1.6 < |y| < 2.4, 3 < p_T < 30$ GeV/c [758]. As visible in Fig. 61, where the nuclear modification factors of the J/ψ and $\psi(2S)$ measured in the rapidity range $2.5 < y < 4$ are shown, the $\psi(2S)$ is significantly more sup-
<figure><img src="image_74.png" /><figcaption>**Fig. 61** (Top) The nuclear modification factor for inclusive $\psi(2S)$ and $J/\psi$, measured in Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV as a function of $N_{part}$ [757], in the range $2.5 < y < 4$. Predictions from the TM-TAMU [725] and SHM [723] models are also shown (Bottom). The $p_T$ dependence of the $\psi(2S) R_{AA}$, measured by ALICE in $2.5 < y < 4$ for 0–90% Pb–Pb collisions [757], compared with corresponding high-$p_T$ results in $|y| < 1.6$ and 0–100% centrality from CMS [759]. Model comparisons are also shown</figcaption></figure>

pressed than the $J/\psi$ over all the centrality and $p_T$ ranges explored by ALICE. Furthermore, a rise of $R_{AA}$ at low-$p_T$ can be observed also for the $\psi(2S)$, suggesting that (re)combination effects can be sizeable also for this state. The calculations of the TM-TAMU model [725] are in good agreement with the experimental results as a function of both $N_{part}$ and $p_T$. The SHM approach [723] tends to overestimate the suppression effects for the $\psi(2S)$ in central events, while a good agreement, as discussed previously, is found for the $J/\psi$. In the region $6.5 < p_T < 12$ GeV/c the results agree with those of CMS [759], obtained at midrapidity, which also explore a higher-$p_T$ range, where the $\psi(2S) R_{AA}$ remains roughly constant.

A similar comparison between the $R_{AA}$ centrality dependence for ground and excited states has been carried out by ALICE also in the bottomonium sector in Pb–Pb collisions

<figure><img src="image_75.png" /><figcaption>**Fig. 62** The nuclear modification factor for inclusive $\Upsilon(2S)$ and $\Upsilon(1S)$ as a function of $N_{part}$ [735], compared to model calculations [736, 738, 740, 744]</figcaption></figure>

at $\sqrt{s_{NN}} = 5.02$ TeV [735]. Bottomonium states are reconstructed at forward rapidity ($2.5 < y < 4$). As for charmonium states, also the $\Upsilon(1S)$ and $\Upsilon(2S)$ show a strong suppression in their production, with a more significant reduction of their yields in the most central collisions. Even if the measurement of the $\Upsilon(2S)$ yields in Pb–Pb is still statistically limited, allowing only a small number of centrality intervals to be defined, a stronger suppression of the excited state, with respect to the $\Upsilon(1S)$, is visible in Fig. 62.

It can be noted that the size of the measured suppression of the $\Upsilon(1S)$ and $\Upsilon(2S)$ production in Pb–Pb with respect to pp collisions is well described by several theory models, whose predictions are shown in Fig. 62, and it is similar to the one measured by the CMS experiment [760], in a complementary rapidity region $|y| < 2.4$ [743]. Some of the models shown in Fig. 62 give explicit predictions for the temperatures needed to describe the relative suppression of the two states. In particular, in the aHydro model [739], the hydrodynamic fireball evolution starts at $\tau_0 = 0.3$ fm/c, with an initial system temperature and shear viscosity to entropy ratio tuned simultaneously to describe the measured charged particle density and anisotropic flow. The most probable values for the above quantities in this model are $T = 641$ MeV and $\eta/s = 1/4\pi$, in agreement with other hydrodynamic calculations which favour the hypothesis of a nearly perfect QGP fluid with a shear viscosity close to the lowest bound conjectured by the AdS/CFT correspondence [107]. Other approaches, as TM-TAMU [737], start the space–time evolution of the medium at a slightly lower temperature, $T \sim 540$ MeV. In any case, more quantitative constraints on the initial QGP temperature and its shear viscosity require a better understanding of all the differential experimental results, especially the rapidity dependence which, as previously shown in Sect. 2.5.2 for the $\Upsilon(1S)$, does not seem to be well reproduced in this model.
### 2.5.5 Discussion of further flow results for quarkonia

A detailed study of the J/ψ anisotropic flow, reported in Refs. [728,761] offers additional insights into the dynamics of the production process. The centrality dependence of the elliptic flow coefficients allows us the study of the sensitivity of the J/ψ production mechanisms to the system size and initial spatial anisotropy. The $v_2$ of light flavoured hadrons was shown to be proportional to the initial state spatial ellipticity [303], as expected from the evolution of a nearly ideal hydrodynamic system, with the proportionality factor being dependent on the medium properties, like the equation of state and viscosity. Since the J/ψ $v_2$ should be determined by the flow of charm quarks in the low- and intermediate-$p_T$ regions where (re)combination processes dominate its production, its centrality dependence provides information on the degree of charm-quark thermalisation/equilibration in the medium as well as on the charmonium formation time.

Figure 63 presents the centrality dependence of the $p_T$-integrated J/ψ $v_2$ in Pb-Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV, measured at forward rapidity in the ranges $p_T < 5$ GeV/c (left) and $5 < p_T < 20$ GeV/c (right) [728]. These results are compared to the $v_2$ of the charged pions measured at midrapidity [290] and at a $p_T$ corresponding to the J/ψ average $p_T$ in the above mentioned intervals. The $v_2$ for low-$p_T$ pions grows fast from the most central towards peripheral collisions, reaching a maximum in the 40–50% centrality interval. For the J/ψ, however, a milder centrality dependence is observed, illustrated quantitatively by the $v_2^\pi/v_2^{J/\psi}$ ratio shown in the bottom panel, which grows linearly from central to peripheral collisions. The significance of the trend is larger than $2\sigma$. Such a trend for the J/ψ may be interpreted as the consequence of an increase towards more central collisions of the fraction of regenerated J/ψ and/or the degree of charm-quark thermalisation, with both favouring the development of anisotropic flow in more central collisions. An incomplete charm thermalisation, arguably affecting more the peripheral collisions, would lead, as the measurements suggest, to a centrality dependent conversion factor of the initial state spatial anisotropy into final state momentum anisotropy. At high-$p_T$, the $v_2$ shows a similar centrality dependence for J/ψ and charged pions, suggesting a possible common origin.

An extension to the picture of the number of constituent quark (NCQ) scaling of the anisotropic flow, first proposed at RHIC [296,762], which assumes that the charm-quark flow magnitude is different than that of the lighter quark species, can be tested as suggested in Refs. [728,763]. Within the NCQ ansatz, the $v_2$ of the charm ($v_2^c$) and of light quarks ($v_2^q$) can be inferred as described in Ref. [728], based on the measured $v_2$ of J/ψ and charged pions, respectively. The $v_2$ of D mesons can be obtained as a superposition of the charm

and light-quark flow

$$v_2^D(p_T^D) = v_2^q(p_T^q) + v_2^c(p_T^c), \qquad (22)$$

where q and c superscripts denote the light and charm quarks, respectively. Unlike for light mesons and charmonia, where the momentum fractions carried by the constituent quarks are assumed equal from symmetry reasons, for D mesons different assumptions need to be made. In coalescence motivated models of hadronisation, constituent quarks have equal velocities and thus carry $p_T$ fractions proportional to their effective masses. For D mesons, a simplistic approach considering the mass of the light quark being one third of the proton mass leads to a light-quark $p_T$ fraction of about 0.2 of the D-meson $p_T$ [763]. Figure 64 shows the $p_T$ dependence of $v_2$ in the 30–50% centrality range measured by ALICE at midrapidity for charged pions [290] and D-mesons [573] (average of prompt D$^0$, D$^+$ and D$^{*+}$ mesons), and at both central and forward rapidity for J/ψ mesons [728]. The red dashed curve is a fit of the J/ψ $v_2$ [728] used to obtain a parameterisation of the charm-quark $v_2$, while the black solid curve corresponds to the D-meson $v_2$ obtained using Eq. 22 and assuming $p_T^q = 0.2p_T^D$. This calculation underestimates the measured D-meson flow at low-$p_T$ and overestimates it at high-$p_T$, which disfavours the constituent quark coalescence inspired assumption. This approach however largely ignores the contribution of gluons to the D meson kinematics. In this naive study, an attempt to compensate for this missing contribution was made by roughly tuning the momentum share carried by the light quark. Surprisingly, the assumption of equal $p_T$ sharing between the light and charm quarks provides a much better description of the data, with the best fit suggesting a value of $p_T^q = 0.4p_T^D$. A very good compatibility with the measurements in different centrality ranges was observed for both $v_2$ and $v_3$ [728], which suggests that the flow of charmonia and open charm mesons can be described assuming a common underlying charm-quark flow different from that of the light quarks.

Finally, Fig. 65 shows a comparison of the $p_T$ dependent $v_2$ measurements for charged pions, prompt D mesons, inclusive J/ψ, electrons from beauty-hadron decays and $\Upsilon(1S)$ [533,573,764]. A rather clear quark flavour hierarchy is observed in the low-$p_T$ range, for both open and hidden heavy flavour hadron species, with the beauty hadrons exhibiting the least amount of flow. Both open and hidden charm hadrons show a significant amount of anisotropic flow, suggesting that charm quarks are at least partly thermalised in the QGP medium. In addition, since the low-$p_T$ J/ψ are not expected to develop a significant collective anisotropic flow on their own, their relatively large observed flow supports the scenario of J/ψ formation via (re)combination during the late stages of the collision, when the charm-quark flow is fully developed. The elliptic flow of the electrons from beauty hadron decays, measured up to the electron $p_T = 6$ GeV/c,
<figure><img src="image_76.png" /><figcaption>**Fig. 63** The centrality dependence of inclusive J/ψ $v_2$ for two transverse momentum regions [728], compared with the corresponding results for charged pions [290]</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>**Fig. 64** The $p_T$ dependence of inclusive J/ψ $v_2$ [728] compared with the corresponding results for prompt D mesons [573] and charged pions [290]. The expected $v_2$ for open charm in the constituent quark approach, and for various values of the ratio $p_T^q/p_T^D$ are also shown</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>**Fig. 65** The $p_T$ dependence of $v_2$ for the centrality range 30–50% in Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. Results are shown for open and hidden charm (D-mesons and J/ψ), for open and hidden beauty (electrons from semi-leptonic beauty-hadron decays and $\Upsilon(1S)$), and for pions [533,573,728,764]</figcaption></figure>

is a convolution of the $v_2$ of beauty hadrons mostly from the $p_T$ range $p_T < 15$ GeV/c and their decay kinematics. Although the $p_T$ dependence of the open beauty hadrons $v_2$ cannot be easily inferred, this measurement clearly shows that open beauty hadrons exhibit flow, even if significantly smaller with respect to open charm hadrons. This observation constitutes also a hint that beauty quarks may participate to some extent to the collective flow of the plasma, although it is difficult to disentangle the contribution of the b and of light quarks in producing the observed b-hadron flow. In the case of hidden beauty hadrons, namely inclusive $\Upsilon(1S)$, the $v_2$ measurements using Run 2 data are compatible with zero elliptic flow and are estimated to be $2.6\sigma$ lower than the J/ψ $v_2$. Within the large uncertainties of these measurements, such an observation is compatible with the expectations of

a negligible contribution from (re)generation in the beauty sector.

### 2.5.6 Conclusions

**Charmonium production.** In Pb–Pb collisions at LHC energies the J/ψ production and the observed anisotropy in its azimuthal distributions is dominated, at low-$p_T$ and for central events, by a (re)generation effect due to the (re)combination of the charm quarks abundantly produced in the collision and at least partially thermalised. This observation constitutes a proof of deconfinement, as it implies that coloured partons can move freely over distances much larger
than the hadronic scale. Results on the weakly bound $\psi(2S)$ state show that this resonance exhibits a stronger suppression compared to $J/\psi$ in all the explored $p_T$ interval. Hints for the presence of (re)combination effects were also detected.

**Bottomonium production.** The bottomonium state $Υ(1S)$ is suppressed by more than 50% in central Pb–Pb collisions, with a nearly flat transverse momentum dependence of the $R_{AA}$ that indicates that (re)generation effects play a small role, as expected because of the much lower multiplicity of b with respect to c quarks. The less strongly bound $Υ(2S)$ is found to be more strongly suppressed (up to about 90% in central Pb–Pb collisions), an observation consistent with a sequential suppression of the bottomonium states.

**Initial state effects.** The influence of gluon shadowing/ saturation is not negligible, with a 30% suppression measured at low-$p_T$ for $J/\psi$ and $Υ(1S)$. Also, feed-down from higher mass charmonium and bottomonium resonances contributes about 20% and 30% at low-$p_T$, respectively. When considering these effects in a semi-quantitative way, it turns out that: (i) they could be responsible for a significant fraction of the measured inclusive $Υ(1S)$ suppression, (ii) that production of low-$p_T$ direct $J/\psi$ is larger than expected from binary scaled pp collisions.

## 2.6 Early electromagnetic fields and novel QCD phenomena

Heavy-ion collisions provide a unique opportunity to study intriguing, novel QCD phenomena that are not directly accessible elsewhere, including local parity (P) and charge conjugation parity symmetry (CP) violating effects in the strong interaction [765–773]. Even though many potential explanations have been proposed [774–780], it is still unclear why P and CP invariances are respected in the strong interaction. This is known as the strong-CP problem [778,781] which still continues to be one of the remaining puzzles of the Standard Model.

The uniqueness of heavy-ion collisions as a tool for the study of the strong P and CP violation effects stems from the presence of the strongest magnetic field in nature, of the order of $10^{19}$ gauss [782–784], produced by the colliding positively charged nucleons, mostly protons, that do not participate in the interaction (so-called spectators) in non-central collisions. This field develops perpendicular to the reaction plane, the plane defined by the impact parameter and the beam axis. The magnetic field, that is rapidly decaying [785,786] with a rate that depends on the electric conductivity of the medium, besides leading to observable effects of the strong P and CP violation, could have direct implications on the kinematics of final state particles. The latter is used for detecting and estimating the strength of the electromagnetic

fields. While so far there exists no unambiguous experimental evidence for the existence of such strong electromagnetic fields, ALICE measurements reveal several important hints, which are discussed below.

The results from ALICE searches for the effects of the early stage electromagnetic fields will be discussed first. This section will be followed by the description of our studies for the discovery of novel QCD phenomena, associated with local P and CP violating effects in the strong interaction, referred to in the literature as the Chiral Magnetic Effect (CME) and the Chiral Magnetic Wave (CMW).

### 2.6.1 Searches for the early stage electromagnetic fields

**Charge-dependent directed flow measurements**

It is proposed in Refs. [787,788] that the strong initial state electromagnetic fields significantly affect the evolution of the produced QGP, thus leaving distinct imprints in the final state particle distributions. Charged particles propagating in the electromagnetic field experience a Lorentz force. Strong electric fields can also be generated due to the decay of the magnetic field with time through Faraday's law. In the centre-of-mass of the collision system the two effects (i.e., the Lorentz force and the force due to the induced electric field) have opposite sign; both effects are also antisymmetric in rapidity (i.e., rapidity-odd) and could be detected experimentally. This would, consequently, lead to the first experimental constrains on the value of the electric conductivity of the QGP.

A way to probe the electromagnetic fields and their rapidity-odd induced effects is the measurement of the charge-dependent anisotropic flow $v_n$. In particular, the charge dependence of the directed flow, $v_1(\eta)$, of the produced particles relative to the spectator plane is directly sensitive to the presence of the electromagnetic fields. The spectator plane is defined by the deflection direction of the collision spectators, measured with the Zero Degree Calorimeter (ZDC), and is strongly correlated with the direction of the magnetic field. As the production times for light (mainly produced through gluon-splitting processes throughout the entire evolution of the system) and heavy quarks (predominantly produced in early stage hard scattering processes) are vastly different, the directed flow of hadrons containing a c-quark is expected to be more sensitive to the early times of the system evolution, of the order of 0.1 fm/c.

Figure 66 presents the differences of the charge-dependent $v_1$, denoted as $\Delta v_1$, in mid-central Pb–Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV as a function of pseudorapidity for charged particles and $D^0$ mesons [789]. The magnitude of the effect is smaller for charged particles than for $D^0$ mesons, as expected from the fact that charm quarks are more sensitive probes of the magnetic field at earlier times, and is thus scaled up by a factor $10^3$ for visibility. The rapidity
<figure><img src="image_79.png" /><figcaption>**Fig. 66** The pseudorapidity dependence of $\Delta v_1$ for charged particles (blue markers) and D mesons (orange markers) in mid-central Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV [789]. The data points for charged particles are scaled by a factor $10^3$. Statistical and systematic uncertainties are represented by the error bars and the boxes around each data point, respectively</figcaption></figure>

slope $d\Delta v_1/d\eta$, extracted with a linear fit function, yields [1.68±0.49 (stat.) ±0.41 (syst.)]×10⁻⁴ for charged hadrons with $p_T > 0.2$ GeV/c for the 5–40% centrality interval and [4.9 ± 1.7 (stat.) ± 0.6 (syst.)] × 10⁻¹ for D mesons with $3 < p_T < 6$ GeV/c in a centrality interval of 10–40%, resulting in a significance of 2.6σ and 2.7σ for having a positive value, respectively. This apparent large slope could point to a stronger effect of the Lorentz force relative to the contribution from the induced electric field and the initial tilt of the source. This measurement constitutes the first experimental hint of the existence of the initial state electromagnetic fields at the LHC. The data samples that will be collected in the upcoming heavy-ion runs at the LHC will give us the potential to make a decisive step in this direction.

## Hyperon polarisation

An independent way of probing and constraining the value of the magnetic field is provided by the measurements of the hyperon global polarisation (see also Sect. 2.2.8).

Fluid vorticity, the main mechanism for the global polarisation effect, leads to the same polarisation between particles and antiparticles. The magnetic field, which is oriented in the same direction as the average vorticity, on the other hand, couples to the magnetic moment ($\mu$) of a particle $q$ and leads to a particle polarisation $P_q^B = \mu_q B/2T = Q_q B/(2m_q T)$, opposite for particles and antiparticles with charge $Q_q$. In the previous formula, $T$ is the temperature of the system at the time of the particle formation and $B$ is the value of the magnetic field. This opens up the possibility to probe the magnetic field with a different, independent observable.

The first observation of the finite global polarisation of $\Lambda$ hyperons in Au–Au collisions at RHIC [359] revealed

<figure><img src="image_80.png" /><figcaption>**Fig. 67** The $\sqrt{s_{NN}}$ dependence of the upper limit for the magnitude of the magnetic field at freeze-out, estimated from the global polarisation results of $\Lambda$ and $\bar{\Lambda}$. The data points are extracted from the results reported in Refs. [359,368]</figcaption></figure>

that the perfect liquid created in such collisions is also the most vortical one. The ALICE measurement of global polarisation at the LHC [368] although presently by large statistical and systematic uncertainties indicates a decrease with increasing collision energy, consistent with expectations based on the correlation between vorticity and the slope of directed flow as a function of centre-of-mass energy [254]. Although no significant difference was found in the polarisation between $\Lambda$ and $\bar{\Lambda}$, one can use the relationships $P_\Lambda \approx 0.5\omega/T + |\mu_\Lambda|B/T$ and $P_\Lambda \approx 0.5\omega/T - |\mu_\Lambda|B/T$, where $\omega$ is the vorticity of the QGP, to constrain the value of the magnetic field at freeze-out by evaluating ($P_\Lambda - P_{\bar{\Lambda}}$). Figure 67 presents the centre-of-mass energy dependence of such an estimate of the value of the magnetic field. At LHC energies, one obtains an upper limit of $eB/m_\pi^2$ of 0.017 (corresponding to $5.7 \times 10^{12}$ T) and 0.044 (or $14.4 \times 10^{12}$ T) at a 95% confidence level for Pb–Pb collisions at $\sqrt{s_{NN}} = 2.76$ and 5.02 TeV, respectively.

## 2.6.2 Searches for chiral anomalies

The ground state of the QCD vacuum consists of a superposition of topologically distinct states [790]. These states, characterised by a topological charge, are separated by a potential barrier and can be connected through tunneling transitions, called instantons [790–795]. However, due to the height of the potential barrier, which is of the order of the QCD scale over the strong coupling constant ($\Lambda_{QCD}/\alpha_s$), instanton transitions are suppressed [791,792]. At high temperatures, like the ones reached when the QGP is created in collisions between two heavy ions, the transitions between two states with different topological quantum numbers can take place by jumping over the barrier, and are known as sphalerons [796–800]. In electroweak theory, transitions
between vacuum states with different topological charges are associated with baryon and lepton number violation [801]. In QCD, they are connected to the non-conservation of chirality. The existence of these transitions is very well motivated by theoretical studies and would lead to the QGP that possessing domains with non-zero net chirality, with a sign that changes from event to event i.e., leading to an event-by-event local P and CP violation.

### Chiral magnetic effect studies

The possibility of observing parity violation in the strong interaction in relativistic heavy-ion collisions was first pointed out in Refs. [765–767] and was further discussed in Refs. [768–773]. A significant development happened in 2004, when it was noticed that in the presence of a chirally restored medium, a system of quarks with non-zero net-chirality immersed in a strong magnetic field leads to an excess of positively-charged particles moving along the magnetic field direction and an excess of negative particles moving in the opposite direction. This introduces a net electromagnetic current and creates an electric dipole moment of QCD matter. This phenomenon is called the Chiral Magnetic Effect (CME) [773], and its existence was recently reported in semi-metals like zirconium pentatelluride (ZrTe$_5$) [802].

The resulting charge separation can be identified by studying the P-odd sine terms in the Fourier decomposition of the particle azimuthal distribution [39,803] according to

$$ \frac{dN}{d\varphi_{\alpha}} \propto 1 + 2 \sum_{n} [v_{n,\alpha} \cos(n\Delta\varphi_{\alpha}) + a_{n,\alpha} \sin(n\Delta\varphi_{\alpha})], \quad (23) $$

where $\Delta\varphi_{\alpha} = \varphi_{\alpha} - \Psi_{RP}$ is the azimuthal angle $\varphi_{\alpha}$ of the particle of type $\alpha$ (either positively or negatively charged) relative to the reaction plane angle $\Psi_{RP}$. The coefficients $v_{n,\alpha}$ are the n-th order Fourier harmonics that characterise the anisotropies in momentum space. The leading order P-odd coefficient $a_{1,\alpha}$ reflects the magnitude of the effects from local parity violation, while higher order terms ($a_{n,\alpha}$ for $n > 1$) describe the specific shape in azimuth. The chiral imbalance that leads to the creation of the CME fluctuates from event to event, and the event average $\langle a_{1,\alpha} \rangle$ is consistent with zero, i.e. compatible with the observation of global parity conservation in strong interactions. Consequently, the effect can be detected only by correlations studies.

In Ref. [803], it was suggested that a suitable way to probe the CME is via two-particle correlations measured relative to the second-order symmetry plane of the form $\gamma_{1,1} = \langle \cos(\varphi_{\alpha} + \varphi_{\beta} - 2\Psi_{2}) \rangle$, where the brackets indicate an average over all events, and $\alpha$ and $\beta$ denote particles with the same or opposite electric charge. The advantage of using this observable is that it probes correlations between two leading order P-odd coefficients $a_{1,\alpha}$ and $a_{1,\beta}$ which do not average to zero over all events, while suppressing correlation that do

not depend on the symmetry plane orientation. In order to independently evaluate the contributions from correlations in- and out-of-plane, one measures at the same time a two-particle correlator of the form $\delta_1 = \langle \cos(\varphi_{\alpha} - \varphi_{\beta}) \rangle$. Both correlators can be generalised according to:

$$ \gamma_{m,n} = \langle \cos(m\varphi_{\alpha} + n\varphi_{\beta} - (m+n)\Psi_{|m+n|}) \rangle, \quad (24) $$

$$ \delta_m = \langle \cos[m(\varphi_\alpha - \varphi_\beta)] \rangle, \quad (25) $$

where m and n are integers.

The observable of Eq. 24 is constructed as the difference between correlations in- and out-of-plane and is thus expected to suppress background effects approximately by a factor of $v_2$. On the other hand, the correlator of Eq. 25, owing to its construction, is affected (if not dominated) by background contributions. This background was recently identified as stemming from an interplay between local charge conservation embedded in an environment that exhibits azimuthal anisotropy [804]. Its isolation and the subsequent quantification of the CME contribution to these observables has been the main focus of this line of research at both RHIC and LHC.

Figure 68 presents the centrality dependence of $\Delta\delta_1$ and $\Delta\gamma_{1,1}$ i.e., the difference of $\delta_1$ and $\gamma_{1,1}$ between opposite- and same-sign pairs in panels (a) and (b), respectively. The data points that correspond to results from the STAR Collaboration in Au–Au collisions at $\sqrt{s_{NN}} = 200$ GeV [108,805] and the ones from Pb–Pb collisions at $\sqrt{s_{NN}} = 2.76$ and 5.02 TeV [806,807] have a characteristic and significant centrality dependence. The correlator $\Delta\delta_1$ is related to the balance function also studied at the LHC [335,808]. The results of $\Delta\delta_1$ are qualitatively consistent with the ones in Refs. [335,808], i.e. oppositely charged particles are more tightly correlated in central events resulting in a narrowing of the balance function width in $\Delta\varphi$ which results in the observed centrality evolution of $\Delta\delta_1$. There is a clear evolution with centre-of-mass energy for $\Delta\delta_1$, consistent with what is reported in Ref. [335]. This energy dependence, however, is less pronounced for $\Delta\gamma_{1,1}$, between the two LHC energies. This lack of a significant energy dependence for $\Delta\gamma_{1,1}$ initially came as a surprise considering the differences in centre-of-mass energy, which between RHIC and LHC is one order of magnitude, and that the magnitude of the magnetic field and the way it evolves is, in principle, different between the two energies. In addition, the particle density is almost three times larger at the LHC compared to RHIC. This might lead to a significant dilution of the CME signal, if any, at this higher energy. Finally, the acceptance used in the two experiments is slightly different, with ALICE using a smaller range in $\eta$. This, in turn, could lead to a more pronounced contribution from both any potential CME signal and background. At the same time, however, there is no significant energy dependence in the effects that constitute the
**Fig. 68** The centrality dependence of $\Delta\delta_1$ (a) and $\Delta\gamma_{1,1}$ (b), measured at RHIC [108,805] and LHC energies [806,807]. Statistical and systematic uncertainties are represented by the error bars and the boxes around each data point, respectively. Comparison of the centrality dependence of $\Delta\delta_1$ (c) and $\Delta\gamma_{1,1}$ (d), measured at the highest LHC energy [807] and estimated with a blast wave [810] and the AVFD [786,811] models (see text for details)

<figure><img src="image_81.png" /></figure>

background to these measurements. In particular, preliminary studies indicate that the correlations between balancing charges, as reflected in the width of the balance function, do not exhibit any significant dependence on collision energy. Furthermore, the values of $v_2$ measured in Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV are around 2% higher than the values at $\sqrt{s_{NN}} = 2.76$ TeV [809].

Figure 68c and d present the comparison between the experimental data points of both $\Delta\delta_1$ and $\Delta\gamma_{1,1}$ in Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV with expectations from a blast wave (BW) model based on a parameterisation from Ref. [810]. Local charge conservation (LCC) is additionally incorporated by generating ensembles of particles with zero net charge produced at the same spatial location, uniformly distributed over the anisotropic particle radiating source. The input parameters of the model are tuned to describe the $p_T$ spectra [402] and the $p_T$-differential $v_2$ values [290] for charged pions, kaons, and protons (antiprotons) measured in Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV (see Ref. [807] for details). In addition, the number of sources that emit balancing pairs is tuned separately for each centrality interval to reproduce the centrality dependence of $\Delta\delta_1$, as presented in the panel (c) of Fig. 68. The tuned model is then used to extract the expectation for the centrality dependence of $\Delta\gamma_{1,1}$, shown in the panel (d) of Fig. 68. The BW model underestimates the measured data points by as much as ≈ 40%, with the disagreement increasing progressively for more peripheral events. A similar study performed in Xe–Xe collisions at $\sqrt{s_{NN}} = 5.44$ TeV [812], revealed that the data points can be described quantitatively by the tuned BW model.

The same figure presents the results from the Anomalous Viscous Fluid Dynamics (AVFD) model [786,811]. This model simulates the evolution of the chiral fermion currents in the QGP on top of the hydrodynamic evolution of a heavy-ion collision as prescribed by VISHNU [460] which couples 2+1 dimensional viscous hydrodynamics (VISH2+1) to a hadron cascade model (UrQMD) [124]. The model allows one to tune the value of the axial current density to entropy ratio ($n_5/s$) that controls the imbalance between left- and right-handed quarks. In addition, it includes local charge conservation effects by emitting particle pairs of opposite charges from the same fluid element. The procedure that was followed relied on finding the proper value of both $n_5/s$ and the percentage of balancing charged particles (LCC) to reproduce simultaneously the centrality dependence of $\Delta\delta_1$ and $\Delta\gamma_{1,1}$. This is illustrated in panels (c) and (d) of Fig. 68, where the AVFD curves are represented by the green bands [813]. The model is able to describe the measurements with a set of input parameters that ranges between 0.03 and 0.06 for $n_5/s$ throughout all centrality intervals while the percentage of charged particles emitted as part of an oppositely charged pair relative to the total multiplicity of an event varies from 30% up to 60% with an increasing trend for more central events. In addition, a background only scenario, i.e., $n_5/s = 0$ for all centralities, is not compatible with the measurements according to this model [813]. Finally, a similar study performed in Xe–Xe collisions at $\sqrt{s_{NN}} = 5.44$ TeV [812] has yielded systematically lower values of $n_5/s$, and a background-only scenario being, at this time, compatible with the data [813].
<figure><img src="image_82.png" /><figcaption>**Fig. 69** (Top) The dependence of $\Delta\gamma_{1,1}$ on $v_2$ for shape selected events together with a linear fit (dashed lines) for various centrality classes. (Bottom) The centrality dependence of the CME fraction, $f_{\text{CME}}$, extracted from the slope parameter of fits to data and MC-Glauber, MC-KLN CGC, and EKRT models, respectively. The dashed lines indicate the physical parameter space of the CME fraction, while the points are slightly shifted along the horizontal axis for better visibility. Figures from Ref. [814]</figcaption></figure>

<figure><img src="image_83.png" /><figcaption>**Fig. 70** Summary of the results for the CME limit obtained from different analyses performed at various LHC energies and colliding systems integrated over centralities (see text for details)</figcaption></figure>

Data-driven methods could offer a more robust approach to constrain and quantify the dominating background while in parallel isolating the signal that comes from the CME. The first measurement in this direction was presented by the ALICE Collaboration in Ref. [814], using an innovative method proposed and developed in [815]. This method, called Event Shape Engineering (ESE), utilises the strong fluctuations of initial geometry (i.e., the position of participating nucleons) even at a fixed impact parameter and allows one to select events with different initial system shapes, e.g., central Pb–Pb collisions with large initial anisotropy. This allows one to select events where the main component of the background, the value of $v_2$, can be varied. The top panel of Fig. 69 illustrates the linear dependence of $\Delta\gamma_{1,1}$ on $v_2$ for different shaped events and centrality ranges, indicative of the dominance of background contributions. This becomes more evident seeing all data points converging on a single line when $\Delta\gamma_{1,1}$ is scaled with the relevant particle density for each centrality interval to account for the dilution of these correlations stemming from the increased multiplicity [814]. The bottom panel of Fig. 69 presents the CME fraction, $f_{\text{CME}}$ which is extracted from the slope parameters of fits to data and MC-Glauber, MC-KLN CGC, and EKRT models, for various centrality intervals of Pb–Pb collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV [814]. All data points corresponding to different initial state models (MC-Glauber, MC-KLN CGC, and EKRT), are consistent with zero within the uncertainties. This study, by combining the data points in the 10–50% centrality interval, set an upper limit of 26–33% at 95% confidence level for the CME signal contribution to $\Delta\gamma_{1,1}$.

Figure 70 presents a summary of the upper limits for the CME signal contribution to $\Delta\gamma_{1,1}$ obtained from different analyses performed at various LHC energies and colliding systems integrated over centralities. These include the data points from studies using the ESE technique [814] for the three initial state models used, analysis of the correlators including higher harmonics at the two different LHC energies [807], and from the comparison of the correlations in Xe–Xe and Pb–Pb collisions in the two-component model [812]. The higher harmonics results are based on measuring charge-dependent correlations relative to the third order symmetry plane of the form $\gamma_{1,2} = \langle\cos(\varphi_{\alpha} + 2\varphi_{\beta} - 3\Psi_3)\rangle$. Considering that the third order symmetry plane ($\Psi_3$) is very weakly correlated with $\Psi_2$ [816], the charge separation effect relative to the third harmonic symmetry plane is expected to be negligible and the relevant correlations are expected to reflect mainly, if not solely, background effects. In this background-only scenario, one could approximate $\Delta\gamma_{1,1}$ and $\Delta\gamma_{1,2}$ according to $\Delta\gamma_{1,1} \propto \kappa_2\Delta\delta_1 v_2$ and $\Delta\gamma_{1,2} \propto \kappa_3\Delta\delta_1 v_3$, with $\kappa_2$ and $\kappa_3$ being proportionality constants [807]. The method relies heavily on the assumption that $\kappa_2 \approx \kappa_3$ if both correlators are studied in the same experimental setup and for the same kinematic region and provides an
upper limit of 15–18% at 95% confidence level for the 0–40% centrality interval. Finally, the analysis of Xe–Xe collisions at $\sqrt{s_{NN}}$ = 5.44 TeV [812] revealed similar results for $\Delta\gamma_{1,1}$ as the ones in Pb–Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV at the same multiplicity. The CME fraction was estimated in this analysis using a two-component model, similar to the one proposed in Ref. [817]. The first component is associated to the CME contribution in the measurement of $\Delta\gamma_{1,1}$ which was considered to be correlated with the value of the magnetic field, modelled with a similar prescription as the one followed in Ref. [814]. On the other hand, the background contribution to $\Delta\gamma_{1,1}$ was associated to the values of $v_2$ in the two systems [279] and was scaled by the corresponding charged particle multiplicity density [176] to account for dilution effects. The study resulted in an upper limit of around 2% and 25% at 95% confidence level for the 0–70% centrality interval in Xe–Xe and Pb–Pb collisions, respectively.

Searches for the Chiral Magnetic Effect have been also pursued by the CMS Collaboration at the LHC and the STAR Collaboration at RHIC. In particular the CMS Collaboration also used the Event Shape Engineering technique, obtaining somewhat lower limits on the CME contribution [818]. Considering that the CMS Collaboration did not account for a possible signal dependence on the elliptic flow value, their upper limits are too low and accounting for this effect might lead to results fully compatible with the ALICE measurements. The CMS Collaboration [818,819] and the STAR Collaboration [820] also attempted to put limits on the CME signal by comparing the measurements in Pb–Pb and Au–Au collisions with corresponding measurements in p–Pb and p(d)–Au collisions, where the CME signal contribution to $\Delta\gamma_{1,1}$ correlator is expected to be negligible due to a very weak correlation between the direction of the magnetic field and the elliptic flow plane. Unfortunately, the measurements of the $\gamma_{1,1}$ correlator in small systems are also plagued by the reaction plane independent background (see Refs. [108,805]), without reliable estimate of which there is no possibility to make a definite conclusion from such measurements. The STAR Collaboration has recently published results [821] from the analysis of the isobar samples collected in 2018. The two isobar nuclei used, $^{96}_{44}$Zr and $^{96}_{40}$Ru, were expected to result into similar flow-driven background contribution due to the same mass number. The difference, however, in the atomic number between the two nuclei would be reflected in a smaller magnitude of the magnetic field in the case of $^{96}_{40}$Ru compared to $^{96}_{44}$Zr, which would translate into a difference in the CME contribution between these two systems. Comparing the outcome of this blind analysis to a set of predefined CME-signal criteria they concluded that the results are consistent with no CME signal. Further work is needed to account for the differences of the background in the two systems, reflected by the corresponding differences

in both the multiplicity and flow harmonics at the same centrality, in order to estimate the final limits.

### Chiral magnetic wave studies

The presence of a net positive electric charge can induce a positive axial current along the direction of the magnetic field i.e., leading to flow of chirality. This is caused by the Chiral Separation Effect (CSE) [822]. The coupling between the CME and the CSE leads to a wave propagation of the electric charge, resulting in an electric charge quadrupole moment of the system, the Chiral Magnetic Wave (CMW) [823–825].

The azimuthal distribution of charged particles due to the presence of the CMW can be written as

$$ \frac{dN^{\pm}}{d\varphi} = N^{\pm}[1 + (2v_2 \mp rA) \cos(2(\varphi - \Psi_2))], \quad (26) $$

where $A = (N^+ - N^-)/(N^+ + N^-)$ is the charge asymmetry, and $r$ is the parameter that encodes the strength of the electric quadrupole due to the CMW. Therefore, one can probe the value of $r$ by measuring the $v_2$ values for different charges as a function of the charge asymmetry. Experimentally, the latter is measured in a specific kinematic region and should be corrected for detector efficiency, leading to the introduction of additional systematic uncertainties. Instead, in Ref. [826], the authors suggested to measure the covariance of $v_\eta$ and $A$ that is a robust observable and does not depend on detector inefficiencies. The first results for this observable measured in Pb–Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV were reported in Ref. [827]. The correlator for the second harmonic, $\langle v_2\{2\}A \rangle - \langle A \rangle\langle v_2\{2\} \rangle$, exhibits a significant centrality dependence compatible with expectations from a CMW signal which becomes more pronounced for more peripheral events where the magnetic field strength increases. On the other hand, local charge conservation effects could also contribute to the development of the measured signal [826].

One direction that could help in identifying whether the CMW contributes significantly to the measured correlations is to perform this study differentially, e.g. as a function of the pseudorapidity gap between the particle for which elliptic flow is measured and another charged particle. Such measurements are expected to differentiate contributions from the CMW signal and the background, especially those accompanied by higher harmonics measurements. The expectation from higher harmonics results is that the CMW contribution will be significantly suppressed due to the symmetry of the CMW effect [827]. The left panel of Fig. 71 presents the results for the second harmonic as a function of $\Delta\eta$. These correlations exhibit a peak with a “typical hadronic width” of about 0.5–1 units of rapidity which qualitatively agrees with possible background contributions from local charge conservation combined with strong radial and elliptic flow [826]. In parallel, noticeable correlations in the third harmonic, presented in the right panel of Fig. 71, indicate a significant contribution from LCC, but the final conclusion requires a
<figure><img src="image_84.png" /><figcaption>**Fig. 71** The pseudorapidity dependence of the results for the second (left) and third harmonic (right) for positive (red squares) and negative (blue circles) particles measured in the 20–60% centrality interval of Pb–Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV [827]</figcaption></figure>

detailed model comparison. Future $p_T$-differential measurements as well as studies employing the ESE technique would also be very helpful to quantify the background and provide a CMW limit. Note that the contribution of the LCC to this correlator and the $\gamma$ correlator used for CME studies are strongly correlated, which can be further used for quantitative characterisation of the background in both measurements.

### 2.6.3 Conclusions

*Electromagnetic fields in heavy-ion collisions.* Measurements of directed flow for oppositely charged particles as well as for $D^0$ and $\bar{D}^0$ provide the first indication that the early stage electromagnetic fields can affect the motion of the final state particles. The differences in the measured global polarisation of $\Lambda$ and $\bar{\Lambda}$ provide an upper limit for the magnitude of the magnetic field at freeze-out of $5.7 \times 10^{12}$ T and $14.4 \times 10^{12}$ T at a 95% confidence level in Pb–Pb collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV, respectively.

*Chiral Magnetic Effect.* Direct studies for the existence of the CME in heavy-ion collisions revealed that background effects are dominating at the LHC. The contribution from the CME to the measurement of charge dependent correlations relative to the second order symmetry plane ($\gamma_{1,1}$) is constrained to an upper limit of 15–33% in Pb–Pb and 2% in Xe–Xe collisions at LHC energies at 95% confidence level. The combination of measurements and model studies in Pb–Pb and Xe–Xe collisions indicates that there is a bigger potential to discover the CME in large collision systems (e.g. Pb–Pb) than in their small collision counterparts, mainly due to the significant differences they exhibit in the value of the magnetic field.

*Chiral Magnetic Wave.* The combination of the second and third harmonic results for CMW studies indicates a significant

background contribution from local charge conservation.

## 2.7 Quantitative characterisation of the QGP at the LHC

### 2.7.1 Macroscopic properties

Macroscopic properties pertain to the global features of the system, such as the temperature, volume and system lifetime. The system formed in heavy-ion collisions undergoes a very rapid radial expansion with velocity up to $v \sim 2/3c$ on a timescale of $10^{-23}$ s. Therefore, each derived macroscopic property will either correspond to some average over this evolution, or to a particular stage. On the other hand, quantities subject to conservation laws such as the total energy or the entropy (assuming dissipative effects are small) will remain constant. In this section, we will review some of these properties that are determined either directly from ALICE measurements, or from theoretical models that are validated by measurements. We will then compare the results to predictions from various theoretical descriptions of the QGP and hadronic stages.

*Temperature evolution of the system.* The deconfined system created in heavy-ion collisions will cool and transition to a hadronic phase as time progresses, with various experimental probes being sensitive to the temperature evolution. In Fig. 72, we show the temperature ranges probed using ALICE measurements from heavy-ion collisions at the LHC, which have varying dependencies on the uncertainties from measurements and theoretical modelling. As shown previously in this chapter, hydrodynamic model frameworks successfully describe many low-$p_T$ observables such as the anisotropic flow coefficients and average $p_T$ of identified particles. These models can be used to determine QGP temperature throughout the system evolution, and therefore, using
constraints from data, they can provide an indirect assess-
ment of the entire temperature range of the QGP created
at the LHC. This range is shown in Fig. 72 for the IP-
Glasma+MUSIC framework [48]. It indicates that at the ear-
liest times, which correspond to the highest temperatures
of the QGP, the temperatures obtained in Pb–Pb collisions
are ~ 0.8 GeV, that is roughly five times higher than the
pseudocritical QGP transition temperature predicted by lat-
tice QCD, T<sub>pc</sub> ≈ 0.155 GeV [119,169]. The hydrodynamic
range is dependent on the choice of initial-state modelling
and the pre-equilibrium descriptions, as well as on the time
the hydrodynamic stage occurs, and the lattice-QCD equa-
tion of state used. The T<sub>RENTO</sub> + VISHNU chain [49] (not
shown in Fig. 72) has a smaller temperature range, with an
upper limit of T = 0.35 GeV. Among other things, this is
related to the time the hydrodynamic evolution starts. This
occurs at a time τ<sub>0</sub> ~ 1 fm/c for T<sub>RENTO</sub> + VISHNU,
and at τ<sub>0</sub> ~ 0.3 fm/c for IP-Glasma+MUSIC. Other mod-
els have also been used to predict the temperatures needed to
describe the relative bottomonia Υ(2S) and Υ(1S) yields (see
Sect. 2.5). The respective model calculations, one of which
uses hydrodynamics [739] and the other is a non-perturbative
transport model [737], use starting temperatures (T<sub>0</sub>) of 0.55–
0.64 GeV. This range is shown in Fig. 72, and can be used as
an additional input to hydrodynamic descriptions of low-p<sub>T</sub>
light-flavour observables to constrain the temperature range
probed by heavy-ion collisions.

Effective temperatures (T<sub>eff</sub>) derived from prompt-photon
and charged-hadron spectra are shown in Fig. 72 as well. The
photon exponential inverse slope is obtained directly from
data (see Sect. 2.1), and due to effects related to the emission
time and blue shift from the expanding source, provides an
effective temperature averaged over the whole system life-
time. The effective charged-hadron temperature was obtained
using a hydrodynamic model [289] to describe ALICE spec-
tra measurements. It can also be interpreted as an effective
temperature averaged over the entire space–time evolution
of the medium. It is lower than the photon effective temper-
ature, which might be a consequence of different sensitiv-
ities to each stage of the medium evolution. The chemical
freeze-out temperature (T<sub>chem</sub>) range obtained from a sta-
tistical model fit to particle yields measured by ALICE (see
Sect. 2.3) is shown in Fig. 72. This temperature is close to the
deconfinement crossover temperature computed from lattice
QCD. Finally, the kinetic freeze-out temperature (T<sub>kin</sub>) range
derived from Blast-Wave fits to spectra and identified particle
v<sub>n</sub> measurements are shown. It is found to be rather highly
model dependent (see Sect. 2.2). The kinetic freeze-out can
either occur around the time of the chemical freeze-out or at
later times that correspond to lower temperatures. We will
investigate this difference by reviewing additional measure-
ments in the next subsection.

<figure><img src="image_85.png" /><figcaption>**Fig. 72** The temperatures ranges probed by central heavy-ion collisions at the LHC derived from ALICE measurements shown throughout the chapter. The region to the right of the red line represents the QGP phase, while the left is the hadronic phase. A lattice QCD calculation for the pseudocritical temperature *T*<sub>pc</sub> and crossover temperature range where the QGP and hadron gas coexist (denoted by the arrows) is also shown [119]</figcaption></figure>

*Space–time extent.* The space–time extent of the system cre-
ated in heavy-ion collisions refers to the properties of the sys-
tem towards the end of its evolution. It can be characterised by
a hypersurface. Such a surface encodes the four-dimensional
properties of the corresponding freeze-out space–time coor-
dinates of the system at the freeze-out temperature. It depends
on the position, time, and velocity of the element that freezes
out. An example of hypersurface profiles from a hydrody-
namic model chain can be found in Ref. [250]. The freeze-
out time on this surface will be longer at freeze-out positions
closer to the centre of the originally created QGP, compared
with positions at the outskirts. The reason for this is that
at the centre, the temperatures are the highest, therefore the
time associated with freeze-out is longer. As the velocity of a
medium element increases, the freeze-out positions relative
to the collision centre decrease, which means higher veloc-
ity elements decouple from the system at earlier times. This
is directly observed from measurements of the femtoscopic
radii in Sect. 2.1, which decrease as the femtoscopic particle
pair velocity (characterised by *k*<sub>T</sub>) increases.

Table 2 shows a summary of the obtained volumes that
corresponds to one unit of rapidity at midrapidity for central
Pb–Pb collisions, at about the freeze-out temperatures shown
in this section. The volume associated with T<sub>chem</sub> is deter-
mined from statistical models of the identified particle yields
described in Sect. 2.3, which corresponds to a static vol-
ume.<fnref n="10" /> As noted in the previous subsection, the correspond-
ing T<sub>chem</sub> temperature is consistent with the deconfinement
temperature. This implies two things: the QGP is born in ther-

<footnote n="10">¹⁰ The four model values shown in Sect. 2.3 are averaged with the
uncertainty determined from the variance.</footnote>
<table><caption>**Table 2** System volumes in central Pb–Pb at $\sqrt{s_{NN}} = 2.76$ TeV collisions for the chemical and kinetic freeze-out temperatures</caption><thead><tr><th colspan="3">Volume: dV/dy at y = 0</th></tr></thead><tbody><tr><td>Static at T<sub>chem</sub> = 156 MeV<br>4730 ± 500 fm<sup>3</sup></td><td>Dynamic at T<sub>kin</sub> = 156 MeV<br>6210 ± 360 fm<sup>3</sup></td><td>Dynamic at T<sub>kin</sub> = 100 MeV<br>7960 ± 460 fm<sup>3</sup></td></tr></tbody></table>

mal equilibrium, and the hadronic inelastic cross sections are too small for the system to maintain chemical equilibrium after the QGP hadronises and the temperature decreases (otherwise $T_{chem} < T_{pc}$). This volume can therefore be interpreted as the volume when the QGP hadronises. The volumes associated with $T_{kin}$ are obtained using measurements of femtoscopic radii shown in Sect. 2.1 at $k_T = 0.25$ GeV/c and for the two temperatures shown. They are determined using $V_{fo} = (2\pi)^{3/2} R_{side}^2 R_{long}$. From the temperature $T_{kin}$, one can determine the mean thermal velocity and the corresponding rapidity interval [828]. In addition, as these pairs are moving, the volume they characterise is reduced compared with the static case. It has been estimated that the transverse expansion leads to a 20–25% reduction in $R_{side}$ [419,829,830]. For the case where $T_{chem} = T_{kin}$, which corresponds to the largest values of $T_{kin}$ from the Blast-Wave fits shown in Sect. 2.2, the dynamic volume is larger than the static volume. Given the reverse is expected if $T_{chem} = T_{kin}$, this then implies $T_{kin} < T_{chem}$. Therefore, based on the available phenomenological models, the hypothesis of $T_{chem} = T_{kin}$, is unlikely, implying the existence of a prolonged time period when the system is still coupled via elastic scatterings. For $T_{kin} = 100$ MeV, which corresponds to the smallest values of $T_{kin}$ from the Blast-Wave fits, these results imply at least a ~70% increase in the kinetic freeze-out volume compared with chemical freeze-out.

An estimate of the lifetime of the system is given by the decoupling (or kinetic freeze-out) time of pions, the most abundant species produced. This decoupling time is in the range 10–13 fm/c in central collisions, based on the measurements of $R_{long}$ discussed in Sect. 2.1, and on the assumption of a hydrodynamic expansion. The hadronic stage of the collision can be further investigated by resonance and femtoscopic measurements, which are summarised in Table 3 and discussed in more detail in Sects. 2.3 and 2.1, respectively. In both cases, the data can be explained by the modelling of the hadronic cascade [71,251] preceded by a hydrodynamic evolution of the system. The lifetimes derived from hadron resonances offer an estimate of the duration of the hadronic stage, assuming negligible regeneration. While all the times are above zero, there are considerable differences. This could reflect a varying competition between the rescattering and regenerating processes for the various resonances – if re-scattering dominates, one would infer shorter lifetimes, however if regeneration contributes significantly, this would imply longer lifetimes due to the delay in the formation of regenerated resonances.

The delay in the emission of $\pi$ and K hadrons was measured with femtoscopic techniques, and was associated with $K^{0*}$ formation and decay via comparison to calculations from a hadronic afterburner, as discussed in Sect. 2.1. We finally note that this observation also requires a reconciliation with the results from the FastReso Blast-Wave model shown in Sect. 2.2, which implies instead the possibility of $T_{chem} = T_{kin}$, namely the absence of such a hadronic stage. While the exact timescale of the hadronic stage remains an open question, evidence based mainly on femtoscopic and resonance measurements suggest indeed that the fireball continues to behave like a coupled system for an extended time after the chemical freeze-out.

### 2.7.2 Microscopic properties

The transport properties of the QGP provide information on how it responds to microscopic excitations, and such excitations can arise from either soft or hard processes. Before their review, it is useful to provide a distinction between the weak and strong coupling regimes that these transport properties probe:

- **Weakly-coupled system.** In this case, the QGP can be considered as a system of well defined quasi-particles (quarks and gluons), whose interactions involve few body processes e.g. $2 \to 2$ scattering or $2 \to 3$ collision-induced gluon radiation, which sometimes are modeled as scatterings off discrete scattering centres. It can be described by perturbative QCD at couplings of $\alpha_S \lesssim 0.3$, and is considered a gas.

- **Strongly-coupled system.** The coupling between QGP constituents is large and therefore dominated by higher-order processes. Such a coupling induces strong correlations between the neighboring constituents, therefore the quasi-particle description is no longer valid, and QGP matter can be treated as a liquid.

It is also important to note that these regimes represent idealised limits - an intermediately coupled system where perturbative and non-perturbative processes compete is of course plausible. A special case of a strongly-coupled system is at infinite coupling. While physically unrealistic, it provides a limiting case for strongly-coupled systems, as will become apparent later in this section. One of the defining features of QCD is the scale dependence e.g. asymptotic freedom. The
<table><caption>**Table 3** Estimates of the kinetic freeze-out time from pion decoupling and of the duration of the hadronic stage implied from resonances and delay in particle emissions from femtoscopic measurements in central Pb–Pb collisions at $\sqrt{s_{NN}} = 2.76$ TeV</caption><thead><tr><th>Kinetic freeze-out time</th><th>Lifetimes of the hadronic stage</th><th></th><th></th><th></th></tr></thead><tbody><tr><td>π decoupling</td><td>K<sup>0*</sup>/K<sup>±</sup></td><td>ρ<sup>0</sup>/π<sup>±</sup></td><td>Λ<sup>*</sup>/Λ</td><td>K − π emission delay</td></tr><tr><td>10–13 fm/c</td><td>~ 3 fm/c</td><td>~ 1 fm/c</td><td>~ 10 fm/c</td><td>1.4 ± 0.6 fm/c</td></tr></tbody></table>

processes associated with a transport property that involve large momentum transfers resolve medium structures at short length scales and will by definition “see” a weakly-coupled QGP, with the opposite being true for processes involving small momentum transfers. In addition, as the temperature of the QGP increases, all transport properties are expected to probe an increasingly-weakly-coupled QGP. Therefore, exploring the temperature dependence of transport properties is critical in establishing whether the QGP remains strongly coupled at the temperatures accessible at the LHC, or whether a weakly-coupled QGP emerges at the highest temperatures possible in the laboratory.

**Viscosity of the QGP.** The success of hydrodynamic model chains in describing a wide variety of soft observables, implies that the system is strongly coupled at momentum scales corresponding to the QGP temperature, as hydrodynamics describes the evolution of the QGP in terms of a liquid. The question then becomes quantitative – how strong is the coupling – and values of the shear ($\eta/s$) and bulk ($\zeta/s$) viscosities per entropy density play a central role in addressing this question. In the strong-coupling picture, both transport properties $\eta/s$ and $\zeta/s$ are proportional to the shear and bulk relaxation times. As the initial state is highly non-uniform, this leads to large non-equilibrium shear and bulk excitations at early times. The shear and bulk relaxation times determine how quickly the system can respond to achieve equilibrium – if the coupling is large, these times will be small, which corresponds to small values of $\eta/s$ and $\zeta/s$, and vice versa.

Measurements of anisotropic flow and hadron spectra by ALICE have allowed for an exploration of the temperature dependence of $\eta/s$ and $\zeta/s$. Figure 73 compares these dependencies for four hydrodynamic model chain calculations constrained by ALICE data. Regarding $\eta/s$, many other fluids exhibit a temperature dependence for $\eta/s$, with a minimum occurring at the phase transition temperature between the gas and liquid phases. Whether such a temperature dependence exists for the QGP phase is an open question for the temperatures probed by the LHC, as demonstrated in Fig. 73. The models assume either an independence (IP-Glasma+MUSIC and EKRT [309]) or a weak dependence ($T_R\text{ENTO}+\text{VISHNU}$). All provide fair descriptions of the anisotropic flow measurements. The correspond-

<figure><img src="image_86.png" /><figcaption>**Fig. 73** The temperature dependence of the shear (top panel) and bulk (bottom) viscosities over entropy density in the QGP phase constrained by the ALICE measurements shown in Sects. 2.1, 2.2, and 2.5 from various hydrodynamic models described in the text. Limits from pQCD [51], AdS/CFT [107], and AdS/Non-Conformal Holographic [833] approaches are also shown. The ranges on the right of the plot represent 90% posterior intervals from the Bayesian analyses</figcaption></figure>

ing shear relaxation time range is $\tau_\pi = 0.15–0.40$ fm/c at $T = 0.5$ GeV, which implies that spatial anisotropies from the initial state are rather quickly smoothed out in the QGP phase. Two coupling limits are shown: an AdS/CFT limit which is calculated for infinite coupling, and a pQCD limit, which is determined in leading order for $\alpha_S = 0.3$. In the weak-coupling picture, $\eta/s$ is inversely related to the coupling strength – therefore the $\alpha_S$ chosen corresponds to the strongest possible coupling (therefore the lowest $\eta/s$) in that scheme. It is clear that the extracted $\eta/s$ values are closer to the infinite-coupling limit. Nonetheless, it should be pointed out that next-to-leading-order corrections for the pQCD determination of $\eta/s$ are very large [831]. We refrain from showing them here, as these corrections largely cancel when ratios of transport parameters are explored, which will be the subject of discussion later. The values of $\eta/s$ from the QGP are roughly four times smaller than for Helium after it transitions to a superfluid [832].

The constraints on $\zeta/s$ from ALICE measurements are shown in the bottom panel of Fig. 73. Predictions from an infinitely-coupled AdS Non-Conformal Holographic approach (AdS/NCH) are also shown for comparison. As the applicability of Conformal Symmetry regarding the strong
potential assumed in the AdS/CFT scheme at temperatures close to the $T_{pc}$ is expected to break down, an alternative approach is needed to determine $\zeta/s$. The breaking of con- formal symmetry leads to $\zeta/s$ rising near $T_{pc}$ (it is zero oth- erwise). This approach was also used to reevaluate the $\eta/s$ in the limit of infinite coupling at all temperatures, and was found to also give $1/4\pi$, which suggests this limit is univer- sal. Its prediction that $\zeta/s$ should depend strongly on the tem- perature in this region is utilised for the $T_{RENTo}+VISHNU$ and IP-Glasma+MUSIC models. The $T_{RENTo} + VISHNU$ model provides the best description of ALICE identified- particle mean-$p_T$ measurements (as shown in Sect. 2.2), to which $\zeta/s$ is sensitive. This suggests that the ranges from both IP-Glasma+MUSIC and EKRT ($\zeta/s = 0$) provide a conservative estimate of the uncertainty for this parameter. The high-temperature pQCD limit for $\zeta/s$ is close to 0, which appears to apply for all the models shown at temperatures above 0.4 GeV. This then implies that bulk excitations in the initial state are washed out in the QGP phase even more quickly than the shear excitations e.g $\tau_\pi < 0.1$ fm/c for IP-Glasma+MUSIC at $T = 0.4$ GeV. We also note that the validity of the initial state models used in each hydrodynamic model chain will be investigated in Sect. 4, using multiplicity and anisotropic-flow measurements that are mainly sensitive to the features of the initial state.

In addition, the posterior distributions for $\eta/s$ and $\zeta/s$ have been evaluated using Bayesian parameter estimation techniques on ALICE data. They have been carried out by the Duke [49],<fnref n="11" /> JETSCAPE [556], Trajectum [834], and Jyväskylä [835] groups. These are shown on the right of Fig. 73, at $T = 0.3$ GeV for $\eta/s$ and $T = 0.2$ GeV for $\zeta/s$. The size of these posterior ranges are influenced by the prior ranges and data-sets included. For example, the JETSCAPE prior ranges were larger than those by the Duke and Tra- jectum groups, and yielded a larger upper limit for $\eta/s$. The Jyväskylä group used an even larger prior range, but included measurements of Symmetric Cumulants, which, as men- tioned in Sect. 2.2, highly constrain the temperature depen- dence of $\eta/s$ and reduce the upper limit. The $\zeta/s$ ranges differ more than the $\eta/s$ ranges. The JETSCAPE group also found that the duration of the pre-equilibrium phase has a strong impact on the extracted viscosity transport parameters. This is clearly demonstrated by the v-USPhydro chain [836]; the hydrodynamic evolution starts without any pre-equilibrium phase, and therefore requires low values of $\eta/s = 0.05$ and $\zeta/s = 0$ in order to develop enough flow to describe the ALICE data [837]. Those values are also shown in Fig. 73.

<figure><img src="image_87.png" /><figcaption>Fig. 74 The temperature dependence of the spatial charm diffusion coefficient $D_s$ in the QGP phase constrained by the ALICE measure- ments shown in Sect. 2.4 from various transport models described in the text. The range on the right of the plot represents a 90% poste- rior interval from the Bayesian analysis at $T = 0.4$ GeV. Limits from pQCD [69], and AdS/CFT [107], are also shown</figcaption></figure>

Charm spatial diffusion coefficient $D_s$. Heavy quark trans- port models attempt to describe how charm and beauty quarks interact with the QGP, given that heavy quarks are pro- duced early and therefore initially out of equilibrium with the system. As described in the Introduction and Sect. 2.4, the key quantity for investigating this processes is the charm relaxation time $\tau_{charm}$, which can be evaluated in the low- momentum limit ($p = 0$). It is proportional to the charm diffusion coefficient $D_s$. In the strong coupling picture, this time will be small, and the charm quark will equilibrate with the QGP and participate in the collective motion, with the reverse being true in the weak coupling picture. In both lim- its, the dimensionless quantity $2\pi T D_s$ is proportional to $\eta/s$, and the proportionality factor depends on whether the system is strongly or weakly coupled as follows [838]:

$$ 
\begin{align*}
2\pi T D_s \\
&= 4\pi \times \eta/s = 1 \quad \text{Infinite Coupling at all} \\
T \text{ and Strong Coupling at } T_{pc} \\
&\approx 1.5 \times 4\pi \times \eta/s \quad \text{Strong Coupling at } T = 400 \text{ MeV} \\
&\approx 2.5 \times 4\pi \times \eta/s \quad \text{Weak Coupling at all } T.
\end{align*} $$

The determination of the heavy-quark spatial diffusion coefficient through model descriptions of experimental mea- surements is a major goal of recent dedicated initiatives [432, 433]. In Fig. 74, we show the temperature dependence of $2\pi T D_s$ from models shown in Sect. 2.4. The model calculations describe the D-meson $v_2$ and $R_{AA}$ measure- ments. As discussed in Sect. 2.4, the models differ regard- ing the treatment of the charm-medium interactions. Some assume a quasi-particle description of QGP constituents e.g.

<footnote n="11">¹¹ The maximum a posterior values were tested with ALICE measure- ments for the $T_{RENTo} + VISHNU$ throughout this chapter, and are shown in Fig. 73.</footnote>
PHSD [581], while others use a hydrodynamic treatment e.g. TAMU [492]. The Catania model [583,584] uses Boltzmann equations for the treatment of charm transport for the predictions shown, while LGR [587] is Langevin-based, and some e.g. MC@sHG+EPOS2 [582] include radiative interactions (all include collisional). Other features such as the hadronisation scheme used play an important role. Figure 74 also shows a leading-order pQCD prediction. As for $\eta/s$, $2\pi T D_s$ is inversely related to the coupling strength in the weakly-coupled limit. The pQCD calculation is consistent with the weak coupling limit regarding $\eta/s$, given the relations shown at the start of this subsection. All of the curves are below the pQCD charm diffusion limit, especially around the deconfinement temperature. Nonetheless, as with $\eta/s$, next-to-leading-order expansions can be large [432], and there are other leading-order calculations [593] that give higher values at high couplings i.e. $2\pi T D_s \sim 30$ at $\alpha_S = 0.4$. An alternative approach [838] (not shown), that attempts to include higher order corrections using the non-perturbative equations with a weak potential (T-matrix), yields predictions at least two times higher than any of curves shown. Therefore, a number of examples of weak-coupling approaches are inconsistent with the $2\pi T D_s$ values constrained by our data.

On the right of Fig. 74, the strong-coupling limit is shown, based on the range of $\eta/s$ values constrained in the soft sector, observed in the previous subsection. The infinitely-coupled AdS/CFT limit is also shown. All of the models are close to the strong-coupling limit around $T_{pc}$. At higher temperatures, the curves start to deviate, with LIDO [839] and MC@sHQ+EPOS2 being the closest to this limit, while the other curves lie significantly above. This observation is consistent with two scenarios: the low-momentum charm quarks either strongly couple with the QGP at all temperatures, or they couple strongly with the QGP at lower temperatures only. The aforementioned model differences play a critical role in addressing these scenarios. The LIDO and MC@sHQ+EPOS2 models include radiative interactions; these are less effective in driving the charm quark to participate in the collective motion of the medium, and therefore require smaller values of $D_s$ (therefore smaller relaxation times) to describe the data. On the other hand, LGR also includes radiative interactions, however unlike MC@sHQ+EPOS2, which uses the Boltzmann equations, it is Langevin-based. The difference between both approaches has been investigated within the Catania model, which led to the finding that the Langevin equations lead to higher values of $D_s$. In addition, using the LIDO model, a Bayesian posterior estimation was used to determine the $2\pi T D_s$ within that scheme [839]. The posterior range at 90% is shown on the right of Fig. 74 at $T = 0.4$ GeV, and is within the range of all of curves at the corresponding temperature. It is also important to point out that the hadronisation mechanisms of charm quarks are one of the main sources of uncertainty

in the estimated values of $D_s$; all models shown in Fig. 74 include a contribution of charm hadronisation from coalescence, but the present data only mildly constrain the hadronisation dynamics [139]. Providing more constraining hadronisation measurements is a major goal of future heavy-flavour programmes.

Finally, the charm relaxation times will be discussed. In Sect. 2.4, these were reported to be in the range of $\tau_{charm} = 3-9$ fm/c at the pseudo-critical temperature of $T_{pc} \approx 0.155$ GeV. At higher temperatures (around the start of the QGP phase), the models in Fig. 74 give $\tau_{charm} = 1.3-2.4$ fm/c at $T = 0.4$ GeV. These times are significantly larger than the shear and bulk relaxation times, which is expected due to the large mass of the charm quark. Nonetheless, despite the uncertainties regarding when low-momentum charm quarks strongly couple to the QGP, all models imply these quarks will readily participate in the collective motion of the QGP after their production in the initial stages.

*Jet transport coefficient* $\hat{q}$. High-energy partons that propagate through the QGP interact with the colour charges in the plasma, leading to changes in both the magnitude and direction of their momenta via elastic interactions, as well as the stimulated radiation of gluons. The resulting modifications of the parton shower and the emerging jets are referred to as jet quenching, and are discussed in Sect. 2.4. The dependence of the gluon-radiation probability and energy distribution on the QGP properties can be encoded in the transport coefficient $\hat{q}$, which is the average of the square of the transverse momentum exchanged with the QGP per unit mean free path (see Eq. 16). The total amount of radiative energy loss ($\Delta E$) of a parton passing through the medium is approximately proportional to $\hat{q}$, see also Sect. 2.4.

The dependence of $\hat{q}$ on the temperature of the QGP and the energy of the propagating parton has been explored in both the weak-coupling regime, using Hard Thermal Loop field theory [840–842] and the strong-coupling regime, using the AdS/CFT correspondence [843]. In both regimes, $\hat{q}$ depends linearly on the temperature cubed ($T^3$) and the jet–medium coupling strength. As a result, the dimensionless ratio $\hat{q}/T^3$ increases with coupling strength and does not saturate to a limiting value at large coupling, unlike $\eta/s$ and $2\pi T D_s$. Specific relations between $\eta/s$ and $\hat{q}/T^3$ have been explored at leading and next-to-leading order for weak- and strong-coupling scenarios [841,842]. Given the soft nature of the interactions in the QGP at the temperatures that are experimentally explored, higher-order contributions, large logarithms and non-perturbative effects may have significant quantitative impact on the energy loss. A full mapping of these effects is actively pursued by the theory community. For example higher-order contributions to $\hat{q}$ have been explored [844], and large logarithms and effects of overlapping formation time are being explored in more
<figure><img src="image_88.png" /><figcaption>**Fig. 75** The temperature dependence of $\hat{q}/T^3$ for gluons at an energy of 20 GeV in the QGP phase from various models that describe the ALICE data (see text for details). The range on the right of the plot represents a 90% posterior interval from the Bayesian analysis at $T = 0.3$ GeV</figcaption></figure>

detail [845,846]. In addition, approaches are being developed to take non-perturbative effects into account by calculating $\hat{q}$ using lattice QCD [847–849].

Figure 75 gives an overview of the temperature dependence of $\hat{q}/T^3$ for gluons with an energy of $E = 20$ GeV, that are used in the radiative+collisional energy loss models, which are compared to the ALICE data in Sect. 2.4. Most of these models describe the $R_{AA}$ and $v_n$ for light and charmed hadrons at high-$p_T$. Only CUJET 2.0 underestimates $v_2$ at high-$p_T$ systematically. In these model calculations, $\hat{q}(T, E)$ is used as input, together with the time-dependent medium density profiles from hydrodynamical calculations. The observables $R_{AA}$ and $v_2$ are calculated by integrating over underlying variables, such as the location of the hard scattering in space, the initial parton momentum, the parton trajectories, as well as intrinsic fluctuations in the energy loss and the fragment distributions. While the different curves in Fig. 75 are relatively close at high temperature, a large spread in $\hat{q}/T^3$ is seen at temperatures closer to the phase transition temperature, indicating that the transport coefficient is not well constrained by current data in that range, possibly because the value of $\hat{q}$ is smaller and the contribution to the total energy loss is small. More systematic comparisons of different measurements with the model calculations using advanced statistical techniques such as Gaussian emulators and Bayesian parameter estimation should allow to constrain the temperature dependence of $\hat{q}$ in more detail and are being pursued, for example by the JETSCAPE collaboration [562]. The result of a first determination of $\hat{q}$ by the JETSCAPE collaboration using $R_{AA}$ measurements at the LHC and RHIC via the LBT model is shown in Fig. 75, with the 90% posterior range at a temperature $T = 0.3$ GeV shown on the right. The curve itself is shown for median

parameter values in that analysis. The earlier results from the JET collaboration [559] are also shown as data points at the average temperatures from RHIC (left) and the LHC (right).

Most of the models for which $\hat{q}/T^3$ is shown in Fig. 75 are based on weak-coupling approaches, except for CUJET 3.1 [604] that incorporates the effect of colour-magnetic monopoles. These monopoles can be interpreted as a strong-coupling mechanism, that only affects $\hat{q}/T^3$ in the temperature range near the phase transition. While the data are consistent with CUJET 3.1, they also agree with the DREENA-A model [603], which uses a purely weak-coupling approach including chromo-magnetic fields, and the JETSCAPE tune based on the MATTER and/or LBT models [562]. The good description of the data by these weak-coupling models and the agreement between the $\hat{q}$ values at high temperature, suggest that strong-coupling scenarios which imply a much larger value of $\hat{q}$ at high-$T$ would lead to a lower $R_{AA}$ (more suppression) and/or higher $v_2$ than seen in the measurements.

To summarise, it can be concluded that there are some indications that the QGP is weakly coupled when probed with high-energy jets at the highest temperatures accessible at the LHC, as all of the models that describe ALICE data are dominated by weakly-coupled interactions in that limit. However, it is an open question whether the same applies around the deconfinement temperature. To further constrain the behaviour of the medium as seen by hard probes in this regime, additional exploration of the temperature dependence of the medium properties is needed, both with theoretical methods and with additional measurements and model comparisons, for example as a function of collision centrality and/or collision energy. At the same time, measurements of the modification of jet structure provide further insight into the microscopic dynamics of parton energy loss.

### 3 High-density QCD effects in proton–proton and proton–nucleus collisions

Our understanding of the conditions necessary for QGP formation was overthrown in 2010 with the arrival of the first LHC data. Until then, small systems such as proton–proton (pp) and proton–lead (p–Pb) collisions were considered as plain references, systems in which the conditions to form the QGP could not be reached. At LHC energies, in a small fraction of these collisions, the number of produced particles is of the same order as that in peripheral nucleus–nucleus collisions. Surprisingly, emblematic signatures of QGP formation, such as a double ridge extended in pseudorapidity [850–853] and strangeness enhancement [508], were observed. This opens two main questions: could the QGP be formed in small systems? What are the mechanisms involved in the initial stages of the collision that might create a sufficiently-high
energy density for the potential phase transition of nuclear matter?

These fundamental issues have received great attention and their investigation can be considered among the main novel aspects of the LHC physics programme, with the related discussions extending outside the community involved in QGP studies. In this chapter, we will describe and discuss the main ALICE contributions to this exciting endeavour, also in relation to complementary results obtained by the other LHC experiments.

Section 3.1 addresses the event classification, a fundamental prerequisite that allows direct comparison of the studied phenomena independently of the considered system, be it (high-multiplicity) pp collisions, p–Pb or Pb–Pb. In Sect. 3.2 a broad discussion of the results on particle production across collision systems, as a function of charged-hadron multiplicity, is reported. Special emphasis is given to the role of strange particles, but the discussion is extended also to results on charmed hadrons relevant in this context. Quantities such as the baryon-to-meson ratios, known to be sensitive to the radial boost in the hot partonic system, are also considered in the presentation of the ALICE results. One of the main avenues for the investigation of collective effects is azimuthal anisotropy studies; they are discussed extensively in Sect. 3.3, in terms of the multiplicity dependence of various flow orders ($n$) calculated with $m$-particle cumulants, $v_n\{m\}$, for different collision systems. Elliptic flow measurements for particles, ranging from inclusive to strange hadrons to $J/\psi$, and showing indications of a mass-ordering effect, are then presented. The last two sections of this chapter report studies of the production of hard probes in p–Pb collisions and the search for effects that may be related to the formation of the QGP. Section 3.4 addresses the quarkonium production and, in particular, the investigation of possible final-state suppression of charmonium- and bottomonium-excited states, which has been established for the $\psi(2S)$ loosely-bound state. Finally, in Sect. 3.5, measurements of observables related to jet quenching are reported, via yield measurements of inclusive hadrons or jets and via coincidence of charged-particle jets recoiling from a high-transverse-momentum trigger hadron.

### 3.1 Event classification for small collision systems

A big success in the exploration of large systems reported in the previous chapter has been the development of the Glauber picture (Sect. 2.1.1), which creates a map between the experimentally-measured multiplicity associated with centrality and key physics quantities such as $N_{part}$, $N_{coll}$ and the initial anisotropies. This approach captures essential features that have simplified the interpretation of several physics measurements in heavy-ion collisions. This simple mapping requires that nuclear geometry dominates over subnucleonic

structure and that fluctuations average out. In small systems, neither of these requirements are fulfilled and, therefore, one becomes very sensitive to the modelling of nucleon-nucleon collisions. Charged-particle multiplicity appears to be the simplest, even if imperfect, event classifier. It is a well defined measurement and allows for a direct comparison among all systems without any model dependence. Multiplicity was largely utilised by the ALICE Collaboration for event classification in small systems and we will discuss our current understanding and the associated biases.

In the case of pp collisions, it is generally accepted in models that initial scatterings give rise to some number of subnucleonic interactions that produce the observed hadrons.<fnref n="12" /> These subnucleonic interactions vary from hard partonic interactions, giving rise to high-$p_T$ jets that produce many particles and to colour exchanges involving soft interactions, which produce mainly longitudinal colour fields that decay into a few hadrons. The number of initial interactions, their hardness, and stochastic hadronisation effects mean that a specific event multiplicity can result from several different combinations of processes. Furthermore, combinations may vary as a function of charged-particle multiplicity leading to non-trivial modifications of measured observables as a function of multiplicity, independently of final state effects.

As the QGP is a state of matter dominated by soft interactions, the event-classification scheme should be sensitive to the number of initial interactions. The dynamics involved in these interactions can be experimentally quantified by comparing the shape of the $p_T$ spectra of charged particles, in particular at high-$p_T$, when performing a certain event selection and when measuring instead all inelastic pp and p–Pb collisions. In the following we focus first on pp collisions before moving on to p–Pb collisions.

Figure 76 (top) shows the ratio of the multiplicity-dependent $p_T$ spectra of charged particles in the midrapidity region to the minimum bias (MB) spectrum. The spectra are obtained by slicing according to the multiplicity at mid (black) and forward (red) rapidity. In this selection, the midrapidity multiplicity estimator refers to particle trajectories within $|\eta| < 0.8$ and forward rapidity refers to the multiplicity in the acceptance of the VOC and V0A scintillators, which is of $-3.7 < \eta < -1.7$ and $2.8 < \eta < 5.1$, respectively. The multiplicity classes are expressed in terms of percentiles of the total cross section of inelastic pp collision with at least one charged particle in $|\eta| < 1$ (INEL > 0). The high-multiplicity ratio for the midrapidity estimator is significantly higher than for the forward estimator. One expects a trivial difference as the $p_T$ spectra are being measured at midrapidity in the same kinematic region where the midrapidity multiplicity selection is done. However, the slope of the

<footnote n="12">¹² These interactions can be between partons, dense colour fields in the CGC, or similar.</footnote>
<figure><img src="image_89.png" /><figcaption>**Fig. 76** Ratios of multiplicity-selected transverse-momentum distributions of charged particles for pp collisions (top) [854] and p–Pb collisions (bottom) [855] to inelastic (INEL) pp collisions, shown in logarithmic scale. For pp collisions, the forward-rapidity event selection is done using the V0M multiplicity estimator ($-3.7 < \eta < -1.7, 2.8 < \eta < 5.1$), while the midrapidity event selection is done using the SPD tracklets in $|\eta| < 0.8$. For p–Pb collisions, the forward-rapidity event selection is done using the A-side of the V0 multiplicity estimator ($2.8 < \eta < 5.1$), corresponding to the Pb-going side, while the midrapidity event selection is done using the first layer of the SPD. The horizontal dashed lines denote the integrated charged-particle multiplicity density ratios as a reference. For this figure, the pp high-multiplicity event classes correspond to 0–1% (V0M-based, forward-rapidity selection) and 0.009–0.088% (SPD-based, midrapidity selection) and the low-multiplicity pp event classes correspond to 70–100% (V0M) and 49.5–100.0%, where percentiles are defined with respect to the total pp INEL > 0 cross section. For all event selections, the p–Pb high-multiplicity and low-multiplicity event classes correspond to 0–5% and 70–100%, respectively, with percentiles defined with respect to the total visible V0A cross section</figcaption></figure>

$p_T$ spectra at high-$p_T$ indicates that the midrapidity estimator selects harder and harder subnucleonic interactions as the multiplicity increases. The ratios obtained with the forward estimator do not show a change in slope at high-$p_T$. Still, the hard high-$p_T$ production is more enhanced than soft low-$p_T$ production in high-multiplicity collisions and vice versa in low-multiplicity collisions. This implies that the scaling with multiplicity of soft and hard processes is fundamentally different in pp compared to nucleus-nucleus collisions.

To gain insight into the origin of the difference between the spectra ratios for the midrapidity and forward estimator, it is interesting to use an event generator in which the initial number of Multiple Parton Interactions (MPIs) utilised in the simulation can be directly related to the charged-particle multiplicity. Figure 77 reports a study with PYTHIA 8 (Monash

2013 tune) which shows that high-multiplicity events are preferably related to large number of MPIs, but are also influenced by other effects, such as the fragmentation of partons into a large number of final state particles. The PYTHIA 8 simulations show that the forward multiplicity estimator has the strongest correlation between the number of MPIs and the multiplicity. For this reason, the forward multiplicity slicing is used for multiplicity selection in the rest of this section unless specifically noted otherwise. As the multiplicity selection is done on charged particles, a second advantage of the forward selection is that it does not create an imbalance between charged and neutral particles at midrapidity, as described in detail in [509].

As discussed in Sect. 2.1.1, and contrary to the MPI description of pp collisions, nucleus-nucleus collisions are traditionally conceived as multiple interactions of individual nucleons. The p–Pb collisions act as a bridge between the two limiting descriptions where both views would need to be reconciled.

In the same way as done for the pp collisions, p–Pb collisions are classified by making use of the V0A detector signal from the pseudorapidity range of $2.8 < \eta < 5.1$, corresponding to the Pb-going side (direction of the fragmentation of the Pb nucleus), where the larger signal amplitudes associated to the nucleus fragmentation region will stabilise short-range multiplicity fluctuations. Figure 76 (bottom) shows the ratio of the $p_T$ spectra obtained in p–Pb, by slicing in multiplicity at mid (black) and forward (red) rapidity, to the MB spectrum obtained in pp collisions at the same $\sqrt{s}$. The midrapidity-estimator ratio has a higher multiplicity dependence than the forward-estimator ratio, although the effect is reduced with respect to the situation in pp collisions, shown in the upper panel. Focusing on the low-multiplicity sample, the ratio is below one in particular for the midrapidity estimator. This poses a fundamental challenge to the Glauber model as one cannot have $N_{\text{part}}$ and $N_{\text{coll}}$ smaller than that of a single pp collision. It is also clear that for the midrapidity estimator, there is no easy fix for this as the scaling for soft and hard processes are different. This indicates that one has to be careful when applying the Glauber model to interpret the results obtained with multiplicity slicing in small collision systems where the effect of particle-production fluctuations for individual nucleon-nucleon collisions, as shown in Fig. 76 (top), can play a big role. Also in peripheral Pb–Pb collisions, the centrality selection based on charged-particle multiplicity, leads to a bias in the determination of the collision geometry parameters [173], which has to be taken into account when interpreting measurements that use these parameters, like the nuclear modification factor (see Sect. 2.4.2).

To circumvent these limitations, the ALICE Collaboration has developed an event-selection method for p–Pb collisions based on the energy deposited from the Pb nuclei in the Zero Degree Calorimeters (ZDCs) [855]. Such an approach has the
<figure><img src="image_90.png" /><figcaption>**Fig. 77** Mean number of Multiple Parton Interactions, $\langle N_{MPI} \rangle$, as a function of mean charged particle multiplicity at midrapidity, $\langle dN_{ch}/d\eta \rangle|_{\eta|<0.5}$, for different event selections, calculated with the PYTHIA 8.2 model (Monash 2013 tune) [856]</figcaption></figure>

ability to classify event activity with an even larger rapidity gap between actual measurement (acceptance of the ALICE central barrel or muon spectrometer) and event selection and a very low sensitivity to fluctuations in particle production introduced by subnucleonic interactions. This method, coupled to an assumption that the total charged-particle multiplicity, dominated by low-$p_T$ particles, scales with $N_{part}$ – a combination henceforth known as the ‘hybrid’ method – has been shown to recover the expected scaling of particle yields at high momenta. This can be observed when calculating the $Q_{pPb}$, a quantity analogous to the nuclear modification factor $R_{AA}$, but calculated in proton–nucleus collisions. While the $Q_{pPb}$ at high momenta deviates from unity in data using typical multiplicity selections, such behaviour is well described by Q-PYTHIA, a model that can be used to do a baseline check including multiplicity fluctuations. When performing ZDC-based event selections and utilising the hybrid method [855], on the other hand, the high-$p_T$ $Q_{pPb}$ values are observed to be compatible with unity. These results are shown in Fig. 78 and the selection based on the ZN signals (neutron ZDCs) will serve as a baseline for the p–Pb studies performed in Sect. 3.5.

Multiplicity was extensively used as an event classifier during the LHC Run 2, although it introduces some difficulties when comparing data and theory. In fact, multiplicity cannot be directly related to the initial energy density of the collision, and the same multiplicity in small and large col-

lision systems can result from completely different initial energy densities.

While the event classification studies performed during Runs 1 and 2 were highly successful, more advanced selection strategies are being developed which will allow for a better determination of the phase space in which several phenomena appear in small collision systems. Notably, the transverse activity estimator $R_T$ [857] and spherocity selections [854] may shed further light on the transition between low and high multiplicities, with some measurements already underway and many more to appear in the near future. In the next sections, we will discuss the results based on the charged-particle multiplicity event classification scheme, since most of these are published at this time.

## 3.2 Dynamics and hadrochemistry of particle production

The ALICE Collaboration studied the production of a large variety of particle species across a broad range of collision energies ($2.76 < \sqrt{s_{NN}} < 13$ TeV) for various hadronic systems (pp, p–Pb, Pb–Pb and Xe–Xe) as a function of the charged-particle multiplicity of the collision at midrapidity in event classes defined by selections on the signal in the forward detectors V0A and V0C ($2.8 < \eta < 5.1$ and $-3.7 < \eta < -1.7$). These results are shown in Figs. 79 and 80.

The first striking observation is that the particle yield ratios with respect to pions increase with multiplicity in pp and p–Pb collisions, while they remain almost unchanged in heavy-ion collisions, as can be seen in Fig. 79. The only exceptions are for resonances such as the $K^{*0}(892)$ and $\Lambda(1520)$ (see Sect. 2.3.3). This evolution is smooth and very similar across collision systems, with measured ratios being consistent within uncertainties for event classes with similar multiplicities in pp, p–Pb, Xe–Xe and Pb–Pb collisions. Moreover, the ratios do not depend on centre-of-mass energy. These are very important observations, which indicate that hadron yields are mostly related to final-state charged-particle multiplicity density rather than collision system or beam energy.

The relative production rates of particles containing strange quarks, measured using the $\Lambda/\pi$, $\Xi/\pi$ and $\Omega/\pi$ ratios and shown in the left panel of Fig. 79, increase faster with multiplicity than those containing u and d quarks only. This effect is commonly referred to as “strangeness enhancement” and it is proportional to the strangeness content of the hadron: largest for triply-strange $\Omega^-$ and progressively less evident for particles with two and one strange quark, being absent for protons as discussed also in Sect. 2.3.1. The $\phi$ meson, a hadron with hidden strangeness $s\bar{s}$, exhibits a relative yield increase that is mid-way between the ones measured for $\Xi^-$ and $\Lambda$. The experimental results concerning strangeness implied not only a major paradigm shift for the heavy-ion community, who had to re-evaluate the assumption that pp was a valid reference system, but were seen with great
interest also by the particle physics community at large. As a consequence, theorists working on QCD-inspired models now consider the description of strongly-interacting matter an integral part of their programme.

A general interpretation of the enhancement could be that the production of (multi-)strange baryons requires sufficiently energetic processes that are not favoured in low-multiplicity collisions. However, this mechanism is solely dependent on available energy at hadronisation, and therefore it should be captured by e.g. the most commonly used QCD-inspired pp event generator PYTHIA. Moreover, PYTHIA itself predicts a different mixture of event processes at low-multiplicity – e.g. diffractive events – and in fact includes no dynamical component able to reproduce the steady increase of particle ratios such as $\Xi/\pi$ and $\Omega/\pi$ as a function of charged-particle multiplicity in pp collisions.

In PYTHIA, hadron production occurs via the incoherent break-up of colour flux tubes called ‘strings’, which exhibit constant energy density, leading to the conclusion that even high-multiplicity events would result in unchanged particle ratios. As a consequence of the observation in [508], PYTHIA modelling had to resort to conceptually different physical mechanisms to reproduce experimental data, such as the inclusion of ‘colour ropes’ – colour flux tubes with increased tension that are created whenever several strings overlap prior to hadronisation in high-multiplicity pp collisions [866]. The predictions from PYTHIA with colour ropes can be seen in Fig. 81 and describe strangeness enhancement in pp collisions within a 10% accuracy in high-multiplicity collisions. However, it is important to note that

the proton-to-pion ratio is not correctly described in this model, which indicates that further theoretical studies are still required for a proper description of hadrochemistry using QCD-inspired, density-dependent modelling. Furthermore, alternative approaches that select on specific strangeness-related string breakup processes in PYTHIA have also been shown to reproduce the correct relative strangeness production increase [867], though a complete dynamical description of hadrochemistry employing this particular mechanism is still to be developed.

In addition, statistical hadronisation models also provide a good description of the relative increase of strangeness production as a function of charged-particle multiplicity density, as can be seen in Fig. 81. A canonical statistical model, incorporating exact conservation of baryon number, electric charge and strangeness, and varying $\gamma_s$ and $T_{\text{chem}}$ using the Thermal FIST package [398,868] is able to capture the evolution of most particle ratios quite well qualitatively across systems, with a maximum of 20% deviation seen to be approximately multiplicity-independent. An alternative canonical calculation restricted to pp and using a fixed $T_{\text{chem}} = 156$ MeV based on [418] is also able to capture the strangeness increase rather well, though also in that case, discrepancies of up to 30% are observed in other particle ratios, indicating that a comprehensive description of all particle ratios within experimental uncertainties remains a theoretical challenge. Further developments in proton-proton yield modelling will also include the extension of the GSI-Heidelberg model employing the S-matrix formalism [434] to pp and are expected to appear in the near future.

<figure><img src="image_91.png" /><figcaption>**Fig. 78** (Left) $Q_{pPb}$ at high transverse momenta calculated using various multiplicity selections (points) and predicted with the Q-PYTHIA model (curves) [855]. (Right) $Q_{pPb}$ calculated with the hybrid method.</figcaption></figure>

<figure><img src="image_92.png" /><figcaption>Spectra are calculated in event classes based on the ZNA energy (neutron ZDC on the A side, Pb-going) and $N_{\text{coll}}$ is obtained assuming $N_{\text{part}}$ scaling at midrapidity. Figure from [855]</figcaption></figure>
<figure><img src="image_93.png" /><figcaption>**Fig. 79** (Left) $p_T$-integrated yield ratios to pions ($\pi^+ + \pi^-$) and $p_T$-integrated yield ratios between resonance and corresponding ground state as a function of $\langle dN_{ch}/d\eta \rangle$ measured in $|\eta| < 0.5$ in pp collisions at $\sqrt{s} = 7$ and 13 TeV [508–511], p–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV [506, 858, 859] and at $\sqrt{s_{NN}} = 8.16$ TeV [860], Xe–Xe collisions at $\sqrt{s_{NN}} = 5.44$ TeV [281] and Pb–Pb collisions at $\sqrt{s_{NN}} =$</figcaption></figure>

2.76 TeV [386, 504, 507] and at $\sqrt{s_{NN}} = 5.02$ TeV [273, 402, 455]. (Right) Ratios involving $\rho$ mesons, $K^*$, $\Sigma^{*\pm}$, $\Lambda(1520)$ and $\Xi^*$. All yields are obtained at midrapidity, i.e. $|y| < 0.5$. The error bars show the statistical uncertainty, whereas the empty and dark-shaded boxes show the total systematic uncertainty and the contribution uncorrelated across multiplicity bins, respectively

As seen in the right panel of Fig. 79, the multiplicity evolution of particle ratios involving resonances and their non-resonant hadronic states, such as $K^{*0}(892)$ and $K^{\pm}$, were also measured by the ALICE Collaboration. As discussed in Sect. 2.3.3, the scope of this study in heavy-ion collisions is to determine the presence and lifetime of the hadronic phase which follows chemical freeze-out in heavy-ion collisions. The same study was performed in small collision systems (pp and p–Pb), showing again a rather smooth trend across multiplicity. Given the statistical and systematic uncertainties in the measurements, there is no evidence of a significant multiplicity evolution for long-lived resonances, while for the $K^{*0}$ a hint of a decreasing trend is observed, pointing to a hadronic stage whose density and energy allow the rescattering of the decay products. The measurement of the $\rho^0$ meson still carries uncertainties that prevent any statement regarding the evolution of the $\rho^0/\pi$ ratio.

The yield of light nuclei also evolves with multiplicity, as shown in the left panel of Fig. 80, increasing smoothly across different systems and energies and reaching the values observed in heavy-ion collisions for the highest multiplici-

ties. The increase is more pronounced for $^3$He nuclei, which contain three nucleons, than for deuterons.

Finally, heavy-flavour yields were also measured from proton–proton to central Pb–Pb collisions, and the corresponding yield ratios can be seen in Fig. 80 (right), which shows clear evidence for an increase in the $(J/\psi)/\pi$ ratio in high-multiplicity proton–proton collisions. However, contrary to what is observed for particles containing only light-flavour quarks, the evolution of the $(J/\psi)/\pi$ is not continuous with charged-particle multiplicity density. This is likely due to the fact that charm quarks are produced in hard scattering processes with a rate that depends on collision system and centre-of-mass energy. While the factorisation approach can be used to calculate single inclusive production in minimum-bias interactions, the multiplicity dependence of charmed-particle production in small collision systems remains a topic in which further theoretical developments are needed for a better understanding of the role of MPIs in the building of the hadronic multiplicity and the interplay between the soft and hard components of the event. More precise measurements in Run 3 will be needed to explore potential dissociation of $J/\psi$ ranging from ground states to excited states through
<figure><img src="image_94.png" /><figcaption>**Fig. 80** $p_T$-integrated yield ratios to protons or charged pions (particles+antiparticles) as a function of $\langle dN_{ch}/d\eta \rangle$ measured in $|\eta| < 0.5$ in pp collisions [511,861,862], p-Pb collisions [858,863,864] and Pb-Pb collisions [387,402,435,865] for nuclei (left) and for $D^0$ and $J/\psi$</figcaption></figure>

mesons (right). All yields are obtained at midrapidity, i.e. $|y| < 0.5$. The error bars show the statistical uncertainty, whereas the empty and dark-shaded boxes show the total systematic uncertainty and the contribution uncorrelated across multiplicity bins, respectively

ratios of relative yields as a function of the charged-particle multiplicity density.

In heavy-ion collisions, the transverse momentum spectra of identified particles are sensitive to the radial flow imparted by the system expansion. This translates into an evolution of the spectra towards higher $\langle p_T \rangle$ for more central A–A collisions, as discussed in Sect. 2.2.2. The same effect has been observed when measuring $p_T$ spectra as a function of the charged-particle multiplicity in pp and p–Pb interactions. The low-momentum part of the spectra (up to around 2.5 GeV/c) becomes progressively more flat as the multiplicity increases, leading to an increase of the $\langle p_T \rangle$ for all particle species. This opens the intriguing question: is the origin of the hardening of the spectra the same in different colliding systems? Is it the QGP in high-multiplicity pp interactions causing the increase in $\langle p_T \rangle$? Can this effect be solely due to radial flow in an expanding system?

To try to answer these questions, the ratio between $\Lambda$ and $K_S^0 p_T$ spectra was measured in pp interactions characterised by different final-state multiplicities [869]. The idea behind this observable is that the radial boost of a collectively-expanding system should impact the heavier hadrons more strongly, giving rise to the observation of an enhanced

baryon-to-meson ratio at intermediate-$p_T$. This enhancement is observed in the $\Lambda/K_S^0$ ratio for small collision systems in a qualitatively very similar way as in heavy-ion interactions. The magnitude of the intermediate-$p_T$ enhancement increases as the multiplicity increases and the peak position moves towards higher values for high-multiplicity collisions, in agreement with the hydrodynamic origin picture. Moreover, the increase at intermediate momenta is accompanied by a corresponding depletion of the ratio at low-momenta, with the integrated $\Lambda/K_S^0$ ratio exhibiting essentially no multiplicity dependence in pp collisions. This observation also holds for proton-to-pion ratios and is qualitatively reminiscent of what is observed in Pb-Pb collisions, as described in Sect. 2.3.2. In order to study the low-$p_T$ depletion and mid-$p_T$ enhancement more quantitatively, specific $p_T$ regions have been selected and the multiplicity dependence of the $\Lambda/K_S^0$ ratio in these intervals is reported in the right panel of Fig. 82.

Similar studies have also been carried out recently for the charmed baryon-to-meson ratio $\Lambda_c/D^0$, where similar behaviour has been observed [862]. The corresponding results are shown in Fig. 82 together with the $\Lambda/K_S^0$ ratio and exhibit a pattern consistent with the non-charmed baryon-
<figure><img src="image_95.png" /><figcaption>**Fig. 81** Selected measurements from the previous figures compared to different theoretical models. (Top left) $p_T$-integrated ratio of $p/\pi$ as a function of $\langle dN_{ch}/d\eta \rangle$ at midrapidity. (Top right) $p_T$-integrated ratio of $\phi/\pi$ as a function of $\langle dN_{ch}/d\eta \rangle$ at midrapidity. (Bottom left) $p_T$-</figcaption></figure>

integrated ratio of $\Omega/\pi$ as a function of $\langle dN_{ch}/d\eta \rangle$ at midrapidity compared to PYTHIA 8.2 (Monash 2013 tune) [856] and PYTHIA colour ropes model [866]. (Bottom right) $p_T$-integrated ratio of $(J/\psi)/\pi$ as a function of $\langle dN_{ch}/d\eta \rangle$ at midrapidity

<figure><img src="image_96.png" /><figcaption>**Fig. 82** (Left) Transverse-momentum dependence of $\Lambda/K_S^0$ at midrapidity for high- and low-multiplicity proton–proton collisions at 13 TeV [869] drawn together with $\Lambda_c/D^0$ ratios [862]. (Right) Mul-</figcaption></figure>

tiplicity dependence of the $\Lambda/K_S^0$ and $\Lambda_c/D^0$ ratios at midrapidity in selected transverse-momentum intervals in pp collisions at 13 TeV
to-meson ratio, although future measurements with smaller uncertainties are still needed to come to a firm conclusion.

Overall, all results discussed here point to the remarkable fact that the hadrochemistry depends solely on the final-state charged-particle multiplicity and is independent of the colliding system and collision energy at the LHC. This observation, together with the specific phenomena that appear to vary smoothly with charged-particle multiplicity, such as strangeness enhancement, constitutes a breakthrough result of the ALICE experiment.

Further understanding regarding particle production rates of specific hadron species is expected to come in the future with studies in which these particles are analysed in association with jets and using other correlation techniques such as net-quantum-number fluctuation measurements.

### 3.3 Collective effects: anisotropic flow

Particles emitted from a collective system created in a collision would exhibit global correlations shared among many particles spanning across a large rapidity range, resulting from the response to an initial anisotropy. This phenomenon can be studied with measurements of azimuthal correlations. A clear illustration can be seen in di-hadron correlations as a function of the difference in azimuthal angle, $\Delta\varphi$, and pseudorapidity, $\Delta\eta$, of the particle pairs, in the form of a double-ridge structure (a cosine modulation in $\Delta\varphi$) extended over a long range in $\Delta\eta$. In particular, the ridge signal in the near-side ($|\Delta\varphi| \approx 0$) is understood as manifestation of a collectively-expanding medium, which is distinctively observed in collisions of heavy ions [321,322].

A significant near-side ridge signal was first reported in high-multiplicity pp collisions by the CMS Collaboration [852]. Shortly after, a double-ridge signal was evidenced in high-multiplicity p–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV by the ALICE Collaboration [851] (Fig. 83) and confirmed by the CMS and ATLAS Collaborations [853,870]. Such observations indicated the possibility that small systems could exhibit a collective behaviour. A similar double-ridge signal, though with smaller magnitude than in p–Pb or Pb–Pb collisions, was also later measured in the high-multiplicity pp collisions [850,871]. Whether the ridge results from multiparticle global correlations having an origin similar to those in heavy-ion collisions will be further discussed in this section using measurements of anisotropic flow $v_n$ obtained from azimuthal particle correlations.

Effects that do not arise from global correlations relative to common symmetry planes, in particular short-range correlations between few particles originating from jets or resonance decays, are called non-flow. This contamination is dominant in low-multiplicity collisions, which makes measurements in small systems more sensitive to non-flow than in heavy-ion collisions. In the case of measurements of two-

<figure><img src="image_97.png" /><figcaption>**Fig. 83** Associated yield per trigger particle in $\Delta\varphi$ and $\Delta\eta$ for pairs of charged particles with $2 < p_{T,\text{trig}} < 4$ GeV/c and $1 < p_{T,\text{assoc}} < 2$ GeV/c in p–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV for the 0–20% multiplicity class, after subtraction of the associated yield obtained in the 60–100% event class under the assumption that the latter multiplicity class is completely dominated by non-flow correlations [851]. A clear double-ridge signal spanning a large range in $\Delta\eta$ can be seen on the near-side $|\Delta\varphi| \approx 0$ and on the away-side $|\Delta\varphi| \approx \pi$, indicative of the possibility of a collective behaviour in p–Pb collisions</figcaption></figure>

particle correlations, non-flow can be reduced by imposing an $\eta$ gap between the correlated particles or between two sections of a phase space called subevents. The contamination can be further suppressed by subtracting correlations measured in low-multiplicity events, representing the non-flow sample. Different methods are used based on the assumption of whether these events contain only short-range correlations, or whether flow is present there too, yielding to slightly different results. The peripheral subtraction method [872] is based on the first assumption and has been widely used in experimental measurements. The results presented in Fig. 83 used the low-multiplicity 60–100% centrality class event sample, selected using the V0A detector, to subtract non-flow from the 0–20% centrality class containing the high-multiplicity p–Pb collisions. The robustness of the subtraction method was studied by projecting the correlations onto the azimuthal axis and comparing the per-trigger yields in low-multiplicity p–Pb collisions and minimum bias pp collisions; no significant difference was found [851]. The latter assumption of an azimuthal correlations even in low-multiplicity events is used in the template fit method developed by the ATLAS Collaboration [850]. Since the main difference between the methods is the way they treat the non-flow in the baseline sample, the resulting flow harmonics differ mainly in the low-multiplicity region, where the measurements of $v_n$ from the template fit method are approximately constant with multiplicity, while the $v_n$ obtained from the peripheral
subtraction method exhibits a decreasing trend with multi-
plicity [873]. In contrast to two-particle correlations, multi-
particle cumulants [874] suppress few-particle correlations
by construction and are thus less prone to non-flow con-
tamination. The remaining non-flow effects in multi-particle
cumulants are expected to be negligible in large systems, but
can remain significant at small multiplicities typical of p–
Pb or pp collisions. Thus, the subevent method consisting
of an η separation between subevents has also been imple-
mented for multi-particle cumulants [875]. Given that more
particles are used to construct the correlation, three or more
subevents may be used for non-flow suppression. While the
two-subevent method largely suppresses short-range correla-
tions e.g. within a jet cone, a signal from correlated particles
between two cones of a dijet separated in η (inter-jet cor-
relations) may still remain. In addition, it is not yet clear
whether this method is affected by possible flow decorrela-
tion effects in the longitudinal direction, which may further
decrease the resulting flow signal when a large η gap is intro-
duced. Therefore, using more subevents further reduces the
non-flow contamination by also suppressing inter-jet corre-
lations and by being less sensitive to flow decorrelations.
Similarly as in the case of the subtraction methods men-
tioned above, the subevent method in multi-particle correla-
tions may also yield differences when comparing results from
narrow (e.g. ALICE) and wide acceptance (e.g. ATLAS and
CMS) measurements due to the varying interplay between
non-flow and possible decorrelation effects. Proper suppres-
sion of non-flow is a crucial element for a correct data-to-
model comparison.

Measurements of $v_n$ in pp and p–Pb collisions as a function of charged-particle multiplicity $N_{ch}(|\eta| < 0.8)$ are shown in Fig. 84a together with measurements from peripheral collisions of Pb nuclei corresponding to centralities higher than 60%. A large pseudorapidity gap $|\Delta\eta| > 1.4$ (1.0) in the case of $v_2$ ($v_3$ and $v_4$) is applied to ensure measurements with maximum possible non-flow suppression. The elliptic flow is found to be compatible between all reported systems at $N_{ch}(|\eta| < 0.8) \le 50$, while at higher multiplicities the rise of $v_2$ for Pb–Pb collisions is in contrast with the weak multiplicity dependence of $v_2$ in pp and p–Pb collisions. On the contrary, the magnitudes of $v_3$ and $v_4$ are similar and without significant dependence on $N_{ch}$ in both small and large systems over the entire measured multiplicity range. These observations support the picture of $v_2$ reflecting the system’s collective response to the initial geometry and higher harmonics originating from its fluctuations, as discussed for heavy-ion collisions in Sect. 2.2. Within this picture, the weak $N_{ch}$ dependence of $v_n$ in small systems shows the dominance of fluctuations, while in Pb–Pb the results of $v_2$ clearly differentiate between the fluctuation-dominated region at low-$N_{ch}$ and the influence of the large geometrical eccentricity of

the nuclear overlap, which is setting in already at $N_{ch} \sim 50$
(roughly corresponding to the centrality interval of 70–80%).

Measurements of *m*-particle cumulants for *m* > 2 can
help to determine whether the measured correlations pre-
sented above are of a collective nature shared among many
particles. The measurement of the four-particle cumulant
$c_n\{4\}$ is of particular interest, since its sign must be nega-
tive to result in a real-valued $v_n$, with $v_n\{4\} = \sqrt[4]{-c_n\{4\}}$.
This condition is fulfilled in large systems [809,877], allow-
ing us to calculate the $v_2\{4\}$, shown in Fig. 84b together
with the same measurements of two-particle cumulants pre-
sented in panel (a). Note that the measurement is reported
with the 3-subevent method in order to be consistent with
the results of small systems shown in the other panels, never-
theless, the method did not yield any significant change with
respect to the default measurement of $v_2\{4\}$ [280]. The rela-
tion $v_2\{2\} > v_2\{4\}$ observed in Fig. 84b is due to the event-
by-event flow fluctuations affecting different orders of cumu-
lants differently, which are discussed in more detail in Sect. 4.
In contrast to heavy-ion collisions, measuring a negative
$c_2\{4\}$ in small systems is not trivial since the dominating non-
flow effects are known to give rise to $c_2\{4\} > 0$ [875,877]. A
negative $c_2\{4\}$ has been observed in high-multiplicity p–Pb
collisions [877]. In pp collisions, this was only possible with
the 3-subevent method that further suppresses dominating
non-flow correlations, and with a specific high-multiplicity
trigger using the V0 detector [280] which acts as a positive
bias towards events with a smaller presence of jets at midra-
pidity (see Sect. 3.1). The resulting $v_2\{4\}_{3-\text{sub}}$ is shown in
Fig. 84c for p–Pb and (d) and (e) for pp collisions, together
with the same measurements of two-particle cumulants pre-
sented in panel (a). The relation $v_2\{2\} > v_2\{4\}$, apparent in
Pb–Pb collisions, is less pronounced in p–Pb collisions, and
the measurements are compatible within uncertainties in pp
collisions. This may indicate different types of flow fluctua-
tions, larger longitudinal decorrelations, or stronger non-flow
effects in the measurements done in small systems. Never-
theless, it should be noted that the region of $N_{ch}$ in which
the measurements presented in Fig. 84 overlap is small. It
should be noted, that a similar comparison was also per-
formed by the ATLAS [873] and CMS [871] Collaborations.
While ATLAS reports that $v_2\{2\}$ is clearly larger than $v_2\{4\}$
in all collision systems at similar multiplicities, this rela-
tion does not hold at low multiplicities in the CMS measure-
ments. Compared to the $v_2\{4\}$ measurement (not shown here,
see Ref. [280]), the subevent method allowed us to obtain a
real-valued $v_2\{4\}_{3-\text{sub}}$ down to lower multiplicities in p–Pb
collisions, suggesting that the presence of collective effects
extends down to regions with multiplicities only few times
larger than that of minimum-bias collisions. In addition, it
has been observed that $v_2\{4\} \approx v_2\{6\}$ within uncertainties in
both p–Pb and pp collisions [280], confirming the collective
nature of small systems.
<figure><img src="image_98.png" /><figcaption>**Fig. 84** Multiplicity dependence of $v_n\{m\}$ in pp, p–Pb and Pb–Pb collisions [280]. Statistical uncertainties are shown as vertical lines and systematic uncertainties as filled boxes. The results of two-particle cumulants $v_n\{2\}$ in all collision systems are shown together in panel (a). The same results together with the four-particle cumulant $v_2\{4\}$ and comparison to models are reported individually in Pb–Pb (b), p–Pb (c) and pp collisions (d). Data are compared with PYTHIA 8.210 (Monash 2013 tune) [856] simulations of pp collisions at $\sqrt{s} = 13$ TeV, and</figcaption></figure>

the same simulations with string shoving mechanism with the string amplitude of the shoving force $g = 10$. Data are further compared to IP-Glasma+MUSIC+UrQMD [48] calculations of pp collisions at $\sqrt{s} = 13$ TeV and p–Pb and Pb–Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV, and iEBE-VISHNU calculations of pp collisions at $\sqrt{s} = 13$ TeV with HIJING initial conditions and parameter set IV [876]. The width of the bands represent the statistical uncertainty of the model

Despite the many similarities with large systems, the origin of the observed collectivity in pp and p–Pb collisions is not clear and further insight can be obtained by comparing the measurements with model calculations. Several approaches are able to explain some features of the experimental data, ranging from purely initial-state gluon momentum correlations [878,879], to final state effects where the particle azimuthal anisotropy arises as the system’s response to initial spatial anisotropy. The latter is described either with a macroscopic model, such as hydrodynamics [876,880,881], or a microscopic approach, such as transport models [882,883], hadronic rescatterings [884,885], or PYTHIA with the string shoving mechanism [886].

Hydrodynamic calculations within the IP-Glasma+MUSIC+UrQMD framework [48] are quantitatively consistent with the $v_n\{m\}$ measurements in Pb–Pb collisions, as shown in Fig. 84b. The same model shows, in panel (c), a good description of the $v_n\{2\}$ and even a qualitative agreement with $v_2\{4\}$ results in p–Pb collisions, which could only be achieved by assuming sub-nucleon fluctuations, i.e. that the proton contains smaller structures modelled as three hot spots [881]. Instead, the data are not reproduced by the model if the pro-

ton is assumed to have a Gaussian profile in the transverse plane. While the contribution from fluctuations at such small scales is relatively mild in heavy-ion collisions where the spatial eccentricity is predominantly defined by the colliding nucleons, they become more important in small collision systems. Despite the successful description of Pb–Pb and p–Pb collisions, the IP-Glasma+MUSIC+UrQMD model [48] does not reproduce the weak dependence of $v_2$ on multiplicity in pp collisions, and it significantly underestimates the $v_3$ (see Fig. 84d). In the same panel the iEBE-VISHNU hydrodynamic model calculations with HIJING initial conditions [876,887] are shown. Only simulations with the parameter set denoted as “IV”, providing the best description of our data, are chosen for comparison. The model tends to follow the weak trend of $v_n$ with multiplicity. While a quantitative agreement is found with measurements of $v_3$ and $v_4$, this was not achieved in the case of $v_2$ with any of the parameter sets used in the simulations. This discrepancy may be due to the fact that the model parameters are tuned to describe the $v_2$ measurements obtained from two-particle correlations with the template fit method used for non-flow subtraction, while our results were performed using two-particle cumu-
lants with a $|\Delta\eta| > 1.4$ gap. Further tuning of the the-
oretical models could improve the comparison, though it
should be noted that non-flow may not be fully removed
when only two-particle correlations are used for the mea-
surement. In addition, the real-valued $v_2\{4\}$ in pp collisions
has not been obtained in any hydrodynamic model so far, pos-
sibly because of the large nonlinear hydrodynamic response
to initial eccentricities causing a positive $c_2\{4\}$ [876,888].

PYTHIA 8 (Monash 2013 tune) model calculations [856, 889], in which there is no thermalised medium, are compared to pp data in Fig. 84e. The calculations cannot fully reproduce all the aspects of the $v_n$ measurements after applying the subevent method and implementing the selection of high-multiplicity events as in data. As shown in Fig. 84e, PYTHIA underestimates the $v_2$ calculated with the two-particle cumulant with an $|\Delta\eta|$ gap. In addition, it should be noted that the PYTHIA results for $v_2\{4\}$ are not shown in the figure due to a positive sign of $c_2\{4\}$, persistent even after the subevent method and the high-multiplicity selection were applied, which had proved to be sufficient in suppressing non-flow effects in data. The higher order harmonics obtained from PYTHIA 8 (Monash 2013 tune) show an opposite trend with respect to the data at high multiplicities, i.e., $v_4 > v_3$. For $N_{ch}(|\eta| < 0.8) \ge 75$ the third harmonic could not be calculated in the model due to corresponding negative two-particle cumulants, while the $v_3$ measured in data exhibits a constant multiplicity dependence. Nevertheless, recent calculations with the additional string shoving mechanism, where overlapping strings repel each other inducing flow-like effects without the presence of a medium, qualitatively reproduce some features of the data. As also shown in Fig. 84e, the shoving mechanism increases the values of $v_2\{2\}$ leading to a better comparison with the data. It further revealed a positive value of $v_3$ at high-multiplicity in contrast to the non-shoving case. Even though the increase of $v_4$ with the shoving mechanism does not agree with the data, the fact that a model without the presence of a medium is able to reproduce some aspects of the measurements shows the potential for an alternative explanation of the collective effects revealed in the smallest collision systems. An interesting finding was shown in [889] where a hint for a negative sign of the four-particle cumulant was revealed after applying the high-multiplicity selection, though further investigations are needed to confirm and fully understand such observations.

Studies of azimuthal anisotropies of identified hadron
species in small systems, in particular their mass depen-
dence, have a large potential to determine the presence of
a partonic medium. In heavy-ion collisions, the mass order-
ing of the anisotropic flow of identified hadrons at low-
$p_T$ is explained as a consequence of a radial expansion of
the medium pushing heavier particles to larger $p_T$, and is
described well by viscous hydrodynamic models (for further

discussion see Sect. 2.2.4). The baryon-meson grouping of $v_n$ at intermediate-$p_T$ is understood as an effect of particle production via quark coalescence, hence pointing to the presence of partonic collectivity in heavy-ion collisions. These effects are thus considered as characteristic flow features, discussed in more detail in Sect. 2.3.2. The analysis of two-particle correlations of identified hadrons in high-multiplicity p–Pb collisions, after the subtraction of correlations from the low-multiplicity data sample, showed the near-side ridge structure [890] as observed previously in charged hadron measurements (see Fig. 83). The Fourier coefficients of this ridge structure, $v_2$, of pions, kaons and protons as a function of $p_T$ shown in Fig. 85 (left) reveal a mass dependence similar to the one caused by a collectively expanding medium in heavy-ion collisions. In particular, at low-$p_T$ they indicate a mass ordering effect where the $v_2$ of protons is shifted to larger $p_T$ with respect to pions and kaons. For comparison, we also show in Fig. 85 (left) measurements of $K_S^0$ and $\Delta$ in p–Pb collisions by the CMS Collaboration [894]. Despite the differences in the selection of multiplicity classes (based on the uncorrected number of tracks at midrapidity), wider pseudorapidity acceptance and larger pseudorapidity difference between correlated particles ($|\Delta\eta| > 1.0$), the measurements tend to follow the trend of the mass dependence of $v_2$ seen by ALICE for $\pi$, K, p hadrons.

The IP-Glasma+MUSIC+UrQMD model calculations [48], which provide a quantitative description of such measurements in Pb–Pb collisions (see Sect. 2.2.3), overestimate the data in p–Pb collisions at low-$p_T$ (see Fig. 85 (right)). While this could point to the need for finding a better estimate of the bulk viscosity in p–Pb collisions, owing to its strong impact on the identified particle flow in the model, it should be noted that in the experimental results the choice of the non-flow reference sample in the subtraction technique may influence the model-to-data comparison on a quantitative level. Besides hydrodynamics, other models based on hadronic rescatterings [884] and on parton transport [895] are able to generate the mass ordering effect. In addition, calculations based purely on initial-state momentum correlations in the IP-Glasma model coupled with Lund string fragmentation [896] revealed this feature too. Thus, the observation of mass ordering alone cannot be taken as evidence of flow-like collectivity in small systems. The measurements of $v_2$ of identified hadrons are further compared to the hydrodynamic model iEBE-VISHNU with the coalescence mechanism [472], allowing for comparison to the data to be extended to the intermediate-$p_T$ region (shown in Fig. 85 (right)). Even though the current statistical uncertainty of ALICE $v_2$ measurements is still considerable, the model provides an excellent description of the $v_2$ of light hadrons over the entire $p_T$ range; such an agreement can not be achieved without invoking coalescence. This indicates the importance of partonic degrees of freedom in the system evolution in
<figure><img src="image_99.png" /><figcaption>**Fig. 85** (Left) Transverse-momentum dependence of $v_2$ of hadrons, pions, kaons, protons, $J/\psi$, electrons from heavy-flavour hadron decays and inclusive muons in the 0–20% multiplicity class of p–Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV and 8.16 TeV after subtraction of the low-multiplicity class [890–893]. The muons at $p_T$ < 1.5 GeV/c are dominated by decays of light hadrons, while they predominantly originate from heavy-flavour hadron decays at $p_T$ > 2 GeV/c. The results of $J/\psi$ are combined results from p–Pb collisions at $\sqrt{s_{NN}}$ = 5.02 and</figcaption></figure>

<figure><img src="image_100.png" /></figure>

8.16 TeV. The data is plotted at the average-$p_T$ for each considered $p_T$ interval and particle species under study. Error bars show statistical uncertainties while boxes denote systematic uncertainties. The bands illustrate the CMS measurements for $K_S^0$ and $\Lambda$ [894]. (Right) Data of light hadrons are compared to the IP-Glasma+MUSIC+UrQMD [48] calculations and iEBE-VISHNU calculations with the coalescence and fragmentation model [472]. The width of the bands represents the statistical uncertainty of the models

high-multiplicity p–Pb collisions at the LHC, which in turn supports the hydrodynamic-like origin of the mass ordering at low-$p_T$ discussed above.

In Fig. 85 (left), $v_2$ measurements of the light hadrons presented above are compared to results of $J/\psi$, electrons from decays of open-heavy-flavour hadrons, and muons [890–893] (the latter are predominantly originating from open-heavy-flavour decays for $p_T$ > 2 GeV/c while being dominated by decays of light hadrons for $p_T$ < 1.5 GeV/c). The early formation time of heavy quarks suggests that their $v_2$ is mainly developed during the early stages of the collision, namely in hard parton scatterings off colour domains in the collision target. In this picture, the colour fields are locally organised into colour domains. Each parton scatters off these domains independently and receives a momentum kick in the process. If more partons that are in the same colour state scatter off the same domain, they will receive a similar kick, which may create anisotropies in the momentum space. As the orientation of the colour fields fluctuates event-by-event independently of the overall geometry, the initial momentum anisotropies are not directly correlated with the global spatial anisotropy [897]. Several model studies tend to support the idea of heavy-flavour hadron flow being driven by the initial gluon momentum anisotropies. These models are rooted in a purely initial-state calculation based on the colour Glass Condensate (CGC) formalism that was able to describe the measurements of flow of heavy-flavour hadrons [898,899], and a study showing that final-state interactions alone cannot generate sufficient $v_2$ of $J/\psi$ [900]. On the other hand, it was found that the momentum anisotropies from the initial state

are uncorrelated with the global spatial anisotropy [901]. Thus, the experimental measurements of the flow of heavy-flavour hadrons performed by correlating them with light hadrons, the flow of which tends to align with the spatial anisotropy (based on the current data and comparison to phenomenological models discussed above), would result in zero flow. This is in contrast to the measurements presented in Fig. 85 (left) exhibiting a finite $v_2$ of heavy-flavour hadrons, which could therefore be interpreted as an indication of final-state induced flow of these particles. It should be noted that it is not yet clear what is the contribution of initial-state correlations to flow of light hadrons, thus more investigations in this area are desirable in the future. While the flow of open heavy-flavour particles is composed of both the heavy and light quarks with poorly known relative contributions in the case of hadronisation via coalescence, the flow of $J/\psi$ may directly reflect the collective behaviour of the charm quark. This is supported by the apparent ordering of $v_2$ in Fig. 85 (left) which tends to be followed by the heavy flavour measurements. It should be noted that different samples of low-multiplicity events were used for non-flow subtraction in the case of the $v_2$ of $J/\psi$, which could potentially result in different remaining non-flow contamination in the final measurements and/or different contribution of $v_2$ from the subtracting peripheral bin. In addition, the rapidity intervals differ between the measurements, and the measured $p_T$ of heavy-flavour decay leptons does not directly correspond to the original $p_T$ of their parents. The data comparison is therefore not straightforward. Future data taking in small collision systems will provide larger samples of the rare particle
species, potentially allowing for their direct reconstruction and thus measurements of correlations of the original parent particles instead of using the decay products carrying partial kinematics information. In addition, it could be possible to perform the measurements using the state-of-the art methods to suppress non-flow, exploited so far only in the light hadron measurements, which would allow for quantitative comparisons.

Similar identified particle $v_n$ measurements in pp collisions are highly desirable to assess whether these features remain for even smaller collision systems. So far, only $v_2$ of $\Lambda$ and $K_S^0$ were measured in pp collisions by the CMS Collaboration [871], with a hint of mass dependence at low-$p_T$. Measurements of more particle species extended to the intermediate-$p_T$ region are therefore necessary to arrive at a conclusive data-model comparison. In addition, a more systematic study of the non-flow subtraction in two-particle correlations, which still results in very large uncertainties on the measurements, should be performed. One of the advanced approaches would be measuring anisotropic flow of identified hadrons using multi-particle cumulants with the subevent method, which has been developed for the ALICE projections of high-energy pp collisions in Run 3 at the LHC [902] but not yet performed in data due to the limited amount of collected data.

It was shown already in Pb–Pb collisions that measurements of $v_n$ alone are not enough to provide detailed constraints on theoretical models. The symmetric cumulants $SC(m, n)$, i.e. correlations between the second moments of different flow harmonics, possess the ability to constrain the initial conditions [307]. The most sensitive correlation is between $v_2^2$ and $v_3^2$, i.e. $SC(3,2)$, while the correlation between the second moments of the second and fourth harmonics, $SC(4,2)$, is sensitive to the transport properties of the medium in heavy-ion collisions (see more details in Sect. 2.2.5). In a similar spirit, the determination of the origin of collectivity in small systems can be addressed with other measurements, sensitive to different aspects of the system evolution and their theoretical description, in addition to $v_n$ coefficients. Results for symmetric cumulants have shown a similar trend with collision multiplicity in both large and small systems [280], indicating that the observed collective correlations share a similar origin, which may vary with multiplicity.

The measurements of azimuthal particle correlations presented in this section revealed the presence of collective effects in small collision systems. Available model comparisons to the results at high multiplicities suggest that the particle correlations originate as a system’s response to the initial spatial anisotropy during final state interactions. In particular, a strong indication for the presence of a partonic system was found in high-multiplicity p–Pb collisions. On the other hand, at multiplicities of the order of the average multiplicity,

the initial-state (gluon momentum) correlations may become significant [896,903–905]. It is very challenging to confirm this with the present experimental results, but the upcoming Run 3 or new observables could promise to provide additional information to address the relative importance of initial correlations to the measured flow.

### 3.4 Charmonium and bottomonium suppression in p–Pb collisions

Studies of quarkonium production in collision systems involving light ions or protons were carried out since the very beginning of QGP investigations. In particular, the study of ion-ion collisions with systems lighter than Pb–Pb or Au–Au were meant to obtain information about the onset of suppression effects or to better understand scaling properties of the suppression as a function of variables related to the geometry of the system (In–In collisions at SPS energy by NA60 [906], Cu–Cu [907] and Cu–Au [908] by PHENIX at RHIC). Studies of p–A (or d–Au at RHIC) collisions were also fundamental for the investigation of the role and the size of cold nuclear matter effects that were found to be responsible for a significant fraction of the observed quarkonium suppression both at SPS and RHIC energies [695,712,714,909–913]. Most of these results concerned the $J/\psi$ state, but also the production of $\psi(2S)$ and $\Upsilon(1S)$ was studied although with larger uncertainties.

At the LHC, ALICE has obtained a result on the inclusive $J/\psi R_{AA}$ in Xe–Xe collisions at $\sqrt{s_{NN}} = 5.44$ TeV [914]. The integrated luminosity available was rather low ($L_{int} \sim 0.34~\mu b^{-1}$), but the observed compatibility between $R_{AA}$ values for Xe–Xe and Pb–Pb at a given $N_{part}$ has shown that the relative contribution of suppression and (re)generation processes is similar for collisions involving a similar number of participant nucleons in different nuclear collision systems.

On the other hand, much more accurate studies of quarkonium production in p–Pb collisions were carried out, at $\sqrt{s_{NN}} = 5.02$ and 8.16 TeV involving $J/\psi$, $\psi(2S)$ and the $\Upsilon(1S, 2S, 3S)$ states [748,749,752,755,915,916]. These investigations were meant to evaluate the influence of initial state effects, and in particular of nuclear shadowing, on the observed quarkonium yields. Remarkably, final-state cold nuclear matter effects, related to the break-up of the resonances via interactions with the nucleons of the colliding ion(s), are negligible at LHC energies due to the extremely small time spent by the heavy-quark pair in the nuclear matter. In fact, $\tau = \langle L \rangle / \beta_z \gamma$ values ranging from $\sim 7 \times 10^{-2}$ to $10^{-4}$ fm/c for low-$p_T$ $c\bar{c}$ pairs are respectively obtained at backward (Pb-going) and forward (p-going) rapidity, with $\langle L \rangle$ being the average length of nuclear matter crossed by the pair, $\beta_z$ the velocity of the $c\bar{c}$ along the beam direction in the nucleus rest frame, and $\gamma = E_{c\bar{c}} / m_{c\bar{c}}$ [917]. In fact, nuclear shadowing was found to be the main effect influenc-
<figure><img src="image_101.png" /><figcaption>**Fig. 86** The nuclear modification factor for J/ψ and ψ(2S) production in p–Pb collisions at $\sqrt{s_{NN}}$ = 8.16 TeV [915]. Data are compared with two model calculations that include final state effects [751,919]</figcaption></figure>

ing the production of ground states (J/ψ, Υ(1S)), as shown in Fig. 59 in Sect. 2.5.3, where it was demonstrated that shadowing calculations based on the EPPS16 set of nPDFs fairly reproduce the observed $R_{pPb}$ values.

However, the situation significantly changes when considering charmonium and bottomonium states with lower binding energy. In the charmonium sector, ALICE studied the rapidity, transverse momentum and centrality dependence of the $\psi(2S)$ $R_{pPb}$ [915–918]. If nuclear shadowing continues to be the main effect at play, similar values of the nuclear modification factors should be expected for J/ψ and ψ(2S), since the Bjorken-x values of the partons involved in the hard scattering are different by at most a few percent, due to the small mass difference between the two states. This comparison is shown in Fig. 86 as a function of rapidity for $\sqrt{s_{NN}}$ = 8.16 TeV. At forward y (p-going), the J/ψ and ψ(2S) $R_{pPb}$ are fully compatible, while at backward y (Pb-going) a stronger suppression is observed for ψ(2S). Considering the correlation of the uncertainties, the significance of this observation is $3.1\sigma$. This effect is found to have no appreciable dependence on collision energy and has no strong dependence as a function of $p_T$ (not shown).

These observations suggest that the suppression seen at backward rapidity for the $\psi(2S)$ must be due to final state effects and more precisely to the dissociation of the weakly bound $c\bar{c}$ pair in the strongly-interacting system created in the collision. The ALICE measurements are compared to calculations from two theory models which include, in addition to initial state effects, also the final state dissociation of the quarkonia. In the comover model [919], the break-up originates from “comovers”, i.e. particles that travel along with the $c\bar{c}$ pair whose density is constrained by the measured $dN_{ch}/dy$, which interact with the quarkonia with

<figure><img src="image_102.png" /><figcaption>**Fig. 87** The double ratio between the $\psi(2S)$ and J/ψ cross sections in p–Pb and pp collisions at $\sqrt{s_{NN}}$ = 5.02 and 8.16 TeV [917,918], shown as a function of $\langle N_{coll} \rangle$ for the backward (Pb-going) rapidity region and compared to the comover model results [919]. The boxes around unity correspond to the global systematic uncertainty on the ratio between the $\psi(2S)$ and J/ψ cross sections in pp collisions</figcaption></figure>

cross sections tuned on lower energy experimental data. The CGC+ICEM model [751] is based on an improvement of the Colour Evaporation Model [920] for the production process and on the Colour Glass Condensate effective theory for the treatment of saturation effects in the small-$x$ regime [921]. In this model, the suppression is here due to parton comovers hadronising on longer time scales than the $c\bar{c}$ pair that have soft colour exchanges with the final-state charmonia. The models are in fair agreement with data (only the comover model can give predictions in the Pb-going rapidity region), advocating that at LHC energies a dense interacting system is produced, able to selectively dissociate some of the quarkonium states.

A suppression effect was already observed by PHENIX when studying the centrality dependence of the ratio of the $\psi(2S)$ and J/ψ nuclear modification factors in d–Au collisions at $\sqrt{s_{NN}}$ = 0.2 TeV, where it was shown that for central collisions ($\langle N_{coll} \rangle \sim 15$) the $\psi(2S)$ suppression was 3 times stronger than the one of J/ψ [910]. The corresponding ALICE result on the centrality dependence of the double-ratio between the $\psi(2S)$ and J/ψ cross sections in p–Pb and pp collisions, shown in Fig. 87 for data taken at $\sqrt{s_{NN}}$ = 5.02 and 8.16 TeV, also indicates an increase of the relative suppression of $\psi(2S)$, which already sets in at smaller values of $\langle N_{coll} \rangle$ [917,918]. The comover model fairly reproduces this suppression [919].

Similar studies were also carried out in the bottomonium sector by comparing the suppression effects between the ground state Υ(1S) and the 2S and 3S states. In particular, a stronger relative suppression of Υ(2S) with respect to Υ(1S) might indicate the presence of final-state effects
on the former resonance which has a binding energy of about 500 MeV, similar to that of the $J/\psi$ and significantly larger than that of the $\psi(2S)$. ALICE measured the double-ratio between the $\Upsilon(2S)$ and $\Upsilon(1S)$ in p-Pb with respect to pp at $\sqrt{s_{NN}} = 8.16$ TeV, obtaining $R_{2S/1S} = 0.85 \pm 0.16$ (stat.) $\pm 0.06$ (syst.) in the Pb-going rapidity region where final-state effects, based on $\psi(2S)$ results, should be stronger [749]. This value is described by comover model calculations, which give $R_{2S/1S} = 0.77 \pm 0.04$ [919]. A stronger experimental indication ($2.8 \sigma$) for a value of $R_{2S/1S}$ smaller than unity was obtained by ATLAS at $\sqrt{s_{NN}} = 5.02$ TeV in the rapidity region $-2 < y_{cms} < 1.5$ ($R_{2S/1S} = 0.765 \pm 0.069$ (stat.) $\pm 0.048$ (syst.)) [922]. Finally, for the even more weakly bound $\Upsilon(3S)$, evidence for a stronger suppression with respect to $\Upsilon(1S)$ was obtained by LHCb ($R_{3S/1S} = 0.44 \pm 0.15$ at $\sqrt{s_{NN}} = 8.16$ TeV for $-4 < y_{cms} < -2.5$ [923]), with the corresponding ALICE results being $R_{3S/1S} = 0.87 \pm 0.27$ (stat.) $\pm 0.07$ (syst.) [749].

### 3.5 Searches for jet quenching in small systems

Measurements of small collision systems explore the question of the limits of QGP formation. Such measurements, at both RHIC and the LHC, have revealed significant effects in the low-$p_T$ sector that are associated with collective dynamics of the QGP in large collision systems (Sects. 3.2 and 3.3). However, interpretations other than QGP formation have been proposed for these small-systems phenomena, based on quantum interference or more conventional hadronic physics (see e.g. [886,924]).

Jet quenching is a necessary consequence of the formation of an extended QGP system. Current theoretical approaches for the magnitude of jet quenching in small systems are based on calculations of jet quenching in large systems. Their predictions vary greatly, however, from significant and readily observable signals [925] to negligible effects [926]. In addition, there may be significant initial-state effects that can mask signatures of jet quenching [927,928]. Current experimental searches for jet quenching in small systems have taken two different approaches: yield measurements of inclusive hadrons or jets, analogous to the measurements of $R_{AA}$ in large systems [565, 566, 623, 929–936] and coincidence measurements [570,937]. To date, no significant jet quenching effects have been observed, within experimental uncertainties.

In this section we discuss ALICE measurements using both approaches to searching for jet quenching in small systems. The focus of this program in ALICE has been on p-Pb collisions, with measurements carried out both inclusively and selected by Event Activity (EA), as measured by forward charged-particle multiplicity in the V0 detectors or by beam-rapidity neutral energy in the ZDCs.

For inclusive observables, jet quenching is measured using $R_{pPb}$, the ratio of the inclusive yield measured in p-Pb collisions to the inclusive production cross section of the same observable in pp collisions scaled by the geometric factor $\langle T_{pPb} \rangle$, where $\langle \dots \rangle$ indicates the average over the EA interval used for event selection. Suppression in inclusive hadron production due to jet quenching would correspond to $R_{pPb} < 1$.

The value of $\langle T_{pPb} \rangle$ is determined using Glauber modeling (Sect. 2.1.1), which is based on the assumption that EA is closely correlated with event geometry or “centrality”. However, as discussed in Sect. 3.1, the correlation between EA and event geometry in small systems is subject to large fluctuations, and the EA distribution is biased by the presence of a hard process in the event.

Figure 88 shows the measurement of this hard-process bias in p-Pb collisions [570]. The blue histogram is the probability distribution in decile bins of the V0A signal (forward charged multiplicity in the Pb-going direction) for the minimum-bias (MB) event population. By construction, this distribution is uniform; indeed, that is the criterion for setting the decile bin boundaries in V0A signal. The black and red histograms show this probability distribution using the same bin boundaries, but for events selected with a high-$p_T$ charged track in the central barrel for two $p_T$ intervals (“TT{12,50}” denotes $12 < p_T < 50$ GeV). In this case the probability is biased towards the low percentile bins, corresponding to higher V0A signal, with little to no dependence on track $p_T$. Such a correlation between soft and hard processes is well-established, and phenomenological approaches have been developed to model it. A recent calculation using the PYTHIA Angantyr model is able to reproduce this measurement, attributing the bias to the increase in the contribution of MPI in the presence of hard processes [938].

Such biases are minimised experimentally by utilising forward neutron energy in the ZDC for EA [855], providing the largest possible phase-space separation between the EA measurement and jet quenching observables in the central barrel. However, this approach requires the model assumption that ZDC signal is correlated with centrality through the emission of “slow neutrons,” where “slow” refers to beam rapidity [939]. Since the modeling of forward slow neutron production is complex, with significant systematic uncertainties, ALICE defines a model-dependent variant of $R_{pPb}$, denoted $Q_{pPb}$, to report jet quenching for inclusive observables using EA-selected p-Pb data.

The irreducible modeling uncertainties for correlating $\langle T_{pPb} \rangle$ with collision geometry can be avoided by carrying out measurements using MB p-Pb data, in which there is no EA-based event selection. For MB pA collisions, hard cross sections scale precisely with the nuclear mass $A$ [714,941] so that $\langle T_{pPb} \rangle = A/\sigma_{pPb}^{\text{inel}}$, where $\sigma_{pPb}^{\text{inel}}$ is the p-Pb inelastic cross section. Therefore, $R_{pPb}$ is constructed as the ratio of the production cross section in p-Pb divided by the produc-
<figure><img src="image_103.png" /><figcaption>**Fig. 88** Probability distribution of EA (VOA) in p–Pb events, in decile bins [570]. Blue: MB events; red and black: events selected by presence of a high-$p_T$ charged track in the central region. See text for details</figcaption></figure>

<figure><img src="image_104.png" /><figcaption>**Fig. 89** Inclusive hadron production in p–Pb collisions. (Top) Inclusive yield ratio $R_{pPb}$ for MB p–Pb collisions, for charged hadrons [594], $\pi^0$ [940], and prompt D mesons [864]. (Bottom) Inclusive yield ratio $Q_{pPb}$ for EA-selected p–Pb collisions (0–10% ZNA), for inclusive charged hadrons [855] and prompt D mesons [864]. Global normalisation uncertainties of 4% for $R_{pPb}$ and 7% for $Q_{pPb}$ are not shown in the figure</figcaption></figure>

<figure><img src="image_105.png" /><figcaption>**Fig. 90** Ratio of semi-inclusive recoil jet distributions, $\Delta_{recoil}$, for high-EA and low-EA p–Pb collisions [570], using V0A. Charged-particle jets are reconstructed with anti-$k_T$, $R = 0.4$. Vertical lines show the statistical error. Grey boxes show the systematic uncertainty of the ratio, which accounts for correlated uncertainty in numerator and denominator. The red line indicates the ratio for a $p_T$-shift of the high-EA distribution of –0.4 GeV, corresponding to the 90% CL limit in medium-induced energy transport out of the jet cone from this measurement</figcaption></figure>

cross section in pp and by the mass number A. While the MB measurement has the advantage of smaller systematic uncertainties, it provides less dynamic range in collision volume than the EA-selected observables, and both should be considered together.

Figure 89 shows ALICE measurements of inclusive hadron production in p–Pb collisions. The upper panel shows $R_{pPb}$ for MB p–Pb collisions, for inclusive charged hadrons, $\pi^0$, and prompt D mesons. The lower panel shows $Q_{pPb}$ for EA-selected p–Pb collisions (0–10% ZNA), likewise for inclusive charged hadrons and prompt D mesons. For $p_T > 5$ GeV, where the yields of all particle species are dominated by hard-scattering processes, $R_{pPb}$ and $Q_{pPb}$ are consistent with unity within the experimental uncertainties for all particle species. Within the sensitivity of these measurements, there is no evidence of jet quenching in p–Pb collisions.

In the second approach to the search for jet quenching, based on coincidence observables, ALICE utilises a semi-inclusive observable to measure the trigger-normalised yield of reconstructed jets recoiling from a high-$p_T$ charged-hadron trigger [568,570]. In this approach, the trigger population is sampled, consistently with the inclusive hadron $p_T$ distribution, and all recoil jets are counted. This ensures that the resulting trigger-normalised recoil jet distribution is equivalent to the ratio of hard cross sections for hadron+jet and inclusive hadron production. Since the numerator and denominator in this ratio both scale with $\langle T_{pPb} \rangle$ for EA-selected data, the geometric factor $\langle T_{pPb} \rangle$ can-
cels exactly [570]. This observable, therefore, enables the measurement of jet quenching in EA-selected data whose interpretation does not require the assumption that EA and collision geometry are correlated, in contrast to $Q_{pp}$. This observable therefore has greater systematic precision for searching for small jet quenching signals in EA-selected data from small systems, and does so based on reconstructed jets rather than hadrons.

Jet measurements in high-energy nuclear collisions must contend with large uncorrelated backgrounds (see Sect. 2.4). In order to subtract such backgrounds in a fully data-driven way, ALICE utilises the distribution of $\Delta_{recoil}$, which is the difference of two semi-inclusive recoil jet distributions with high-$p_T$ TT and low-$p_T$ TT [568,570]. Figure 90 shows the ratio of $\Delta_{recoil}$ distributions for high-EA p-Pb (0–20% V0A) and low-EA p-Pb (50–100% V0A), for charged-particle jets with $R = 0.4$. The ratio is compatible with unity, again indicating no evidence of jet quenching in high-EA p-Pb events, consistent with the inclusive measurements based on $R_{ppb}$ and $Q_{ppb}$.

In this case, however, the experimental uncertainty is dominated by the statistical error, and it is possible to measure a limit on jet quenching effects. The $\Delta_{recoil}$ distribution as a function of $p_T$ is well-approximated by an exponential function. Assuming that jet quenching corresponds to medium-induced energy transport out of the jet cone whose magnitude is independent of $p_T^{jet}$, such an energy loss corresponds to a uniform $p_T$-shift in the $\Delta_{recoil}$ spectrum, and uniform suppression in the ratio shown in the figure. The red line in the figure corresponds to a $p_T$-shift of –0.4 GeV of the high-EA spectrum, which is the one-sided 90% CL limit on medium-induced energy transport out of the jet cone $R = 0.4$ in high-EA p-Pb collisions. This limit should be compared to the result of the equivalent measurement for central Pb-Pb collisions, in which significant recoil yield suppression is observed, corresponding to a recoil jet spectrum shift of (–8±2) GeV [568]. The limit of –0.4 GeV shown in Fig. 90 for high-EA p-Pb is therefore a factor 20 smaller than the magnitude of jet quenching effects observed in central Pb-Pb collisions. While this measurement provides the most precise limit on jet quenching in small systems currently available, theoretical calculations are needed to determine whether it is still compatible with the mechanisms thought to underlie jet quenching in large systems.

Projections for Runs 3 and 4 at the LHC indicate that the limits on jet quenching from such a semi-inclusive approach will be in the vicinity of 0.1–0.2 GeV for an array of small systems, including high-EA pp high-EA p-Pb, and EA-selected O–O collisions. Additional studies are underway to elucidate and quantify the correlations between forward EA and high-$Q^2$ processes in the central region.

## 3.6 Conclusions

**Event selection.** The study of event and multiplicity selection in small collision systems has been proven to be crucial for a proper understanding of the physics observables across different colliding nuclei. These studies indicate that further measurements in the field of relativistic heavy-ion collisions must take proper care to understand such selections carefully, especially in high-multiplicity proton–proton and proton-nucleus collisions.

**Identified particle production.** A comprehensive set of identified particle measurements from low-multiplicity proton–proton to highest multiplicity Pb–Pb collisions has conclusively demonstrated that identified particle ratios such as $\Xi/\pi$ and $\Omega/\pi$ evolve continuously as a function of charged-particle multiplicity density, regardless of collision system or beam energy. Furthermore, specific final-state effects such as the depletion/excess of baryons with respect to mesons at low-/intermediate- transverse momentum have been observed to be present throughout all collisions systems in different magnitudes. Recent measurements have even extended this pattern to the $\Lambda_c/D^0$ ratio. A successful description of all these phenomena for all collision systems remains a theoretical challenge of the field.

**Collective behaviour.** Measurements of azimuthal correlations and anisotropic flow in small collision systems exhibit features of a collectively expanding system, similar to those observed in heavy-ion collisions, where they are believed to originate from the presence of QGP medium. The origin of the collective effects depends on particle multiplicity of a collision. While hydrodynamic-like description seems to be favored by data especially at high multiplicities, the effects of initial state effects from initial gluon momentum correlations may play an important role at low-$N_{ch}$.

**Charmonium and bottomonium.** Charmonium measurements in p–Pb collisions indicate the $\psi(2S)$ suppression in the backward rapidity region (Pb-going) to be larger with respect to that of the $J/\psi$. This is not expected from cold nuclear matter effects and is suggestive of final state interactions between the weakly bound $c\bar{c}$ and a strongly interacting partonic or hadronic medium, leading to its dissociation. Similar suppression patterns are observed also in the bottomonium sector for $\Upsilon(2S)$ and $\Upsilon(3S)$.

**Searches for jet quenching.** Measurements in small systems have not observed significant jet quenching effects thus far within the sensitivity of the measurements. These measurements however provide the first quantitative limit on jet quenching in small systems, corresponding to medium-induced energy transport out of a jet cone with radius 0.4 less
than 400 MeV at 90% confidence level in high-event-activity p–Pb collisions.

## 4 The initial state of the collision

An accurate description of the initial state of nucleus–nucleus collisions is crucial for the interpretation of all measurements carried out in heavy-ion collisions at the LHC. In particular, the description of hard processes strongly depends on the longitudinal momentum distribution of partons inside the nucleus, while the spatial transverse profile of partons and its fluctuations are important to properly interpret observables related to azimuthal anisotropies and correlations.

Parton distribution functions (PDFs) of protons, describing the probability to find a parton in the proton with a longitudinal momentum fraction $x$ at resolution scale $Q^2$, show a rapid growth towards small values of $x$. This growth is expected to slow down and saturate due to unitarity constraints [942]. However, no compelling evidence for saturation has been found so far. Besides, PDFs are significantly modified inside the nucleus with respect to the free proton PDFs. In particular, nuclear PDFs at $x < 10^{-2}$ demonstrate a clear suppression compared to proton PDFs, a phenomenon referred to as nuclear shadowing [943]. The most recent parameterisations of nuclear PDFs, such as nCTEQ15 [46], EPPS16 [47], or nNNPDF [944], are based on global QCD fits to available data samples of nuclear deep-inelastic scattering including p–Pb data on dijets and electroweak bosons. These parameterisations are affected by large uncertainties due to the limited kinematic coverage of the available data and the largely indirect determination of the gluon distributions. Latest results on the heavy-flavour and charmonium nuclear modification factors in p–Pb collisions at the LHC provide new stringent constraints on nuclear PDFs [945–947]. However, the role of final state effects in these measurements, such as energy loss and radial flow, still needs to be clarified. On the contrary, electroweak-boson measurements are not affected by final-state effects and can serve as particularly clean probes of light quark PDFs in nuclei at large scales $Q^2 \sim M_{W.Z.}^2$. ALICE results on $Z^0$ and $W^{\pm}$ boson production in p–Pb and Pb–Pb collisions are summarised in Sect. 4.1.

ALICE measurements of heavy-vector-meson photoproduction in ultra-peripheral collisions, described in Sect. 4.2, are sensitive to the gluon distribution in nuclei at hard scales corresponding to the heavy-quark mass [948], while $\rho^0$ photoproduction also allows for the study of the approach to the black-disk limit in QCD at semi-hard scales [949]. The rapidity dependence of the coherent cross section of heavy vector mesons gives information on the energy evolution of the gluon distribution, providing one of the cleanest probes of shadowing at small $x$ down to $10^{-5}$. So

far, these results have not been included in the global parton analyses to constrain the nuclear gluon distributions in heavy nuclei, the reason being uncertainties in the choice of scale and implementation of NLO effects. Intense theoretical work, however, is ongoing to resolve these issues [950,951]. Furthermore, the transverse-momentum distribution of the coherently-produced vector meson encodes information on the average distribution of colour fields in the impact-parameter plane [952], while the incoherent production of vector mesons is sensitive to fluctuations of the gluon distribution, particularly to the different possible geometries of the nuclear initial state in the impact-parameter plane [953].

The characterisation of the spatial distribution of nuclear matter in the initial state is also a crucial task in the study of relativistic heavy-ion collisions. Not only it allows mapping the structure of the nucleus at high-energy, but this knowledge of the initial state also helps constrain calculations of key QGP medium properties, which were discussed in the previous chapters. The initial state models generate energy or density profiles prior to the time of the formation of the QGP in the overlap region of a collision, and this can be accomplished in several ways. Almost all implementations take advantage of the assumptions associated with the Monte Carlo (MC) Glauber approach [42]. This approach assumes the nucleons in the nucleus are positioned randomly according to the Woods–Saxon distribution, the nucleons travel in an unperturbed trajectory irrespective of whether they interact with other nucleons, and the criterion for a nucleon–nucleon interaction depends on the inelastic cross section (which can be inferred from experimental measurements in pp collisions). In the simplest MC Glauber nucleon model, the density of the overlap region in a collision is the sum of the densities of the participating nuclear matter. The $T_{RENTO}$ model [462] also uses nucleons as the relevant degrees of freedom, however in its optimal configuration, it assumes the overlap densities are the square root of the product of the individual nuclear densities. The MC-KLN [954,955] and IP-Glasma [44,956] models take a different approach, and assume that the relevant degrees of freedom are gluons in the nucleons. They use the Colour Glass Condensate effective theory coupled with saturation equations to determine the density profiles of gluons in the overlap region. Another initial state model, EKRT [309,957], determines parton densities from next-to-leading-order perturbative QCD using the saturation conjecture. Finally, the MC Glauber constituent quark model [958] is an extension of the MC Glauber nucleon model, however the relevant degrees of freedom are constituent quarks in the nucleons.

Particle multiplicity distributions and the anisotropic flow measurements can be extremely valuable in the understanding of the initial state. After two decades of developments, viscous hydrodynamic calculations using various initial-state models can often describe the multiplicity distributions,