Phenomenological Review on Quark-Gluon Plasma: Concepts vs. Observations

In this review, we present an up-to-date phenomenological summary of research developments in the physics of the Quark--Gluon Plasma (QGP). A short historical perspective and theoretical motivation for this rapidly developing field of contemporary particle physics is provided. In addition, we introduce and discuss the role of the quantum chromodynamics (QCD) ground state, non-perturbative and lattice QCD results on the QGP properties, as well as the transport models used to make a connection between theory and experiment. The experimental part presents the selected results on bulk observables, hard and penetrating probes obtained in the ultra-relativistic heavy-ion experiments carried out at the Brookhaven National Laboratory Relativistic Heavy Ion Collider (BNL RHIC) and CERN Super Proton Synchrotron (SPS) and Large Hadron Collider (LHC) accelerators. We also give a brief overview of new developments related to the ongoing searches of the QCD critical point and to the collectivity in small ($p+p$ and $p+A$) systems.

Experimental attack on producing QGP under laboratory conditions started in the world-leading particle physics facilities at CERN (Conseil Européen pour la Recherche Nucléaire, Geneva, Switzerland) and Brookhaven National Laboratory (BNL, New York, NY, USA) in the late 1980s [10][11][12].
In the year 2000, after finishing the main part of its heavy-ion program at the Super Proton Synchrotron (SPS) accelerator, CERN announced circumstantial evidence for the creation of a new state of matter in Pb + Pb collisions [13]. The real discovery of QGP took place in 2005, when four international collaborations studying Au + Au collisions at the Relativistic Heavy Ion Collider (RHIC) at BNL announced the results of their first five years of measurements [14][15][16][17]. Surprisingly, the properties of the new state of matter [18] differed markedly from predictions made over many years before its discovery 2 .
This review aims at giving a short historic introduction into the vast research field of QGP physics and the underlined phenomenological aspects with a comprehensive list of corresponding references. Effects of the hot/dense medium such as the nuclear suppression, initial-state interactions, in-medium energy loss, color screening and saturation are important for a proper understanding of the collective phenomena in heavy-ion collisions, and are included in the scope of this review. Besides, we have qualitatively overviewed and confronted with existing observations such fundamental theoretical concepts as the QCD vacuum and phase diagram, equation of state (EoS) of deconfined QCD matter, initial-state effects, collectivity, flow, hydrodynamic properties of the QGP, and associated electromagnetic effects.
The paper is organized as follows. Section 2 presents a short history of the theoretical understanding of extreme states matter. The phase diagram of QCD is discussed in Section 3. Section 4 gives a historical perspective on experiments operating with with collisions of heavy ions. Section 5 describes the basic signatures of QGP production, while current developments in QGP research are provided in Section 6. The concluding remarks are given in Section 7. For further reading on fundamental concepts and the latest studies of QGP dynamics, we recommend several textbooks [9,[20][21][22][23][24] and reviews [25][26][27][28][29][30] published in recent years.

Matter under Extreme Conditions
The properties of matter under extreme conditions at high values of state variables have always attracted the curiosity of scientists, owing to the possibility of advancing to new domains of the phase diagram and producing the exotic states of matter in laboratory [31,32]. The first attempts to discuss the properties of matter at densities well above the normal nuclear density ρ 0 = 2 × 10 14 g·cm −3 (≈0.16 GeV·fm −3 ) date back to the seminal Oppenheimer-Volkoff paper from 1939 [33]. A study of the gravitational stability of a new phase of neutron matter suggested a few years earlier by Landau [34] led them to carry out their computations to several tens of ρ 0 before smoothly extrapolating the results to the black-hole singularity. In 1962, when discussing relativistic limitations of the equation of state (EoS) of matter made of baryons interacting via a massive vector field, Zeldovich used the density twenty times exceeding ρ 0 [35]. In 1976, the same value of density was shown by Baym and Chin [36] to be energetically favourable for the neutron matter-quark phase transition.
In cosmology, a very dense matter with ρ ∝ 10 6 g·cm −3 (≈1 eV·fm −3 ) was first studied in 1946 by Gamow [37] when discussing the relative abundances of elements in the universe. The key discovery of the cosmic microwave background radiation by Penzias and Wilson in 1965 [38] not only provided a strong basis for the hot universe scenario which was used by Gamow, but also motivated Sakharov [39] to push it ad extremum. Considering the properties of hot matter at densities when gravitational interaction between photons becomes significant, he established that the absolute maximum of the temperature of any substance in equilibrium with radiation is on the order of Planck temperature 2 For a representative collection of papers tracing the development of theoretical ideas on the QCD deconfining phase transition before nineties , see Ref. [19]. Current knowledge on the QCD phase diagram is summarized on the left panel of Figure 1.
The arrows indicate the expected crossing through the deconfinement transition during the expansion phase in heavy-ion collisions at different accelerators. The red and black full circles denote the critical endpoints of the chiral and nuclear liquid-gas phase transitions, respectively. The (dashed) freeze-out curve indicates where hadro-chemical equilibrium is attained at the final stage of the collision. The nuclear matter ground-state at T = 0 and µ B = 0.93 GeV, and the approximate position of the QCD critical point at µ B ∼ 0.4 GeV are also indicated. The dashed line is the chiral pseudo-critical line associated with the crossover transition at low temperatures. Comparing this diagram to the phase diagram of water shown on the right panel, one notices that (at least theoretically) the complexity of the former approaches the latter.

The Role of the QCD Ground State
The quantum ground state of QCD plays an immensely important role in both particle physics and cosmology. In particular, the quark-gluon condensate is responsible for spontaneous chiral symmetry breaking, color confinement, and hadron mass generation (for a comprehensive review on the QCD vacuum, see Refs. [73][74][75][76][77] and references therein). It determines the properties (and possibly the generation mechanism) of the quark-gluon plasma, and dynamics of the phase transitions and hadronization. The latter phenomena are the most critical QCD phenomena taking place beyond the Perturbation Theory (PT), and thus are very difficult to explore by means of the well-known approaches. This strongly motivates further even deeper studies in this direction.
Let us start with the classical Yang-Mills (YM) gauge theory in the SU(N c ) (N c = 3 for QCD) determined by the gauge-invariant Lagrangian where F a µν = ∂ µ A a ν − ∂ ν A a µ + g s e abc A b µ A c ν is the gluon field stress tensor with SU(N c ) adjoint a, b, c = 1, . . . N 2 c − 1 and Lorentz µ, ν = 0, 1, 2, 3 indices, and with the strong coupling constant g s . The generating functional of such a classical theory is given by the Euclidean functional integral which is dominated by minima of the classical action S cl [A] corresponding to the classical vacuum state with F a µν = 0 unaltered by quantum corrections. The field excitations about the classical YM vacuum are referred to as instantons [78,79].
In fact, the classical YM equations of motion corresponding to Equation (1) are form-noninvariant with respect to small quantum fluctuations which break the conformal invariance of the gauge theory [80]-the effect known as the conformal (or trace) anomaly. Indeed, there is no threshold for the vacuum polarisation of a massless quantum gluon field by its classical component such that the solutions of the classical YM equations are unstable w.r.t the radiative corrections and cannot be used in physical applications. The conformal anomaly in QCD has notable implications; for example, in Cosmology, leading to an appearance of the Lorentz-invariant negative-valued contribution to the cosmological constant, where the one-loop expression for the QCD β-function β = −bα s /(4π), b = b eff = 9 accounting for three light flavours u, d, s (for pure gluodynamics, b = b g = 11) is typically used. Besides the wrong sign, the QCD vacuum density QCD is over forty orders of magnitude larger in absolute value than the positive cosmological constant observed in astrophysical measurements, The nonperturbative QCD vacuum effect is expected to be dynamically cancelled at macroscopically large distances in the course of cosmological expansion (see Ref. [81] and references therein). A dynamical mechanism of such a cancellation of vacua terms is yet unknown (for the existing scenarios discussed in the literature, see Refs. [81][82][83][84].
Consider a consistent effective Lagrangian formulation of the YM theory incorporating the conformal anomaly. In the corresponding variational technique, the strong coupling g s is treated as an operator depending on operators of quantum gluon fields by means of the RG equations in the operator form. Namely, the gauge field operator A a µ is considered as a variational variable which-together with the corresponding stress tensor operator-are related to those in the standard normalisation as follows The effective action and Lagrangian operators of the quantum gauge theory is given in terms of the gauge-invariant operator of the least dimension J by [85] respectively, whose variation w.r.t A a µ leads to the energy-momentum tensor of the gauge theory T ν,g In Equation (6), as a normalisation point, one can choose, for example, the strong coupling in the minimum of the effective action g 2 s, * = g 2 s (J = J * ). One distinguishes chromomagnetic B 2 > E 2 and chromoelectric B 2 < E 2 condensates, such that one or both of the corresponding ground-state solutions (minima of the effective action) should be stable in order to contribute to the physical QCD vacuum.
The strong coupling dependence on J is determined by the RG evolution equation The effective action (6) can be considered as an effective classical model [85] which possesses well-known properties of the full quantum theory such as (i) local gauge invariance; (ii) RG evolution and asymptotic freedom; (iii) correct quantum vacuum configurations; and (iv) trace anomaly.
In Ref. [84], it was noticed that the effective YM equation of motion in an expanding universe with the conformal metric g µν = a 2 diag(1, −1, −1, −1) (g ≡ det(g µν )) has a partial nonperturbative solution β(g 2 s, * ) = 2, where g 2 s, * = g 2 s (J * ) is the solution of the RG Equation (8) indicates that the QCD vacuum QCD < 0 indeed has a chromomagnetic nature B 2 > E 2 for g 2 s, * > 0 in the deeply non-perturbative domain. This means that the corresponding solution is stable; namely, any small perturbation around the vacuum state effectively vanishes at the typical QCD time scale ∆t ∼ 1/Λ QCD .
In asymptotically free gauge theories like QCD, the quantum vacuum configurations are controlled by the strong coupling regime. Performing an analysis in Euclidean spacetime, in Ref. [85] it was shown that the vacuum value of the gauge invariant J in a strongly-coupled quantum gauge theory does not vanish as it does in the classical gauge theory, and the corresponding functional integral is not dominated by the minima of the classical action (2). Moreover, it was shown that there are no instanton solutions to the effective action (6), such that the ground state of the quantum YM theory does not contain the classical instanton configurations. Instead, the quantum vacuum can be understood as a state with ferromagnetic properties which undergoes spontaneous magnetisation, providing a consistent description of the nonperturbative QCD vacuum and confinement alternative to the conventional instanton model.
How to understand the smallness of the observed cosmological constant within the effective QCD action approach? One way elaborated in Ref. [84] is to assume that such a compensation happens due to the presence of an additional QCD-like dynamics-Mirror QCD-with a confinement scale Λ mQCD Λ QCD . The corresponding nonperturbative Mirror QCD vacuum contribution may have an opposite sign to that in QCD, and can therefore compensate the QCD one at a certain time scale in the course of cosmological expansion. Another interesting possibility explored in Ref. [82] is to assume that the QCD vacuum is degenerate itself, and at a given time scale consists of two opposite-sign (quantum-topological and quantum-wave) contributions. Then, as soon as such a compensation occurs, the observable small cosmological constant can be generated by means of weak gravitational interactions in the QCD vacuum. Both possibilities, however, require a fine tuning of vacuum parameters in order to provide an exact compensation of bare (zeroth-order in gravitational interactions) QCD contributions to the ground state density.
Can one avoid such a major fine tuning problem? The stable ground-state solution J * > 0 can actually be both chromomagnetic (when g 2 s, * > 0, B 2 > E 2 , and − < 0) and chromoelectric (corresponding to g 2 s, * < 0, B 2 < E 2 , and + > 0). Indeed, the standard argument in favor of the positive definiteness of g 2 s is given in a classical YM theory, where in Minkowski space, such that g 2 s < 0 would lead to infinitely fast growth of the field A a i and action S cl = L cl d 4 x would not have a minimum. In the quantum case, however, g 2 s is a function of J, and can take negative values as long as the effective action S eff has a minimum for g 2 s < 0. Besides, in close vicinity of the ground-state solution β(g 2 s, * ) = 2, the corresponding solution of the RG Equation (8) takes a linear behaviour if both contributions coexist in the QCD vacuum, thus canceling each other beyond the confinement radius or after the QCD phase transition epoch in the cosmological history of the universe. The latter is an important example of conformal anomalies' cancelation in the classical limit of a YM theory without any fine tuning. In the deconfined (QGP) phase (i.e., at temperatures T T QCD ), the chromoelectric contribution + (T) should quickly vanish such that QCD − (T = T QCD ), providing a consistency with hadron physics phenomenology. This effect becomes plausible as long as + (T) is attributed to the ground state of hadronic degrees of freedom, which indeed becomes relevant only as soon as QGP is hadronised [82]. What is the possible role of the QCD vacuum in the generation of QGP? As was demonstrated in Ref. [86] within a semi-classical analysis, an interacting YM system of homogeneous gluon condensate and inhomogeneous wave modes evolves in real time in such a way that the amplitude of waves parametrically grows at the expense of decaying gluon condensate. The corresponding effect of the energy "swap" from the condensate to the gluon plasma waves is illustrated in Figure 2. Together with the growth of the plasma wave amplitudes, a vacuum average of their bi-linear products do not vanish and grow as well, inducing the positive-valued component + (T) of the ground-state density. This effect-if it holds in the full quantum formulation-can then be the basis for the QGP generation and reheating mechanism, both in heavy ion collisions and in Cosmology. A similar mechanism may be responsible for reheating of the cosmological plasma and particle production in the end of Cosmic Inflation due to decay of the dominant spatially-homogeneous chromoelectric condensate (inflaton) driving the inflationary epoch.  Can such parameteric growth of plasma modes due to a decay of the gluon condensate be studied in particle (e.g., heavy-ion) collisions? Can this effect be detectable, and if so, be used as a tool for the investigation of the dynamical evolution of the quantum ground state in QCD? Answers to these and other related questions are big unknowns, and very little has been done so far in this direction. An interesting insight into the problem of QCD ground state can be offered by the low-p ⊥ (<200 MeV) spectra of pions measured at the LHC, which show up to ∼30%-50% enhancement compared to the hydrodynamic models (see Refs. [87][88][89]). A possible interpretation could be found in the framework of a hypothesis about hadronization and freeze-out processes in chemical non-equilibrium [90]. Among the possible reasons for the non-equilibrium dynamics are the QGP supercooling [91,92] and gluon condensation [93] phenomena. A particularly interesting possibility has been proposed in Refs. [94][95][96], where it was shown that the Bose-Einstein pion condensate at the level of 5% can account for the missing low-p ⊥ charged pion yields coming from a coherent source in Pb + Pb collisions ( √ s = 2.76 TeV) at various centralities. Moreover, if there is such a condensate, there must be large fluctuations of pions, which should be seen starting from the fourth moment of the multiplicity distribution [97]. Further studies of the non-equilibrium QCD dynamics accounting for the ground state are certainly required from both theoretical and experimental standpoints.

Strongly Interacting Quark-Gluon Plasma
The surprising fact that the deconfined matter found at RHIC [14][15][16][17] does not behave as a gas of almost-free quarks and gluons, but as a strongly interacting liquid [18,26] was anticipated by only a few [98][99][100]. The fact that QGP close to the critical temperature T c is a strongly interacting system was used in Ref. [101,102] to exploit its analogy with strongly-coupled classical non-relativistic plasmas [62] in order to understand experimental observations and to interpret the lattice QCD results.
By definition, plasma is a state of matter in which charged particles interact via long-range (massless) gauge fields [26]. This distinguishes it from neutral gases, liquids, or solids in which the inter-particle interaction is of short range. So, plasmas themselves can be gases, liquids, or solids, depending on the value of the plasma parameter Γ, which is the ratio of interaction energy to kinetic energy of the particles forming the plasma [62].
A non-relativistic electromagnetic plasma is called strongly-coupled if the interaction energy (Coulomb energy) between the particles is larger than the thermal energy of the plasma particles; i.e., if the Coulomb coupling parameter Γ EM = q 2 /(aT 2 ) > 1, where q is the particle charge, a is the interparticle distance, and T is the plasma temperature (in the system of units whereh = c = k B = 1). Let us note that the strongly-coupled classical electromagnetic plasmas are not exotic objects at all [62]. For example, table salt (NaCl) can be considered as a crystalline plasma made of permanently charged ions Na + and Cl − [26]. At T ≈ 10 3 K (still too small to ionize non-valence electrons), it transforms into a molten salt, which is a liquid plasma with Γ ≈ 60. An estimate of the plasma parameter for QGP was considered in Ref. [101] where it was found that Γ QGP = 2Cα s /(aT 2 ), where-depending on the type of plasma-C = 4/3 or C = 3 is the Casimir invariant for quarks or gluons and a is the interparton distance a = 0.5 fm. For QGP at temperature only slightly above the critical de-confining temperature (i.e., T = 200 MeV), the corresponding coupling constant α s = 0.3-0.5, and Γ QGP = 1.5-6, the plasma can be considered as a strongly interacting one.
The strongly interacting plasmas Γ ≥ 1 are also a special case of strongly correlated systems where correlated behavior means a deviation from the trivial ideal gas behavior [104]. Prominent properties of all strongly correlated systems (see Figure 3) can be quantified by a few dimensionless parameters: the coupling parameter Γ, the degeneracy parameter χ = nλ 3 th , and the Brueckner parameter r s = a/a B , where n is the number density of the particles, λ th = 2π/(mT) is the thermal de Broglie wavelength, a is the average interparticle distance, m is the particle mass, and a B = 1/(me 2 ) is the Bohr radius.
Strongly interacting plasmas that can be studied in laboratory are ultracold atomic Fermi gases [105]; in particular, strongly-coupled 6 Li atoms [106,107]. A distinctive property of these plasmas is that-similarly to the strongly-coupled QGP (sQGP)-their shear viscosity-to-entropy density ratio η/s (see Sec. 4.3.1 for definition) characterizing how close the fluid is to a perfect liquid [108] is effectively negligible [26,103,106]. Cold atomic gases are produced in optical or magneto-optical traps containing typically 10 5 − −10 6 atoms [109]. The hydrodynamic behaviour is observed when the trapping potential is modified, or if the local density or energy density is modified using laser beams [106]. In this way, the scattering length a (and hence the interaction strength between the atoms) can be made almost infinite [26]. This is also the case of data point 6 Li (a = ∞) shown in Table 1, where the thermodynamical parameters for several other substances of interest are summarized. For H 2 O and 4 He, two points are displayed. First are the data at atmospheric pressure and temperatures just below the boiling point and the λ transition, respectively. These data points roughly correspond to the minimum of η/n at atmospheric pressure. Second are the data near the critical point, which roughly corresponds to the global minimum of η/s.

QCD at High Temperatures and Vanishing Chemical Potentials
The grand canonical partition function in SU(N c ) gauge theory (such as QCD) with N f fermion flavours having a common chemical potential µ reads where M is the Dirac operator, S g is the gauge part of the QCD action which depends on temperature T through boundary conditions. In the Hamiltonian formulation, whereN is the number operator, andĤ is the Hamiltonian. It is needless to mention that the analytic properties of the free energy F(T, µ) as a function of general complex µ are known to be useful for studying the phase structure of QCD on the lattice [110]. Let us consider QCD thermodynamics at high temperatures and zero chemical potentials-the region relevant for the LHC and also partly for RHIC-by recalling how the basic bulk thermodynamic observables can be obtained from the grand canonical partition function with vanishing quark chemical potentials, Z(T, V) ≡ Z(T, µ)| µ→0 [55]. The grand canonical potential, Ω(T, V), normalized in such a way that it vanishes at zero temperature, can be used to obtain the thermal part of the pressure (p) and energy density ( ) both vanishing at small temperature, by construction. Using these relations, one can express the difference between and 3p; i.e., the thermal contribution to the trace of the energy-momentum tensor Θ µµ (T) (also called the trace anomaly or the interaction measure) in terms of a derivative of the pressure with respect to temperature: In fact, it is Θ µµ (T) which is the basic thermodynamic quantity conveniently calculated on the lattice as the total derivative of ln Z with respect to the lattice spacing a [111]: Before moving to the results of lattice calculations, it is useful, for comparison, to recall a description of the strongly interacting matter below deconfinement temperature T c . Here, all thermodynamic quantities are expected to be well-described by the hadron resonance gas (HRG) model consisting of non-interacting hadrons as proposed by Hagedorn in the mid 1960s [41] (see also Ref. [42]). The trace anomaly in the HRG model is given by where K 1 ( km i T ) is a modified Bessel function, the different particle species of mass m i have degeneracy factors d i and η i = −1(+1) for bosons (fermions), and the sum runs over all known hadrons up to the resonance mass of m max = 2.5 GeV.
The results on the temperature dependence of the trace anomaly and suitably normalized pressure, energy density, and entropy density from lattice calculations, together with the HRG predictions, are shown on the left and right panels of Figure 4, respectively. The vertical band in the right panel marks the crossover region, T c = (154 ± 9) MeV. The horizontal line at 95π 2 /60≈15.6 corresponds to the ideal Stefan-Boltzmann (SB) gas limit for the energy density of relativistic massless gas consisting of N f = 3 quark flavours and gluons with N c = 3 colors having altogether g degrees of freedom: The fact that even at T ∝ 400 MeV, the pressure, energy density, and entropy density of the QGP are far from their ideal gas values indicates substantial remaining interactions among the quarks and gluons in the deconfined phase. It is interesting to compare (at least qualitatively) this behaviour with the gaseous two-component plasma of particles with charge q = ±ze. The pressure normalized to that of the ideal gas can be deduced from the standard textbook formula (e.g., Ref. [112]) which is valid for n T 3 /q 6 . The non-ideal behaviour of gaseous two-component plasma thus increases very rapidly with charge q of the plasma particles, but much more slowly with their density n, and decreases rather quickly with the plasma temperature T.  An important property of the QGP phase transition is the presence of a local minimum in the ratio of pressure to energy density p/ as a function [113,114]. The possible existence of this softest point in the QCD EoS is distinguishable by a very small sound velocity of the deconfined medium, and has thus been suggested as a signal of the first-order phase transition. This also becomes evident in second-order derivatives of the QCD partition function with respect to temperature. The speed of sound, c s , is related to the inverse of the specific heat, C V = d /dT, The quantity Td( /T 4 )/dT can be calculated directly from the trace anomaly and its derivative with respect to temperature, In Figure 5 (left panel), we show the speed of sound as a function of temperature. The softest point of the EoS predicted in Ref. [114] at T (145 − 150) MeV (i.e., at the minimum of the speed of sound) lies on the low temperature side of the crossover region. At this point, the speed of sound is only slightly below the corresponding HRG value. Furthermore, the value c 2 s 0.15 is roughly half-way between zero, the value expected at a second-order phase transition with a divergent specific heat, and the value for an ideal massless gas, c 2 s = 1/3 [111]. At the high temperature end, T ∼ 350 MeV, it reaches within 10% of the ideal gas value. The softest point of the EoS is of interest for the phenomenology of heavy ion collisions, as it characterizes the temperature and energy density range in which the expansion and cooling of matter slows down. The system spends a longer time in this temperature range, and one expects to observe characteristic signatures from this regime. The quantity c 2 s as a function of the energy density is shown in Figure 5 (left). At the softest point, the energy density is only slightly above that of normal nuclear matter, ρ 0 = 160 MeV/fm 3 . In the crossover region, T c = (154 ± 9) MeV, the energy density varies from 180 MeV/fm 3 at the lower edge to 500 MeV/fm 3 at the upper edge-slightly above the energy density inside the proton proton = 450 MeV/fm 3 .
Thus, the QCD crossover region starts at or close to the softest point of the EoS, and the entire crossover region corresponds to relatively small values of the energy density, (1.2-3.1) nuclear . This value is about a factor of four smaller than that of an ideal quark-gluon gas in this temperature range.

Testing the Properties of the Medium with Infinitely Heavy Static Test Charges
An important property of the QGP medium is the color screening: the range of interaction between heavy quarks becomes inversely proportional to the temperature. This effect also forms the basis of the most common description of dynamics of quarkonia (mesons consisting of heavy QQ), produced in heavy ion collisions-the potential between the heavy quarks cc or bb becomes screened by deconfined quarks and gluons, and the heavy quarks separate from each other, leading to a suppression of quarkonia yields [115].
On the lattice, this phenomenon is studied using (infinitely) heavy static test charges [116,117]. The color screening effect is estimated from the spatial correlation function G(r, T) of a static quark and anti-quark, which propagate in Euclidean time from τ = r = 0 to τ = r = 1/T, where T is the temperature. The free energy of static quark pair QQ is then calculated as the logarithm of the correlator F(r, T) = −T ln G(r, T) [116,118]. In the zero temperature limit, the singlet free energy coincides with the zero temperature potential calculated on the lattice [119]. However, as argued in Ref. [117], using the free energies instead of potentials is preferable, since the latter are not gauge invariant. On the other hand, the gauge-invariant static quark-antiquark pair free energy is a non-perturbatively well-defined quantity that carries information about the deconfinement properties of the QGP.
The heavy-quark free energies for different temperatures T of the medium from two different lattice calculations are presented in Figure 6 (left). The solid black line on the upper plot is a parameterisation of the zero temperature potential. One can see that with increasing temperature the free energy-and hence, the spatial correlations between Q andQ-get more and more diluted.  [116] and (bottom) continuum-extrapolated values of the static QQ free energy for different temperatures [117]. The solid black line on the right top plot is a parameterisation of the zero temperature potential; Right: The S-wave (upper) charmonium and (lower) bottomonium spectral functions calculated in potential models. Insets: correlators compared to lattice data. The dotted curves are the free spectral functions. Reproduced from Ref. [116].
The correlation functions in time variable are related to the spectral functions σ(ω, T) by While a stable QQ quarkonium state in the vacuum contributes a δ-function-like peak at the value of its mass m H to the spectral function, in the medium, it gives a quasi-particle-like smeared peak with the width being the thermal width. As the temperature increases, the width increases, and at sufficiently high temperatures, the contribution from the meson state in the spectral function becomes sufficiently broad so that it is no longer meaningful to speak of it as a well-defined state ( Figure 6, right). The effect is more prominent for the lighter mesons like charmonia, consisting of cc pairs, and much weaker for the bottomonia-bb mesons.

QCD at High Temperatures and Non-Zero Chemical Potentials
As direct lattice QCD calculations at non-zero µ B are not yet possible, one has to analyze the EoS using Taylor expansion in quark chemical potentials µ u , µ d , and µ s : [52,120] p where µ u , µ d , and µ s are related to the chemical potentials corresponding to the baryon number B, electric charge Q, and strangeness S of hadrons, as follows: The EoS at non-zero µ B,Q,S can thus be obtained from the coefficients χ BQS ijk of the Taylor expansion in hadronic chemical potentials expressed via χ uds ijk [52]. Here, we report the result [120] for µ Q = µ S = 0, which sufficiently illustrates the relative importance of higher-order corrections in different temperature and µ B regions. The Taylor series for the pressure is given by The leading-order correction to the pressure at non-vanishing µ B is proportional to the quadratic fluctuations of the net baryon number. The next-to-leading order corrections are proportional to the quartic fluctuations. In Figure 7 B ) correction rapidly loses importance relative to the leading O(µ 2 B ) term. Moreover, the results for the µ B -dependent contribution to the total pressure evaluated for different values of µ B /T [120] suggest that the EoS given by Equation (30) works well for all values of the chemical potential below µ B /T = 2, corresponding to the region of nuclear collisions at energies √ s NN ≥ 20 GeV. . The expansion coefficients of the pressure at the non-zero baryon chemical potential adopted from Ref. [120]. Left: The leading-order correction; Right: The relative contribution of the next-to-leading order lcorrection. BNL: Brookhaven National Laboratory.
Let us note that χ BQS ijk are interesting in their own right, as they are related to the fluctuations and correlations of conserved charges. The latter are sensitive to the underlying degrees of freedom, which could be hadronic or partonic, and so they are used as sensitive probes of deconfinement. While the off-diagonal expansion coefficients are related to correlations among conserved charges (e.g., , the diagonal ones describe their second-and higher order fluctuations

Heavy Ion Accelerators
The basic hopes and goals associated with investigations of very hot and dense nuclear matter in laboratory were first formulated in mid-seventies [121][122][123]. It was the experience with astrophysical objects like supernovae and neutron stars, and with thermonuclear ignition which led the authors to the idea that the nuclear matter shock compression [31] of about five-fold normal nuclear density should be accomplished in violent head-on collisions of heavy nuclei [11]. The goal was to find out the response of the nuclear medium under compression by pressure resisting that compression; i.e., to study the nuclear matter EoS. The original question was: is such a bulk nuclear matter EoS accessible within the dynamics of relativistic heavy ion collisions? [124,125]. The prospect of observing a phase transition in highly compressed nuclear matter [126] was lurking behind.
The interest in collisions of high-energy nuclei as a possible route to a new state of nuclear matter was substantially strengthened with the arrival of QCD as the microscopic theory of strong interactions. The particle physics community began to adapt existing high-energy proton accelerators to provide heavy-ion nuclear beams in the mid-1970s. The Berkeley Bevalac and JINR Synchrophasotron started to accelerate nuclei to kinetic energies from few hundreds of MeV to several GeV per nucleon [11,125]. By the mid-1980s, the first ultra-relativistic nuclear beams became available. Silicon and gold ions were accelerated to 10 GeV/nucleon at Brookhaven's Alternating Gradient Synchrotron (AGS) [10]. The first nuclear collisions took place at CERN in the early 1980s when alpha particles were accelerated to center-of-mass energy per nucleon-nucleon pair √ s NN = 64 GeV at the ISR collider.
The new era of research began at CERN in fall 1986, when oxygen-and later on in the summer of 1990, sulphur ions-were injected into the SPS and accelerated up to an energy of 200 GeV/nucleon ( √ s NN = 19.6 GeV) [10][11][12]. However, the genuine heavy ion program only started in 1994, after the CERN accelerator complex was upgraded with a new lead ion source which was linked to pre-existing interconnected accelerators, the Proton Synchrotron (PS) and the SPS. Seven large experiments involved (NA44, NA45/CERES, NA49, NA50, NA52, WA97/NA57, and WA98) studied different aspects of Pb + Pb and Pb+Au collisions at √ s NN = 17.3 GeV and √ s NN = 8.6 GeV [11,12].
In the meantime, at the Brookhaven National Laboratory (BNL), the Relativistic Heavy Ion Collider (RHIC) [127] rose up from the ashes of the ISABELLE/CBApp collider project abandoned in 1983 by particle physicists. In 1984, the first proposal for a dedicated nucleus-nucleus machine accelerating gold nuclei up to √ s NN = 200 GeV was submitted. Funding to proceed with the construction was received in 1991, and on 12 June 2000, the first Au + Au collisions at √ s NN = 130 GeV were recorded by the BRAHMS, PHENIX, PHOBOS and STAR experiments [14][15][16][17]. The idea of the Large Hadron Collider (LHC) [128] dates even further back, to the early 1980s. Although CERN's Large Electron Positron Collider (LEP, which ran from 1989 to 2000) was not yet built, scientists considered re-using the 27-kilometer LEP ring for an even more powerful pp machine running at the highest possible collision energies ( √ s = 14 TeV) and intensities. The ion option ( √ s NN = 5.4 TeV per nucleon-nucleon pair for Pb + Pb collisions) was considered since the beginning.
The LHC was approved in December 1994, and its official inauguration took place on 21 October 2008. The first proton-proton collisions occurred on 23 November 2009, and the first Pb + Pb collisions on 7 November 2010. The ALICE, ATLAS, and CMS experiments are currently involved in the heavy-ion program at the LHC [25,29].
Many years ago, American accelerator physicist M. Stanley Livingston noted that advances in accelerator technology increase the energy records achieved by new machines by a factor of ten every six years. This trend is illustrated on the left panel of Figure 8, summarising the worldwide advances in high-energy accelerators in the period of 1960-2008. One can see that the increase in the energy is even faster for the ion accelerators than for the proton accelerators. However, even in the most central nucleus-nucleus collisions, not all of the energy could be converted into a thermalised form of energy needed for the phase transition from hadronic into a QGP state of matter to occur. The experimentally accessible quantity measuring this transformation is the density of transverse energy E T per unit of pseudorapidity η. The latter can be used to estimate the initial energy density using the Bjorken formula [129]: where S ⊥ is the transverse overlap area of the nuclei and τ is the time scale for the evolution of the initial non-equilibrium state of matter into a (locally) thermalized system. The dependence of BJ τ on √ s NN for most central Au + Au and Pb + Pb collisions is presented on the right panel of Figure 8.
Even for a rather pessimistic value of the eqilibration time τ = 1 fm/c [17], the achieved energy density increases from 1.4 to 14 GeV/fm 3 . It is thus not only much higher than than the normal nuclear density, but 3-30 times bigger than the energy density inside the nucleon N = 0.45 GeV/fm 3 , and is definitely higher than the 1 GeV/fm 3 scale required for the QCD deconfining transition from the lattice calculations that were discussed in Section 3.3. Available center-of-mass energy √ s-2m N versus time for (anti)proton (blue triangles) and ion (magenta circles) accelerators, adapted from Ref. [71]; Right: The Bjorken estimate of the initial energy density ε BJ (Equation 32) multiplied by τ calculated from the data on transverse energy distributions in 5% most central Au + Au [130] and Pb + Pb [131] collisions as a function of c.m.s. energy per one nucleon-nucleon pair √ s NN . The red line corresponds to a power law fit. From Ref. [132].

Heavy Ion Collisions as a Source of Strong Electromagnetic Fields
The important quantity determining the lifetime of heavy ion beams inside the accelerator tube (and hence their potential to produce an adequate number of nuclear collisions) is the loss of ions in the bunch. Its decay rate λ T is given by the formula: in which k is the number of bunches, N is the number of particles per bunch, N IR is the number of interaction regions, σ T is the total cross-section, and L is the available luminosity: which-in addition to k and N-also depends on the revolution frequency f , the Lorentz factor of the beam γ, the emittance , the beta function at the collision point β * , and the geometric luminosity reduction factor F due to the crossing angle at the interaction point [128]. Since the event rate for a certain process is given by its cross-section times the luminosity N evt = Lσ evt , the possibility of studying rare phenomena depends on the maximum luminosity accessible. At the same time, the rate of background processes (which in general have large cross-sections) will increase L, at some point reaching the maximum event rate that the experiment can handle. Secondary beams created by these background processes can limit the collider heavy ion luminosity, since they have a different charge-to-mass ratio than the primary beam, and can be lost in cryogenically cooled magnets.
Quite surprisingly, at RHIC and the LHC, these background processes are not part of the strong nuclear interaction cross-section σ R determined primarily by the nuclear geometry , but are solely accounted for by the coherent action of all electric charges in colliding nuclei. The cross-section of electromagnetic processes-primarily due to the creation of e + e − pairs with subsequent e − atomic shell capture and electro-magnetic dissociation-becomes important for ions with Z > 30. It is as large as hundreds of barns; i.e., about 30 (60) times larger than σ R for Au + Au collisions at RHIC (Pb + Pb collisions at the LHC) [133].
A classical description of the electromagnetic action of a fast-moving charged particle on another one based on the equivalence between the perturbative action of its field and the flux of electromagnetic radiation dates back to Fermi [134], Weizsäcker [135], and Williams [136]. This equivalence is true as far as the effects caused by different spectral components add up incoherently (i.e., a perturbation caused by the fields is small enough). The solution for the time-dependent electromagnetic fields mutually seen by the two incident ions can be found, for example, in the textbook on "Classical Electrodynamics" [137]. The longitudinal ( ) and transversal (⊥) field components induced by a heavy ion I passing a target I I at distance b and with velocity β are given by the following formulas: Let us note that for γ 1, these fields act on a very short time scale of order ∆t ∝ b/γ. During this time, fields E ⊥ (t) and B ⊥ (t) are equivalent to a linearly polarized pulse of radiation incident on a target in the beam direction. Thus, according to the equivalent photon method, the strong and rapidly time-varying field of the point charge Z I is seen by a passing charge as a flux of virtual (nearly real) photons with intensity where E(ω), B(ω), and E ⊥ (ω) are the Fourier components of the fields E, B, and E ⊥ . The energy spectrum of these photons falls as ∝ 1/E γ , up to a maximum energy E max The interaction between the colliding nuclei becomes dominantly electromagnetic for impact parameters b exceeding the size of the radii of colliding nuclei b > b min = R I + R I I = √ σ R /π. Interactions between ultra-relativistic nuclei taking place at b > b min are called the ultra-peripheral collisions. By taking advantage of the photon fields carried by relativistic nuclei, they are used to study photoproduction and two-photon physics at hadron colliders. This field of ultra-relativistic heavy ion collisions is sometimes called "non-QGP" physics, and is thus outside the scope of this article. We refer the interested reader to reviews [138][139][140][141] where more detailed information on these aspects can be found.

Quark-Gluon Plasma in a Strong Magnetic Field
More important from the point of view of QGP physics are the strong magnetic fields accompanying ultra-relativistic heavy ion collisions [142]. Consider the collision of two identical nuclei of radius R with electric charge Ze, and use the Biot-Savart law to estimate the magnitude of the perpendicular magnetic field they create in the center-of-mass frame Here, γ = √ s NN /2m N is the Lorentz factor. At RHIC, heavy ions are collided at √ s NN = 200 GeV per nucleon, hence γ = 100. Using Z = 79 for Gold and b ∼ R A ≈ 7 fm, we estimate eB ≈ m 2 π ∼ 10 18 G. At the LHC at √ s NN = 5.02 TeV and Z = 82, this value is even 30 times larger. To appreciate how strong this field is, compare it with the magnetic field of a neutron star 10 10 -10 13 G [69], or that of its slowly rotating magnetic variant, the magnetar, 10 15 G [143]. It is very likely the strongest magnetic field in nature, though existing only for a minute period of time.
Calculation with the realistic distribution of protons in a nucleus shows that magnetic field rapidly decreases as a power of time, and after the first 3 fm/c drops from its maximal value (37) by more than three orders of magnitude [144]. However, different estimates to be discussed in the next sections indicate that a strongly interacting thermalised medium is formed as early as 0.5 fm/c. Therefore, a more realistic calculation going beyond the above field in the vacuum calculation has to include the response of the medium determined by its electrical conductivity. It has been found by lattice calculations [145] that the gluon contribution to the electrical conductivity of static quark-gluon plasma is This result was confirmed and further extended by more elaborate lattice simulations with 2 + 1 dynamical flavors for temperatures T = (120-350) MeV [146,147]. The calculations have shown that σT already starts to deviate from zero for T < T c (i.e., in the confined phase), and increases towards the QGP value (38) and further on. The non-zero electrical conductivity in the QGP and (probably also) in the hadronic phase when taken at its face value would inevitably lead to a substantially prolonged lifetime of the magnetic field inside the medium, and might thus even influence the hadron decay widths [148].
A plethora of novel non-dissipative transport phenomena related to the interplay of quantum anomalies with the magnetic field and vorticity in systems with chiral fermions, including the QGP, is reviewed in Ref. [149]. The most direct effect of magnetic field B on the QGP is the induction of electric currents carried by the charged quarks and antiquarks in the plasma, and later, by the charged hadrons. In Ref. [150], it was suggested that it may leave its imprint on the azimuthal distributions and correlations of the produced charged hadrons. Charged particles moving along the magnetic field direction y are not influenced by the magnetic Lorentz force, while those moving in the xz-plane (i.e., in the reaction plane) are affected the most. The result is an azimuthally anisotropic flow of expanding plasma in the xy-plane, even when the initial plasma geometry is completely spherically symmetric.
Another effect is related to the chiral symmetry restoration. In such a state, within a localized region of space-time, gluon fields can generate nontrivial topological charge configurations that lead to parity violation in strong interactions [126]. In ultra-relativistic heavy-ion collisions, interactions between quarks and these gluonic states can lead to an imbalance in left-and right-handed quarks which violates parity symmetry [151]. The presence of a strong magnetic field induced by the spectator protons transforms this chirality imbalance into an electromagnetic current perpendicular to the reaction plane. This interesting phenomenon stemming from the interplay of chirality, magnetic field, and the chiral anomaly is called the Chiral Magnetic Effect (CME) [149].
Several manifestations of the phenomena related to strong magnetic fields produced in Au + Au or Pb + Pb collisions have been reported by RHIC [152][153][154] and the LHC [155] experiments. Some doubts on the prevailing interpretation were cast by the recent observation of charge-dependent azimuthal correlations also in p + Pb collisions at the LHC [156]. Moreover, in the presence of elliptic flow (for its definition, see Ref. 5.1), practically all conventional two-particle correlations like the local charge conservation [157] may contribute to the reaction-plane-dependent correlation function used to quantify the CME [158]. Obviously, more investigations are needed. The program of varying the magnetic field by a controlled amount while keeping all else fixed by using nuclear isobars (pairs of nuclei with the same mass number A but different charge Z) is now under consideration at RHIC. The most attractive isobars are Zr+Zr and Ru+Ru, or some other combinations having charge differences of four like Sn 124 /Xe 124 , Te 130 /Ba 130 , and Xe 136 /Ce 136 [158].
An interesting suggestion addressing the simultaneous effects of the huge vorticity of nearly-perfect fluid and the strong magnetic field generated in non-central heavy-ion collisions was made in Ref. [159]. The authors suggest the measurement of a global polarization of the final hadrons in order to estimate the thermal vorticity due to the large orbital momentum of colliding nuclei as well as the electromagnetic field developed in the plasma stage of the collision.

Transport Models
One of the main tasks of the theory is to link experimental observables to different phases and manifestations of the QCD matter. To achieve this goal, a detailed understanding of the dynamics of heavy-ion reactions is essential. This is facilitated by transport theory, which helps to interpret or predict the quantitative features of heavy-ion reactions. It is particularly well-suited for a non-equilibrium situation, finite size effects, non-homogeneity, N-body phase space, particle/resonance production and freeze-out, as well as for collective dynamics. Microscopic [160][161][162][163][164][165][166][167][168], macroscopic (hydrodynamical) [169][170][171][172], or hybrid [173][174][175] transport models attempt to describe the full time evolution from the initial state of a heavy-ion reaction up to the freeze-out of all initial and produced particles after the reaction. This is illustrated in Figure 9, where a comparison of the data from heavy-ion collisions to the microscopic and hydrodynamical models is presented. Data (red circles) are compared to two microscopic models: EPOS LHC [168], DPMJET [167], to the hydrodynamics calculation [170], and to the Blast-Wave fit with formula Equation (53). Reproduced from Ref. [176]; Right: Comparison of the experimental p T -spectra of π + , K + , and p + from Au + Au collisions at (top) √ s NN = 200 GeV and (bottom) (for the case of π + ) of their centrality dependence [177] with the hydrodynamical model calculations. Reproduced from Ref. [178].
The hadronic cascade models (some with mean-field interactions) have succeeded in reproducing the gross and many-detailed features of the nuclear reactions measured at SIS, AGS, and SPS [160,161,163,164]. They have become indispensable for experimentalists who wish to identify interesting features in their data or to make predictions to plan new experiments. The main strength of the models based on superposition of pp collisions, relativistic geometry, and final-state hadronic rescattering is not that it gives a precise agreement with experiment for individual observables in particular kinematic regions, but in its ability to give an overall qualitative description of a range of observables in a wide kinematic region. The price to be paid for this simplicity is to assume that either hadrons or hadron-like objects can exist at the earliest stage of the heavy-ion collision just after the two nuclei pass through each other; i.e., that the hadronization time in the frame of the particle is short and insensitive to the environment in which it finds itself. The general success of these models at lower energies can nonetheless easily lead to misconceptions at higher energies. The main concern is the relevance of these models at high particle densities, which are so characteristic for collisions of heavy systems. Here, all the models based on hadronic dynamics are fundamentally inconsistent [179]. Studying the size of the fraction of the energy contained in known hadrons and that temporarily stored in more elusive objects (such as pre-hadronized strings), it was found [163] that up to a time of 8 fm/c, most of the energy density resides in strings and other high-mass continuum states that have not fully decayed. The physical properties of these objects are poorly known, even when they occur in isolation [180]-not to speak about their interactions (or even their existence) in a dense environment. The application of these models to the early phase of collision of two ultra-relativistic heavy nuclei is therefore ill-founded [179].
A complementarity between the microscopic and macroscopic descriptions becomes obvious for the case of strongly-interacting plasmas. The fact that neither Boltzmann equation nor cascades can be used for liquids stems from the fact that particles are strongly correlated with several neighbours at all times. The very idea of "scattering" and cross-section involves particles coming from and going to infinity: it is appropriate for dilute gases, but not for condensed matter where interparticle distances do not exceed the range of the forces at any time [26].
The idea to use the laws of ideal hydrodynamics to describe the expansion of the strongly interacting matter formed in high energy hadronic collisions was first formulated by Landau in 1953 [181]. Later on, Bjorken [129] discovered a simple scaling solution that provides a natural starting point for more elaborate solutions in the ultra-relativistic domain. The phenomenological success of the Landau model was a big challenge for high energy physics for decades [182]. First, because hydrodynamics is a classical theory, and second because it assumes local equilibrium. Both of these assumptions imply many degrees of freedom, and it is by no means clear that the highly excited but still small systems produced in violent nuclear collisions satisfy the criteria justifying treatment in terms of a macroscopic theory [183]. Therefore, the Landau model (and other statistical models of strong interactions) were considered up to the mid-seventies as exotic approaches, lying outside of mainstream physics [182]. Then, the authors of Refs. [121,122] realized that the exploitation of hydrodynamics in an interpretation of data is the only chance of proving the existence of a new state of matter in laboratory. This is a trivial corollary of the well-known fact that a state of matter is defined by its EoS, and there is no other way to get information about the EoS than by using hydrodynamics [182][183][184].

Elements of Relativistic Hydrodynamics
Let us now briefly recall some of the basic results from relativistic hydrodynamics [9,20,22,[185][186][187] upon which the contemporary models are based. The basic hydrodynamical equations describe the energy-momentum and the current conservation where j µ i , i = B, S, Q is the conserved current. Both quantities can be decomposed into time-like and space-like components using natural projection operators, the local flow four-velocity u µ , and the second-rank tensor perpendicular to it ∆ µν = g µν − u µ u ν : where = u µ T µν u ν is the energy density, p = p s + Π = − 1 3 ∆ µν T µν is the hydrostatic + bulk pressure, i is the charge current, and π µν = T µν is the shear stress tensor. The angular brackets in the definition of the shear stress tensor π µν stand for the following operation: To further simplify our discussion, we restrict ourselves in the following to only the one conserved charge, and denote the corresponding baryon current as j µ = j µ B . The various terms appearing in the decompositions (40) and (41) can then be grouped into ideal and dissipative parts Neglecting the dissipative parts, the energy-momentum conservation and the current conservation (39) define ideal hydrodynamics. In this case (and for a single conserved charge), a solution of the hydrodynamical Equation (39) for a given initial condition describes the space-time evolution of the six variables-three state variables (x), p(x), n(x), and three space components of the flow velocity u µ . However, since (39) constitute only five independent equations, the sixth equation relating p and -the EoS-has to be added by hand to solve them. For this, one can either use the relativistic non-interacting massless gas EoS or its generalization to the case of a non-zero interacting measure Θ µµ (T) = − 3p. In addition to many different phenomenological parameterizations of Θ µµ , one can exploit the relation (20) to obtain the EoS directly from the lattice QCD simulations. Examples of this approach were given in Section 3.3.1 (see particularly Equations 24 and 25), and are illustrated in Figure 5.
Two definitions of flow can be found in the literature [9,20,185]; one related to the flow of energy (Landau) [181] reads while the other relating to the flow of conserved charge (Eckart) [188] as follows: Let us note that W µ = 0 (V µ = 0) in the Landau (Eckart) frame. In the case of vanishing dissipative currents, both definitions represent a common flow. In ultra-relativistic heavy-ion collisions, the Landau definition is more suitable when describing the evolution of matter in the region with a small or zero baryon number deposition (i.e., when j = j B = 0) like the mid-rapidity region at the LHC and at the top RHIC energy; see Figure 1. In this case, the heat conduction effects can be neglected.
In order to solve the hydrodynamic equations with the dissipative terms, it is customary to introduce the following two phenomenological definitions (so-called constitutive equations) for the shear stress tensor π µν and the bulk pressure Π [185], where the coefficients η and ζ are called the shear viscosity and bulk viscosity, respectively. For the boost-invariant Bjorken flow [129] which is also called the one-dimensional Hubble flow, since velocity in the z direction, v z , is proportional to z, where τ = √ t 2 − z 2 is the proper time, one obtains the following equation of motion [185]: Neglecting the last two terms in Equation (49), one obtains the famous Bjorken solution of ideal hydrodynamics [129]. The last two terms on the right-hand side in Equation (49) describe a compression of the energy density due to viscous corrections. The first one is due to the shear viscosity in compressible fluids, while the second comes from the bulk viscosity. Two dimensionless coefficients in the viscous correction, η/s and ζ/s, reflect the intrinsic properties of the fluids; see Table 1 and Figure 15 (left panel). The value η/s = 1/4π has been obtained in the framework of N = 4 SUSY Yang-Mills theory [108]. The conformal nature of this theory gives ζ/s = 0 automatically. Moreover, η/s = O(0.1 − 1) for gluonic matter is obtained from the lattice calculations of pure SU(3) gauge theory [189], while the bulk viscosity ζ has a prominent peak around T c resulting from the trace anomaly of QCD [190].
Hydrodynamics provides an effective description of a system that is in local thermal equilibrium, and can be derived from the underlying kinetic description through Taylor expansion of the entropy four-current S µ = su µ in gradients of the local thermodynamic variables. In zeroth order in gradients, one obtains ideal fluid dynamics. Then, the higher orders describe dissipative effects due to irreversible thermodynamic processes such as the frictional energy dissipation between the fluid elements that are in relative motion, or heat exchange of the fluid element with its surroundings on its way to approach thermal equilibrium with the whole fluid. The relativistic Navier-Stokes description given above in Equation (47)-which accounts only for terms that are linear in velocity gradient-leads, unfortunately, to severe problems. Qualitatively, it can be seen from Equation (47) observing what happens if one of the thermodynamic forces ∇ µ u ν or ∇ µ u µ is suddenly switched off/on [191]. In this case, the corresponding thermodynamic flux π µν or Π which is a purely local function of the velocity gradient also instantaneously vanishes/appears. The linear proportionality between dissipative fluxes and forces causes an instantaneous (acausal) influence on the dissipative currents, leading to numerical instabilities [192].
Any numerical implementation of relativistic dissipative fluid dynamics thus requires the inclusion of terms that are second order in gradients [186]. The resulting equations for the dissipative fluxes π µν and Π are relaxation-type equations [187] with microscopic relaxation times τ π ≡ 2ηα and τ Π ≡ ζ β, which encode the time delay between the appearance of thermodynamic gradients that drive the system out of local equilibrium and the associated build-up of dissipative flows in response to these gradients, thereby restoring causality [187,193]. Since the relaxation times must be positive, the Taylor expansion coefficients α and β must all be larger than zero. Accounting for non-zero relaxation times at all stages of the evolution constrains departures from local equilibrium, thereby both stabilizing the theory and improving its quantitative precision.

Blast Wave Parametrization
Interpretation of the results of hydrodynamical calculations or of the experimental data in terms of the collective flow of matter [87,194] is greatly facilitated by the use of the analytical so-called Blast Wave (BW) parametrization [195][196][197]. Within the boost-invariant scenario of Bjorken [129], and for the full azimuthal symmetry which is valid in central collisions of two nuclei, the velocity field of expanding matter is given by where ρ = tanh −1 β T and η are transverse and longitudinal rapidities, respectively, and e r is the unit vector in the transverse plane. The transverse velocity distribution β T (r) of the thermalized matter in the region 0 ≤ r ≤ R is described by a self-similar profile where β s is the surface velocity, and parameter k is usually given the value k = 2 to resemble the solutions of hydrodynamics [196]. The spectrum of locally thermalized matter is constructed as a superposition of the individual thermal components [198]: where σ is the hypersurface defining a borderline between the hydrodynamical behaviour and free-streaming particles (the so-called freeze-out hypersurface), and T kin is the temperature of the kinetic freeze-out. Boosting each component with the transverse rapidity ρ = tanh −1 β T , one obtains the transverse momentum spectra of particles from the collective radial flow of expanding matter: where m T = (m 2 + p T 2 ), I 0 (x), and K 1 (x) are the Bessel functions. Formulas for the case of non-central collisions when the transverse shape (51) is controlled not by one but two parameters, R x and R y , can be found in Ref. [197]. In full generality, there are eight parameters describing the blast wave parametrization: T, ρ 0 , ρ 2 , R x , R y , a s , τ 0 , and ∆τ. Here, T is the temperature, ρ 0 and ρ 2 describe the strength of the zero-and second-order oscillation of the transverse rapidity, the parameter a s corresponds to a surface diffuseness of the emission source, and τ 0 and ∆τ are the mean and width of a Gaussian longitudinal proper time τ = √ t 2 − z 2 freeze-out distribution. In Figures 9 and 10, the examples of the BW fit analysis are presented. As can be seen from Figure 10, the kinetic freeze-out temperature T kin that determines the shape of the p T -spectra of particles is strongly anti-correlated with the radial flow velocity β : the higher T kin is, the lower β is, and vice versa. Nevertheless, the radial flow reveals itself as a shoulder structure at small transverse momenta in the p T -spectra of Λs, protons, and kaons; see Figure 9. For the pions, there is almost no sensitivity to distinguish between the two cases-a reduction of temperature is almost compensated by the radial flow.

Initial State Description of Nuclear and Hadronic Interactions
An indispensable part of the full description of the experimental data from heavy-ion collisions comes not only from the understanding of its dynamics starting from the moment of thermalization, but also at earlier times. In particular, the question of where the observed (local) thermalization of deconfined matter comes from is still quite open [201][202][203][204]. The importance of event-by-event initial state fluctuations on anisotropic collective flow and other final state observables is also worth mentioning [205][206][207][208]. Since these topics currently remain a significant source of uncertainty in predicting the final state observables, in the next two paragraphs we will provide two alternative ways of describing the initial state of the collision.
In high-energy nucleus-nucleus (A + B) interactions, the de Broglie wavelength of the nucleons (N) of the incoming nucleus is much smaller than the inter-nucleon distances inside the partner nucleus. To each incoming nucleon, the positions of the nucleons within the partner nucleus appear to be frozen. After a single elementary NN (elastic or inelastic) collision, both participating nucleons acquire a transverse momentum, which in the majority of cases is very small compared to their longitudinal one, and so the longitudinal momenta before and after the collision are very close to each other p z ≈ p z . High incident energies and small scattering angles mean that the scattering is dominated by a large orbital momentum , and so it is convenient to replace the partial-wave expansion of the scattering amplitude by an impact parameter b = (1 + )/p representation. The A + B collision can thus be described using a semi-classical approach due to Glauber [209][210][211][212], which treats the nuclear collision as multiple NN interactions [213,214]. The nucleons which have suffered at least one NN collision are called the participants, and those who have avoided it are called the spectators (see Figure 11). The total number of spectators and participants thus adds up to N spec +N part = A + B. On the other hand, the total number of collisions suffered by all participants fulfils inequality N coll ≥ N part /2.

Glauber Model
Using the nuclear mass number density dzd 2 s ρ A,B (z, s) = A,B and the inelastic NN cross-section σ in NN , we can express N part and N coll analytically [9,22,28,215]: where T A (b) and T AB (b) are the nuclear thickness and the nuclear overlap functions, respectively: The Glauber model calculations are also often carried out via Monte Carlo [215][216][217]. Nucleons inside the colliding nuclei are distributed randomly according to a nuclear density profile. At a given impact parameter, b, the impact parameter s of all pairs of nucleons is calculated. Interaction occurs when πs 2 < σ in NN ; see Figure 11. The calculated N part (b) and N coll (b) are then used to make a contact with the measured bulk observables, like the multiplicity of charged particles N ch measured by the tracking detectors at midrapidity (Figure 12, left panel), or the energy left by the spectator nucleons in the Zero Degree Calorimeter E ZDC (Figure 12, right panel). The left panel also explains how different bins in multiplicity of charged particles N ch can be transformed into the bins in collision centrality. The left panel of Figure 13 illustrates that even for quite different nuclear systems (Au + Au, Cu + Au, Cu + Cu), the number of participants N part selects the collisions with almost the same energy density ε BJ . The right panel of Figure 13 shows that centrality dependence of the yields of EW-interacting particles (like direct photons, W ± , or Z-bosons) is completely determined by the N coll .  A generalisation of the Glauber model beyond its original non-relativistic potential description of the scattering process was first formulated by Gribov [220,221], who used the effective field theory to describe multiple interactions proceeding via the Pomeron exchange. The interference terms appearing naturally in this treatment of multiple scattering automatically assure the unitarity of the theory [222]. With the advent of the QCD, the parton-based multiple scattering models became popular [223][224][225][226][227]. Such models allow the treatment of A + B and NN collisions on more equal footing. Their participants could thus be not only nucleons, but also the pQCD partons, the valence quarks [223,228], or any effective sub-nucleon degrees of freedom [229].

QCD Scattering in the Dipole Picture: Initial-State Energy Loss and Shadowing
The color dipole formalism [230,231] provides phenomenological means for the universal treatment of inclusive and diffractive p + p, p + A, and A+A collisions at high energies. In reactions involving a hard scale such as the Deep Inelastic Scattering (DIS), Drell-Yan (DY), heavy quark production, etc., this formalism effectively accounts for the higher-order QCD corrections and enables the quantification of the initial state interaction, saturation, gluon shadowing, and nuclear coherence effects [230,[232][233][234][235][236][237][238]. At small x, the eigenstates of interaction correspond to color dipoles with a definite transverse separation, the main ingredients of the dipole picture. Provided that the initial projectile state can be considered as a collection of dipoles, any production process in the target rest frame is then viewed as a superposition of universal dipole-target partial amplitudes f el (b, r; x) at different dipole separations r, impact parameters b, and Bjorken x. Integrating the amplitude over b, one arrives at the universal dipole cross-section σq q = σq q (r, x), which is determined phenomenologically. For example, the DIS process in the target rest frame is viewed as a scattering of the "frozen" qq dipole of size r∼1/Q, originating as a fluctuation of the virtual photon γ * → qq off the target nucleon. The DY pair production is considered as a bremsstrahlung of massive γ * (and Z 0 boson) off a projectile quark before after and the quark scatters off the target [233]. The projectile qq dipole probes the dense gluonic field in the target at high energies when nonlinear (e.g., saturation) effects due to multiple scatterings can become relevant. The latter, however, cannot be reliably predicted due to unknown soft and non-perturbative QCD dynamics.
In the limit of point-like color-neutral dipole r → 0, the corresponding dipole cross-section vanishes quadratically σq q ∝ r 2 due to color screening, which is known as the color transparency [231,239,240]. A naive saturated shape of the dipole cross-section known as the Golec-Biernat-Wüsthoff ansatz originally fitted to the DIS data [241] with the energy-dependent saturation scale Q s denoting a characteristic boundary between the linear and non-linear QCD regimes. Such a simple ansatz has proven to successfully capture the major features of a vast scope of observables in p + p and p + A, while the DGLAP evolution of the gluon density in the target affects the saturation scale at large scales µ 2 such that a QCD-corrected ansatz Q 2 s ∝ α s (µ 2 ) xg(x, µ 2 ) is often used in practice [242][243][244][245].
The saturation effects are amplified in nuclear collisions, such that the nuclear saturation scale Q A s is expected to be enhanced w.r.t. the nucleon one by a factor of A 1/3 . Then, the dipole-nucleus cross-section σ A qq can be represented in terms of the forward partial dipole-nucleus amplitude f A el analogically to Equation (57). Evolution of this amplitude in rapidity Y = ln(1/x) is predicted, for example, in the Color Glass Condensate formalism (see below) ; while no first-principle calculation capable of predicting the b-dependence exists, only phenomenological fits inspired by saturation physics are implemented (aside from what is mentioned above, see Refs. [246][247][248][249][250][251][252]). In heavy-ion collisions at high energies, due to the absence of final state interactions, the DY pair production process on nuclear targets in the dipole picture is often considered to be an excellent probe accessing the impact parameter dependence of the Initial State Interaction (ISI) effects as well as nuclear shadowing and nuclear broadening-the crucial information that cannot be derived within the parton model (e.g., [253][254][255] and references therein). The ISI effects emerging due to multiple rescattering of projectile partons in a medium prior to a hard scattering are only relevant close to kinematic boundaries due to energy conservation significantly suppressing the nuclear cross-sections [256,257]; e.g., when the Feynman variable x L ≡ x F = 2p L / √ s → 1 and x T = 2p T / √ s → 1. In order to account for the ISI-induced energy loss in the Glauber approximation, one sums up over initial state interactions in a p + A collision at a given impact parameter b, leading to a nuclear ISI-improved projectile quark PDF [256,257]: where T A (b) is the nuclear thickness function defined in Equation (56), C v is fixed by the Gottfried sum rule, and σ eff = 20 mb is the effective cross-section controlling the multiple ISI. The ISI-induced energy loss can induce a significant suppression at large M ll , p T , and forward rapidities. In fact, it has been noticed earlier in π 0 production in central dAu collisions [258], and in direct photon production in central AuAu collisions [259,260] at midrapidity and large transverse momenta (where no shadowing is expected) by the PHENIX Collaboration.
In the course of propagation, a projectile parton experiences multiple rescatterings in the color medium of a target nucleus. These cause a notable suppression of the gluon radiation and energy loss by the parton-a QCD analog of the Landau-Pomeranchuk-Migdal (LPM) effect in QED [261][262][263] (for more details, see Ref. [264] and references therein). In heavy ion collisions, the LPM effect was studied in Refs. [265][266][267].
In the nuclear DY reaction, for example, the shadowing is controlled by the coherence length corresponding to the mean lifetime of γ * -quark fluctuations given in terms of the quark mass m f , nucleon mass m N , dilepton invariant mass M ll , its transverse momentum p T , and light-cone momentum fractions of dilepton α and the target gluon In Ref. [254], it was noticed that the long coherence length (LCL) compared to the nuclear radius R A (namely, l c R A ) is practically useful in kinematic domains of RHIC and the LHC experiments. The quark shadowing in the LCL limit is automatically accounted for by scattering of the leading Fock components in the dipole formula which represents a process-dependent convolution of the light-cone wave function for a given Fock state with a superposition of partial dipole amplitudes. In p + A collisions, the latter are given by the eikonal Glauber-Gribov form [231]: which resums high-energy elastic scatterings of the dipole in a nucleus. Such an eikonalisation is justified in the LCL regime where the transverse separation of partons in the projectile multiparton Fock state is "frozen" in the course of propagation of the dipoles through the nuclear matter such that it becomes an eigenvalue of the scattering operator. At large scales and at small collision energies ( √ s 200 GeV), one recovers l c 1 fm such that the LCL approximation fails and an alternative description is necessary. In Refs. [268][269][270], a universal path-integral approach to multiple dipole interactions with a nuclear target which can be used for any value of l c has been developed. It is also known as the Green function approach, and consistently incorporates the color transparency and quantum coherence effects responsible for nuclear shadowing. Recently, the Green functions technique accounting for the nuclear attenuation of the colorless dipole in the medium has been employed in Ref. [255] in analysis of the nuclear coherence effects in DY reaction in p + A collisions. In Ref. [271], it was demonstrated that the Green function technique is equivalent to the diagrammatic formulation developed by the Baier-Dokshitzer-Mueller-Peigne-Schiff Collaboration (e.g., Refs. [272][273][274]). In Ref. [275], it was understood that the basic reason for the observed suppression of high-p T hadron production in heavy ion collisions is the attenuation of early produced colorless dipoles ("pre-hadrons") propagating through a dense absorptive matter, but is not an energy loss because of a much shorter production length than was typically assumed.
In the case of DY, the lowest Fock state |qγ * scattering is not enough in the LHC kinematics domain, where the higher Fock components containing gluons (e.g., |qγ * g , |qγ * gg etc.) also become noticeable, causing an additional suppression known as gluon shadowing (GS). These components are heavier than |qγ * , and thus have a shorter coherence length. In the high energy limit, the additional suppression factor R G induced by the GS corrections has been computed in Ref. [276][277][278] using the Green function technique, going beyond the LCL limit as in terms of the inelastic correction ∆σ to the total nuclear DY cross-section σ γ * A tot , related to the formation of a next-to-leading |qq g Fock component. The GS is a leading-twist phenomenon effective at small x 2 1 such that R G → 0 very slowly at Q 2 → ∞.

Color Glass Condensate
In addition to the sQGP matter, another instance where quarks and gluons cannot be treated as independent degrees of freedom is the case of parton coherence. A generalization of pQCD to hard collisions of small-x (x 1) partons (called also a semi-hard regime) was first discussed by Gribov, Levin, and Ryskin [279]. The basic failure of the standard DGLAP approach [280][281][282] is that it predicts too-fast increase of small-x parton density with the scale Q 2 . Consequently, the growth of hadronic cross-sections proceeds at rate which would sooner or later violate unitarity. The proposed solution-parton recombination and saturation-is at variance with the standard assumption that the partons themselves can be considered as independent free particles. The parameter determining the probability of parton-parton recombination is the ratio of the parton-parton cross-section to the square of the average distance between partons. The fact that the cross-section of such a semi-hard process (which now complies with the unitarity) increases rapidly with incident energy gives rise to expectations that (at least asymptotically) bulk particle production in hadron-hadron collisions can be described via pQCD [18].
The modern implementation of the above ideas is the Color Glass Condensate (CGC) formalism [283][284][285][286]-a natural generalization of pQCD to dense partonic systems. When applied to heavy nuclei, it predicts strong color fields in the initial stage of the collision. The strength of the fields is due to the condensation of low-x gluons into single macroscopic (i.e., classical) field state called the CGC. Since the characteristic scale of the parton saturation grows as Q s ∝ A 1/3 [18,284], it is enhanced on nuclear targets. According to the CGC-motivated phenomenology, the saturation phenomena are expected to show up (if not already) in p(d) + Au collisions at RHIC and , for sure, in nuclear collisions at the LHC. For example, due to the gluon saturation, the growth of the inelastic nucleon-nucleon cross-section σ in NN with increasing collision energy √ s may result in a broadening of the nucleon density distribution in position space. This in turn leads to a natural smoothing of the initial energy density distribution in the transverse plane of the matter created near midrapidity in heavy-ion collisions [207].
The CGC is described by an effective field theory that separates two kinds of degrees of freedom-fast frozen color sources and slow dynamical color fields. The basic evolution equation of such an effective field theory is an RG equation known as the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation [287][288][289][290][291][292], which reflects the independence of physical quantities with respect to variations of the cutoff separating these degrees of freedom (for more details, see Ref. [285] and references therein).
A supporting argument for the CGC as a possible state of QCD matter comes from successful analysis of HERA data in terms of geometrical scaling [293]. The geometrical scaling is the statement that the total γ * p cross-section depending a priori on two independent variables-the photon virtuality Q 2 and the Bjorken variable x-is a function of a single variable τ = Q 2 /Q 2 s , where the so-called saturation scale Q 2 s depends nontrivially on x, with dimensions given by a fixed reference scale Q 2 0 . However, calculations of Ref. [294] show that the standard linear leading-order DGLAP perturbative evolution is able to explain the geometric scaling. The situation with CGC applicability at current energies is thus unsettled (see also Refs. [253,254]). The experimental data from RHIC and the LHC, as well as exploitation of non-CGC-based models [228] are needed to resolve this problem.
There are two popular representative models of the initial state which are based on the CGC-the KLN model [295] and the IP Glasma model [296]. A Monte Carlo implementation of KLN CGC initial state [174,207,297] is based on the number distribution of gluons produced in the transverse plane given by the k T -factorisation formula [295]: Here, p ⊥ and y denote the transverse momentum and rapidity of the produced gluons, and x 1,2 = p ⊥ exp(±y)/ √ s NN are the light-cone momentum fractions of the colliding gluons. The running The gluon distribution function is given by An overall normalisation factor κ is chosen to fit the multiplicity data in most central Au + Au collisions at RHIC. In the MC-KLN model [295], the saturation momentum is parameterized by assuming that the saturation momentum squared is 2 GeV 2 at x = 0.01 in Au + Au collisions at b = 0 fm at RHIC, where ρ part = 3.06 fm −2 ; i.e., Here, λ is a free parameter which is expected to be in the range of 0.2 < λ < 0.3 from the global analysis of e + p scattering for x < 0.01 [241,293].
The IP-Glasma model [296] solves the classical Yang-Mills equations in which initial charge distributions of two colliding nuclei are sampled from a Gaussian distribution with the impact parameter and Bjorken x-dependent color charge distributions. A parameterization of x and impact parameter dependence of the saturation scale is taken from the IP-Sat (Impact Parameter Saturation) model [242,251]. Fluctuations in the IP-Glasma model have a length scale on the order of the inverse of the saturation scale Q −1 s (x ⊥ ) ∼ 0.1-0.2 fm. A comparison of the initial energy density distribution among the IP-Glasma, MC-KLN, and MC-Glauber models is shown in Figure 14.   [103]. The shading corresponds to the freeze-out temperature. The freeze-out occurs when the viscous terms become large compared to the ideal terms. Note that hydrodynamics breaks down not only at late, but also at early times (see the curve η/s = 0.4 in Figure 15). The right panel displays the centrality dependence of the elliptic flow coefficient v 2 (Equation 66) for two models for the initial density in the transverse plane-one is motivated by the parton saturation (CGC), and the other exploits nucleons only (Glauber). The calculations [298] were done within a hybrid model, where the expansion of the QGP starting at τ 0 = 0.6 fm/c is described by ideal hydrodynamics with a state-of-the-art lattice QCD EoS, and the subsequent evolution of hadronic matter below switching temperature T sw = 155 MeV is described using a hadronic cascade model. This nicely illustrates the strength of hydrodynamics-either the viscosity of QGP from RHIC to the LHC increases, or the CGC initial condition is ruled out [298].

Other Initial State Models
For completeness, let us here briefly list some other approaches to initial state description used in the analysis of ultra-relativistic heavy ion collisions. Fluctuating initial conditions given by a multi-phase transport (AMPT) model [165] were applied in an event-by-event partonic transport plus hydrodynamics hybrid approach [166,300,301] to study collective flow. In Refs. [302,303], NLO pQCD together with saturation-like suppression of low-energy partons were used to calculate the initial energy densities and formation times which have been used further on in 3D hydro calculations of space-time evolution of the QCD matter with dissipative fluid dynamics, event by event. A new initial conditions model for high-energy p + p, p + A, and A+A collisions that deposit entropy proportional to the generalized mean of nuclear overlap density was introduced in Ref. [304]. The model assumes that N one-on-one nucleon collisions produce the same amount of entropy as a single N-on-N collision.

Experimental Signatures of Deconfined QCD Matter
Evolution of the high energy nucleus-nucleus collision is schematically depicted in Figure 16. Two Lorentz-contracted pancakes of nuclear matter collide, thermalize, and form a deconfined QGP medium which expands, cools down, and hadronises to final state hadrons. Experimentally we do not observe each stage separately, but only through the time-integrated final state quantities-the momentum spectra of hadrons, photons, or leptons, particle multiplicities, energy flow, etc. Nevertheless, some time ordering of different processes giving rise to the final state observables exists. At very early collision times when colliding matter thermalizes, the entropy is produced which later-after (almost) isotropic expansion-transforms into particle multiplicities [28,181,205]. Early collision times also favour production of high p T partons [29,305] or heavy quarks (c, b) [306]. The formation of QGP reveals itself in many ways, including radiation of low momentum direct or virtual photons serving as a thermometers, enhanced production of hadrons containing strange (s) quarks [21,307,308], and melting of cc or bb mesons [115,116,306] called quarkonia. The subsequent rapid expansion of deconfined matter having more than ten times the degrees of freedom than the hadronic matter (see Equation 22)-and therefore also much higher internal pressure-produces a strong radial flow which leaves its imprint on the spectra of final state particles and their yields [28,183,205,309,310].
In the following, we present several examples of observables related to different stages of dynamics of nucleus-nucleus collisions at high-energies.

Bulk Observables
Traditionally, the very first measurements of heavy-ion collisions at a new energy regime comprise the charged-particle density at midrapidity dN ch /dη η=0 , also including its centrality dependence. Its collision-energy dependence for the 5% (6%) most central heavy-ion collisions-normalized per participant pair (i.e., N part /2)-is presented in Figure 17 (left panel). The right panel of Figure 17 shows that the normalized charged-particle density is rising with centrality, which means that the particle multiplicity at mid-rapidity increases faster than N part , presumably due to the contribution of hard processes to the particle production [29]. However, this increase is very similar to that observed at the top RHIC energy. Figure 17. Left: The charged particle rapidity density per participaiting pair, dN ch /dη/(0.5N part ), in central Au + Au, Cu + Cu, and Pb + Pb collisions from SPS to LHC energies. The star denotes an extrapolation to Pb + Pb collisions at 5.5 TeV. The IP -saturation model calculation [312] is illustrated by a dashed curve. Adapted from Ref. [30]; Right: dN ch /dη/(0.5N part ) vs. N part in Pb + Pb and Au + Au collisions at the LHC and RHIC, respectively. The RHIC data are multiplied by 2.15. The inset shows the N part < 60 region in more detail. Reproduced from Ref. [313].
One of the most celebrated predictions of the collective behaviour of matter created in non-central collisions of ultra-relativistic nuclei concerns its evolution in the transverse plane, which results from the pressure gradients due to spatial anisotropy of the initial density profile [183,314] (see Figure 18). The azimuthal anisotropy is usually quantified by the Fourier coefficients [315]: where φ is the azimuthal angle of the particle, Ψ n is the angle of the initial state spatial plane of symmetry, and n is the order of the harmonic. In a non-central heavy ion collision, the beam axis and the impact parameter define the reaction plane azimuth Ψ RP . For a smooth matter distribution in the colliding nuclei, the plane of symmetry is the reaction plane Ψ n = Ψ RP , and the odd Fourier coefficients are zero by symmetry. However, due to fluctuations in the matter distribution (including contributions from fluctuations in the positions of the participating nucleons in the nuclei-see Figure 11), the plane of symmetry fluctuates event-by-event around the reaction plane. This plane of symmetry is determined by the participating nucleons, and is therefore called the participant plane Ψ PP [316].
Since the planes of symmetry Ψ n are not known experimentally, the anisotropic flow coefficients are estimated from measured correlations between the observed particles [299,317]. This initial spatial anisotropy with respect to the reaction plane translates via pressure gradients into (c) a momentum anisotropy of the produced particles. Reproduced from Ref. [30].
In the following, we shall restrict ourselves to the properties of the Fourier coefficients v n with n = 2 and n = 3, which provide the dominant contributions to the observed azimuthal elliptic and triangular asymmetry, respectively. The sensitivity of v 2 to initial condition is illustrated on Figure 15 (right panel), where the centrality dependence of the elliptic flow in Pb + Pb collisions at √ s NN =2.76 TeV is shown. For more details on the corresponding initial state models, see Section 4.4.3.
The left panel of Figure 19 shows the measured energy dependence of the integrated elliptic flow coefficient v 2 in one centrality bin. Starting from √ s NN ≈ 5 GeV, there is a continuous increase of v 2 .
Below this energy, two phenomena occur. At very low energies, due to the rotation of the compound system generated in the collision, the emission is in-plane (v 2 > 0). At the laboratory kinetic energy around 100 MeV/nucleon, the preferred emission turns into out-of-plane, and v 2 becomes negative. The slowly moving spectator matter prevents the in-plane emission of participating nucleons or produced pions, which appear to be sqeezed-out of the reaction zone [318]. As the spectators move faster, their shadowing disappears, changing the pattern back to the in-plane emission.  Figure 19. Left: A compilation of data on the dependence of the integrated elliptic flow, v 2 , on the beam energy. The data correspond to (∼20-30%) Au + Au or Pb + Pb most central collisions. Reproduced from Ref. [193]; Right: The differential elliptic flow of charged particles. The Pb + Pb collisions at √ s NN = 2.76 TeV (colored symbols) are compared to Au + Au collisions at √ s NN = 200 GeV (grey lines).
Let us note that at RHIC, for the first time, the magnitude of the elliptic flow ( Figure 19) was found to be consistent with the EoS expected from the QGP [16,183]. The integrated value of v 2 for the produced particles increases by 70% from the top SPS energy to the top RHIC energy (see left panel of Figure 19), and it appears to do so smoothly. In comparison to the elliptic flow measurements in Au + Au collisions at √ s NN = 200 GeV, at the LHC, v 2 increases by about 30% at √ s NN = 2.76 TeV.
However, this increase is not seen in the differential elliptic flow of charged particles shown on the right panel of Figure 19. Thus, the bulk medium produced at RHIC and LHC has similar properties, and the 30% increase of v 2 between the two energies is due to an enlarged available phase space, resulting in the same increase of the average transverse momentum of particles <p T > between the RHIC and LHC energies.
As was first noted in Ref. [314], at high energies, only the interactions among the constituents of matter formed in the initially spatially deformed overlap can produce v 2 > 0. A transfer of this spatial deformation into momentum space provides a unique signature for re-interactions in the fireball, and proves that the matter has undergone significant nontrivial dynamics between its creation and its freeze-out [183]. The rapid degradation of the initial spatial deformation due to re-scattering causes the "self-quenching" of elliptic flow: if the elliptic flow does not develop early (when the collision fireball was still spatially deformed), it does not develop at all [183]. In particular, the transformation of the rapidly expanding ideal gas of non-interacting quarks and gluons into strongly interacting hadrons is unable to produce a sufficient elliptic flow. The elliptic flow thus reflects the pressure due to re-scattering-the induced expansion and stiffness of the EoS during the earliest collision stages. Its continuous rise with the energy up to its highest value at the LHC indicates that the early pressure also increases.
The energy dependence of the integrated triangular flow coefficient v 2 3 {2} of charged hadrons is shown on the left panel of Figure 20 in four bins of centrality, 0-5%, 10-20%, 30-40%, and 50-60%. As v 2 3 {2} is sensitive to the fluctuations in the initial matter distribution, it is interesting to observe that at √ s NN = 7.7 and 11.5 GeV, values of v 2 3 {2} for 50-60% centrality become consistent with zero. For more central collisions, however, v 2 3 {2} is finite-even at the lowest energies-and changes very little from 7.7 GeV to 19.6 GeV. Above that, it begins to increase more quickly, and roughly linearly with log( √ s NN ). Generally, one would expect that higher energy collisions producing more particles should be more effective at converting the initial state geometry fluctuations into v 2 3 {2}. Deviations from that expectation could indicate interesting physics, such as a softening of the EoS [113,114] discussed already in Section 3.3. This can be investigated by scaling v 2 3 {2} by the charged particle rapidity density per participating NN pair, n ch,PP = dN ch /dη/(0.5N part ), see the left panel of Figure 20. A local minimum of v 2 3 {2}/n ch,PP in the region near 15-20 GeV observed in the centrality range 0-50% and absent in the more peripheral events could indicate an interesting trend in the pressure developed inside the system.  As was already briefly mentioned in Section 4.4.4, recent years have witnessed a growth of interest in studying collective phenomena on event-by-event basis. In addition to the key publications [165,166,[205][206][207][208][300][301][302][303], it is worth mentioning the study [320] where the anisotropic flow coefficients v 1 -v 5 were computed by combining the IP-Glasma flow with the subsequent relativistic viscous hydrodynamic evolution of matter through the quark-gluon plasma and hadron gas phases. The event-by-event geometric fluctuations in nucleon positions and intrinsic color charge fluctuations at the sub-nucleon scale are expected to result in experimentally measurable event-by-event anisotropic flow coefficients [321]. Let us note that fluctuating initial profiles observed in over-many-events integrated triangular and higher odd flow coefficients are also revealed in difference between various participant plane angles Ψ n introduced in Equation (66). Correlations between the event plane angles Ψ n of different harmonic order can not only yield valuable additional insights into the initial conditions [322], but are also experimentally measurable [323]. The same is also true for the correlations between different flow harmonics [324][325][326][327].
Enhanced production of hadrons with the quantum numbers not present in colliding matter is one of the oldest signals of the deconfined QGP medium [307,308]. Measurements of the yields of multistrange baryons were carried out at CERN SPS by WA85, and later on by WA97/NA57 collaborations since the mid-eighties. After 2000, more data came from RHIC, and starting from 2010 also from the LHC. The current status is summarized in the five panels of Figure 21. In the top left and middle panels (Figure 21a,b), a compilation of the results from SPS, RHIC, and the LHC in terms of strangeness enhancement defined as normalized (to p + p or p + Be) yield per participants is presented. On the top right (Figure 21c), the hyperon-to-pion ratios as functions of N part for Pb + Pb, Au + Au, and p + p collisions at the LHC and RHIC energies are displayed. The normalized yields are larger than unity for all the particles, and increase with their strangeness content. This behaviour is consistent with the picture of enhanced ss pair production in a hot and dense QGP medium [307,308]. Boxes on the dashed line at unity indicate the statistical and systematic uncertainties on the p + p or p + Be reference; (c) Hyperon-to-pion ratios as a function of N part for Pb + Pb, Au + Au, and p + p collisions at LHC and RHIC energies. The lines are the thermal model [309] (full line) and [328] (dashed line) predictions. Reproduced from Ref. [329]; Bottom: (left) (Ξ − + Ξ + )/(π + +π − ) and (right) (Ω − +Ω + )/(π + +π − ) ratios as functions of dN ch /dη for p + p, p + Pb, and Pb + Pb collisions at the LHC. The Pb + Pb data points [329] represent, from left to right, the 60-80%, 40-60%, 20-40%, 10-20%, and 0-10% centrality classes. The chemical equilibrium predictions by the GSI-Heidelberg [309] and the THERMUS 2.3 [328] models are represented by the horizontal lines. Reproduced from Ref. [330].
The two bottom plots represent a comparison between the hyperon-to-pion ratios from p + p, p + Pb, and Pb + Pb collisions. Interestingly, the ratios in p + Pb collisions increase with multiplicity from the values measured in p + p to those observed in Pb + Pb. The rate of increase is more pronounced for particles with higher strangeness content. Let us note that the Grand canonical statistical description of Pb + Pb data shown as full and dashed lines in Figure 21 may not be appropriate in small multiplicity environments such as those produced in the p + Pb case. It appears that for the latter case, the evolution of hyperon-to-pion ratios with the event multiplicity is qualitatively well described by the Strangeness Canonical model implemented in THERMUS 2.3 [328]. In this case, a local conservation law is applied to the strangeness quantum number within a correlation volume V c , while treating the baryon and charge quantum numbers grand-canonically within the whole fireball volume V [330].

Hard Probes
Heavy quarks, quarkonia, and jets-commonly referred to as hard probes-are created in the first moments after the collision, and are therefore considered as key probes of the deconfined QCD medium. Production of these high transverse momentum (p T Λ QCD ) objects occurs over very short time scales (τ≈1/p T ≈0.1 f m/c), and can thus probe the evolution of the medium. Since the production cross-sections of these energetic particles are calculable using pQCD, they have been long recognised as particularly useful "tomographic" probes of the QGP [331][332][333].
Let us start our discussion with the results on the inclusive production of high-p T hadrons. The latter are interesting on their own because it was there where for the first time the suppression pattern was observed [14][15][16][17]. In an inclusive regime, the comparison between d 2 N/dp T dη (the differential yield of high-p T hadrons or jets per event in A+B collisions) to that in p + p collisions is usually quantified by introducing the nuclear modification factor For collisions of two nuclei behaving as a simple superposition of N coll nucleon-nucleon collisions, the nuclear modification factor would be R AB = 1. The data of Figure 22 reveal a very different behaviour. The left panel shows a compilation of R AA from Au + Au and Pb + Pb collisions, the right panel the result of R pPb from three LHC experiments at the same energy √ s NN = 5.02 TeV. In the R AA case, the suppression pattern of high-p T (>2-3 GeV/c) hadrons in the deconfined medium-predicted many years ago [331][332][333] as a jet quenching effect-is clearly visible at RHIC and the LHC. However, for proton-nucleus collisions ( Figure 22, right panel), no suppression is seen, even at the highest LHC energy. Moreover, R AA in the 5% most central Pb + Pb collisions at the LHC shows a maximal suppression by a factor of 7-8 in the p T region of 6-9 GeV. This dip is followed by an increase, which continues up to the highest p T measured at √ s NN = 5.02 TeV, and approaches unity in the vicinity of p T = 200 GeV [334].

High-p T Hadrons and Jets
The suppression of high-p T hadrons in the deconfined medium was thoroughly studied at RHIC using azimuthal correlations between the trigger particle and associated particle; see Figure 23. Near-side peaks in central (0-5%) Au + Au collisions present in all panels of Figure 23 (left) indicate that the correlation is dominated by jet fragmentation. An away-side peak emerges as p trig T is increased. The narrow back-to-back peaks are indicative of the azimuthally back-to-back nature of dijets observed in an elementary parton-parton collision. Contrary to the latter, the transverse-momentum imbalance of particles from the jet fragmentation due to different path lengths of two hard partons in the medium is apparent. The azimuthal angle difference ∆φ for the highest p trig T range (8 < p trig T < 15 GeV/c) for mid-central (20-40%) and central Au+Au collisions-as well as for d+Au collisions-is presented in Figure 23 (right panel). We observe narrow correlation peaks in all three p assoc T ranges. For each p assoc T , the nearside peak shows a similar correlation strength above background for the three systems, while the away-side correlation strength decreases from d + Au to central Au + Au. For the d + Au case, the yield of particles on the opposite side ∆φ = π prevails over the same side. Moreover, for Au + Au collisions, the nearside yields obtained after subtraction of the background contribution due to the elliptic flow show a little centrality dependence, while the away-side yields decrease with increasing centrality [335]. Unfortunately, the advantage of the large yield of dijets is offset by a loss of information about the initial properties of the probes (i.e., prior to their interactions with the medium). Correlating two probes that both undergo an energy loss also induces a selection bias towards scatterings occurring at-and oriented tangential to-the surface of the medium. It is thus interesting to study correlations when one of the particles does not interact strongly with the medium. Triggering on the high-p T isolated photon (i.e., not from π 0 → 2γ decays) would do the job. While in p + p collisions an emerging quark jet should balance its transverse momentum with the photon, in the heavy-ion collisions, much of its momentum is thermalized while the quark traverses the plasma. This is illustrated in Figure 24 (left panel), where a single hard photon with p T = 402 GeV emerges unhindered from the de-confined medium produced in Pb+Pb collisions at the LHC. The accompanying quark jet produced via the QCD Compton scattering qg → qγ loses 1/3 of its energy (≈140 GeV!) inside the hot and dense matter.
The measurement presented on the middle and right panels of Figure 24 shows that for more central Pb + Pb collisions, a significant decrease in the ratio of jet transverse momentum to photon transverse momentum-x Jγ -relative to the PYTHIA reference [200] is observed. Furthermore, significantly more photons with p T > 60 GeV/c in Pb + Pb are observed to not have an associated jet with p T > 30 GeV/c jet, compared to the reference. However, no significant broadening of the photon + jet azimuthal correlation has been observed. An important progress in the theoretical understanding of the suppression of energetic partons traversing a deconfined matter was the introduction of the diffusion coefficientq relevant for the transverse momentum broadening and collisional energy loss of partons (jets) [338]. This quantity, which is commonly referred to as the jet quenching parameter, can be determined either via weak coupling techniques [339][340][341], a combination of lattice simulations and dimensionally-reduced effective theory [342], or from the gauge/gravity duality [343]. Typical estimates for this quantity at RHIC and LHC energies range between 5 and 10 GeV 2 /fm, demonstrating the currently still sizable uncertainties in these calculations.

Quarkonia
Melting of the quarkonia-bound states of heavy quark and anti-quark qq where q = c, b-due to a colour screening in the deconfined hot and dense medium was proposed thirty years ago as a clear and unambiguous signature of the deconfinement [115]. However, shortly after that, it was noticed that not only diffusion of the heavy quarks from melted quarkonium, but also the drag which charm quarks experience when propagating through the plasma is important [344].
The latter might lead to an enhancement instead of a suppression. This is in variance with the original proposal that the heavy quarks, once screened, simply fly apart. With the advent of the strongly-interacting QGP, the Langevin equation model of quarkonium production was formulated, where the charm quark-antiquark pairs evolve on top of a hydrodynamically expanding fireball [345]. A heavy quark and anti-quark interact with each other according to the screened Cornell potential and interact, independently, with the surrounding medium, experiencing both drag and rapidly decorrelating random forces. An extension of this approach to bottomonium production [346] shows that a large fraction of bb pairs that were located sufficiently close together during the initial hard production will remain correlated in the hot medium for a typical lifetime of the system created in heavy-ion collisions. The distribution of the correlated bb pair in relative distance is such that it will dominantly form 1S bottomonium. A study of quarkonia production in heavy-ion collisions thus provides an interesting window not only into static, but also into dynamical properties of the hot, dense, and rapidly expanding medium [116,306].
On the left panels of Figure 25, the invariant-mass distributions of µ + µ − pairs (di-muons) produced in the p + p (a) and Pb + Pb (b) collisions at the LHC are presented. A prominent peak due to the production of the heavy quarkonium state, the bottomonium Υ(1S), can be clearly seen in both p + p and Pb + Pb data. Peaks from the higher excited states of Υ, Υ(2S), and Υ(3S), although discernible in the p + p case, are barely visible in the Pb + Pb data. More quantitative information on this effect can be found in the right panels of Figure 25, where the centrality dependence of the double ratio [Υ(2S)/Υ(1S)] PbPb / [Υ(2S)/Υ(1S)] pp (top) and of the nuclear modification factors R AA of Υ(1S) and Υ(2S) (bottom) are displayed. Let us note that the observed suppression of the relative yield is in agreement with the expectations that different quarkonium states will dissociate at different temperatures with a suppression pattern ordered sequentially with the binding energy; i.e., the difference between the mass of a given quarkonium and twice the mass of the lightest meson containing the corresponding heavy quark [347]. Moreover, the observed pattern is now also confirmed in Pb + Pb collisions at √ s NN = 5.02 TeV [348]. The double ratio is significantly below unity at all centralities, and no variation with kinematics is observed, confirming a strong Υ suppression in heavy-ion collisions at the LHC.

Penetrating Probes
Electromagnetic probes like photons [352,353] and di-leptons [354] (for recent developments, see Refs. [253][254][255]) have long been expected to provide crucial information on the properties of QGP. The absence of strong final-state interactions makes them an ideal penetrating probe of strongly-interacting matter [355]. In collisions of ultra-relativistic nuclei, the photons and leptons can be produced either in the initial hard collisions between partons of the incident nuclei (e.g., qg → qγ, qq → γg, or qq → ¯ ), or radiated from the thermally equilibrated partons and hadrons or via hadronic decays. The direct photons are defined to be all produced photons, except those from hadron decays in the last stage of the collision. The high-p T isolated photon can be used to estimate the momentum of the associated parton, allowing a characterisation of the in-medium parton energy loss (see Figure 24). The prompt photons also carry information about the initial state and its possible modifications in nuclei, and should thus be one of the best probes of the gluon saturation. The thermal photons emitted from the produced matter in nuclear collisions carry information on the temperature of QGP.
The first observation of direct photons in ultra-relativistic heavy-ion collisions has been made by the CERN SPS experiment WA98 [356]. In 10% most central Pb + Pb collisions at √ s NN = 17.2 GeV, they observed a clear excess of direct photons in the range of p T > 1.5 GeV/c which was not present in more peripheral collisions. The extraction of the direct photon signal (which is extremely difficult) was described in-depth in Ref. [357]. The data from RHIC and the LHC are presented in Figure 26 on the left and right panels, respectively. The invariant cross-section (p + p) and invariant yield (Au + Au) of direct photons as functions of p T . The three curves on the p + p data represent the NLO pQCD calculations, and the dashed curves show a modified power-law fit to the p + p data, scaled by T AA = N coll /σ in NN . The dashed (black) curves are the same, but with the exponential plus scaled p + p data. The dotted (red) curve near the 0-20% centrality data is a theory calculation. Reproduced from Ref. [350]; Right: The direct photon spectra in Pb + Pb collisions at √ s NN = 2.76 TeV for the 0-20% (scaled by a factor of 100), the 20-40% (scaled by a factor of 10), and the 40-80% centrality classes compared to the NLO pQCD predictions for the direct photon yield in pp collisions at the same energy, scaled by a number of binary nucleon collisions for each centrality class. Reproduced from Ref. [351].

Search for the Critical Point of QCD Phase Diagram
The search for the (tri)critical point (CEP) in the T − µ B phase diagram-where the phase transition between the QGP and hadron matter changes from the first-to the second-order one-represents one of the most active fields of contemporary ultra-relativistic heavy-ion physics, both experimentally [59,[360][361][362] and theoretically [58,63,363,364]. In order to gain more insight into the CEP location, quite advanced techniques from condensed matter physics, such as Finite-Size Scaling (FSS) analysis of data [365] or thermal fluctuations characterized by the appropriate cumulants of the partition function [363,364,366] are being exploited.
The search for the CEP exploiting the potential of the RHIC accelerator complex was mounted by the STAR and PHENIX collaborations in 2010 within the Beam Energy Scan (BES) program [59]. Going down from the RHIC maximum energy √ s NN = 200 GeV, they have scanned the available phase space down to √ s NN = 7.7 GeV. Some of the results from that scan were already mentioned in previous sections, and can be found in Figures 8,10,17,19,20. In the following, we therefore restrict ourselves to the measurements which provide a direct link to the lattice results discussed in Section 3.4. As already mentioned there, the study of ratios of the Taylor expansion coefficients given by Equation (31)-which are also known also as susceptibilities-seems to be very attractive, since both the temperature and volume dependences drop out. In particular, the ratio χ B 4 /χ B 2 calculated from the moments of the net-baryon multiplicity N B has different values for the hadronic and partonic phases [366]. For HRG, it equals unity, but is expected to deviate from unity near the CEP. Other interesting ratios of the net-baryon charge moments which can be expressed using mean (M B ), variance (σ 2 B ), skewness (S B ), and kurtosis (k B ) of the net baryon number distributions read Experimentally, the net-baryon number N B fluctuations and their cumulants are not accessible, and so one has to resort to measurements of the cumulants of the net-proton number N P fluctuations [358,367]. On the other hand, the electric charge fluctuations are experimentally accessible [359]. This is illustrated in Figure 27, where the measurements from the STAR experiment at RHIC are shown. Generally, study of the p T and rapidity acceptance dependence for the moments of the net-proton distributions shows that the larger the acceptance is, the larger are the deviations from unity. Preliminary results from the STAR collaboration [368] reveal that increasing the transverse momentum acceptance of protons and anti-protons from 0.4 < p T < 0.8 GeV/c to 0.4 < p T < 2.0 GeV/c leads to a pronounced non-monotonic structure in the energy dependence of κσ 2 of net-proton distributions from 5% most central Au + Au collisions. At energies above 39 GeV, the values of κσ 2 are close to unity, while for energies below 39 GeV, it shows a significant deviation below unity around 19.6 and 27 GeV and a large increase above unity observed at 7.7 GeV.

Collectivity in Small Systems
Recent years have witnessed a surprising development in multiparticle dynamics of high multiplicity p + p [369][370][371] and p + A [372,373,[378][379][380][381] collisions. It all started in 2010 with the observation of ridge-like structures in p + p collisions by the CMS experiment at the LHC [369]. The surprise was due to the fact that a very similar effect was found just a few years before in heavy-ion collisions: first in Au + Au collisions at √ s NN = 200 GeV at RHIC [382,383], and later on also in Pb + Pb collisions at the LHC [384] and in Cu + Cu collisions at √ s NN = 62.4 GeV and 200 GeV at RHIC [385].
In heavy-ion collisions, it was found that pairs of particles are preferentially emitted with small relative azimuthal angles (∆φ = φ 1 − φ 2 ≈ 0). Surprisingly, this preference persists even when the particles are separated by large pseudo-rapidity (η) gaps (−4 < |∆η| < 2). These long-range correlations-known as the ridge-have been traced to the conversion of density anisotropies in the initial overlap of the two nuclei into momentum space correlations through subsequent interactions in the expansion [386].
In p + p minimum bias collisions at the LHC, the peak in the correlation function of particles with p T > 0.1 GeV/c observed at small angular differences (∆η, ∆φ≈ 0, see Figure 28a) is due to several effects: resonance decays, Bose-Einstein correlations, and near-side jet fragmentation. The fragmentation due to back-to-back jets is visible as a broad elongated ridge around ∆φ≈ π. The pattern does not change much even when selecting the events with very high multiplicity N ≥ 110 (see Figure 28c). The cut on the multiplicity enhances the relative contribution of high p T jets, which fragment into a large number of particles, and therefore, has a qualitatively similar effect on the shape as the particle cut 1 < p T < 3 GeV/c on minimum bias events (see Figure 28b). However, using now the particle cut 1 < p T < 3 GeV/c in conjunction with a high multiplicity cut changes the picture dramatically (see Figure 28d). A novel feature never seen before in p + p collisions at lower energies shows up-a clear and significant ridge-like structure at ∆φ≈ 0 extending to |∆η| of at least four units [369].
Let us note that for two particles with approximately the same energy E 1 ≈ E 2 ≈ E, the correlations at ∆η≈ ∆y are by the uncertainty relation ∆x ≈ 1/∆p = 1/(E∆y) connected to the correlations in coordinate space. While for pions with ∆η≈ 1 and p T ≈ 0.1 GeV/c, we have ∆x ≈ 1 fm, for ∆η≈ 4 and p T ≈ 1 GeV/c, one gets ∆x ≈ 0.02 fm. It is obvious that at such small inter-parton distances, one enters the realm of the initial state description of nuclear collisions when the density of matter even inside the proton fluctuates-see Section 4.4.3. The CGC scenario was recently exploited in Ref. [387] to predict long-range photon-jet correlations in p + p and p + A collisions at near-side for low transverse momenta of the produced photon and jet in high-multiplicity events.
The relevance of the saturation approach is further supported by the observation of the same ridge-like structure in p + p collisions at √ s =13 TeV [370,371], see Figure 29. For the associated yield of long-range near-side correlations for high-multiplicity events (N >110) peaks in the region 1 < p T < 2 GeV/c, see Figure 29a-the yield reaches a maximum around p T ≈1 GeV/c and decreases with increasing p T . No center-of-mass energy dependence is visible. The multiplicity dependence of the associated yield for 1 < p T < 2 GeV/c particle pairs is shown in Figure 29b. For low-multiplicity events, the associated yield is consistent with zero. At higher multiplicity, the ridge-like correlation emerges, with an approximately linear rise of the associated yield with multiplicity for N ≥ 40. Let us note that within the CGC models, the observation that the integrated near-side yield as a function of multiplicity is independent of collision energy is a natural consequence of the fact that multiparticle production is driven by a single semi-hard saturation scale [388].  Figure 28. The 2-D two-particle correlation functions for p + p at √ s = 7 TeV: (a) minimum bias events with p T > 0.1 GeV/c; (b) minimum bias events with 1 < p T < 3 GeV/c; (c) high multiplicity (N offline trk ≥ 110) events with p T > 0.1 GeV/c; and (d) high multiplicity (N offline trk ≥ 110) events with 1< p T < 3 GeV/c. The sharp near-side peak from jet correlations is cut off in order to better illustrate the structure outside that region. Reproduced from Ref. [369].  . The associated yield for the near-side of the correlation function averaged over 2 < |∆η| < 4 for p + p data at √ s = 13 TeV (filled circles) and 7 TeV (open circles) [370]. Panel (a) shows the associated yield as a function of p T for events with N offline trk ≥ 105. In panel (b), the associated yield for 1 < p T < 2 GeV/c is shown as a function of multiplicity-N offline trk . The p T selection applies to both particles in each pair. Curves represent the predictions of the gluon saturation model [388].
Another interesting phenomenon observed during recent years is a flow-like pattern in super-central p + p and p + Pb collisions at the LHC [373][374][375][376][377][378] and in d + Au, 3 He + Au, and p(d) + Au collisions at RHIC [380,381]. These collisions not only reveal a similar elliptic flow v 2 , but in some cases, also v 3 anisotropy previously observed only in collisions of large nuclei. The measurement of higher-order cumulants of the azimuthal distributions has strengthened the collective nature interpretation of the anisotropy seen in p + Pb collisions. Moreover, the collective radial flow analysis of p + p data enabled authors of Ref. [389] to claim that in this case the fireball explosions start with a very small initial size, well below 1 fm. This raises questions about whether the perfect liquid sQGP is also formed in these much smaller systems [390][391][392][393][394][395]. Are these flow-like structures only similar in appearance to what one observes in heavy-ion collisions, or do they have the same physical origin [396][397][398]? Obviously, answering these questions may help us to understand the emergence of collective phenomena in strongly-interacting systems in general. To make further progress on these fundamental questions, more analyses-both experimental and theoretical-are really needed [399].

Conclusions
In this review, we have tried to present an up-to-date phenomenological summary of a relatively new and rapidly developing field of contemporary physics-the physics of Quark Gluon Plasma. We have also explored its broader ramifications when discussing matter under extreme conditions, strongly interacting plasmas, physics of strong electromagnetic fields, history of ultra-relativistic heavy-ion collisions, relativistic hydrodynamics, the role of the QCD ground state, and QCD saturation phenomena. In the experimental part, we were overwhelmed by a huge amount of results, and so, in order to keep this review of tolerable length, we had to skip several quite important topics. To name just a few-flow and suppression of identified particles [30,400], femtoscopy [401,402], or identified particle yields and chemical freeze-out conditions [29,30,309]. We hope that the interested reader will find in the given references enough information to pursue a deeper study of these important subjects. However, in skipping these topics, we hope that we have not given up our main goal-to give the reader a possibility to see the QGP landscape at large. At last, we would like to say that it is never enough to stress how important this field is also for other branches of physics, and so we finish with yet another argument. The study of ultra-relativistic heavy-ion collisions appears so far to be our only way of studying the phase transitions in non-Abelian gauge theories (most likely taken place in the early universe) under laboratory conditions.