The exploration of hot and dense nuclear matter: introduction to relativistic heavy-ion physics

The approaches described in the previous section allow theorists to relate various fundamental properties of the quark-gluon plasma to specific experimental observables:

• (i)  
  The equation of state of the quark-gluon plasma encodes the dependence of the pressure p on the two intensive variables energy density ε and net baryon density $n={n}_{{\rm{b}}}-{n}_{\bar{{\rm{b}}}}$, where n_b (${n}_{\bar{{\rm{b}}}}$) denotes the density of baryons (antibaryons). It is reflected in the final particle spectra, collective flow properties, and the propagation of initial-state density fluctuations into the final state. The speed of sound, ${c}_{s}=\sqrt{\partial p/\partial \varepsilon }$, determines how the pressure gradients propagate during the expansion of the system.
• (ii)  
  The viscosities η (shear) and ζ (bulk), which may be expressed as autocorrelation functions of the components of the stress tensor T_ij [28],

  $$\begin{eqnarray}&&\eta =\frac{1}{T}\int {{\rm{d}}}^{4}x\,\langle {T}_{{xy}}(x){T}_{{xy}}(0)\rangle ,\end{eqnarray} \tag{ 1 }$$

  $$\begin{eqnarray}&&\zeta =\frac{1}{9T}\displaystyle \sum _{i,j}\int {{\rm{d}}}^{4}x\,\langle {T}_{{ii}}(x){T}_{{jj}}(0)\rangle .\end{eqnarray} \tag{ 2 }$$

  influence the collective transverse flow pattern, quantified by Fourier decomposition of the azimuthal distribution of particles in momentum space. The anisotropy of the transverse flow is especially sensitive to the dimensionless ratio η/s, where s denotes the entropy density.
• (iii)  
  Momentum transport coefficients, including the so-called jet quenching parameter $\hat{q}$, are expressed as correlation functions of the color force experienced by a fast parton in the quark-gluon plasma:

  $$\begin{eqnarray}&&\hat{q}\propto \int {\rm{d}}{x}^{-}\langle \mathrm{Tr}\,[{U}_{0,x}{F}^{+i}(x){U}_{x,0}{F}_{i}^{+}(0)]\rangle .\end{eqnarray} \tag{ 3 }$$

  Here F denotes the field strength tensor and U are so-called gauge links that ensure the gauge invariance of $\hat{q}$. The minus index in x⁻ indicates that the correlation is sampled along the light-cone. Together with the coefficient $\hat{e}$, which quantifies the rate of longitudinal momentum loss due to collisions with medium constituents [29], $\hat{q}$ determines the rate of energy loss of a fast parton traversing the quark-gluon plasma. Experimentally, these coefficients can be determined by measuring the suppression of hadrons emitted with high momentum transverse to the beam direction.
• (iv)  
  The Debye screening length λ_D of the color force between a heavy quark-antiquark pair in the quark-gluon plasma, which is determined by the correlator of the static gauge potential between two points separated by a spatial distance x:

  $$\begin{eqnarray}&&{\lambda }_{{\rm{D}}}^{-1}=-\mathop{\mathrm{lim}}\limits_{| x| \to \infty }\displaystyle \frac{\langle \mathrm{Tr}\,[{U}_{0,x}{A}^{0}(x){U}_{x,0}{A}^{0}(0)]\rangle }{| x| },\end{eqnarray} \tag{ 4 }$$

  controls the ability of a pair of heavy quarks to form a bound state. When λ_D is less than the size of the bound state, the state dissolves in the medium. Experimentally, this phenomenon is expected to be revealed by a strong suppression of charmonium ($c\overline{c}$) and bottomonium ($b\overline{b}$) states.
• (v)  
  The electromagnetic current-current correlation function

  $$\begin{eqnarray}&&{C}_{\mathrm{em}}^{\mu \nu }(q)=\int {{\rm{d}}}^{4}x\,{{\rm{e}}}^{{\rm{i}}q\cdot x}\langle {j}^{\mu }(x){j}_{\nu }(0)\rangle \end{eqnarray} \tag{ 5 }$$

  encodes the response of the quark-gluon plasma (or hadronic medium) to electromagnetic fields. It is directly only sensitive to quarks, but indirectly also to gluons, because quarks can be off-shell due to their interactions with thermal gluons. Observables carrying information about ${C}^{\mu \nu }{\left(q\right)}_{\mathrm{em}}$ are photons, sensitive to q² = 0, and dileptons, sensitive to time-like q² > 4m² with m being the lepton mass.

These well-defined quantities summarize the main properties of the hot and dense QCD matter. Their linkage to experimental observables lays the ground for robust theoretical and experimental research programs in heavy ion physics. In making this linkage it is important to keep in mind that only the momentum space properties of the matter are accessible to experiments, because the fireball is too small and short-lived to allow direct spatial of temporal resolution. It is possible to infer some information about the space- and time-dependence of the collision from two-particle correlations, which will be discussed in section 2.4 under the topic 'HBT Femtoscopic Radii' ⁶ , but even these indirect measurements are limited to the last stage of the collision and require detailed modeling of the reaction.

From a practical viewpoint, on the experimental side, it is important that precision measurements of the relevant observables are made over a wide kinematic range so that theoretical models of the collision can be constrained. On the theoretical side, precise and reliable calculations of the transport coefficients under given conditions are needed in combination with realistic simulations of the dynamical evolution that connect the matter properties to the observables.

The following sections introduce the basic concepts underpinning the most commonly considered observables and point out their generally accepted connections to the properties of hot and dense QCD matter.

Single-particle observables are all the observables that can be measured by a detector that records only individual particles. Below we discuss a few specific observables of this type that are of interest in heavy ion physics.

Particle yields and multiplicities

The most basic analysis in a heavy ion collision is to count the number of particles of a certain species that are produced. One can either restrict the measurement of the yield to a region in momentum space around mid-rapidity or, if detector acceptance allows, extrapolate the number to the full multiplicity including the whole phase-space (4π yields). By fitting the ratios of particle multiplicities to a grand canonical thermal distribution one can infer information about the conditions at the time of chemical freeze-out, after which the particle yields do not change anymore (see section 5.1).

Rapidity distributions

Rapidity distributions dN/dy are sensitive to the initial energy deposition and thus to the stopping power of the colliding nuclei. Extreme limits of the charged particle rapidity distribution are associated with Landau (Gaussian shape) or Bjorken (flat shape) hydrodynamic behavior. The Landau limit corresponds to full stopping followed by explosive expansion. The Bjorken model corresponds to complete transparency for the valence quarks resulting in longitudinal boost invariance—this approximation is only strictly valid in the infinite collision energy limit.

Transverse momentum spectra

Counting the particles in bins of transverse momentum p_T or transverse mass ${m}_{T}=\sqrt{{m}^{2}+{p}_{T}^{2}}$ conveys information about the kinetic decoupling temperature and the flow profile

$$\begin{eqnarray}&&\frac{1}{{m}_{T}}\frac{{\rm{d}}N}{{\rm{d}}{m}_{T}}\propto {\exp }^{-{m}_{T}/{T}_{\mathrm{eff}}}.\end{eqnarray} \tag{ 6 }$$

By fitting the particle spectra to a thermal distribution one can extract an effective 'temperature' T_eff, which is correctly called the slope parameter. If the emitted hadrons originate from a locally equilibrated and collectively flowing momentum distribution, the slope parameter represents a blue-shifted temperature. The outward radial flow that the fireball develops during the dynamical evolution generally results in a flatter shape of the spectra of higher mass particles, e.g. baryons, at low transverse momentum (see section 5.2).

Nuclear modification factor R_AA

The nuclear modification factor R_AA is the ratio of the yield at a given transverse momentum p_T and rapidity y measured in a heavy ion collision (AA) normalized by the yield measured in the corresponding number of independent nucleon–nucleon collisions (p+p)

$$\begin{eqnarray}&&{R}_{{AA}}({p}_{T},y;\langle b\rangle )=\frac{{{\rm{d}}}^{2}N/{\rm{d}}{p}_{T}{\rm{d}}{y}_{{AA}}}{\langle {N}_{{coll}}\rangle {{\rm{d}}}^{2}N/{\rm{d}}{p}_{T}{\rm{d}}{y}_{p+p}}.\end{eqnarray} \tag{ 7 }$$

Here 〈N_coll〉 denotes the average number of binary nucleon–nucleon collisions in heavy ion collisions within a certain centrality window; 〈b〉 is the average impact parameter associated with this window. In any given collision N_coll reflects the geometric overlap of the colliding nuclei but fluctuates around the average value by virtue of quantum fluctuations in the nuclear densities and the nucleon–nucleon scattering probability. Sometimes good reference data from p+p collisions do not exist. In such cases, one often takes peripheral heavy ion collisions as a reference point and uses the similarly defined binary collision-weighted ratio of yields measured in central and peripheral collisions, denoted as R_CP.

If there are final-state interactions with the medium that degrade the transverse momentum of energetic partons, R_AA falls below unity, which explains the alternative name of this observable: nuclear suppression factor. ⁷ This phenomenon is referred to as jet quenching. For light-heavy collision systems, such as p + A or d + A, the R_p(d)A is sensitive to initial-state effects that can be attributed to the behavior of cold nuclear matter relative to the reference of proton–proton collisions. For heavy-heavy collision systems the nuclear modification factor is mainly sensitive to final-state effects. The first determination of the energy loss parameter $\hat{q}$ is based on measurements of the nuclear modification factor (see section 5.9). R_AA is also used to quantify the nuclear suppression or enhancement (with respect to p–p collisions) of many other observables in heavy ion collisions, such as the yields of heavy quarkonia, heavy-flavor hadron, photons, etc.

Photon spectra

Direct photon spectra reflect the cumulative emission of thermal and prompt photons from all stages of a heavy ion reaction. The definition of a 'direct photon' is that it does not originate from a particle decay, such as π⁰ → γ γ. In the initial non-equilibrium evolution high-p_T photons are produced in hard processes. Thermal emission from the quark-gluon plasma and the subsequent hadron gas phase are added during later collision stages. Since the mean free path of photons is larger than the fireball size, once photons are produced they end up in the detector and, different from strongly interacting probes, are minimally distorted by rescattering processes. By detailed dynamical modeling and matching of the various contributions to the photon spectrum, the goal is to extract information about the initial temperature and the lifetime of the QGP fireball from this observable. Photon spectra are also sensitive to the equation of state, since the emission depends strongly on the lifetime of the system that is substantially longer when a first-order phase transition takes place (see section 5.11).

Dilepton spectra

Another electromagnetic probe are lepton pairs, commonly referred to as dileptons. These originate from virtual photons that decay into a e⁺ e⁻ or μ⁺ μ⁻ pair. The theoretical advantage is that leptons do not interact strongly and therefore can be measured from all stages of the heavy ion reaction. On the other hand, the cross-sections for leptonic pair production are orders of magnitude smaller than for hadronic emission. Since some hadron decays (e.g. π⁰ → γ e⁺ e⁻, D⁰ → K⁻ μ⁺ ν) produce leptons which generates a large combinatorial background of uncorrelated lepton pairs, dilepton spectra constitute challenging measurements. Vector mesons, such as the ρ and J/ψ mesons, can decay into lepton pairs and show up as peaks in the dilepton invariant mass spectrum. By carefully analyzing the height and width of these peaks, one hopes to learn something about the influence of the medium on the resonance properties. Changes in the spectral function of the ρ meson are associated with chiral symmetry restoration and deconfinement (see section 5.11).

Charmonium suppression

One telltale signature for the formation of quark-gluon plasma is the suppression of charmonia (J/Ψ, ϒ, etc). While the number of heavy quark pairs created in initial hard interactions is deemed to be independent of the final-state nuclear medium by virtue of QCD factorization (see section 4.2), they are predicted to be unable to form bound meson states in a quark-gluon plasma due to color screening. Lattice calculations for the heavy quark potential and the spectral function predict a sequential melting of charmonia states with increasing medium temperature. Therefore, the measurement of the ratio of charm spectra in nucleus-nucleus collisions divided by the spectrum expected from scaled proton–proton collisions is thought to serve as a direct way to extract the deconfinement temperature. This simplistic picture is complicated due to recombination of heavy quarks during hadronization and other effects (see section 5.10).

Heavy flavor

The mass dependence of medium interactions can be probed by measuring the spectra of hadrons containing unpaired c or b-quarks (open charm or beauty hadrons) or those of single leptons (electrons or muons) that originate from the decay of these hadrons. High statistics measurements of D-mesons identified by their hadronic decays have become the most precise tool for this purpose. Due to the dead-cone effect the expectation is that energetic charm and bottom quarks lose less energy in the quark-gluon plasma than light quarks. Measurements of the azimuthal angular distribution of open heavy-flavor mesons can provide insight into the question whether heavy quarks thermalize similar to light quarks (see section 5.10).

Two- and few-particle correlations are coincidence measurements of the probability of finding a particle with a certain property given that the event contains one particle with another (or the same) property.

Photon-hadron correlations

Photon-hadron correlations use a direct photon as the trigger particle for the correlation function. This has the advantage that the energy of the original parton in the hard process is constrained more tightly by the photon energy and, therefore, allows to access information about the fragmentation function of energetic partons, which measures the fraction z of the primary parton momentum carried by a secondary hadron. This observable is quite clean for theoretical calculations, but poses a big challenge to experiments due to the contamination of the photon yield by decay photons (mainly from π⁰) (see section 5.9).

Di-hadron correlations

Hard di-hadron correlations are a sensitive tool for studying the path length dependence of jet energy loss. The trigger particle is a high-p_T particle (usually with p_T > 5 GeV/c), whereas the associate particle is chosen in a lower p_T range. Di-hadron correlations are called 'hard', if both p_T ranges are 3 GeV/c and higher, so that a perturbative QCD (pQCD) description is applicable. These correlations are often quantified as a distribution in the difference of the azimuthal angle between trigger and associated particle. One can then separate the near-side and the away-side structures, by an angle difference of roughly 180° that originate from the back-to-back hard parton scattering. Given a near-side high-p_T trigger particle the ratio of the yield of particles on the away-side in heavy ion collisions to the one in proton–proton collisions is called I_AA in analogy to R_AA. Since I_AA is a conditional yield, it potentially contains more detailed information about the energy loss, but also is subject to the so-called 'trigger bias' (see section 5.9).

Triggered correlations

In addition to looking for correlations between two high-p_T particles, jet-medium interactions can be investigated by studying correlations between only one high-p_T particle and an associated particle of lower p_T . The high-p_T trigger particle preferentially emerges from the surface of the medium, having lost little energy on its way out, and therefore the backwards emitted parton experiences above average modification by the medium. Also, the medium itself might be modified by the traversing high energy parton resulting in characteristic excitation patterns like a Mach cone. To measure the effect of high p_T partons on the medium, multi-particle correlations over various p_T ranges and in longitudinal (Δη) and azimuthal (Δϕ) phase space have been studied.

Untriggered correlations

The last option for two- or three-particle angular correlations is to measure so-called untriggered correlations in minimum bias events, where one does not require one of the hadrons to have high transverse momentum. These d² N/dΔϕΔη correlations have revealed interesting structure: there is an elongated enhancement in pseudorapidity on the near-side Δϕ ≈ 0) that is often referred to as the 'ridge'. This structure is associated with higher Fourier coefficients in the azimuthal distribution with respect to the reaction plane (see below under collective flow) and will be discussed further in section 5.5.

Charge asymmetry correlations

Another 3-particle measurement is the charge asymmetry with respect to the reaction plane. Because up and down quarks carry different electric charges (+2e/3 and −e/3, respectively) they could be separated by the very strong transient magnetic field that is generated by the colliding nuclei. The chiral magnetic effect (see section 5.13 for more details) is predicted to result in a charge asymmetry of the yield of charged particles with respect to the reaction plane. This asymmetry can be detected by observables of the type:

$$\begin{eqnarray}&&{\gamma }^{\pm ,\pm }=\langle \cos ({\phi }_{\alpha }^{\pm }+{\phi }_{\beta }^{\pm }-2{{\rm{\Psi }}}_{\mathrm{RP}})\rangle ,\end{eqnarray} \tag{ 8 }$$

where ϕ_α , ϕ_β denotes the azimuthal emission angles of the two particles, Ψ_RP denotes the reaction plane defined by the impact parameter vector, and ± indicates the sign of the electric charge of the detected hadron. The observable (8) requires the detection of at least three particles, one of which is used to identify the reaction plane. A similar asymmetry of charged particle production with respect to the reaction plane can also be generated by a combination of balance functions and elliptic flow (see section 5.14), which means that the observable has a large non-specific background.

This section contains all of the observables that are based on multi-particle measurements and potentially reflect the collective behavior of the system.

Anisotropic flow

One of the most important observables in heavy ion reactions is collective flow of the particles. The Fourier coefficients v_n of the azimuthal distribution of the final particles in momentum space are used to quantify anisotropic flow.

$$\begin{eqnarray}&&E\frac{{{\rm{d}}}^{3}N}{{\rm{d}}{p}^{3}}=\frac{1}{2\pi }\frac{{{\rm{d}}}^{2}N}{{p}_{t}{\rm{d}}{p}_{t}{\rm{d}}y}+\displaystyle \sum _{n=1}^{\infty }2{v}_{n}\cos [n(\phi -{{\rm{\Psi }}}_{n})],\end{eqnarray} \tag{ 9 }$$

whee ϕ denotes the azimuthal angle around the beam axis, and Ψ_n stands for the angle defined by the maximum of the nth Fourier component of the angular distribution in a given event.

The most prominent and well-studied of these coefficients is v₂, commonly called elliptic flow, which arises due to the different pressure gradients in the transverse plane in non-central collisions. The initial almond shape in the coordinate space is translated by hydrodynamics into a momentum space anisotropy. As it is driven by pressure gradients, elliptic flow is highly sensitive to the shear viscosity and the equation of state of the expanding medium. Its collision energy dependence can, therefore, serve as a signature of the transition from hadronic matter to the quark-gluon plasma. Whereas v₂ can be measured by an event average, higher-order flow coefficients are non-zero require event-by-event measurements. They can be used to constrain the initial state profile of the fireball. The v_n values can be determined by different analysis techniques like the event-plane method, the cumulant method, the flow vector method, or by means of Lee-Yang zeros (see section 5.4).

HBT femtoscopic radii

Hanbury-Brown–Twiss (HBT) correlations of identical particles allow to infer the spatial volume, lifetime, and outward flow velocity of the fireball. The HBT correlations are caused by quantum (Bose or Fermi) interference between the two identical particles. The measurements are sometimes referred to as density interferometry. The measured correlation functions are usually fitted with a Gaussian that allows to extract three different radii

$$\begin{eqnarray}&&C({\bf{q}})=1+\lambda {{\rm{e}}}^{-({q}_{o}^{2}{R}_{\mathrm{out}}^{2}+{q}_{s}^{2}{R}_{\mathrm{side}}^{2}+{q}_{l}^{2}{R}_{\mathrm{long}}^{2})},\end{eqnarray} \tag{ 10 }$$

where R_long is aligned with the beam direction, R_out points along the center-of-momentum of the particle pair, and R_side is perpendicular to both. The vector q denotes the difference between the momenta of the two identical particles, and q_o , q_s , q_l are its components with respect to these three orthogonal directions. The HBT radii for different particle species contain information about the coherence region of the emission at kinetic freeze-out. The ratio of R_out/R_side is indicative of the lifetime of the fireball and therefore sensitive to certain aspects of the equation of state, such as the nature of the QCD phase transition. Results and more details are discussed in section 5.3.

Balance functions

Balance functions are a tool to study charge correlations in heavy ion reactions. They are similar to observables used to investigate hadronization in jets produced in $p\bar{p}$ or e⁺ e⁻ collisions. The balance function describes the conditional probability that a particle in the momentum bin p₁ will be accompanied by a particle of opposite charge in the momentum bin p₂

$$\begin{eqnarray}\begin{array}{rcl}{B}_{{ab}}({p}_{2}| {p}_{1}) & = & \displaystyle \frac{1}{2}\left(\rho (b,{p}_{2}| a,{p}_{1})-\rho (b,{p}_{2}| b,{p}_{1})\right.\\ & & \left.+\rho (a,{p}_{2}| b,{p}_{1})-\rho (a,{p}_{2}| a,{p}_{1}\right),\end{array}\end{eqnarray} \tag{ 11 }$$

where ρ(b, p₂∣a, p₁) is the conditional probability of observing a particle of type b in bin p₂ given the existence of a particle of type a in bin p₁. Balance functions can be defined for any conserved quantum number.

Reconstructed jets

Another observable sensitive to parton propagation in the quark-gluon plasma are fully reconstructed jets. For p+p collisions with usually only a single hard interaction the outgoing hard partons fragment into hadrons, and the energies of all fragment particles are measured in calorimeters. Applying sophisticated clustering algorithms the full jet can be reconstructed. Events containing a jet are identified online using a dedicated jet trigger. In the high-multiplicity environment of a heavy ion collision these measurements are much more challenging due to the underlying event background. Reconstructed jets allow comprehensive studies of jet energy loss and modifications of the internal jet structure (fragmentation functions, jet shape, etc) by the comparison with predictions from jet shower Monte-Carlo algorithms tuned to jets in p+p collisions (see section 5.9).

Energy flow and track functions

Calorimeters measure energy flow carried by hadrons (hadronic calorimeter) or by photons and electrons (electromagnetic calorimeter) in a given direction. The most general realization of this concept is the energy flow operator [32]

$$\begin{eqnarray}&&{ \mathcal E }(\vec{n})=\mathop{\mathrm{lim}}\limits_{r\to \infty }{\int }_{0}^{\infty }{\rm{d}}t\,{r}^{2}{n}^{i}{T}_{i}^{0}(t,r\vec{n}),\end{eqnarray} \tag{ 12 }$$

which measures energy flow into the direction $\vec{n}$. The T⁰ _i are the components of the energy-momentum tensor ${T}_{\nu }^{\mu }$ related to momentum or energy flow. (For photons these are the components of the Poynting vector.) Observables of interest that provide a measure of the energy correlators of the form $\langle { \mathcal E }({\vec{n}}_{1})\cdots { \mathcal E }({\vec{n}}_{N})\rangle$, where one looks at the correlation of energy flow within a jet into different directional domains with an opening solid angle ΔR [33].

A related concept is that of track functions T_q→h (x), [34, 35] which measure the probability that the total (light-cone) momentum fraction x of a quark-initiated jet is carried by a certain type of hadron h. The track functions are similar to fragmentation functions, but instead of considering the momentum fraction of a single hadron, they sum over all hadrons of a given type within the jet. For a recent analysis of the QCD evolution of track function moments and a comparison with LHC data, see [36].

Event-by-event fluctuations

Event-by-event fluctuations of conserved quantum numbers, such as net electrical charge Q, net baryon number B, or net strangeness S, are a prominent signal for the phase transition to the quark-gluon plasma. Since quarks and gluons have different elementary units of the charges the predictions for event-by-event fluctuations are very different based on the assumption that a quark-gluon plasma or a hadron gas is present. The fluctuations of the mean transverse momentum can be regarded as the analog to the temperature fluctuations in the cosmic microwave background, which are vestiges of quantum fluctuations in the initial state. The susceptibilities that are directly related to the fluctuations of conserved quantum numbers or their higher moments such as the skewness and the kurtosis are also calculable on the lattice and offer a direct connection between data and QCD predictions (see section 5.7).

The quantitative ab initio calculation of the thermodynamic properties of QCD matter requires the model-independent evaluation of the functional integral that defines the quantum field theory. The only known rigorous method that allows to do this starts with the discretization of the QCD Lagrangian on a large space-time lattice and then evaluates the very high, but finite-dimensional functional integral by Monte-Carlo methods. In order to study the properties of QCD at non-zero temperature, periodic boundary conditions (anti-periodic for quarks) are imposed on the imaginary time coordinate with period ℏ/T. The functional integral defining the quantum field theory is evaluated by means of Monte-Carlo sampling techniques that are used to generate a sufficient number (hundreds or thousands) of statistically independent, representative field configurations. Observables of interest are then evaluated by averaging over these stored configurations. The numerically most expensive part of the calculation is the integration over the fermion fields, which requires the evaluation of the determinant of a very large matrix. Much of the progress that has occurred over the past decade consists of finding improved ways to represent the QCD Lagrangian on a discrete lattice and speed up the evaluation of the fermion determinant. Other improvements in some simulations concern the implementation of manifest chiral symmetry by means of the domain wall or overlap fermion algorithms [37, 38].

Lattice QCD simulations have come a long way since the first calculations that demonstrated quark liberation in the high-temperature phase of pure SU(2) lattice gauge theory [39, 40]. State-of-the-art lattice calculations include both, SU(3) gauge fields and dynamical u, d, and s quarks (sometimes even c-quarks) with physical masses, employ improved lattice actions, extrapolate to the continuum limit (lattice distance a → 0), and investigate the convergence of the results as a function of the number N_t of Euclidean time slices. Fully converged lattice QCD simulations for baryon-symmetric QCD matter (μ_B = 0) have been made for the QCD equation of state and various other thermodynamic properties [41–43].

Rigorous algorithms for the evaluation of the functional integral for thermal QCD are known only for μ_B = 0, because the integrand is not positive definite for non-zero real values of μ_B . Nevertheless, attempts have been made to explore the QCD phase diagram away from the μ_B = 0 axis by means of lattice simulations using a number of indirect approaches. One method evaluates a certain number of derivatives of an observable with respect to μ_B at μ_B = 0 and then reconstructs the value at μ_B ≠ 0 from the truncated Taylor series [44]. Another method evaluates the functional integral for imaginary values of μ_B , where the integrand is well-behaved, and then attempts an analytic continuation to real values of μ_B [45]. A third method uses field configurations generated for μ_B = 0, but reweights them to simulate the μ_B -dependence of the lattice action [46]. The 'holy grail' of such efforts is to locate a possible critical point in the QCD phase diagram. In spite of considerable progress in obtaining reliable results for μ_B /T < 2, none of these approaches currently gives a positive indication for the possible location of such a critical point. Even the existence of a critical point in the T − μ_B plane, although extremely plausible on the basis of physical arguments and predicted by a multitude of QCD models, has not been established with absolute certainty [47].

3.1.1. Pseudocritical temperature

It has now been firmly established that the transition in baryon-antibaryon symmetric (μ_B = 0) QCD matter is a smooth cross-over [3]. As a consequence, the value of the pseudocritical transition temperature T_c cannot be defined unambiguously. Usually, one identifies T_c as the location of the peak or inflection point of a relevant thermodynamic quantity, but even this definition is ambiguous, because the value of T_c depends upon the choice of this quantity. The transition temperature defined by the inflection point of the renormalized quark condensate Δ_l,s [48] is ${T}_{c}^{(\chi )}\approx 158\,\mathrm{MeV}$ [49] with an uncertainty of less than 1 MeV (see figure 2).

Zoom In Zoom Out Reset image size

On the other hand, the transition region of quantities that are more direct measures of quark liberation, such as the expectation value of the Polyakov loop 〈L〉 or the strange quark number susceptibility is considerably broader. The expectation value of the Polyakov loop measures the interaction energy E_Q of an isolated heavy, static quark with the thermal bath of gluons via the relation ${E}_{Q}=-T\mathrm{ln}\langle L\rangle$. In the absence of light, dynamical quarks, it vanishes in the confined phase of QCD, which implies that an isolated quark would have infinite energy. At high temperature, in the deconfined phase, the Polyakov loop approaches unity because the interaction energy is small compared with the temperature. This behavior is illustrated in figure 3, which shows the temperature dependence of the strange quark susceptibility (left panel) and of the renormalized Polyakov loop (right panel). Both observables yield values of the pseudocritical temperature, defined by the location of the inflection point, of ${T}_{c}(L)\approx {T}_{c}^{(s)}\approx 175\,\mathrm{MeV}$. (The higher value may be simply due to the greater width of the transition for these quantities.) The transition occurs over a range of about 20 MeV for the strange quark condensate and an even wider range for the Polyakov loop, which indicates that quark deconfinement is occurring gradually as the quark-gluon plasma is heated beyond T_c .

Zoom In Zoom Out Reset image size

The transition is seen first, and over the narrowest interval, in the chiral quark condensate. This suggests that for the physical values of the parameters of QCD (the confinement scale Λ_QCD and the quark masses m_u , m_d , m_s ) the transition is driven by the restoration of approximate chiral symmetry for the light quarks. The broad range of the transition seen in the strange quark susceptibility indicates that strange quark degrees of freedom thaw more gradually and at slightly higher temperature. The even more gradual transition of the renormalized Polyakov loop may be an indication that remnants of quark confinement in terms of correlations among quarks and antiquarks persist over a larger temperature range. This behavior is also evident in the QCD equation of state, which we discuss next.

3.1.2. Equation of state

The QCD equation of state at μ_B = 0 exhibits a rapid rise in the ratio of all thermodynamic quantities to the Stefan–Boltzmann (SB) limit for massless quanta over the temperature range 100 MeV < T < 500 MeV, as shown in figure 4 [48]. Perhaps the clearest manifestation of this behavior is in the effective number of degrees of freedom d_eff defined via the formula for the entropy density s(T) of a massless, noninteracting gas of particles:

$$\begin{eqnarray}&&s(T)=\displaystyle \frac{45}{2{\pi }^{2}}\,{d}_{\mathrm{eff}}{T}^{3}.\end{eqnarray} \tag{ 13 }$$

As seen in the right panel of figure 4, the effective number of degrees of freedom in hot QCD matter rises steadily over this temperature range from d_eff ≈ 8.5 at T = 140 MeV to d_eff ≈ 38 at T = 500 MeV, which corresponds to 80% of the Stefan–Boltzmann limit. The remaining deviation from the SB limit can be understood as the effect of perturbative color interactions in the quark-gluon plasma.

Zoom In Zoom Out Reset image size

Effective descriptions of the low-temperature and high-temperature domains are indicated in the left panel of figure 4. At low temperatures, the hadron resonance gas (HRG) provides a rather good description, although in detail the equation of state indicates the existence of hadron resonances beyond those listed in the Particle Data Book [50, 51]. The HRG gives a good description for most quantities up to T ≈ 150 MeV. Some especially sensitive quantities may not only depend on the hadron mass spectrum, but also on the widths of resonances [52].

At high temperature (T > 350 MeV), hard-thermal loop resummed perturbation theory gives a good description of the equation of state and related observables. As indicated in the left panel of figure 4 the next-to-next-to-leading order (NNLO) calculations in this scheme [53, 54] match the lattice results for the pressure quite well, albeit with a large uncertainty deriving from the choice of the renormalization scale μ. The three dash-dotted lines shown in the figure correspond to μ = π T (top), 2π T (center), 4π T (bottom); the central curve obviously matches best.

An especially interesting feature is the broad peak in the QCD interaction measure, defined as I(T) = (ε − 3P)/T⁴, shown in the left panel of figure 5, which measures the manifest violation of scale invariance in the QCD equation of state. The nonperturbative 'interaction measure' I(T) is distinct from the perturbative violation of scale invariance in the light quark and gauge sector of QCD that derives from the temperature dependence of the running coupling constant α_s (T). There have been multiple attempts to interpret this nonperturbative feature of hot QCD in terms of some quasi-particle picture [55, 56], but the strong coupling among the degrees of freedom has made it difficult to draw definite conclusions. The scale anomaly remains large throughout the entire temperature domain (T ≤ 400 MeV) that is probed in relativistic heavy ion collisions.

Zoom In Zoom Out Reset image size

For T < T_c , the growth of I(T) can be understood in terms of the hadron resonance gas model, which predicts that more massive states are being excited with increasing temperature. As the figure shows, the lattice result slightly exceeds the resonance gas model, which accounts for all experimentally known resonances, below 175 MeV. A possible explanation for this excess is the excitation of a large number of unknown resonances [50]. Although the microscopic structure of these states is not known, many of them must contain internal gluon excitations, either as hybrid states (hadrons that contain both valence quarks and valence gluons) or as glueballs (hadrons composed solely of valence gluons). If this is so, the question of what happens to the gluons when the quark-gluon plasma hadronizes, has a simple resolution: the thermally excited gluons become internal gluonic excitations of short-lived hadrons. As the hadron gas expands, these internal excitations will quickly decay into light hadrons, mostly pions, which absorb the additional entropy.

Lattice gauge theory currently does not permit to perform reliable computations of real-time processes. For this reason, insight into transport processes in the hot QCD matter presently relies either on thermal perturbation theory or on QCD-inspired models of strongly coupled gauge theories with a holographic gravity dual. Holographic models will be discussed in the next subsection; here we focus on thermal perturbation theory.

Thermal perturbation theory is usually formulated in terms of the hard thermal loop (HTL) effective theory [57–59]. The HTL formalism resumes the leading thermal contributions to the gluon and quark self-energies into momentum-dependent effective masses. The effective theory relies on the separation of the scales T ≫ gT ≫ g² T, where g ≡ g(T) is the QCD coupling constant at the thermal scale. The scale separation requires weak coupling (formally g ≪ 1 corresponding α_s = g²/4π ≪ 0.1). QCD is not so weakly coupled at any physically relevant temperature, but it turns out that the power counting implied by the HTL scheme appears to work even when the scale ordering is not strictly realized, as it often is the case for effective field theories.

Gauge invariance dictates that the thermal modifications of the propagators are balanced by vertex corrections that satisfy the Ward-Takahashi identities. In the HTL effective theory the gauge field develops a collective longitudinal mode, the plasmon, with mass

$$\begin{eqnarray}&&{m}_{\mathrm{pl}}^{2}=\displaystyle \frac{2{N}_{c}+{N}_{f}}{6}{\left({gT}\right)}^{2}+O({g}^{4}{T}^{2}),\end{eqnarray} \tag{ 14 }$$

for a plasmon at rest in the medium. Here N_c and N_f denote the number of colors and light quark flavors, respectively. m_pl represents the additional mass scale that enters into the HTL effective Lagrangian and controls the thermal properties of the gauge theory. For example, the HTL formalism predicts that the static chromo-electric field around a heavy color charge in a QGP with N_f light quark flavors is screened with a Debye mass ${m}_{{\rm{D}}}=\sqrt{3}{m}_{\mathrm{pl}}={gT}\sqrt{1+{N}_{f}/6}$. Furthermore, all HTL gauge field modes are Landau damped in the space-like domain at the same scale.

An estimate of the domain of applicability of HTL perturbation theory can be obtained by calculating and comparing quantities that can be computed reliably on the lattice. An example is the QCD equation of state, which has been calculated up to three-loop (NNLO) order for both, μ_B = 0 [53, 54] and μ_B ≠ 0 [60], in the resummed HTL formalism and found to agree well with the lattice results for T > 350 − 400 MeV (see figure 4). It is also clear from these results that the HTL effective theory has very limited applicability in the temperature range T < 2T_c that is most relevant to the heavy ion experiments.

Many properties of the quark-gluon plasma have been calculated in the HTL effective theory, mostly at leading order, but NLO results have increasingly become available. Examples include the shear viscosity η [61–63], the so-called jet quenching parameter $\hat{q}$ [64–66], the rate of collisional energy loss of a fast parton $\hat{e}$ [67], the heavy-quark diffusion constant D [68] and the rates of thermal photon [69, 70] and lepton-pair [71–73] production.

The discovery of the Anti-de Sitter/Conformal Field Theory (AdS/CFT) duality by Maldacena in 1997 fundamentally changed theorists' ability to perform exact calculations for strongly coupled gauge theories [24, 74–76]. The AdS/CFT duality is based on the notion that conformal field theories in (d + 1)-dimensional space-time can be mapped onto quantum gravity in (d + 2)-dimensional Anti-de Sitter space (AdS_d+2), which has (d + 1)-dimensional space-time as a boundary at infinity. The AdS space is known as the 'bulk'. Because the mapping is between theories defined on space-times with a different number of dimensions, one speaks of a holographic duality.

The most ubiquitous case for relativistic heavy ion physics is the original duality between ${ \mathcal N }=4$ super-Yang-Mills theory in (3 + 1) dimensions, which can be considered as a 'cousin' of QCD, and superstring theory on the space AdS₅ × S₅. In the limit of an infinite number of colors N_c , expressed as the limit of strong 't Hooft coupling (λ = g² N_c → ∞), the holographic dual is reduced to (super-)gravity on AdS₅. While QCD is not conformally symmetric—as we already mentioned, the violation of scale invariance is quite large in the temperature range of greatest interest to relativistic heavy ion physics—there are a number of modifications of the original holographic model that incorporate scale symmetry breaking in a fundamental (see [77] for an overview) or phenomenological way (see [78, 79]).

An especially attractive feature of AdS/CFT duality is that it enables rigorous calculations of dynamical processes at strong coupling, which is not possible in the framework of Euclidean lattice gauge theory [80]. In particular, it is possible to study thermalization at strong coupling by mapping thermalization onto the process of formation of a black hole in the bulk [81, 82], as illustrated in the left panel of figure 6.

Zoom In Zoom Out Reset image size

The dynamics of space-time near the event horizon of black holes have long been known to be governed by viscous relativistic hydrodynamics, a phenomenon known as the membrane paradigm [83]. In the context of AdS space with a black hole, this translates into properties of viscous hydrodynamics in the gauge theory on the AdS boundary, some of which are deemed universal for any strongly coupled nonabelian gauge theory [84, 85]. The best known of these properties is the ratio of the shear viscosity to the entropy density, η/s = 1/(4π) characteristic of a 'perfect' fluid.

Other experimentally relevant properties of strongly coupled gauge theories have also been calculated using holographic techniques, e.g. the energy loss of heavy quarks [87] and light quarks [88, 89]. Quite generally, holographic approaches lend themselves most easily to the study of energy and momentum transport as these are directly mapped onto the dynamics of the gravitational field in the bulk. Modeling of other processes associated with quark flavor or spin is possible, but requires adding additional fields or geometric objects, such as D-branes, to the model thereby rendering the conclusions less universally valid.

The holographic approach makes it possible to model the approach to hydrodynamics from arbitrary initial conditions without arbitrary approximations. Under a wide range of initial conditions, boost invariant (Bjorken) hydrodynamics with η/s = 1/(4π) emerges as asymptotic description of the dynamical evolution [88, 90–92]. The complete process of energy deposition, thermalization, and hydrodynamic expansion can be modeled in terms of collisions of two gravitational shock waves. Such shock wave collisions in AdS space can be rigorously solved using numerical techniques [93–95]. By varying the height and width of the shock waves the whole range from transparency to full stopping [96] including asymmetric collisions [97, 98] can be explored. By adding a scalar field to AdS gravity, nonconformal models can be studied, which leads to hydrodynamic behavior including bulk viscosity [99]. When one combines holographic modeling with other descriptions that are more appropriate for the late stages of a nuclear collision, a seamless end-to-end description of heavy ion collisions can be implemented [100].

The dynamical evolution of heavy-ion collisions can be divided into several stages (see figure 7) that are described in detail in the following subsections. Initially, there are two Lorentz contracted nuclei approaching each other. Right after the impact there is a period characterized by non-equilibrium evolution. Eventually relativistic viscous fluid dynamics becomes applicable for the hot and dense stage. When the system dilutes enough, the fluid is converted to hadrons, which decay and interact with each other until the freeze-out.

Zoom In Zoom Out Reset image size

4.1.1. Initial conditions

The initial conditions are an indispensable component of any dynamical description of heavy ion reactions. We have rather extensive knowledge of the structure of the colliding nuclei, but we lack a first-principle treatment on the basis of QCD of the evolution from the colliding nuclear ground states to the initial conditions that are used as input for hydrodynamics.

The simplest picture for the initial geometry of the overlap region in the transverse plane is based on the Glauber model. While the original Glauber model provides a general formalism for multiple scattering, its application to relativistic heavy ion physics commonly makes use of the eikonal limit, where all particles are assumed to travel along the straight line trajectories. (For a simplified introduction to the formalism, see [101]; for an overview of its various uses in relativistic heavy ion physics, see [102].) In this eikonalized Glauber model, the two nuclei are described by thickness functions ${T}_{A}(\vec{s})$, which represent Woods-Saxon density profiles denoted by ρ_A that are integrated over the longitudinal direction: ${\hat{T}}_{A}(\vec{s})=\int {\rho }_{A}(\vec{s},z){\rm{d}}z$. $\vec{s}$ is the transverse position with respect to the center-of-mass of the nucleus A. The joint probability for the overlap of nucleus A with nucleus B colliding at an impact parameter $\vec{b}$ can be expressed as

$$\begin{eqnarray}&&{\hat{T}}_{{AB}}(\vec{b})=\int {\hat{T}}_{A}(\vec{s}){\hat{T}}_{B}(\vec{s}-\vec{b}){{\rm{d}}}^{2}s.\end{eqnarray} \tag{ 15 }$$

Participant or 'wounded' nucleons are defined as the nucleons that interact with at least one nucleon of the opposite nucleus, where the interaction probability is obtained from the inelastic nucleon–nucleon cross section (σ_NN) at a given energy:

$$\begin{eqnarray}&&{N}_{\mathrm{part}}(\vec{b})=A\int {\hat{T}}_{A}(\vec{s})\left\{1-{\left[1-{\hat{T}}_{B}(\vec{s}-\vec{b}){\sigma }_{\mathrm{NN}}\right]}^{B}\right\}{{\rm{d}}}^{2}s\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\,+B\int {\hat{T}}_{B}(\vec{s}-\vec{b})\left\{1-{\left[1-{\hat{T}}_{A}(\vec{s}){\sigma }_{\mathrm{NN}}\right]}^{A}\right\}{{\rm{d}}}^{2}s,\end{eqnarray} \tag{ 17 }$$

where A and B denote the nuclear mass numbers. The nucleons outside the interaction region are called spectators and are not considered further. The number of binary nucleon–nucleon collisions (N_coll) can be calculated as well.

$$\begin{eqnarray}&&{N}_{\mathrm{coll}}(\vec{b})=A\,B\,{\hat{T}}_{{AB}}(\vec{b}){\sigma }_{\mathrm{NN}}.\end{eqnarray} \tag{ 18 }$$

In some heuristic models [103] of the centrality dependence of the charged particle multiplicity binary collisions are assumed to contribute a fraction α < 1 and while a fraction (1 − α) is attributed the participants:

$$\begin{eqnarray}&&{N}_{{AA}}(\vec{b})=\left[\alpha {N}_{\mathrm{coll}}(\vec{b})+\displaystyle \frac{1-\alpha }{2}{N}_{\mathrm{part}}(\vec{b})\right]{N}_{{pp}},\end{eqnarray} \tag{ 19 }$$

where N_pp is the number of particles produced in a proton–proton collision at the same energy. The Glauber model is usually implemented as a Monte-Carlo process and therefore provides different initial condition profiles for individual events. Figure 8 shows the distribution of participant nucleons (solid circles) and spectator nucleons (dotted circles) in a randomly chosen Au+Au collision event at top RHIC energy. Many improvements of this simplistic scheme have been studied over the years. These include treatments of nucleon–nucleon correlations in the initial state [104], the finite extent and internal structure of nucleons, fluctuations in the energy deposition by individual nucleon–nucleon collision, and much more (see [105] for a recent review).

Zoom In Zoom Out Reset image size

The Color Glass Condensate (CGC) model [107] aims to describe the colliding nuclei directly in terms of their quark-gluon content and is thus limited to collisions at sufficiently high-energy. The CGC model assumes that energy deposition at midrapidity is dominated by the liberation of gluons from the colliding nuclei, and that the relevant gluon content of the nuclei can be described by semiclassical gauge fields. These fields are generated by sources (primarily valence quarks) that move at near-projectile rapidity and whose dynamics during the collision can be ignored because of time dilatation. In the simplest version of the CGC, the McLerran-Venugopalan model [108, 109], the colliding nuclei are modeled as Gaussian distribution of color currents ${J}^{a\pm }(x)={\rho }^{a}({\vec{x}}_{\perp })\delta ({x}^{\mp }){n}^{\pm }$ along the light-cone. The distribution of these color sources is given by

$$\begin{eqnarray}&&W[{\rho }^{a}]=N\exp \left(-\frac{1}{2{g}^{2}{\mu }^{2}}\int {{\rm{d}}}^{2}{x}_{\perp }{\left({\rho }^{a}({\vec{x}}_{\perp })\right)}^{2}\right),\end{eqnarray} \tag{ 20 }$$

where μ is a scale setting parameter that is related to the saturation scale Q_s ≈ 0.6g² μ [110]. The semiclassical color fields are obtained as solutions of the classical Yang-Mills equation ${D}_{\mu }^{{ab}}{F}^{a\mu \nu }={{gJ}}^{a\nu }$, where D^ab is the gauge covariant derivative and J^{a μ} is the color current due to quarks.

The classical approximation is only warranted when the occupation number of gauge field modes is of the order of ${\alpha }_{s}{\left({Q}_{s}\right)}^{-1}\gg 1$. This condition requires that Q_s ≫ Λ_QCD. Because the saturation scale Q_s scales as ${Q}_{s}^{2}\sim {A}^{1/3}{x}^{-\lambda }$ with λ ≈ 0.3, the condition can be fulfilled either for a large nuclear mass number A ≫ 1 or for very small momentum fraction x ≪ 1. Figure 9 gives a graphical representation of this scaling. Small values of x are only accessible at high collision energies, because the transverse momentum p_T of a produced hard probe is given by the relation ${p}_{T}^{2}={x}_{1}{x}_{2}{s}_{\mathrm{NN}}$, where x₁ and x₂ are the momentum fractions of the participating partons in the two colliding nuclei and $\sqrt{{s}_{\mathrm{NN}}}$ is the center-of-mass collision energy per nucleon pair.

Zoom In Zoom Out Reset image size

The dependence of the color sources ρ^a (x_⊥) on the momentum fraction x is governed by renormalization group (RG) evolution. At a given value of x₀ all sources of color fields residing at x > x₀ contribute (valence quarks for all x₀ and gluons and sea quarks for small values of x₀). The RG evolution can either be expressed as a differential equation for the functional W[ρ^a ], the so-called JIMWLK equation [111], or as a nonlinear differential equation for the so-called dipole density N_xy , the Balitsky-Kovchegov (BK) equation [112]. Whereas the scaling with A and x can be calculated from geometry and QCD renormalization group arguments, the overall scale is nonperturbative. This means that the range of applicability of the CGC model is still under debate and awaits precise data from the future electron-ion collider (EIC).

The color fields generated in this way are virtual and remain an integral part of the nuclei until an interaction occurs. During the collision of two heavy ion the CGC fields of the nuclei interact, and a large fraction of the virtual gluons are scattered on-shell, as sketched in figure 10. The resulting nonequilibrium state is known as glasma. Because the CGC color fields are transversely polarized, the resulting energy-momentum tensor is far off equilibrium and cannot be directly used as an initial condition for the hydrodynamic evolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to describe individual heavy ion collision events, it is necessary to not only model the average energy deposition, but also the event-by-event fluctuations. There are two main sources of such initial-state fluctuations. One are the quantum fluctuations in the distribution of nucleons in the colliding nuclei or, if one instead looks at the collision at the parton level, in the distribution of valence quarks in the nuclei. The other source of fluctuations are the quantum fluctuations of particle production that are encoded in the S-matrix of a nucleon–nucleon or quark-quark collision. In the standard Glauber model, whether it is implemented at the nucleon or the valence quark level, the second type of fluctuations is commonly parametrized by a negative binomial distribution with the generating function [113]

$$\begin{eqnarray}&&F(t;k)=\frac{{k}^{k}}{{\rm{\Gamma }}(k)}{t}^{k-1}{{\rm{e}}}^{-{kt}},\end{eqnarray} \tag{ 21 }$$

where t is a continuous variable and the parameter k (with 1 < k < 4) depends on the energy. The probability distribution of the number of produced particles is given by

$$\begin{eqnarray}&&P(n)={\int }_{0}^{\infty }{\rm{d}}t\,F(t;k){{\rm{e}}}^{-\bar{n}t}\frac{{\left(\bar{n}t\right)}^{n}}{n!},\end{eqnarray} \tag{ 22 }$$

where $\bar{n}$ is the average number of produced particles.

For the CGC model, one can obtain fluctuating initial conditions using a formalism similar to the Glauber model, where the positions of nucleons are chosen randomly from the nuclear density distribution and the color sources are distributed within the individual nucleons. This implementation of the CGC model is known as impact parameter saturation (IP-sat) model [114]. When applied to nuclear collisions, the IP-sat model generates fluctuating initial conditions for the glasma (IP-glasma), which naturally result in a negative binomial distribution for the initial energy density [115, 116]. The energy densities obtained with these two alternative approaches are depicted in figure 11, where the nucleon-based Glauber model (left panel) is seen to result in a much smoother initial distribution than the parton-based IP-sat model (right panel).

The initial pattern of energy or entropy deposition in these models does not lend itself directly to be used as initial conditions for viscous hydrodynamics, because they do not correspond to an energy-momentum tensor near thermal equilibrium and the spatial gradients are too large in the case of event-by-event fluctuations. In the case of the standard Glauber model, this problem is commonly resolved by a short period of free streaming evolution, which smoothes out the spatial gradients. Alternatively, one widely used version of a general Glauber-type initial condition generator, the TRENTO model [117], uses the entropy density to map the initial conditions generated by the model directly onto the initial conditions for viscous hydrodynamics.

For CGC initial conditions, the initial time evolution of the glasma is described by the nonlinear classical Yang-Mills equation; the resulting energy-momentum tensor is then mapped onto initial conditions for hydrodynamics [116]. The microscopic thermalization process is thought to follow a scenario known as 'bottom-up' thermalization [118]. In this scenario semi-hard gluons contained in the colliding nuclei radiate soft gluons which form a thermal bath that drains energy from the semi-hard gluons via elastic scattering and thereby thermalizes them. An effective kinetic description of this thermalization process has recently been developed and applied to the initial condition of hydrodynamics in relativistic heavy ion collisions [119]. There are many additional details and possible refinements of this process, which are under investigation, but no generally accepted approach has emerged yet.

4.1.2. The QGP phase

Due to asymptotic freedom, the quark-gluon plasma was originally expected to behave as a weakly interacting gas with almost massless degrees of freedom, the quarks and gluons, for which collective flow would be difficult to establish on the length scale of the fireball. However, when the first measurements of anisotropic flow in Au+Au collisions at $\sqrt{{s}_{\mathrm{NN}}}=200\,\mathrm{GeV}$ at RHIC were published, a surprisingly good agreement with predictions from ideal fluid dynamic calculations was observed [120, 121], if the initial conditions were matched to fit the overall particle production yields and transverse momentum spectra. The elliptic flow is then a natural consequence of collective behavior for a geometrically anisotropic initial state.

Fluid dynamics is based on the equations for the conservation laws

$$\begin{eqnarray}&&{\partial }_{\mu }{T}^{\mu \nu }=0,\quad {\rm{and}}\quad {\partial }_{\mu }{N}_{i}^{\mu }=0\end{eqnarray} \tag{ 23 }$$

of the energy-momentum tensor T^{μ ν} and the current density ${N}_{i}^{\mu }$ of any conserved quantum number (e.g. B, S, Q). If the fluid is in local equilibrium, the energy-momentum tensor of an ideal fluid takes the form

$$\begin{eqnarray}&&{T}^{\mu \nu }=(\epsilon +P){u}^{\mu }{u}^{\nu }+{{Pg}}^{\mu \nu },\end{eqnarray} \tag{ 24 }$$

with being the energy density, P the pressure and u^μ denoting the flow 4-velocity. g^{μ ν} is the metric tensor and typically chosen as (+, − , − , − ) in special relativistic applications. The quantum number current density has the simple form N^μ = nu^μ , where n is the density measured in the local rest frame.

In viscous hydrodynamics, the energy-momentum tensor is decomposed as

$$\begin{eqnarray}&&{T}^{\mu \nu }=\varepsilon {u}^{\mu }{u}^{\nu }-{{\rm{\Delta }}}^{\mu \nu }(p+{\rm{\Pi }})+{\pi }^{\mu \nu },\end{eqnarray} \tag{ 25 }$$

with Π denoting the bulk viscous contribution and π^{μ ν} the shear viscous tensor. First order viscous hydrodynamics is not causal, since it permits signal propagation that is faster than the speed of light. Therefore, second order viscous fluid dynamics must be applied with the following relaxation equations for shear and bulk viscosity including the coupling between the two

$$\begin{eqnarray}&&D{\rm{\Pi }}=\displaystyle \frac{-\zeta \theta -{\rm{\Pi }}}{{\tau }_{{\rm{\Pi }}}}-\displaystyle \frac{{\delta }_{{\rm{\Pi }}{\rm{\Pi }}}}{{\tau }_{{\rm{\Pi }}}}{\rm{\Pi }}\theta +\displaystyle \frac{{\lambda }_{{\rm{\Pi }}\pi }}{{\tau }_{{\rm{\Pi }}}}{\pi }^{\mu \nu }{\sigma }_{\mu \nu }\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}\begin{array}{rcl}D{\pi }^{\langle \mu \nu \rangle } & = & \displaystyle \frac{2\eta {\sigma }^{\mu \nu }-{\pi }^{\mu \nu }}{{\tau }_{\pi }}-\displaystyle \frac{{\delta }_{\pi \pi }}{{\tau }_{\pi }}{\pi }^{\mu \nu }\theta +\displaystyle \frac{{\phi }_{7}}{{\tau }_{\pi }}{\pi }_{\alpha }^{\langle \mu }{\pi }^{\nu \rangle \alpha }\\ & & -\displaystyle \frac{{\tau }_{\pi \pi }}{{\tau }_{\pi }}{\pi }_{\alpha }^{\langle \mu }{\sigma }^{\nu \rangle \alpha }+\displaystyle \frac{{\lambda }_{\pi {\rm{\Pi }}}}{{\tau }_{\pi }}{\rm{\Pi }}{\sigma }^{\mu \nu }.\end{array}\end{eqnarray} \tag{ 27 }$$

θ = ∂_μ u^μ is the expansion scalar, σ^{μ ν} is the strain tensor, τ_π and τ_Π are the relaxation times for the shear and bulk viscous corrections and λ_Ππ , λ_πΠ, δ_{π π} , τ_{π π} and ϕ₇ are higher-order couplings, whose forms are taken from [122].

There are two typical choices for the local rest frame, the Landau frame which moves with the energy density and the Eckart frame which moves with the conserved density, if one exists. More details on viscous relativistic hydrodynamics and its application to heavy ion physics can be found in [123–125].

Since the set of hydrodynamics equations listed above is insufficient to determine all independent components of the quantities of interest, one needs an additional equation as an input. This is the equation of state, which encodes the properties of the fluid under consideration. In the context of heavy ion physics this has the major benefit that hydrodynamics provides a controlled handle on the phase transition between a hadron gas and the quark-gluon plasma phase. One can study different equations of state to explore its influence on the dynamics of the system and various experimentally measurable observables.

The simplest equation of state is that of a relativistic ideal gas composed of massless particles:

$$\begin{eqnarray}&&P=\displaystyle \frac{\varepsilon }{3}.\end{eqnarray} \tag{ 28 }$$

At vanishing baryon chemical potential the equation of state is known from lattice QCD calculations and can be parametrized as originally done in [126]. Currently, there are many other possibilities including the extension to multiple conserved currents (see e.g. equations of state by the BEST collaboration [127], the BNL group [128] or from holography [129]). At finite density and vanishing temperatures there are constraints from the mass-radius relation of neutron stars. Also any equation of state should take the constraints based on stable nuclear matter into account.

At high collision energies, hydrodynamics predicts the evolution of the fireball at midrapidity to be approximately boost invariant [130]. The hydrodynamical variables then become functions of the transverse coordinates and the local proper time τ. The evolution of the fireball can be visualized by contour plots of the energy density, as shown in figure 12 for a Au+Au collision at top RHIC energy. At lower beam energies the framework has to be extended to include finite net baryon charges and three spatial dimensions, since the assumption of a boost invariant longitudinal expansion breaks down.

Zoom In Zoom Out Reset image size

The second input into hydrodynamics calculations, which contains information about the microscopic properties of the medium, are the transport coefficients. Main consideration has been given to the shear and bulk viscosity (η and ζ), but recently calculations have begun to include quantum number diffusion coefficients and even their cross-talk [132]. The most interesting dimensionless quantities are η/s and ζ/s, where s is the entropy density. It is extremely challenging to compute transport coefficients from fundamental lattice QCD, therefore one has to rely on models and approximations as well as limits for the high temperature–high density region (perturbative QCD) and the low temperature–low density region (hadron gas models). One of the major goals of the heavy ion program is to infer the values of the QCD transport coefficients from experimental data by detailed model-to-data comparison as explained in section 5.

Fluid dynamics is only applicable when the system is sufficiently close to local equilibrium. Very shortly after the initial collision of the two nuclei, which originally only have longitudinal momentum, this condition is not satisfied. There is some transition time required for the transverse momenta of the medium constituents to become of similar magnitude as the longitudinal momenta. The initial state and associated thermalization process has been discussed in section 4.1.1. In the later stages of the reaction, the increasingly dilute medium will drop out of equilibrium, when the expansion rate is faster than the collision rate. This stage is the subject of the next section 4.1.3. The boundary conditions for the differential equations that govern the fluid dynamic behavior, are very important for a complete understanding of the dynamics.

Conservation equations can be solved using a variety of algorithms. It is important to ensure that they also apply to relativistic hydrodynamics and event-by-event calculations with potentially spiky initial conditions without creating numerical artefacts. For this purpose a suite of tests of several scenarios, for which an analytic solution is known, was developed within the TEC-HQM collaboration. ⁸ Some widely used implementations of viscous fluid dynamics for heavy ion collisions are: MUSIC [133], vHLLE [134], VISH2 + 1 [135], and CLVisc [136].

On the theoretical side, there are activities aimed at rooting relativistic viscous hydrodynamics in kinetic theory (see e.g. [137]). Formalisms extending the applicability of fluid dynamics are also being developed. One example is anisotropic fluid dynamics that allows for larger differences between the longitudinal and transverse pressure in the system and therefore applies at earlier times than second-order viscous hydrodynamics [138, 139]. Owing to the interest in the chiral magnetic effect and polarization observables, attention has been given to developments of fluid dynamic formulations with spin and or chiral currents (see e.g. [140–142]).

Fluid dynamics is not the only option to describe the hot and dense stage of heavy ion reactions, where the quark-gluon plasma is formed. Another possibility are microscopic transport codes that include partons as their degrees of freedom. Quarks and gluons are either treated as nearly massless particles interacting according to cross sections derived from perturbative QCD [143, 144] or as dynamical quasiparticles [145]. These microscopic non-equilibrium calculations have the advantage of providing detailed information about the complete phase-space over the entire dynamical evolution of the fireball at the expense of making more detailed assumptions that may be difficult to test.

4.1.3. Hadron transport and freeze-out

When the matter created in a heavy ion collision becomes more dilute again, the quark-gluon plasma hadronizes and the hadrons subsequently decouple from the system. The moment in time when hadron abundancies are fixed and inelastic interactions no longer occur is called chemical freeze-out, while the moment when also elastic scattering ceases is called kinetic freeze-out. The standard approach to the dynamical evolution at high beam energies involves hybrid models, based on (viscous) fluid dynamics as discussed in the previous section 4.1.2 for the hot and dense stage and Boltzmann-type hadronic transport equations for the late dilute stage of the reaction (see [146] for a review).

The transition from partons to hadrons happens within the fluid dynamics calculation via a change of degrees of freedom in the equation of state. When a chosen switching criterion, usually either a certain temperature or a certain energy density, is reached, a three-dimensional hypersurface is constructed [147] within the four-dimensional space-time that serves as the basis for the particlization process, where fluid quantities are mapped into particle properties [148]. An example of such a Cooper–Frye freeze-out hypersurface is shown as the dashed curve in figure 12. The equation of state needs to be identical on both sides of the Cooper–Frye hypersurface and the transition happens according to

$$\begin{eqnarray}&&\frac{{\rm{d}}{N}_{i}}{{{\rm{d}}}^{3}p}=\frac{{d}_{i}}{{\left(2\pi \right)}^{3}}\ {\int }_{{\rm{\Sigma }}}\left[{f}_{0}(x,\vec{p})+\delta {f}_{\mathrm{shear}}(x,\vec{p})+\delta {f}_{\mathrm{bulk}}(x,\vec{p})\right]\ \frac{{p}^{\mu }\ {\rm{d}}{{\rm{\Sigma }}}_{\mu }}{{E}_{\vec{p}}}.\end{eqnarray} \tag{ 29 }$$

Here Σ is the hypersurface, f₀ the equilibrium distribution function, and δ f denotes the viscous corrections to the distribution function. The specific form of the viscous corrections is one of the major uncertainties in the current hybrid approaches [149]. Once the particles are sampled according to their distributions on the hypersurface one can feed them into hadronic transport approaches like UrQMD [150, 151] or SMASH [152].

One generally assumes that the majority of particles escapes outward of this hypersurface but, in principle, there can be also inward directed contributions. These are particles that reenter the hydrodynamic evolution and are therefore very hard to deal within practice [154]. While all other parts of the evolution satisfy exact conservation laws event by event, in the sampling process they are only fulfilled for an ensemble of events. Therefore, very often many such samples and hadronic evolutions are calculated for each hydrodynamic event (oversampling method). To study correlation and fluctuation observables one might need a more careful consideration of conservation laws at the Cooper–Frye transition [155]. For calculations of dilepton emission one needs to include the finite-width spectral functions of the resonances in the sampling process. For this purpose, it is necessary to match the degrees of freedom and their properties in the hadronic transport approach that is chosen for the final rescattering and decays.

The advantage of the microscopic non-equilibrium treatment of the late hadronic stage is that the chemical and kinetic freeze-out are happening automatically, whenever the corresponding interactions cease for each particle species individually. An alternative option is to keep running hydrodynamics and introducing partial chemical equilibrium. In this approach one encounters a kinetic freeze-out hypersurface, an example of which is shown in figure 13. The explicit particle degrees of freedom have the added advantage that the outcome is very similar to what is measured in experiment. When using unified output formats, like OSCAR or the HepMC library ⁹ , the same analysis can then be applied to the simulation results and to the experimental data, which allows for more accurate comparisons. Recently, the RIVET library ¹⁰ is also getting extended to include heavy ion observables, which could be an interesting option in the future.

Zoom In Zoom Out Reset image size

The same information on hadronic degrees of freedom and their interactions that one needs for final-stage hadron transport calculations in high-energy heavy ion collisions is employed for the full evolution in low energy heavy ion reactions. Hadronic transport approaches are based on the relativistic Boltzmann equation where, appropriate for low net baryon density, the mean-field term is omitted and only collisions are taken into account. The input for cross-sections and processes is taken from the Particle Data Book, experimental data for elementary reactions, or estimated using effective models. For low energy binary scatterings, the reactions are modeled via resonance excitation and decay, while at higher energies string excitation and fragmentation is used to model inelastic processes. Such microscopic approaches have the advantage that one has access to the complete phase-space information of all particles at all times, and they often serve as event generators for simulations of experimental results.

Hard probes of hot and dense QCD matter are defined as those that are, in important parts, perturbatively calculable. This applies to strongly interacting (QCD) probes at high virtuality and to electromagnetic probes. In most instances, hard probes rely on input that is nonperturbative and must be derived from other experimental data. Parton distribution functions (PDFs) f(x) and fragmentation functions D(z) are salient examples.

Generally, hard probes in p+p or A+A collisions rely on the principle of factorization, illustrated in figure 14: A scattering process is factorized into a nonperturbative initial-state matrix element, a perturbatively calculable hard scattering S-matrix element, and a nonperturbative final-state matrix element:

$$\begin{eqnarray}&&{\rm{d}}\sigma ={f}_{{\rm{A}}}({x}_{1},{Q}^{2})\otimes {f}_{B}({x}_{2},{Q}^{2})\otimes {\rm{d}}\hat{\sigma }({Q}^{2})\otimes {D}_{h/c}(z,{Q}^{2}).\end{eqnarray} \tag{ 30 }$$

The process at the core of this scheme, the hard scattering cross section ${\rm{d}}\hat{\sigma }({Q}^{2})$, is universal and only depends on the hard scale Q², but the initial- and final-state matrix elements that give rise to the parton distribution functions f(x, Q²) and the fragmentation function D(z, Q²) can be modified by nuclear effects—cold nuclear matter in the case of the initial state and hot QCD matter in the case of the final state.

Zoom In Zoom Out Reset image size

The deceptively simple sketch in figure 14 hides many of the complexities that arise as a result of the nonperturbative interactions occurring in the initial and final states. This is illustrated in figure 15. Their sensitivity to the nuclear medium is what makes them of core interest to relativistic heavy ion physics. The challenge is to encapsulate these modifications in physics-motivated descriptions that can be used to probe the properties of the nuclear medium. Examples are parton energy loss, diffusion coefficients, saturation scales, and so on.

Zoom In Zoom Out Reset image size

We have already covered models of the initial state in section 4.1.1. We now discuss models for the final-state modifications. The output of the hard scattering process are two or more energetic partons, which then can interact with the medium created in the collision, often referred to as the underlying event. These interactions can result in energy or momentum loss [67], deflection [156], emission of radiation [157], or transmutation of the parton into another particle [158]. A sequence of such interactions leads to modifications of the jet shower evolution, which manifests itself in various observables (see section 5.9) and can be modeled using Monte-Carlo methods (see e.g. [159, 160]).

The mechanisms causing an energetic parton to lose energy during its propagation through matter can be divided into two categories: elastic and inelastic. Elastic energy loss occurs in two-body scattering, where energy is transferred to the scattering partner. Inelastic energy loss occurs predominantly by radiation of a gluon following the interaction with a scattering center in the matter. At high energies, inelastic energy loss dominates; at low energies elastic energy loss dominates. The threshold at which inelastic energy loss begins to dominate increases with the mass of the energetic parton, because radiation from heavy particles is suppressed by the dead-cone effect [161].

Elastic energy loss is one aspect of momentum change in the medium caused by scattering. For heavy quarks with mass M ≫ T multiple scatterings are required to appreciably change the momentum. The process can then be considered as a diffusion process and described by a Langevin equation with a fluctuation constant κ and drag coefficient η. ¹¹ In a thermal medium κ and η are related to the spatial diffusion constant D and to each other by the Einstein relation

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{2{T}^{2}}{D},\quad \eta =\displaystyle \frac{T}{{MD}}=\displaystyle \frac{\kappa }{2{MT}}.\end{eqnarray} \tag{ 31 }$$

The diffusion constant can be calculated in thermal perturbation theory [68], on the lattice [162], or by holographic techniques [163]. The Langevin framework can be easily simulated within a fluid dynamical description of the quark-gluon plasma. Rare scatterings with larger momentum transfer are not described by the Langevin approach. These can be treated in the framework of linearized Boltzmann transport in a thermalized, hydrodynamically flowing medium [164]. Both descriptions can be merged into a unified framework that effectively describes all aspects of heavy quark transport [165].

Inelastic energy loss, i.e. energy loss by gluon radiation, is the dominant mode of energy loss for light quarks and gluons. Medium-induced radiation requires off-mass shell scattering of the energetic parton in the medium, followed by gluon emission. In a dense medium, multiple scattering events contribute coherently to a single radiation event. The resulting suppression of the radiation yield is known as Landau–Pomeranchuk–Migdal (LPM) effect. The scattering power of the medium (momentum transfer squared per unit length) is expressed in the transport coefficient $\hat{q}$ defined in (3). Several different formalisms have been developed to describe the radiative energy loss of quarks or gluons. These correspond to different approximations or truncations of the multiple scattering process. For an overview and comparison of existing approaches, see [166]. General overviews of energy loss and jet quenching can be found in [29].

A fraction of the energy lost by an energetic parton traversing the quark-gluon plasma is deposited into the thermal medium. This applies to all energy lost by elastic collisions and to the soft component of the radiated energy which thermalizes. The response of the medium to the injection of energy and momentum along the trajectory of the energetic parton is treated by adding a source term to the hydrodynamic equations. The linear response of the medium can be understood in terms of sound propagation and results in the possible formation of a Mach cone [167]. A more microscopic approach treats the individual interactions within the evolution of a full jet as sources in a numerical hydrodynamics code [168]. For the effect of the medium response on jet observables, see [169].

Electromagnetic (EM) probes, i.e. photons and dileptons (often called real and virtual photons), are special in that they do not suffer from final state interactions. (They may, however, suffer from large backgrounds.) Their production can be related to the photon spectral function $\rho (\omega ,\vec{p})$ in the medium [58]. In the vacuum, ρ is a function of the invariant mass M² = ω² − p² only; in a thermal medium, ρ(ω, M) is a function of, both, M and the frequency measured in the fluid rest frame ω = −u_μ p^μ . The rate of photon emission can be expressed as [170]

$$\begin{eqnarray}&&\omega \frac{{\rm{d}}R}{{{\rm{d}}}^{3}p}=\frac{{n}_{{\rm{B}}}(\omega )\rho (\omega ,0)}{{\left(2\pi \right)}^{3}}\end{eqnarray} \tag{ 32 }$$

with ${n}_{{\rm{B}}}(\omega )={\left({{\rm{e}}}^{\omega /T}-1\right)}^{-1}$. The dilepton emission rate is given by

$$\begin{eqnarray}&&\frac{{\rm{d}}R}{{{\rm{d}}}^{4}p}=\alpha \frac{{n}_{{\rm{B}}}(\omega )\rho (\omega ,M)}{12{\pi }^{4}{M}^{2}}\left(1+\frac{2{m}^{2}}{{M}^{2}}\right)\sqrt{1-\frac{4{m}^{2}}{{M}^{2}}}\theta ({M}^{2}-4{m}^{2}),\end{eqnarray} \tag{ 33 }$$

where m denotes the lepton mass. Apart from the fact that the measured yields involve an integration over the space-time evolution of the emitting matter, EM probes thus provide direct information about the photon spectral function. A broad overview of electromagnetic probes of relativistic heavy ion collisions can be found in [171]; concise up-to-date reviews are contained in [172, 173]. For a comprehensive recent review of direct photon emission see [174].

The photon spectral function is predicted to change from one that is dominated by vector mesons at T < T_c to a broad continuum corresponding to free quark-antiquark pair excitations at T > T_c (see figure 16). It can be estimated using lattice QCD [176, 177], thermal perturbation theory [72], or using QCD sum rules [178, 179]. In the sum rule approach, the transition of the photon spectral function is closely related to the phenomenon of chiral symmetry restoration. At low frequencies and long wavelengths, the spectral function is related to the electrical conductivity σ of the QGP by

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{1}{2}\mathop{\mathrm{lim}}\limits_{\omega \to 0}\displaystyle \frac{1}{\omega }{\left.\rho (\omega ,p)\right|}_{{\omega }^{2}={p}^{2}}=\displaystyle \frac{1}{3}\mathop{\mathrm{lim}}\limits_{\omega \to 0}\displaystyle \frac{1}{\omega }{\left.\rho (\omega ,p)\right|}_{p=0}.\end{eqnarray} \tag{ 34 }$$

Zoom In Zoom Out Reset image size

It is a challenge for the future to combine the soft bulk evolution with the hard probes consistently. Many calculations use state-of-the-art descriptions of the bulk matter evolution for the evaluation of hard probes, but including the medium response and developing a fully coupled picture over the different kinematic regions remain challenges for the future.

Different kinematic regions of the particles produced in heavy ion collisions are governed by different regimes of quantum chromodynamics. The soft region below ∼3 GeV is typically described by the bulk dynamic evolution (see section 4.1), while the hard region above ∼6 GeV is covered by perturbative QCD (see section 4.2). Experimentally, it is only possible to detect final-state hadrons, while we would like to learn about the hot and dense quark-gluon plasma stage of the reaction. Therefore, the hadronization process plays an important role in our interpretation of the data from heavy ion reactions.

Figure 17 shows an example of the nuclear modification factor R_AA as a function of transverse momentum for several particle species from the ALICE collaboration in PbPb collisions at $\sqrt{{s}_{\mathrm{NN}}}=2.76$ TeV. This picture serves nicely as a basis to introduce the hadronization models that are applied in different regions.

Zoom In Zoom Out Reset image size

Below ∼2 GeV the particle production is described by the viscous hydrodynamic evolution, and hadronization happens via the change of degrees of freedom in the equation of state. On the Cooper–Frye surface the hadronic fluid is converted to individual particles and those are then propagated by a hadronic cascade. This is one of the big advantages of hydrodynamics, that the transition from partonic to hadronic degrees of freedom is under nice theoretical control as long as local equilibrium can be assumed. In this region the nuclear modification factor shows the typical mass dependence expected from collective behavior. For a more detailed discussion see section 5.2.

Above a transverse momentum of p_T ∼ 8 GeV one observes a rather universal behavior, this is the region where hadronization is described by fragmentation of color flux tubes. The strings break apart when the leading quark-antiquark (or quark-diquark) pairs are pulled apart and the color field becomes strong enough to create particle-antiparticle pairs out of the vacuum. Subsequently, the newly produced partons connect with each other to form hadrons with well-defined quantum numbers.

The intermediate region p_T = 2−8 GeV is governed by the recombination and coalescence models for hadronization and the interaction of hard fragmentation and soft hadrons. Recombination means that the partons from the quark-gluon plasma are clustered into hadrons when they are close in phase space (see section 5.8 for more details). In this region, unusually large baryon-to-meson ratios can be observed in the nuclear modification factor.

Understanding hadronization on a microscopic level poses a challenge for the future. It is important to advance the understanding on how hadrons emerge from the collective partonic state that is formed in ultra-relativistic heavy ion collisions. This involves understanding QCD on multiple different scales connecting the different regions outlined above in a smooth fashion.

The first question to be answered when studying heavy ion reactions is which temperature and density is reached in the collision. While full dynamical models as described in section 4 provide (model-dependent) answers to this question, the possibility of defining a temperature for the system has more fundamental implications. To be able to assign a single temperature and net baryon chemical potential to the system created in such a highly dynamical process hints at the fact that an equilibrated plasma has formed. Another interpretation involves the idea that hadronization follows statistical behavior and, therefore, the numbers of produced hadrons follow thermal expectations. The second interpretation is supported by the finding that even in particle production from e⁺ e⁻ collisions thermal particle yields are measured [181].

The basic idea to connect final yields at midrapidity or 4π multiplicities to a temperature relies on the formation of a fireball in global equilibrium, that emits all particles at the same instant in time. The thermal model actually has no notion of time evolution or spatial variations, therefore it is limited to a small set of observables, namely the particle abundances. The grand canonical partition function includes all mesons and baryons of the non-interacting hadron resonance gas. The fireball is treated as a grand canonical ensemble with Bose and Fermi statistics depending on the particle spin. The distribution function for particle species i is given by

$$\begin{eqnarray}&&\frac{{\rm{d}}{N}_{i}}{{\rm{d}}{p}^{3}{\rm{d}}{x}^{3}}=\frac{{g}_{i}}{{\left(2\pi \right)}^{3}}\,\frac{1}{{\exp }^{(E-{\mu }_{i})/T}\pm {\alpha }_{i}},\end{eqnarray} \tag{ 35 }$$

where g_i is the degeneracy factor, T is the temperature, μ_i = μ_B B_i + μ_S S_i the particle-specific chemical potential and α_i = ±1 for fermions and bosons, respectively. Integrating over the entire phase space and assuming Boltzmann statistics—a good approximation for all hadrons except pions when ∣μ_i ∣/T is small—one obtains the following well-known expression for thermal particle production

$$\begin{eqnarray}&&{N}_{i}=\frac{{g}_{i}{{VT}}^{3}}{2{\pi }^{2}}\,{{\rm{e}}}^{{\mu }_{i}/T}\,\frac{{m}_{i}^{2}}{{T}^{2}}\,{K}_{2}\left(\frac{{m}_{i}}{T}\right),\end{eqnarray} \tag{ 36 }$$

where N_i is the number of particles of species i produced in the fireball of volume V at temperature T and chemical potential μ_i . m_i denotes the mass of the particle and K₂ is the modified Bessel function of the second kind. The main quantity determining the particle yields is the mass of the particle species. One expects higher mass particles to be produced with smaller probabilities than lighter particles. To remove sensitivity to the volume, the usual strategy is to look at particle yield ratios that are assumed to be emitted from the same volume. By fitting ratios of several different particle species one can then obtain temperatures and chemical potentials.

It is impressive how well a vast amount of data on particle production in heavy ion collisions can be understood with such simple assumptions. One recent example is shown in figure 18 (left) where particle ratios at midrapidity from Pb+Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}=2.76$ TeV are fit within the statistical model. The bars denote the model fit and the symbols show experimental data from the ALICE collaboration. Even the light nuclei follow the trend of thermal production at a unique global temperature and chemical potential. A widely adopted interpretation is that these thermodynamic properties reflect the conditions prevailing at the chemical freeze-out, the moment at which the inelastic reactions cease and the abundances are frozen. (Elastic rescattering is still possible until the kinetic freeze-out which will be discussed in more detail in the next section 5.2.) It is worth pointing out, however, that this interpretation is not logically consistent, because the model also describes the yields of particles that cannot be formed under the chemical freeze-out conditions, such as light nuclei. A physically credible interpretation of the observed yields of such states can be made either on the basis of the dynamical coalescence model [183] or by invoking general properties of strongly coupled, highly excited quantum systems [184].

Zoom In Zoom Out Reset image size

The right-hand part of figure 18 shows the result of thermal fits to particle yields and multiplicities over a large range of beam energies. As expected, the temperature rises as a function of beam energy, while the net baryon chemical potential decreases. This can be understood by the stopping dynamics of the two nuclei passing through each other. At lower energy the nuclear passing time is longer and more baryons are stopped at midrapidity, while at high beam energies the nuclei fly through each other so fast that the baryon number remains localized at forward rapidities (transparency). The parametrizations shown in figure 18 (right) are [182]:

$$\begin{eqnarray}&&{T}_{\mathrm{CF}}=\displaystyle \frac{{T}_{\mathrm{CF}}^{\mathrm{lim}}}{1+\exp (2.60-(\mathrm{ln}\sqrt{{s}_{\mathrm{NN}}})/0.45)},\qquad {\mu }_{B}=\displaystyle \frac{a}{1+0.288\sqrt{{s}_{\mathrm{NN}}}},\end{eqnarray} \tag{ 37 }$$

with $\sqrt{{s}_{\mathrm{NN}}}$ in GeV, ${T}_{\mathrm{CF}}^{\mathrm{lim}}=158.4\,\mathrm{MeV}$ and a = 1307.5 MeV.

When plotting the values for temperatures and net baryon chemical potential obtained in thermal model fits on the phase diagram one obtains a curve as specified in figure 19. Since those values correspond to a hadron resonance gas, the quark-gluon plasma formation can only happen beyond this line. These values define a lower bound of temperatures and densities that are reached in the corresponding heavy ion collisions, since the chemical freeze-out happens at the end of the evolution. The shape of the line is consistent with the hypothesis that the chemical freeze-out is associated with a constant energy per baryon in the system [185].

Zoom In Zoom Out Reset image size

At lower beam energies and for smaller system sizes, the assumptions of a grand canonical ensemble may no longer be valid, especially for particles carrying a conserved quantum number such as strangeness. Since s and $\overline{s}$ quarks must be produced in pairs, the conservation of net strangeness in the volume must be taken into consideration. Using a thermal model that includes these ideas, one can even obtain decent fits to particle ratios from low energy heavy ion reactions as shown in the left panel of figure 20 by the HADES collaboration. While it is remarkable how well the extracted temperatures and net baryon chemical potentials fit into the world data (see right panel of figure 20), some discrepancies are also appearing. For example, the high yield of Ξ baryons needs a different explanation.

Zoom In Zoom Out Reset image size

At high beam energies, there are some indications that strange particles have a higher chemical freeze-out temperature than hadrons consisting of light quarks only [187]. This could be understood by the smaller cross-sections of strange hadrons with the surrounding matter. Those studies go beyond fitting yields by investigating the relevant susceptibilities through particle number fluctuations of various orders. There are a lot of caveats when comparing these observables directly to grand canonical lattice QCD calculations but, if properly taken into account, allow to extract valuable information on the thermodynamic properties of the system.

If strangeness is found to be in thermal equilibrium with the rest of the produced particle species, it is an indication of quark-gluon plasma formation since the strangeness production is enhanced through partonic production channels compared to reactions at the hadron level [188]. There have been even claims that the peak in the kaon-to-pion ratio as a function of beam energy around E_lab = 40A GeV is a signature for the phase transition to the quark-gluon plasma [189]. Also the 'step' in the mean transverse momentum of kaons as a function of beam energy has been attributed to the latent heat of the phase transition. However, in dynamical model calculations both phenomena can be attributed to more conventional effects, namely the transition from a baryon-dominated to a meson-dominated system, as well as to the softening of the equation of state due to the excitation of more active degrees of freedom associated with hadron resonances [190, 191].

Besides canonical suppression, there are other more advanced variants of the thermal model including, for example, van der Waals interactions [192] or excluded volume effects that can vary for different species [193]. The limitations of the thermal model are that it is a static picture. It is not possible to address observables of more dynamical origin (e.g. anisotropic flow, HBT radii,..) within this approach. Also, resonances cannot be described since the rescattering of daughter particles restricts the ability of experimental collaborations to reconstruct all resonances produced at chemical freeze-out. At LHC energies, baryon annihilation during the late stage of the fireball expansion plays a significant role and affects the proton yields, which also leads to challenges in a fully thermal description [194].

In addition to looking globally at particle production at midrapidity or integrated yields, one can investigate the variation of particle yields as a function of rapidity. For example, the rapidity distribution of net protons gives valuable insight into the baryon stopping dynamics [195]. In the context of the present section, it is important to realize that by studying more forward rapidities higher net baryon densities are accessible also at colliders like RHIC and LHC [196]. Of course, the high rapidity region is typically not equipped with as many detector elements capable of particle identification, but the high energies and high multiplicities might still provide an advantage over low beam energies, especially in fixed target experiments as they have been conducted in the STAR detector at RHIC and at the CERN-SPS.

The basic dynamic observables related to transverse flow are transverse momentum distributions of (identified) particles. A more integral, commonly reported measure is provided by the mean transverse momentum 〈p_T 〉. These observables carry information about the radial expansion of the system in addition to thermal information. While the yields of particles are fixed at chemical freeze-out when inelastic reactions cease, the slopes of the transverse momentum spectra contain information about the kinetic freeze-out of the different particle species. Kinetic freeze-out is the moment when the particles no longer suffer elastic collisions and the momenta become fixed. Because elastic cross sections differ from particle to particle the effective temperatures extracted from spectra of different particles do not necessarily agree. Note that the decay of a resonance into the same (final) state from which it was formed, such as the reaction π N → Δ → N π, is understood as a pseudoelastic process, since it does not alter the chemical composition.

Employing the ideas introduced in the previous section 5.1 one can calculate the thermal expectation for particle spectra differentially along transverse momentum under the assumption of a Boltzmann distribution and obtains that

$$\begin{eqnarray}&&\frac{1}{{p}_{T}}\frac{{\rm{d}}N}{{\rm{d}}{p}_{T}}\propto {{\rm{e}}}^{-{p}_{T}/{T}_{\mathrm{eff}}}.\end{eqnarray} \tag{ 38 }$$

where p_T is the transverse momentum, and T_eff is the slope of the distribution in a semi-logarithmic plot. In figure 21 the transverse momentum spectra for pions, kaons and protons in heavy ion collisions at LHC energies are shown. The shoulder at low p_T indicates a deviation from the purely thermal expectation that becomes more pronounced with increasing hadron mass, in particular, for protons. To obtain exponential distributions one can apply a variable transformation and fit transverse mass spectra ${m}_{T}^{-1}({\rm{d}}N/{\rm{d}}{m}_{T})$ instead of transverse momentum spectra. It is interesting that most of the information on particle production can be summarized by the integrated yield at midrapidity and the mean transverse mass. If a theoretical calculation reproduces both, the full spectra are typically described rather well.

Zoom In Zoom Out Reset image size

By measuring transverse momentum distributions in different centrality windows and rapidity bins, one can obtain the yields as a function of centrality and rapidity by integrating over these spectra. Since the acceptance of experimental detectors is often limited in the very low momentum region, one needs adequate fit functions. The particles with very low momenta are bent strongly in the magnetic fields and therefore do not reach the tracking detectors. To estimate systematic uncertainties, different functional forms are employed, e.g. pure exponential functions in p_T or m_T , functions including quantum statistics or Tsallis-Levy distributions.

To extract kinetic freeze-out temperatures and average flow velocities from transverse momentum spectra of several particle species the blast-wave model is employed. The assumption here is that all particles are emitted from a boosted fireball that follows a common radial expansion profile. The different behavior of the different species is then mainly attributed to the different masses that result in different mean transverse momenta even for the same underlying velocity profile. To obtain integrated yields and mean transverse momenta, the spectra shown in figure 21 have been fitted individually for each species with a blast-wave function [199]:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{p}_{T}}\displaystyle \frac{{\rm{d}}N}{{\rm{d}}{p}_{T}}\propto {\int }_{0}^{R}r{\rm{d}}r\,{m}_{{\rm{T}}}\,{I}_{0}\left(\displaystyle \frac{{p}_{{\rm{T}}}\sinh \rho }{{T}_{{kin}}}\right){K}_{1}\left(\displaystyle \frac{{m}_{{\rm{T}}}\cosh \rho }{{T}_{{kin}}}\right),\end{eqnarray} \tag{ 39 }$$

where the velocity profile ρ is described by

$$\begin{eqnarray}&&\rho ={\tanh }^{-1}{\beta }_{{\rm{T}}}(r)={\tanh }^{-1}\left({\beta }_{{\rm{s}}}{\left(r/R\right)}^{n}\right).\end{eqnarray} \tag{ 40 }$$

Here, ${m}_{T}=\sqrt{{p}_{T}^{2}+{m}^{2}}$ is the transverse mass, I₀ and K₁ are modified Bessel functions, r is the radial distance in the transverse plane, R is the transverse radius of the fireball, β_T(r) is the transverse expansion velocity and β_s is the transverse expansion velocity at the surface.

Combining the spectra of different particle species in a joint fit makes it possible to extract kinetic freeze-out temperatures and average flow velocities for different beam energies and centralities as shown in figure 22. The results are rather similar at high RHIC and LHC energies. As a function of centrality, one can see a trend that more central collisions result in lower temperatures and higher transverse velocities, while more peripheral collisions end up at higher temperatures and with less explosive radial expansion. The apparently counter-intuitive ordering of the deduced temperatures indicates that the fireball in central collisions produces more particles and is 'cooking' longer, such that the particles decouple at lower temperatures when the system had more time to expand. One caveat is that the results are rather sensitive to the fit ranges, which have to be carefully chosen (see the black bars in figure 21).

Zoom In Zoom Out Reset image size

The world data shown in the right panel of figure 22 reveal the expected behavior of rising radial expansion velocities with increasing collision energy. The extracted temperatures show a more interesting pattern. At low beam energies the chemical freeze-out overlaps with the kinetic freeze-out, while above $\sqrt{{s}_{\mathrm{NN}}}\sim 8\,\mathrm{GeV}$ the kinetic freeze-out temperatures are lower than the chemical ones and stay constant or even decrease at higher beam energies. One way to understand this invokes the same argument as above for central collisions, that the fireball has a longer lifetime at high beam energies and, therefore, the particles decouple later in the evolution.

During the final expansion of the hadronic fireball both, the momenta of the hadrons and the chemical composition of the fireball, change. The latter can occur in two ways: by decay of unstable, excited hadrons, and by inelastic collisions among hadrons. The latter are generally thought to diminish rapidly as the fireball cools, since most low-energy interactions among ground-state hadrons are elastic. A notable exception to this rule are baryon-antibaryon annihilation reactions, which have large cross-sections at low energies. Hadron momenta, of course, can be changed by elastic scattering. This applies especially to the most abundant hadrons at high beam energies, pions and nucleons, which scatter with a very large cross section due to the Δ-resonance.

One important consideration when comparing model calculations to experimental data is to 'match apples to apples'. For transverse mass spectra this seems straightforward, but there can be complications like feed-down corrections that have to be taken into account. In many models the Λ hyperon is regarded as stable particle, because it is stable with respect to the strong interaction. In practice, some of the produced Λ hyperons undergo a weak decay before they reach the detector, which is at a macroscopic distance from the collision vertex. Therefore, the reported number of protons differs by about 40 % depending on whether the feed-down from Λ hyperons is being considered or not. Similarly, the Σ⁰ decays to almost 100% into Λ and needs to be added in model calculations before comparing to experimental data.

Transverse momentum spectra and particle yields are also the basic observables that a hydrodynamic evolution needs to describe. In figure 23 calculations within the VISHNU hybrid approach based on the VISH2 + 1 dimensional viscous hydrodynamic code and hadronic rescattering by UrQMD are compared to pion and proton spectra in Au+Au collisions at the highest RHIC energy. One can see that different levels of initial state fluctuations described by Monte-Carlo Glauber or a saturation-based MC-KLN approach do not affect the spectra much. Also, varying the constant shear viscosity-over-entropy density ratio η/s employed in the hydrodynamic evolution from 0 to 0.24 does not have a large effect. (See below for the dependence on bulk viscosity.) Of course, smaller variations might become visible when one compares results on a linear scale. In general, transverse momentum spectra can be used to gauge the basic hydrodynamic evolution parameters. The anisotropic flow discussed in the next section 5.4 has a higher sensitivity to the initial state fluctuations and to the properties of the quark-gluon plasma and its transport coefficients.

Zoom In Zoom Out Reset image size

Figure 24 (left) shows the mean transverse momentum of pions, kaons and protons in heavy ion collisions at LHC energies, measured by the ALICE collaboration, in comparison with a hybrid calculation relying on IP-Glasma fluctuating initial conditions, the 3 + 1 dimensional viscous hydrodynamic implementation MUSIC and hadronic rescattering by UrQMD. As a function of centrality, the mean transverse momenta are found to be very sensitive to the inclusion of bulk viscosity in the evolution. A temperature-dependent bulk viscosity with a peak around the transition temperature was employed in this study. While there is the caveat that this calculation still used an equation of state with a transition temperature around T_c = 190 MeV, the qualitative differences are expected to remain also with an equation of state fitted to state-of-the-art lattice calculations (see e.g. the recent study of deuterons and their sensitivity to bulk viscosity in [203]). The right panel of figure 24 shows the effect of hadronic rescattering. While pions and kaons are only mildly affected, the transverse momentum of protons increases by about 30% during the late stage due to rescattering. This can be attributed to the large cross sections of nucleons with fast-moving pions; the phenomenon is sometimes referred to as 'pion wind'.

Zoom In Zoom Out Reset image size

Measurements of correlation functions between particle pairs provide information about the geometry of the emission region in heavy-ion collisions. The inferred size of the so-called 'region of homogeneity' reflects the kinetic freeze-out, since these measurements are performed on final state hadrons. This technique is originally applied in the size estimation of stars by Hanbury-Brown and Twiss (HBT). The shape of the correlation function is typically assumed to be Gaussian with the following parametrization:

$$\begin{eqnarray}&&\begin{array}{l}C\left(\vec{q}\right)=\left(1-\lambda \right)+\lambda \,{K}_{\mathrm{Coul}}({q}_{\mathrm{inv}})\\ \,\times \,\exp \left(-{q}_{{\rm{o}}}^{2}{R}_{{\rm{o}}}^{2}-{q}_{{\rm{s}}}^{2}{R}_{{\rm{s}}}^{2}-{q}_{{\rm{l}}}^{2}{R}_{{\rm{l}}}^{2}-2{q}_{{\rm{o}}}{q}_{{\rm{s}}}{R}_{\mathrm{os}}^{2}-2{q}_{{\rm{o}}}{q}_{{\rm{l}}}{R}_{\mathrm{ol}}^{2}\right),\end{array}\end{eqnarray} \tag{ 41 }$$

where a value of the chaoticity parameter λ < 1 accounts for dilution of the HBT signal due to particle correlations from resonance decays, and the different radii (R_o, R_s and R_l) refer to different directions. The out-direction (R_o) points along the direction of the average momentum of the particle pair, the long-direction (R_l) along the beam direction and the side-direction R_s denotes the third orthogonal direction. q_inv denotes the relative momentum of the two particles measured in the rest frame of the particle pair; modifications of the relative wave function of the pair is most conveniently expressed in terms of this quantity. The contribution proportional to λ accounts for quantum correlations. When the data are insufficient for a three-dimensional analysis, the correlation function is often assumed to be spherically symmetric,

$$\begin{eqnarray}&&C\left(\vec{q}\right)=\left(1-\lambda \right)+\lambda \,{K}_{\mathrm{Coul}}({q}_{\mathrm{inv}})\exp \left(-| \vec{q}{| }^{2}{R}_{\mathrm{inv}}^{2}\right)\end{eqnarray} \tag{ 42 }$$

with a single Gaussian source radius R_inv. A review on the topic can be found in [31]. As a function of mean transverse mass of the particle pairs, there is a global decreasing trend [204, 205] This can be easily understood due to the fact that faster particles are emitted earlier in the evolution and therefore from smaller source sizes than slower ones.

Identical particle correlations have also been measured for neutral and charged kaons and for protons [206]. The results are generally consistent with those obtained for pions, but slightly smaller source sizes have been measured for kaons perhaps pointing to their earlier kinetic freeze-out. Information about source sizes can also be gleaned from unlike charged particle correlations, where the information resides in the factor K_Coul(q) that accounts for final-state Coulomb and other interactions between the emitted particles (see e.g. [207]). In the future, looking at HBT correlations for dileptons and photons can provide more information about the geometry of the system during the evolution since electromagnetic probes leave the fireball from all stages without disturbance from rescattering.

The excitation function of the radius parameters as well as the λ parameter and the R_o/R_s ratio is summarized in the left panel of figure 25. The ratio of the two transverse radii is predicted to be sensitive to the lifetime of the system. Interestingly, this ratio stays rather constant over a large range of beam energies. The fraction of correlated particles λ is higher at lower beam energies and saturates at higher energies to a lower constant value, which points to an increasing fraction of particles originating from longer-lived resonances in the sample. R_out and R_side stay more or less constant as a function of beam energy while R_long is rising significantly indicating a longer lifetime of the overall system evolution.

Zoom In Zoom Out Reset image size

The right panel of figure 25 shows the dependence of the HBT radii on the charged particle rapidity density [205]. If the particle freeze-out always occurs at the same density, the radii should be linear functions of ${\left({\rm{d}}{N}_{\mathrm{ch}}/{\rm{d}}\eta \right)}^{1/3}$. The data, especially those for R_l, are roughly consistent with this hypothesis. By combining the transverse HBT radii one can infer a measure of the fireball volume at kinetic freeze-out as shown in figure 26. These volumes can be compared to the one required by the thermodynamic description in section 5.1 and the qualitative behavior as a function of beam energy is consistent. In general, the analyzes of HBT radii and their comparison to hydrodynamic calculations allow one to connect coordinate and momentum space descriptions that are otherwise hard to achieve, see e.g. [153]. The complexity of this connection was analyzed in [208] where it was shown that multiple aspects of the dynamical evolution influence the size of the HBT radii and that it is non-trivial to get the bulk evolution in agreement with the correlation measurements.

Zoom In Zoom Out Reset image size

A large value of the R_o/R_s ratio was proposed as one of the signatures of a first-order phase transition as a long duration of particle emission increases the outward-directed correlation length R_o. When the system undergoes a phase transition and the pressure is reduced significantly, the lifetime is expected to be significantly extended [210, 211]. In such a scenario one expects an increase in the R_o/R_s ratio as a function of beam energy. Figure 27 shows a calculation within a hybrid approach based on UrQMD initial conditions, (3+1)-dimensional ideal hydrodynamic evolution, and final hadronic transport within UrQMD. Different switching energy density criteria as well as the equation of state influence the lifetime. The noticeable result is that the only curve that qualitatively follows the slight peak as a function of beam energy that the experimental data from the NA49 collaboration suggests is the one with a first-order phase transition. Detailed correlation observables have the potential to hint at a first-order phase transition between hadron gas and quark-gluon plasma although, when the dynamical evolution is modeled more realistically, the effect is much smaller than originally estimated.

Zoom In Zoom Out Reset image size

The anisotropic flow coefficients as defined in section 2.4 contain interesting information about the properties of the quark-gluon plasma as well as the initial conditions. The angular modulations of collective flow in the plane transverse to the beam axis are analyzed in a Fourier decomposition. Non-vanishing values of v_n have been measured to high precision at RHIC and LHC as a function of transverse momentum, particle species and centrality. When hydrodynamic calculations are tuned to describe the yield and transverse momentum spectra for a system, the anisotropic flow can be very well predicted [120, 121]. This central finding that lies at the foundation of our 'standard model' of heavy ion collisions is based on a ≤10% effect (for the average values) on the background radial expansion discussed in the previous section 5.2. The left panel of figure 28 shows one such example of a (3+1)-dimensional viscous hydrodynamic calculation for Au+Au collisions at $\sqrt{{s}_{\mathrm{NN}}}=200\,\mathrm{GeV}$ that fits the charged particle anisotropic flow as a function of transverse momentum.

Zoom In Zoom Out Reset image size

The agreement of hydrodynamic calculations with anisotropic flow measurements forms the basis for the paradigm that the quark-gluon plasma exhibits an extremely low ratio of shear viscosity to entropy density and is therefore one of the most 'perfect' fluids in nature. Interestingly, this feature is shared with ultra-cold matter that can be probed in the laboratory, where atoms are confined in a trap in an almond-shaped configuration. Once released the particle stream along the pressure lines and the coordinate space anisotropy is transformed into a momentum space anisotropy [212]. The complex response function of the nonlinear hydrodynamic evolution can be dissected into eccentricities of different order and combinations of eccentricities that lead to the final-state flow coefficients (see e.g. [213]).

Higher order flow coefficients probe smaller-scale structures and are more sensitive to the properties of the plasma. Figure 28 (right) demonstrates the effect of increasing the shear viscosity-over-entropy density ratio from 0 to 0.08 and 0.16. Higher viscosities lead to smaller flow coefficients since the initial state structures are diluted more quickly. Without initial state fluctuations the odd flow coefficients would be zero by symmetry. But since 2010 it has been recognized that triangular flow (and higher coefficients) are actually non-zero resulting from small scale structures in the initial state [106]. Figure 29 shows the first measurement of triangular flow as a function of the number of participants in Au+Au collisions at $\sqrt{{s}_{\mathrm{NN}}}=200\,\mathrm{GeV}$. While all calculations agree for elliptic flow in central collisions, the v₃ measurement is highly sensitive to the initial state structures (see [215] for a review). On the one hand, this poses the challenge that collective flow is not only sensitive to the medium properties but also to the details of the initial state. On the other hand, this also offers the opportunity to learn something about the initial state created by two nuclei colliding close to the speed of light.

Zoom In Zoom Out Reset image size

By systematic improvements of theory and experiment it has been confirmed that all observed collective flow data are consistent with a quark-gluon plasma that exhibits a very low specific shear viscosity that is close to the lower bound predicted by AdS/CFT for a strongly coupled gauge plasma of $\eta /s={\left(4\pi \right)}^{-1}=0.08$. The averaged flow coefficients as a function of centrality are the most sensitive observables for this purpose. The full event-by-event distributions of the flow coefficients are mainly sensitive to the initial state and its fluctuations, while the viscosity mainly influences the average. Nowadays, this complex inverse problem involving many parameters and many observables is attacked with Bayesian analysis techniques as discussed in section 5.6.

The flow data for identified particles are sensitive probes of the hadronization mechanism and the dynamics of hadronic rescattering in the late stage that can be described by transitioning to a transport treatment for the fireball evolution (see [146] for a review). Additional information on the details of the initial conditions and the hydrodynamic response can be inferred from a vast amount of data on correlations of flow coefficients of different order, the correlation of the flow coefficients with transverse momentum and similar more advanced observables.

Even at the highest currently available collision energies the hot matter produced in the collision is not entirely boost invariant. One way by which this shows up is the gradual decorrelation of the event plane angles Ψ_n (9) measured in different pseudorapidity windows Δη. The correlation coefficient r_n (η) is defined in terms of the complex flow vectors ${{\bf{q}}}_{n}={\sum }_{i}{w}_{i}{{\rm{e}}}^{{\rm{i}}n{\phi }_{i}}/({\sum }_{i}{w}_{i})$ measured in each event, where the sum runs over all particles in a selected pseudorapidity window and the w_i account for detector efficiencies [220], as

$$\begin{eqnarray}&&{r}_{n}(\eta )=\displaystyle \frac{\langle {{\bf{q}}}_{n}(-\eta ){{\bf{q}}}_{n}^{* }({\eta }_{0})+{{\bf{q}}}_{n}(\eta ){{\bf{q}}}_{n}^{* }(-{\eta }_{0})\rangle }{\langle {{\bf{q}}}_{n}(\eta ){{\bf{q}}}_{n}^{* }({\eta }_{0})+{{\bf{q}}}_{n}(-\eta ){{\bf{q}}}_{n}^{* }(-{\eta }_{0})\rangle }.\end{eqnarray} \tag{ 43 }$$

Here brackets indicate an event average, the asterisk indicates the complex conjugate, and η₀ denotes a reference pseudorapidity window usually chosen at a far-forward or far-backward pseudorapidity. Figure 30 shows the longitudinal decorrelation of the event planes Ψ_n for n = 2, 3, 4 in Xe+Xe collisions at the highest LHC energy. The n = 2 event plane is seen to decorrelate less rapidly than the n = 3, 4 planes. This phenomenon can be attributed to the predominantly geometric origin of elliptic flow, whereas the higher anisotropic flow coefficients are more sensitive to initial density fluctuations, for odd n exclusively so. Results for different energies at RHIC, shown in figure 31, show that the decorrelation occurs much more rapidly at lower energies, in agreement with the expectation of larger deviations from boost invariance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The current theoretical understanding is that the decorrelation can be traced back to rapidity-dependent fluctuations in the gluon densities of the colliding nuclei and to the gluon-gluon interactions that seed the initial state of the fireball. In addition, hydrodynamic fluctuations during the expansion of the quark-gluon plasma probably also contribute to the decorrelation. Theorists have successfully modeled these fluctuations and found rough agreement with the measured decorrelation (see e.g. [223]).

In general, collective flow is among the main observables that are sensitive to the phase transition between the hadron gas and the quark-gluon plasma. This idea can be easily understood: since the equation of state is encoded in the pressure as a function of the energy density (and possibly net baryon density), the anisotropic flow will react differently depending on the pressure profiles that are probed by heavy ion collisions as a function of time. In other words, the flow coefficients that are primarily built up in a quark-gluon plasma with its very low specific viscosity should eventually vanish when the beam energy is reduced and the plasma is no longer formed. Figure 32 shows the integrated elliptic flow of charged particles as a function of the beam energy in Au+Au collisions. Interestingly, at low collision energies the elliptic flow is generated entirely by hadronic transport while at high energies the fraction of the flow built up during the hydrodynamic stage is more than 60%. At lower collision energies, the so-called directed flow v₁ is an important observable as well. v₁ = 〈p_x /p_⊥〉 quantifies the 'bounce-off' of particles in the reaction plane. This rapidity-odd flow coefficient can be summarized by fitting the slope dv₁/dy around midrapidity. The beam energy dependence of dv₁/dy for protons is expected to show non-monotonic behavior in case of a first-order phase transition [225]. In fact, understanding the collective flow and, in particular, the directed flow poses a challenge, since it is also very sensitive to the treatment of the interactions with the spectators and the interface between hydrodynamics and transport in hybrid calculations [226].

Zoom In Zoom Out Reset image size

Long-range correlations among emitted particles in rapidity, akin to the phenomenon of anisotropic collective flow, are also found in high-multiplicity p+p and p+Pb collisions at the LHC [227–230]. Figure 33 shows two-body correlations of charged particles with transverse momenta in the range 1 GeV/c < p_T < 3 GeV/c versus azimuthal angle difference Δϕ and pseudorapidity difference Δη in high-multiplicity (N_ch > 110) events. The left panel is for p+p collisions, the right panel is for p+Pb collisions measured by CMS. The feature of interest in the near-side 'ridge' visible at Δϕ ≈ 0 and extending over the full pseudorapidity acceptance ∣Δη∣ ≤ 4. The ridge is clearly more pronounced in p+Pb collisions than in p+p collisions, which may indicate a larger degree of collectivity.

Zoom In Zoom Out Reset image size

A crucial test whether the ridge phenomenon observed in these small collision systems is caused by geometry-driven collective flow was carried out by the PHENIX collaboration, which compared p+Au, d+Au, and ³He+Au collisions at RHIC. Because the geometric structure of the light ions (p, d, ³He) is distinctively different (see left panel of figure 34), one expects significant, geometry-driven differences in the elliptic and triangular anisotropic flow coefficients v₂ and v₃. This is, indeed, what was observed by the experiment, as seen in the right panel of figure 34. In v₂ the p+Au system stands out as less elliptically deformed; in v₃ the ³He+Au system stands out as having a significantly larger triangular deformation. Signatures of collective flow are also seen in the particle-specific elliptic flow v₂(p_T ) in high-multiplicity p+Pb events at $\sqrt{{s}_{\mathrm{NN}}}=8.16$ TeV (see figure 35). The v₂ of K⁰ _s and D⁰ show a pronounced mass splitting at low p_T , but become approximately equal for p_T > 5 GeV/c. The v₂(p_T ) for identified baryons and mesons exhibit the usual valence-quark number splitting.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Additional evidence for medium-like behavior in small systems comes from the gradual approach to full thermal (grand canonical) equilibrium of multi-strange baryon production [233]. Figure 36 shows the ratio of Ω and $\overline{{\rm{\Omega }}}$ baryons to charged pions in p+p, p+Pb, and Pb+Pb collisions at various LHC energies as a function of ${\rm{d}}{N}_{\mathrm{ch}}/{\rm{d}}\eta$. The steady increase with multiplicity is interpreted as evidence that the volume and lifetime of the hot quark-gluon plasma are growing with multiplicity. It is difficult to imagine that this steady growth is not accompanied by collective flow in response to the thermal pressure inside the fireball.

Zoom In Zoom Out Reset image size

One often presented counter-indication against formation of a collectively behaving QGP in p+A collisions is that jet quenching has not been observed. Recent data from ALICE for $15\,\mathrm{GeV}/c\lt {p}_{T}^{\mathrm{ch}}\lt 50\,\mathrm{GeV}/c$ jets in high-multiplicity p+Pb collisions limit the energy loss for a R = 0.4 jet cone to less than 0.4 GeV/c, more than an order of magnitude less than the energy loss measured for R = 0.4 jets in central Pb+Pb collisions [235]. This can likely be understood as a consequence of the fact that hard scattered partons need time to evolve down to the virtuality scale Q² at which the surrounding medium interacts with them [236]. For a parton with energy p_T , the kinematic virtuality within a transverse distance L from the scattering vertex is Q² > p_T /L. On the other hand, the virtuality scale associated with rescattering in the medium is

$$\begin{eqnarray}&&{Q}_{\mathrm{med}}^{2}\leqslant {\int }_{{\tau }_{0}}^{L}{\rm{d}}\tau \hat{q}(\tau )\sim (\hat{q}/{T}^{3}){T}_{0}^{3}{\tau }_{0}\mathrm{ln}(L/{\tau }_{0})),\end{eqnarray} \tag{ 44 }$$

where the subscript '0' indicates values at the initial moment of hydrodynamical evolution. For L ≤ 3 fm in p+Pb one estimates ${Q}^{2}\gt 2{Q}_{\mathrm{med}}^{2}$, which means that the parton has left the medium before it reaches the virtuality scale at which the medium can begin to influence its QCD evolution.

The large flow gradients in the initial phase of the heavy ion reaction, especially for smaller collision system, severely strain the applicability of viscous hydrodynamics, which relies on a gradient expansion of the collective flow field. A measure of the applicability of hydrodynamics is the Knudsen number Kn, which describes the ratio of the mean free path λ to the characteristic system size L:

$$\begin{eqnarray}&&\mathrm{Kn}=\lambda /L.\end{eqnarray} \tag{ 45 }$$

Here L could either be the size of the system or the typical length scale on which the fluid properties vary. A small Knudsen number, Kn ≪ 1 implies that viscous hydrodynamics should provide a good description of the system. At early times during the collision and for small collision systems, such as p+A, this condition is violated.

Numerical simulations have shown, however, that viscous hydrodynamics provides a remarkably robust description even under conditions where the Knudsen number is not small. This behavior has been traced to the existence of hydrodynamic 'attractors', collective flow patterns that emerge as asymptotic solutions of the transport equations for many different initial conditions and even different microscopic transport dynamics (see [237] for a review). The existence of these attractors ensures remarkable insensitivity of the expansion dynamics against many uncertainties about the initial state and against the not well understood dynamics of the fluid before equilibration.

Figure 37 shows the hydrodynamic attractors for three different exactly solvable microscopic transport dynamics in the boost invariant (Bjorken) scenario. Flows for widely different initial conditions are seen to rapidly converge on the attractor within τ < 1/T, well before the hydrodynamic gradient expansion becomes reliable.

Zoom In Zoom Out Reset image size

Attractors not only apply to hydrodynamical flow patterns but also emerge in the framework of kinetic theory. In this setting, attractors have been found to smoothly interpolate between the very early stage of the collision where free streaming provides for a good approximation and the late stage when viscous hydrodynamics applies. This behavior is illustrated in figure 38, where the evolution of the energy density is shown for two different values of the 't Hooft coupling constant λ = g² N_c in the kinetic theory, which merge smoothly into the early (free streaming) and late (hydrodynamic) evolution. Altogether, these discoveries have strengthened our confidence in the robustness of the existing framework for modeling relativistic heavy ion collisions over the entire duration of the collision.

Zoom In Zoom Out Reset image size

Heavy ion reactions follow a rather complex dynamical evolution that involves switching between non-equilibrium descriptions for the initial (parton gas) and final (hadron gas) stage and viscous hydrodynamic evolution during the intermediate stage (QGP fluid). All the pieces of the theoretical framework come with their own parameters controlling, for example, the amount of initial state fluctuations, the transition criteria between descriptions, as well as the relevant physics properties, such as the equation of state and the transport coefficients. As discussed in section 5.4 the observables are often sensitive to more than one parameter in a non-trivial interplay. Note that there is a difference between the model parameters and the parameters of the system under consideration. The beam energy, centrality and other choices are the controllable settings that have to be matched between theory and experiment for meaningful comparisons.

One option to cope with such a multi-parameter, multi-observable problem is to apply the methodology of Bayesian multi-parameter analysis. The main advantage of such an approach is that one obtains a multi-dimensional sensitivity analysis as a by-product and statistically meaningful parameter extractions with quantified uncertainty. Important inputs are the statistical and systematic uncertainties of the measurements and of the theoretical models, which are sometimes difficult to assess. The pioneering studies that have applied Bayesian analysis for the first time in the field of heavy ion physics were carried out by the MADAI ¹² collaboration [241]. By now, many groups have adopted the procedures and there are further collaborative efforts between statisticians and heavy ion groups (e.g. JETSCAPE ¹³ , MUSES ¹⁴ and BAND ¹⁵ collaborations). The application of such advanced statistical techniques only makes sense, when the underlying theoretical description is established, otherwise no meaningful results can be obtained. If major physics components are missing ('unknown unknowns'), a Bayesian analysis will not produce valid results.

The following steps are usually taken:

• (i)  
  Select the theoretical models and define the parameters to be varied.
• (ii)  
  Define the prior ranges for the parameters as large as practically feasible.
• (iii)  
  Run the full evolution model at enough points in the multi-parameter space, typically determined by Latin hypercube sampling.
• (iv)  
  Select the observables to be employed in the study.
• (v)  
  Employ principal component analysis (PCA) to determine orthogonal components of the observables.
• (vi)  
  Construct an emulator for each observable trained on the full model runs.
• (vii)  
  Use Markov Chain Monte-Carlo (MCMC) techniques to determine the posterior distribution of the parameters.
• (viii)  
  Perform a closure test for verification that the posterior distribution covers results of full model runs that have not been used in the training of the emulator.

While explaining all the details of such a multi-parameter Bayesian analysis goes beyond the purpose of this review, the most important concepts will be covered here. Bayes' theorem is formulated as follows

$$\begin{eqnarray}&&P(A| B)=\displaystyle \frac{P(B| A)P(A)}{P(B)}.\end{eqnarray} \tag{ 46 }$$

The posterior distribution P(A∣B) containing the information of how likely parameter set A is under the set of observations B is given by the inverse probability of observables B given parameters A (P(B∣A)) with a prior distribution P(A) and normalized by the integral P(B) = ∫dA P(B∣A)P(A). In this manner, Bayes' relation permits to invert the conditional probabilities to one that is easier to determine. In other words, we can formulate the question which parameters fit the observations best based on the outcome of simulations for a certain set of parameters. Since the integral in the denominator is usually hard to calculate, one generally ignores the denominator and restricts the analysis to relative probabilities instead of absolute ones.

Since evaluating the theoretical model for a multitude of parameter combinations becomes quickly too expensive even for current high-performance computing facilities, surrogate models are important to enable a wide parameter scan. Gaussian process emulators are constructed in the following way. Any arbitrary functional form describing the dependence of the observables on the parameters can be approximated with a Gaussian process GP that encodes the mean and the covariance

$$\begin{eqnarray}&&f(x)\unicode{x00303}\mathrm{GP}(\mathrm{mean}(x),\mathrm{cov}(x,{x}^{{\prime} })).\end{eqnarray} \tag{ 47 }$$

In general, one can envision such a Gaussian process emulator as a way to fit the training data points with associated quantified uncertainties for the values between training points. Therefore, a crucial input is the chosen kernel function

$$\begin{eqnarray}&&k(\vec{{x}_{p}},\vec{{x}_{q}})={C}^{2}\,\exp \left(-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{s}\displaystyle \frac{| {x}_{p,i}-{x}_{q,i}{| }^{2}}{{l}_{i}^{2}}\right),\end{eqnarray} \tag{ 48 }$$

for parameter values ${\vec{x}}_{p}$ and PCA values of observables ${\vec{x}}_{q}$. The parameters C and l_i have to be chosen with care not to overfit the training data, but also to not leave too much freedom in the construction of the emulator. Figure 39 shows the behavior of unconditioned draws (left) as well as the result after proper training (right). The 95% confidence interval indicated by the gray bands is fully constrained at the training points and increases further away from them as expected. Uncertainties at the training points can be incorporated by an additional white noise kernel into the Gaussian process.

Zoom In Zoom Out Reset image size

In figure 40 the first application of Bayesian methods to heavy ion collisions is shown. A (2+1)-dimensional hybrid approach was applied to Au+Au collisions at the highest RHIC energy and to Pb+Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}=2.76$ TeV. The results were compared to data for bulk observables including particle spectra, anisotropic flow and HBT correlations. When drawing parameters weighted with their posterior distribution, it is apparent that the experimental data prefer an equation of state that is consistent with lattice QCD calculations. This is a very nice confirmation that our dynamical description actually prefers an equation of state similar to lattice QCD results based on observables at two different collider energies. This is an important finding that confirms the predicted properties of hot and dense QCD matter from experimental data.

Zoom In Zoom Out Reset image size

Figure 41 shows the posterior distributions for shear and bulk viscosity obtained from a hybrid calculation compared to experimental data. Within this work, 14 parameters were varied within the Trento+VISH2 + 1+UrQMD hybrid approach and the bulk properties of hadrons measured in Pb+Pb collisions at 2 LHC energies were taken into account. The Bayesian analysis quantifies constraints on the temperature dependence of transport coefficients and properties of the initial state at the same time. Other analyzes have focused on small systems [245] or employed different implementations for the hydrodynamic and hadronic transport evolution [149].

Zoom In Zoom Out Reset image size

Extending the list of observables to transverse momentum fluctuations and anisotropic flow correlations has been the subject of recent studies (see e.g. [246, 247]). Obtaining comprehensive quantitative conclusions on the QGP properties from all these efforts will be a major task in the future. While the different analyzes agree on major features, details like the preferred size of hot spots in the initial state or the maximum value of the bulk viscosity differ substantially between different analyzes. A very interesting further application of Bayesian analysis is demonstrated in [248], where the sensitivity of the parameters to potential future measurements in oxygen–oxygen collisions at LHC has been assessed based on 'best fit' parameters from prior studies.

While most Bayesian analysis have concentrated on the soft sector at RHIC and LHC energies, where the well established 'standard model' for the dynamical evolution of heavy ion reactions is applicable, there are a few studies targeting hard probes and lower beam energies. In [249] the jet quenching transport coefficient $\hat{q}$ governing the transverse momentum transport of a hard parton traveling through a medium has been quantified within the JETSCAPE framework [160]. Employing the MATTER and LBT energy loss modules for high and low virtuality partons, respectively, the temperature and energy dependence of $\hat{q}$ can be constrained by a comparison to data from RHIC and LHC. The results are compatible with prior constraints from the JET collaboration [166].

In [250] a Bayesian analysis for the charm diffusion coefficient was reported, and in [251] results for the beam energy dependence of the shear viscosity-over-entropy density were presented. Both results are consistent with prior works and expectations for the qualitative behavior of the transport coefficients. As the range of Bayesian model-data comparison efforts grows, one must be careful to avoid applying Bayesian analyzes in regimes where the theoretical model is under insufficient control and account properly for the variability in modeling choices (see e.g. [149]). In the future, it will be rewarding to see more and more observables from the hard and soft sector confronted with a unified theoretical description. A straightforward extension of previous work will be to do a Bayesian analysis for hard probes on a well-calibrated (by Bayesian methods) soft background.

Modern experiments with heavy ion collisions typically record many millions of events under identical conditions. Due to quantum fluctuations there will always be event-by-event fluctuations, even if the species, beam energy, and impact parameter selection is restricted. In heavy ion reactions there are many sources of fluctuations, some of them trivial (like statistical fluctuations), some of them far from trivial (like the dynamical fluctuations associated with a critical endpoint). Here, we will concentrate on fluctuations of conserved quantities at high beam energies as well as fluctuations associated with the quark-gluon plasma phase transition at lower beam energies. The event-by-event fluctuations in the initial state resulting in higher order flow coefficients have been addressed in section 5.4.

As shown in section 5.1 the system formed in heavy ion collisions can be regarded to first approximation as being in thermal equilibrium. In that case, the fluctuations of conserved quantum numbers follow the expectations from the grand canonical ensemble as known from statistical mechanics (see [252] for a review). There is one caveat: when the charge under investigation is only produced in small quantities or the volume under consideration is close to the entire system, then exact conservation laws have to be taken into account in a canonical, instead of grand canonical, approach.

The idea to measure fluctuation observables originates from their association with properties of the system created in heavy ion collisions. For example, the fluctuations of mean transverse momentum are expected to reflect temperature fluctuations. If the system is hotter the particles are emitted with larger transverse momenta, while they obtain less transverse momentum in a colder fireball. For thermal fluctuations of conserved quantum charges, such as net baryon number B, strangeness S, and electric charge Q, the expected size is sensitive to the degrees of freedom that are active in the system. Figure 42 depicts the expected differences in fluctuations of net baryon number and net charge as a function of beam energy for a quark-gluon plasma and a hadron gas. Since the quarks carry fractions of baryon number and electric charge the corresponding fluctuations are smaller. Experimental data are found to be mainly consistent with the fluctuations expected from a hadron resonance gas. One possible explanation is that the hadronization process washes the partonic fluctuations out, and finally only hadronic fluctuations are observed. This has been demonstrated in a dynamic coalescence approach, where the hadronization process was modeled microscopically in an expanding system (see the right panel of figure 42). Nowadays the interest in the mean number of pairwise produced conserved charges has shifted to correlation observables, such as balance functions for charged particles [255, 256].

Zoom In Zoom Out Reset image size

Another very intriguing application of fluctuation measurements is their direct comparison to lattice QCD calculations (see [257, 258] for a review and lecture notes on this topic). Fluctuations of conserved charges cannot be created by rapidly occurring local processes, only by much slower diffusion of particles between different regions of the fireball [253, 259]. Because diffusion proceeds much more rapidly in the QGP stage than in the hadronic gas phase, these fluctuations are expected to reflect the conditions just before hadronization. The thermal δ Q_i are related to susceptibilities ${\chi }_{n}^{{Q}_{i}}$ by

$$\begin{eqnarray}&&\langle {\left(\delta {Q}_{i}\right)}^{n}\rangle ={T}^{n}\displaystyle \frac{{\partial }^{n}}{\partial {\mu }_{i}^{n}}\mathrm{ln}Z(T,{\mu }_{i})={{VT}}^{3}{\chi }_{n}^{{Q}_{i}},\end{eqnarray} \tag{ 49 }$$

where Q_i is the conserved charge of interest, T the temperature, μ_i the corresponding chemical potential and Z the partition function. Experimentally, the moments of the conserved charge distribution can be measured by standard statistical quantities:

$$\begin{eqnarray}&&M\,(\mathrm{mean}),\quad {\sigma }^{2}\,(\mathrm{variance}),\quad S\,(\mathrm{skewness}),\quad \kappa \,(\mathrm{kurtosis}).\end{eqnarray} \tag{ 50 }$$

By associating each beam energy $\sqrt{s}$ with pairs of temperature and baryon chemical potential (T, μ_B ), these moments correspond to certain susceptibilities as stated in equation (49). To remove the volume dependence, one usually considers ratios

$$\begin{eqnarray}\begin{array}{rcl}\,S\sigma & = & {\chi }_{3}/{\chi }_{2}\quad ;\quad \kappa {\sigma }^{2}={\chi }_{4}/{\chi }_{2}\\ M/{\sigma }^{2} & = & {\chi }_{1}/{\chi }_{2}\quad ;\quad S{\sigma }^{3}/M={\chi }_{3}/{\chi }_{1},\end{array}\end{eqnarray} \tag{ 51 }$$

where the symbols refer to quantities defined in (49) and (50). Studies of the influence of volume fluctuations on such comparisons between experiment and lattice QCD calculations can be found in [260, 261].

The left panel of figure 43 shows the lattice-QCD results for the quartic baryon number susceptibility ${\chi }_{4}^{B}$ as a function of temperature and compares them with the predictions from several effective approaches. At low temperature the lattice results merge smoothly to the hadron resonance gas model (HRG), while they approach perturbative calculations (HTLpt, partial 4-loop DR) at high temperature, eventually approaching the free parton gas limit.

Zoom In Zoom Out Reset image size

The comparison of fluctuation measurements to susceptibilities of conserved quantum numbers calculated from lattice QCD is an alternative method to determine the chemical freeze-out conditions in terms of temperature and chemical potential. The results for net baryon number and net charge fluctuations are in rough agreement with the findings according to the thermal model (see 5.1). This agreement, which is illustrated in the right panel of figure 43, is expected since the lattice results in that regime correspond to a hadron-resonance gas.

In [187] it was suggested that comparisons between lattice calculations and strangeness fluctuations indicate a somewhat higher decoupling temperature than the one for light hadrons. This is in line with the expectations from microscopic models, where strange particles have a smaller cross section with other hadrons than the non-strange particles, most prominently protons and pions.

While it is exciting to directly compare ab initio lattice calculations to experimental data, one has to be aware of the limitations: depending on the kinematic cuts of the measurement, the comparison to a grand canonical ensemble calculation may be appropriate or not: typically only net proton fluctuations are measured and the mapping to net baryon number fluctuations carries uncertainties; also final state interactions and non-equilibrium effects may affect fluctuations.

The second important application of fluctuation observables is related to their expected sensitivity to the QCD critical endpoint. Finding signatures of the critical endpoint of the first-order phase transition between quark-gluon plasma and hadron gas is one of the main motivations for the heavy ion program at finite densities. The theory of phase transition dynamics predicts that the correlation length increases when the system passes through a critical region. In particular, higher moments of the distributions of conserved quantum numbers are related to a higher power of the increased correlation length. In the idealized equilibrium scenario, the correlation length as well as the higher moments diverge. In heavy ion collisions this divergence is prevented by finite-size and finite-lifetime effects [265].

The fluctuation moments (51) can be determined experimentally as the cumulants C_n of the measured particle distributions [266]. For example, the kurtosis κ is obtained by evaluating the event average

$$\begin{eqnarray}&&\kappa {\sigma }^{4}={C}_{4}=\langle {\left({\rm{\Delta }}N\right)}^{4}\rangle -3{\left(\langle {\left({\rm{\Delta }}N\right)}^{2}\rangle \right)}^{2}\end{eqnarray} \tag{ 52 }$$

with the variance ${\sigma }^{2}={C}_{2}=\langle {\left({\rm{\Delta }}N\right)}^{2}\rangle$. Figure 44 (left) shows the scaled kurtosis κ σ² = C₄/C₂ in the phase diagram of nuclear matter. There are interesting structures and sign changes expected in the region around the critical endpoint of the QCD phase transition, which are analogous to those known in standard liquid–gas phase transitions. If one follows a typical freeze-out line in the phase diagram, figure 44 (right) depicts the expected beam energy dependence of the kurtosis involving a peak and then a dip structure when going from high to low beam energies.

Zoom In Zoom Out Reset image size

Experimentally, measurements of the excitation function of the kurtosis of the net proton distribution have been carried out within the (first) RHIC beam energy scan program by the STAR collaboration and by the HADES collaboration at GSI (see figure 45 for a compilation of results). These are extremely challenging measurements since the impact of efficiencies and kinematic cuts for a multi-particle correlation measurement has to be controlled. As can be seen in figure 45 the data indicate a non-trivial behavior as a function of collision energy that is not compatible with the prediction of a transport (UrQMD) calculation which does not contain information about a critical point. The global net baryon number conservation influences the results, but not enough to bring them into agreement with experimental data. NA61 has looked at lower-order fluctuation observables as a function of beam energy and system size but could, so far, find no sign of critical behavior. In the future, higher precision data from the Beam Energy Scan II at RHIC, CBM at FAIR, MPD at NICA and other heavy ion physics programs at lower beam energies will complement the existing measurements.

Zoom In Zoom Out Reset image size

On the theory side, there have been many developments targeting a dynamic non-equilibrium evolution through a critical endpoint (see [272] for a recent summary). All of these efforts are based on extending the fluid dynamics description to include thermal fluctuations. One approach involves adding thermal noise to the relativistic hydrodynamic evolution, which is numerically challenging. A complementary approach propagates the two-particle correlations on top of a hydrodynamic background ('Hydro+' formalism [273]). Both calculations agree in their finding that the non-equilibrium evolution has significant effects on the magnitude and behavior of the kurtosis, even though they are carried out in simplified settings. A full (3+1)-dimensional non-equilibrium evolution including the critical dynamics still remains a future challenge for a quantitative understanding of the experimental measurements.

Event-by-event fluctuations in heavy ion physics can also be used to select events of interest. The 'event shape engineering' technique groups, for example, events according to certain features like the magnitude of radial flow or certain anisotropic flow coefficients to provide further handles beyond centrality and beam energy. Machine learning techniques may make it possible to access interesting information by feeding information from single events into an artificial intelligence (AI) system. Of course this has to be done with great care, since the machine is not smarter than the best theoretical model on the market that was used to train the neural network. Systematic uncertainties inherent in such approaches are difficult to assess.

Hadron production from heavy ion collision in the transverse momentum region below a few GeV/c exhibits two striking features: (1) The baryon-to-meson ratio, for both protons and antiprotons, in central Au+Au grows steadily with p_T for p_T ≤ 3 GeV/c reaching a value three times as large as in peripheral collisions as shown in figure 46 [274]. (2) The elliptic flow coefficient v₂(p_T ) for mesons and baryons shows a distinctly different behavior, with the baryon v₂ saturating at larger values than the meson v₂ [275], as shown in figure 47 for several species of identified hadrons [276]. Especially noteworthy is the observation that ϕ-mesons (green symbols), which have approximately the same mass as nucleons, follow the proton data points at low p_T indicating hadronic collective flow, but follow the data points for other (much lighter) mesons at p_T > 2 GeV/c indicating that the collective flow is established at the quark level. This finding is consistent with the prediction that the formation of low-p_T hadrons occurs mainly during the bulk QGP hadronization, whereas higher-p_T hadrons are preferentially emitted from the QGP surface where quark recombination dominates (see figure 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both phenomena can be explained by the mechanism of hadron formation by quark recombination (or coalescence) from the quark-gluon plasma, in which each valence quark inherits the collective flow properties of the QGP fluid. Since baryons contain three valence quarks whereas mesons contain only two, baryons experience a stronger push from collective flow towards higher-pT than mesons [277–279]. For the elliptic flow v₂(p_T ) this mechanism implies a scaling law with valence quark number n_q [280]:

$$\begin{eqnarray}&&{v}_{2}({p}_{T})={n}_{q}{v}_{2}^{({\rm{q}})}({p}_{T}/{n}_{q}),\end{eqnarray} \tag{ 53 }$$

where ${v}_{2}^{({\rm{q}})}({p}_{T})$ denotes the elliptic flow coefficient for (anti-)quarks in the QGP.

Strictly speaking, the scaling with p_T /n_q can only be justified in the kinematic domain where hadron masses can either be neglected or described additively by constituent quark masses. In order to apply the scaling law heuristically over a wider momentum range, especially down to small momenta p_T , it is customary to compare the elliptic flow of different hadrons as a function of the scaled transverse mass (m_T − m₀)/n_q , where m₀ is the rest mass of the hadron and ${m}_{T}({p}_{T})=\sqrt{{p}_{T}^{2}+{m}_{0}^{2}}$. This version of the valence quark scaling law has been found to be remarkably well obeyed by a large number of hadron species over a wide collision energy range. Two examples from the recent RHIC beam energy scan are shown in figure 48. A review of theoretical and experimental aspects of the quark recombination mechanism can be found in [281].

Zoom In Zoom Out Reset image size

The valence quark scaling (53) has also been observed in identified particle emission patterns at the LHC, where the scaling is observed to hold even for the higher flow anisotropy coefficients v₃ and v₄ [276]. Figure 49 shows valence quark-number scaled flow coefficients v_n (p_T /n_q )/n_q for n = 2, 3, 4 for several identified hadrons in Pb+Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}=5.02$ TeV.

Zoom In Zoom Out Reset image size

The scaled flow coefficients in figure 49 exhibit broad peaks around p_T /n_q ≈ 1.5 GeV/c. In the fragmentation-recombination scenario, this peak corresponds to a gradual transition to the fragmentation dominated regime. In the transition region it is plausible that quarks from a parton shower recombine with thermal, collectively flowing partons [283]. An implementation of this idea in a dynamical model (EPOS) describes hadron formation at intermediate values of p_T as string (color flux-tube) fragmentation in the presence of a thermal parton fluid [284, 285]. Figure 50 compares ALICE data [286] for the p_T -dependence of the hyperon-to-kaon ratio N_Λ/N_K for different centrality windows in Pb+Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}=2.76$ TeV with results from the EPOS model that accounts for recombination of the leading parton with thermal partons as well as fragmentation by quark-pair production from the vacuum [285].

Zoom In Zoom Out Reset image size

At low momenta (p_T < 1.5 GeV/c) hadronic rescattering affects baryon and meson flow differently and amplifies the mass splitting observed in the unscaled elliptic flow v₂(p_T ) [287]. Pions move much faster than baryons and push them out to larger p_T ('pion wind', see also the end of section 5.2), while heavy baryons have the opposite effect on pions and other mesons. Note that this mechanism does not work for ϕ-mesons, as these do not have large cross sections with pions or nucleons [288]. The increase of the mass splitting in v₂(p_T ) is clearly visible in figure 51, which shows v₂(p_T ) for pions, kaons, and protons calculated with and without hadronic rescattering in comparison with PHENIX data from Au+Au collisions at $\sqrt{{s}_{\mathrm{NN}}}=200\,\mathrm{GeV}$ [289].

Zoom In Zoom Out Reset image size

The emission of hadrons at high transverse momentum (p_T > 6 GeV/c) in relativistic heavy ion collisions is suppressed because high-momentum quarks and gluons lose energy when they propagate through the quark-gluon plasma [290, 291]. For light quarks or gluons the dominant energy mechanism is gluon radiation in association with scattering off a virtual gluon in the QGP. In the BDMPS-Z formalism of multiple scattering the scattering power of the QGP is encoded in the parameter $\hat{q}$, which describes the average squared momentum exchange with the medium per unit path length, $\hat{q}=d\langle {q}_{T}^{2}\rangle /{\rm{d}}x$ [292–294]. The BDMPS-Z approach is well suited to describe energy loss in a thick QGP. The GLV formalism, which is based on an opacity expansion, is better suited for a thin QGP. The higher-twist formalism [157] aims at the description of the full virtuality evolution of a jet created by an energetic quark or gluon and is expected to apply to both, short and long path lengths (see also [295]).

Experimental evidence that high-p_T hadron suppression is a final-state effect comes most directly from the comparison of the nuclear modification factors R_AA (7) for hadrons with that for photons [296]. Figure 52 shows that R_AA ≈ 1 for photons in Au+Au collisions at the top RHIC energy, whereas R_AA ≈ 0.25 for pions and η-mesons. Similar results have been obtained for Z⁰ or W^± bosons in p+Pb and Pb+Pb collisions at LHC [297–301], where the data are consistent with nuclear modification of the initial-state parton distributions but show no final-state suppression of the boson yield when identified by their leptonic decay modes. However, as Z⁰ bosons decay within 0.1 fm/c, the hadron showers created in the hadronic decay mode ($q\bar{q}$) may provide for a calibrated probe of the jet-medium interactions [302].

Zoom In Zoom Out Reset image size

Schematically, the hadron spectrum can be expressed as a convolution of the parton distribution functions ${f}_{i}^{(A)}(x)$ in the colliding nuclei with the hard QCD scattering cross section and a fragmentation function ${D}_{i\to h}^{(\mathrm{med})}(z)$ that is modified by the medium

$$\begin{eqnarray}&&\frac{{\rm{d}}{N}_{h}}{{\rm{d}}{p}_{T}}({p}_{T})={f}_{i}^{(A)}({x}_{1})\otimes {f}_{j}^{(A)}({x}_{2})\otimes \frac{{\rm{d}}{\sigma }_{{ij}}}{{\rm{d}}{p}_{T}}({p}_{T}/z)\otimes {D}_{i,j\to h}^{(\mathrm{med})}(z).\end{eqnarray} \tag{ 54 }$$

It is usually assumed that the hadronization of the quarks and gluons contained in the evolving parton cascade occurs in vacuum outside the medium. For hadrons deriving from the leading parton in this cascade, the medium modification can then be expressed as a shift in the momentum fraction in the vacuum fragmentation function:

$$\begin{eqnarray}&&{D}_{i\to h}^{(\mathrm{med})}(z)\approx \left\langle {D}_{i\to h}^{(0)}\left(\displaystyle \frac{z}{1-z{\rm{\Delta }}{p}_{T}/{p}_{T}}\right)\right\rangle ,\end{eqnarray} \tag{ 55 }$$

where Δp_T is the momentum loss of the leading parton and 〈⋯〉 indicates an average over the position and orientation of the hard scattering event. The magnitude of the suppression thus depends not only on the amount of energy loss but also on the steepness of the unmodified hadron spectrum. As the spectrum becomes flatter at high collision energy, the same energy loss causes less suppression. In spite of this effect, experimental data from RHIC and LHC confirms that the suppression increases with collision energy, as shown in figure 53 for neutral pions with p_T > 6 GeV/c, implying a strong increase in the energy loss. This observation agrees with expectations, as the energy loss parameter $\hat{q}$ grows rapidly with temperature: $\hat{q}/{T}^{3}\approx 2-5$ (see figure 54).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The increase with N_part reflects the strong path-length dependence of the radiative energy loss which, for a medium of constant density, is approximately given by [294]:

$$\begin{eqnarray}&&{\rm{\Delta }}E=-\displaystyle \frac{1}{2}{C}_{2}{\alpha }_{s}\hat{q}{L}^{2},\end{eqnarray} \tag{ 56 }$$

where C₂ denotes the SU(3) Casimir operator for the energetic parton. Figure 55 shows an attempt to relate the fractional momentum loss S_loss of energetic partons to the average squared path-length L² of the parton inside the medium. S_loss is obtained by comparing the spectrum of hadrons measured in $\sqrt{{s}_{\mathrm{NN}}}=200\,\mathrm{GeV}$ Au+Au collisions with the binary collision-scaled hadron spectrum in p+p collisions in different centrality windows. The data are consistent with a linear dependence on L² as predicted by (56).

Zoom In Zoom Out Reset image size

An independent assessment of the path-length dependence of parton energy loss can be obtained by measuring the azimuthal anisotropy of the nuclear suppression with respect to the collision plane. This anisotropy can be expressed by the Fourier coefficient v₂(p_T ), defined as

$$\begin{eqnarray}&&\frac{{\rm{d}}N({p}_{T})}{{\rm{d}}\phi }\propto 1+2{v}_{2}({p}_{T})\cos (2\phi )\end{eqnarray} \tag{ 57 }$$

as the 'elliptic flow' coefficient. The average difference in path length dL for partons emitted perpendicular to the collision plane compared with those emitted along the plane is quite large as shown in the left panel of figure 56, implying that a significant dependence on the emission angle relative to the collision plane is to be expected. This expectation is confirmed by data from Au+Au collisions at RHIC, see the right panel of figure 56 which shows a sizable value of v₂ for hadrons up to 10 GeV/c momentum. Data from Pb+Pb collisions at LHC for much higher-pT, shown in figure 57, reveal a strong correlation with the elliptic flow coefficient v₂ measured in the low-p_T region and indicate that the anisotropy at high p_T has the same geometric origin.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Unless they are absorbed by the QGP, the gluons that are radiated by a fast parton in the medium, remain part of the full jet. What fraction of the radiated energy is recovered in measurements of the full jet energy depends on the opening angle R (the 'jet radius') that is used to define the jet. Typical values used for such studies are R = 0.2 − 0.6. Smaller jet radii imply a larger energy loss. This trend is readily apparent in figure 58, which shows the relative magnitude of R_CP(p_T ) for jets with total transverse momentum p_T and radius 0.2 ≤ R ≤ 0.5. (See the text below equation (7) for an explanation of R_CP.)

Zoom In Zoom Out Reset image size

Full jet quenching exhibits many of the same phenomena as the suppression of single hadrons at high p_T , except that all results depend quantitatively on the jet radius R. This means that full jet quenching measurements are not only sensitive to longitudinal energy loss, but also to the angular redistribution of the energy within the jet [309]. Measurements of full jet suppression thus enable more differential measurements, e.g. the study of how quenching modifies the jet shape in terms of the longitudinal momentum fraction of a hadron within the jet, ${z}_{h}={p}_{T}^{h}/{p}_{T}^{\mathrm{jet}}$, and the relative angle r < R of the momentum of a hadron with respect to the jet axis (see [310] for a review on jet measurements).

Often the jet shape ρ(ξ, r) is expressed in terms of the variables $\xi =\mathrm{ln}(1/z)$ and r. Examples of the modification of the jet shape in central Pb+Pb collisions compared with p+p collisions are shown in figure 59. The data show 'softening' of the shape of the jet in terms of a redistribution of the energy in the jet to smaller z and larger angles r. This is precisely the pattern expected from in-medium gluon radiation, which involves lower parton virtuality than vacuum radiation and does not exhibit the same angular ordering that suppresses low-z, large-angle radiation.

Zoom In Zoom Out Reset image size

Jets produced back-to-back with an isolated high-energy photon are preferentially initiated by a quark that has been scattered to large angles, e.g. by the Compton-like process gq → γ q. Simulations with an event generator, such as Pythia, Herwig or Sherpa, show that photon tagging enhances the fraction of jets initiated by hard scattered quarks over those initiated by gluons, from 35%–50% for inclusive jets to 70%–80% for photon-tagged jets at LHC energies over the jet energy range 60 GeV <E_T < 200 GeV [313]. A comparison between photon-tagged jets and inclusive jets thus allows to probe the color charge dependence of parton energy loss expressed by the dependence of the energy loss (56) on the color-SU(3) Casimir operator (C₂ = 4/3 for quarks and C₂ = 3 for gluons). Data from ATLAS shown in figure 60 confirm the expectation that jets initiated by quarks are less suppressed than those initiated by gluons, manifested in the stronger suppression of inclusive jets.

Zoom In Zoom Out Reset image size

Jets generally occur in pairs (di-jets) where one jet balances the transverse momentum of the other. This means that the relative di-jet distribution is strongly peaked at 180° in azimuthal angle. While highly correlated in azimuthal angle ϕ, di-jets are not strongly correlated in pseudorapidity η but separated by a variable gap ${\rm{\Delta }}\eta \sim \mathrm{ln}({x}_{2}/{x}_{1})$, where x_i are the momentum fractions of the colliding partons that produce the di-jet. In order to localize both partners of the di-jet, one employs two trigger particles, one for each jet, or two calorimeter-based triggers. The particle distribution in each jet is then measured relative to the (η,ϕ) coordinates of the respective trigger and denoted as [314]

$$\begin{eqnarray}&&\frac{1}{{N}_{\mathrm{trig}}}\frac{{{\rm{d}}}^{2}N}{{\rm{d}}({\rm{\Delta }}\eta ){\rm{d}}({\rm{\Delta }}\phi )}.\end{eqnarray} \tag{ 58 }$$

In order to isolate di-jets in heavy ion collisions one typically also imposes a lower p_T -cutoff on the included particles and subtracts the randomized background from mixed minimum-bias events [314].

There are two main observables that have been studied for di-jets in A+A collisions. One is the additional nuclear suppression of high-p_T hadron pairs (di-hadrons) relative to the suppression of single inclusive hadrons. Such a suppression is to be expected, because both di-jet precursor partons propagate through the quark-gluon plasma and lose energy. The additional suppression for inclusive hadron pairs is expressed in terms of the quantity

$$\begin{eqnarray}&&{I}_{\mathrm{AA}}=\displaystyle \frac{{R}_{\mathrm{AA}}^{\mathrm{di}-\mathrm{jet}\,\mathrm{triggers}}}{{R}_{\mathrm{AA}}^{\mathrm{single}\,\mathrm{triggers}}}.\end{eqnarray} \tag{ 59 }$$

The left panel of figure 61 shows the I_AA in Au+Au collisions at $\sqrt{{s}_{\mathrm{NN}}}=200\,\mathrm{GeV}$, relative to the baseline from d+Au collisions at the same energy. As one can see, the additional suppression in central collisions is comparable to the single suppression factor R_AA shown in figure 53. This is to be expected as the average path lengths of both scattered partons in the medium are comparable, and thus both partons suffer similar energy loss. The right panel of figure 61 shows the analogous di-jet suppression factor ${I}_{\mathrm{AA}}^{\mathrm{dijet}}$ for the complete jets measured in Pb+Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}=5.02$ TeV. This quantity is to be compared with the single jet ${R}_{{AA}}^{\mathrm{jet}}$ shown in figure 60. The additional suppression grows with centrality, but is generally less severe than the suppression observed for single jets.

Zoom In Zoom Out Reset image size

The other observable is the di-jet asymmetry A_J or, equivalently, the di-jet imbalance ratio x_J, defined as

$$\begin{eqnarray}&&{a}_{{\rm{J}}}=\displaystyle \frac{{p}_{T,1}-{p}_{T,2}}{{p}_{T,1}+{p}_{T,2}},\qquad {x}_{{\rm{J}}}=\displaystyle \frac{{p}_{T,2}}{{p}_{T,1}},\end{eqnarray} \tag{ 60 }$$

where p_T,1 and p_T,2 denote the transverse momenta of the leading and sub-leading jet, respectively. The di-jet balance ratio x_J in central Pb+Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}=2.76$ TeV is shown in figure 62 in comparison with the same ratio in p+p collisions. The trigger conditions were either p_T,1 > 100 GeV (left panel) or p_T,1 > 200 GeV (right panel) and p_T,2 > 25 GeV for jets within ∣η∣ < 2.1. For the lower trigger energy (left panel) the di-jet balance ratio distribution in Pb+Pb collisions is peaked around x_J ≈ 0.5 and differs strongly from the distribution observed in p+p collisions.

Zoom In Zoom Out Reset image size

This behavior can be interpreted as follows: when the jet is produced well outside the center of the fireball, one of the jets traverses a substantially shorter distance through the medium than the other. This causes a larger energy loss, which is reflected in a ratio x_J significantly smaller than unity. Interestingly, the difference between Pb+Pb and p+p collisions shrinks with increasing trigger threshold until the distributions are statistically indistinguishable for p_T,1 > 200 GeV. The same trend is found when one goes from central to peripheral collisions [316, 317].

Hadrons containing heavy quarks (c or b-quarks) are of interest both theoretically and experimentally. One distinguishes hadrons with open heavy flavor, such as D- and B-mesons or Λ_c baryons, and those with hidden heavy flavor, such as charmonium ($c\overline{c}$) and bottomonium ($b\overline{b}$). Open heavy flavor hadrons in their ground states decay weakly and live long enough so that their decay can be identified by micro-vertex detectors. Hadrons with hidden heavy flavor can decay either electromagnetically or by gluon-mediated strong interaction. For states with quantum number J^PC = 1⁻⁻, e.g. J/ψ and ϒ, the strong decay channel (J/ψ(ϒ) → ggg) is suppressed by a factor $({\pi }^{2}-9){\alpha }_{s}^{3}$ involving an algebraic near cancellation and a higher power of alpha_s . Accordingly, the decay into lepton pairs has a rather large branching ratio, which makes these states readily detectable. An extensive survey of heavy quark physics in relativistic heavy ion collisions (as of 2015) can be found in [318, 319].

On the theoretical side, heavy quarks are interesting because they are almost exclusively produced during the initial stage of the reaction by hard QCD processes, mainly $g+g\to Q+\overline{Q}$. They may or may not subsequently thermalize, which is an interesting question by itself, but their number remains essentially conserved from then until hadronic freeze-out. Inside the QGP, when light quarks are deconfined, hadrons containing both heavy and light quarks cannot exist. During hadronization, such hadrons are created by recombination of deconfined heavy quarks with light quarks or antiquarks. However, hadrons containing solely heavy quarks may survive under conditions not too far above the deconfinement threshold because their binding radii are small and their binding energies are large compared with the temperature.

The question, above which temperature ${T}_{{\rm{m}}}^{(H)}$ a specific heavy heavy quark bound state H 'melts', has been studied in great detail using lattice QCD. Initial investigations focused on static color screening studies [320] but more recently the focus has shifted to the investigation of dynamic properties encoded in the spectral functions [321, 322], which include non-static effects, such as ionization by thermal gluons. More generally, the heavy quark mass provides for a large momentum scale M_Q ≫ Λ_QCD that enables various effective field theory approaches to QCD, known as heavy quark effective theory (HQET) or nonrelativistic QCD (NRQCD). In combination with HTL perturbation theory techniques, these approaches form a rigorous theoretical framework for the calculation of transport properties of heavy quarks in the QGP, including the formation and destruction of quarkonia [323].

Originally, the inclusive study of the transport properties of heavy quarks in the QGP mainly relied on the measurement of leptons (e, μ) emitted in their semi-leptonic weak decays [324]. The discrimination between leptons from b-decays and those of c-decays further requires the identification of the decay vertex where one uses, on a statistical basis, the property that hadrons containing single b-quarks (open beauty hadrons) have a longer average lifetime than those containing single c-quarks (open charm hadrons). The comparison of the inclusive lepton spectrum measured in A+A collision with the binary collision-scaled spectrum measured in p + p collisions yields information about the transport of heavy quarks in the QGP. This information is usually presented as nuclear modification factor R_AA plotted as a function of p_T . In order to relate to nuclear modification of the heavy quark spectrum, the lepton spectrum requires unfolding with the decay spectrum of the parent hadrons in their rest frame, which is an ill-defined procedure. One therefore usually compares the data with calculations that include the weak decays of open heavy flavor hadrons.

Figure 63 shows that lepton spectra from heavy flavor decays exhibit similar nuclear modification features as those of light hadrons. The left panel, which shows flavor separated R_AA for leptons from c versus b decays in Au+Au collisions at RHIC [325, 327], provides evidence that c-quarks experience strong rescattering in the QGP, resulting in a suppression which is comparable to that of light quarks. The suppression effect for leptons from b-quark decays is significantly smaller. This is expected, as the energy loss of a b-quark in collisions with thermal partons is reduced by a factor m_c /m_b when compared with that of c-quarks. In addition, radiative energy loss by heavy quarks exhibits a dead-cone effect [161, 328] that increases with quark mass. Although medium-induced gluon radiation is predicted to partially fill the radiation dead-cone, a mass dependent reduction of the radiative energy loss is predicted [329].

Zoom In Zoom Out Reset image size

The right panel of figure 63 shows the muon R_AA for unseparated b- and c-decays in central $\sqrt{{s}_{\mathrm{NN}}}=2.76,5.02$ TeV Pb+Pb collisions at the LHC [326]. Again, the nuclear modification exhibits similar features as that measured for light charged hadrons with a minimum around p_T ∼ 10 GeV/c. A similar plot for $\sqrt{{s}_{\mathrm{NN}}}=5.02$ TeV [330] shows that the observed trend is generally well reproduced by theoretical calculations that include collisional and radiative energy loss.

With the advent of high precision tracking and vertexing of charged particles in heavy ion collisions, high statistics spectra of identified open charm hadrons have become available that do not require unfolding. Open charm hadrons that exhibit a reconstructable hadronic weak decay mode include D⁰, D^±, and Λ_c . The best data are available from ALICE for identified D-mesons for Pb+Pb collisions at LHC [331]. The left panel of figure 64, showing R_AA(p_T ) in central collisions, confirms that all species of D-mesons are equally suppressed, in agreement with the hypothesis that the suppression mainly reflects the energy loss of c-quarks in the QGP. The right panel of figure 64 shows that the suppression strongly depends on centrality as expected from the path-length and energy density dependence of the energy loss.

Zoom In Zoom Out Reset image size

The strong suppression of D-mesons, similar to the pattern observed for light hadrons, raises the question whether c-quarks thermalize in the QGP and participate in the collective flow and if so, to what degree. A partial answer to this question is afforded by the measurement of the elliptic flow coefficient v₂ for D-mesons. Results for v₂(p_T ) of identified D-mesons measured by ALICE are shown in figure 65 in comparison with v₂(p_T ) for charged pions in two centrality windows. Except at the lowest measured p_T the D-mesons exhibit almost the same amount of elliptic flow as pions, which indicates that they participate in the overall collective flow of the QGP. The reduced v₂ at low p_T is expected because D-mesons have almost twice the mass of a proton and thus should show an even stronger kinematic reduction of v₂ at low p_T than protons.

Zoom In Zoom Out Reset image size

Bound states of heavy quarks (charmonium, Upsilon, and B_c ) are sensitive to the color screening length in the QGP [334]. When the color screening length r_D , also called the color Debye length, is shorter than the radius of the heavy quark bound state, the bound state dissolves ('melts') in the QGP and becomes a broad resonance. For each bound state there is a characteristic threshold temperature [320], also called the Mott temperature T_M [335]. In addition to color screening, the other important contribution to this process is ionization by thermal gluons in the QGP.

The two determinants of the Mott temperature are the radius of the bound state and its binding energy. In the non-relativistic limit and using natural units, these are given by ${R}_{Q\bar{Q}}={N}^{2}{\left({\alpha }_{s}{m}_{Q\bar{Q}}\right)}^{-1}$ and ${B}_{Q\bar{Q}}={N}^{2}{\alpha }_{s}^{2}{m}_{Q\bar{Q}}$, where ${m}_{Q\bar{Q}}$ is the reduced mass and N ≥ 1 is the principal quantum number of the Coulombic bound state. For charmonium the 1s and 2s states are bound in the vacuum (J/ψ and $\psi ^{\prime}$); for bottomonium the 3s state is also bound (ϒ, ${\rm{\Upsilon }}^{\prime}$, and ϒ''.) One thus expects that when the temperature is raised above the critical temperature T_c , first the $\psi ^{\prime}$ and ϒ'' states melt (close to T_c ), then the J/ψ and ${\rm{\Upsilon }}^{\prime}$ states (around 1.5T_x ) and eventually the ϒ ground-state (slightly above 2T_c). This predicted phenomenon is known as 'sequential melting'.

Over the past three decades the suppression of quarkonium production in heavy ion collisions, compared with scaled proton–proton collisions, has been measured in great detail. The suppression is commonly expressed in terms of the ratio R_AA, similar to the suppression of jet production. Results for charmonium (J/ψ) suppression in Au+Au collisions at RHIC and in Pb+Pb collisions at LHC are shown in figure 66. The suppression of bottonium (ϒ(ns)) production in Pb+Pb collisions at LHC is shown in figure 67.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The interpretation of the R_AA data is complicated by several effects. The two most important ones are:

• The primary production process for heavy quark pairs, ${gg}\to Q\bar{Q}$, is suppressed in nuclear collisions because the nuclear gluon distribution at small Bjorken-x is screened ('shadowed'). This effect, which can be studied in p(d) + A collisions, is mainly observed at low transverse momentum p_T , as is visible in figure 68 (left panel).
• In Pb+Pb collisions at LHC energies, $c\bar{c}$ pairs are copiously produced and are thought to be thermalized in the QGP. These pairs can coalesce into J/ψ and $\psi ^{\prime}$ mesons when the QGP hadronizes. This regeneration mechanism leads to a striking difference in the centrality dependence of J/ψ suppression at LHC (right panel of figure 66), which is constant over a wide centrality range, compared with the centrality dependence at RHIC (left panel of figure 66), which shows a continued drop toward central collisions. The different behavior is caused by a relative enhancement of charmonium production at low p_T at LHC, which is shown in the right panel of figure 68. The same mechanism, albeit less pronounced, may already be at work in Au+Au collisions at RHIC at central rapidity (∣y∣ < 0.35) where J/ψ production is found to be less suppressed than at forward rapidity (∣y∣ > 1.2). The stronger suppression of $\psi ^{\prime}$ compared with J/ψ seen in this figure is also evidence for the sequential melting concept.

Zoom In Zoom Out Reset image size

Electromagnetic probes are theoretically interesting, because they do not interact strongly and therefore are emitted without final-state effects from all stages of the reaction. Many hadrons and other particles are also produced during the whole evolution of the heavy ion reaction, but the main advantage of electromagnetic probes is that they reach the detector undisturbed. A dilepton pair or a photon escapes even the hot and dense reaction zone without further interaction. However, their relatively weak interaction makes them experimentally challenging, because their yields are small, and the experimental signal often has large backgrounds from weak or electromagnetic decays of abundant hadrons that create photons or leptons in the final state (e.g. π⁰ → γ γ, D⁰ → K⁻ μ⁺ ν_μ ). Any interpretation of experimental results relies therefore on theoretical input on the origin of the contributions from the different sources and stages of the reaction.

There are two main physics questions that can be addressed with dilepton and photon measurements. By investigating the invariant mass spectrum of dileptons emitted from vector mesons one can get insight about the properties of the spectral functions of resonances in the medium (see [346] for a recent review). The ρ-meson is of special interest, since the idea is to study chiral symmetry restoration by observing the spectral functions of the ρ and its chiral partner a₁ become degenerate. More recently, the focus has shifted to signatures of chiral mixing, since the measurement of the a₁ spectral function is out of reach. The second main topic is the idea of extracting a temperature of the quark-gluon plasma from the thermal radiation. For this purpose, photons might seem more straightforward, but one has to account for a blue shift in the photon spectrum due to radial flow. The slope of dilepton spectra in the invariant mass region between the ϕ-meson and the J/ψ, 1.1GeV < m_inv < 3 GeV provides a more direct measure of the thermal radiation from the plasma, because m_inv is not affected by flow.

Figure 69 depicts the most precise dilepton measurement in a heavy ion environment to date. The invariant mass spectrum of dimuons was recorded for In–In collisions at E_lab = 158A GeV at the CERN-SPS. After contribution from the 'cocktail' of known hadronic sources has been subtracted the remaining excess yield (see the left panel of figure 69) provides a measure of the spectral function of the ρ-meson inside the hot and dense medium. From these results, it can be concluded that the ρ-meson is strongly modified in the medium and mainly broadened, while a mass shift of the pole position has not been observed. The differential dimoun measurements make it possible to fit transverse momentum spectra in different invariant mass bins and extract the effective temperature as the inverse slope shown in the right panel of figure 69. In the low mass region (LMR) the extracted values agree with the ones from the hadronic spectra, consistent with a hadronic origin, while above 1 GeV in invariant mass the thermal slope saturates and suggests emission from an equilibrated quark-gluon plasma. In this mass region, there is a correlated background from heavy quark decays that needs to be carefully subtracted.

Zoom In Zoom Out Reset image size

Moving to higher beam energies the STAR collaboration has measured dielectron invariant mass spectra as shown in figure 70. The excess yield above the hadronic cocktail emission indicates that the spectral function of the ρ-meson is also broadened at these higher beam energies. Theoretical calculations include thermal dilepton rates from effective field theory folded with a fireball model [348, 349], coarse-grained UrQMD transport calculations with the same thermal dilepton rates [350, 351] and fully microscopic non-equilibrium calculations by within the PHSD approach [145, 352]. In the future, ALICE is expected to be also able to measure precise, background subtracted dilepton spectra in Pb+Pb collisions at LHC energies.

Zoom In Zoom Out Reset image size

At GSI, the HADES experiment is dedicated to investigating dilepton emission from elementary and heavy ion reactions. In this baryon-dominated environment the ρ-meson is mainly modified due to its interactions with the baryonic resonances. Figure 71 shows the extracted thermal emission from Au+Au collisions at E_lab = 1.23A GeV. There are clear indications of medium modifications in comparison with the vacuum environment indicated by the agreement with coarse-grained (CG) transport calculations including medium-modified spectral functions for vector mesons. The calculation (see right panel of figure 71) shows this difference explicitly for the radiation from the ρ- and ω-mesons, the main vector mesons contributing in this mass range. In the future, HADES and the CBM experiment will measure the excitation function of thermal dilepton emission with the goal to identify signatures of a first-order phase transition.

Zoom In Zoom Out Reset image size

Photon production from heavy ion collisions is dominated by sources from hadronic decays, most prominently the π⁰ that decays to photons with an almost 100% branching ratio. Therefore, it is very challenging to experimentally extract the primary photons, usually called 'direct photons'. The first measurement of a direct photon spectrum was accomplished by the WA98 collaboration at SPS [364]. More recently, PHENIX and ALICE have published transverse momentum spectra and flow measurements for direct photons. A compilation of the spectra is shown in figure 72. In order to compare spectra measured in different centrality windows they have been divided by a normalization factor ${\left({\rm{d}}{N}_{\mathrm{ch}}/{\rm{d}}\eta \right)}^{1.25}$. This factor is motivated by the fact that the charge particle multiplicity $({\rm{d}}{N}_{\mathrm{ch}}/{\rm{d}}\eta )$ scales with the number of participant nucleons, N_part, whereas the direct photon yield is expected to scale with the number of binary nucleon–nucleon collisions ${N}_{\mathrm{coll}}\sim {N}_{\mathrm{part}}^{4/3}$. At higher transverse momenta the spectra match the expectations from perturbative QCD calculations [365], while at lower momenta an exponential behavior is observed. By fitting the slope of the transverse momentum spectra of direct photons, one can infer an effective temperature of the quark-gluon plasma as depicted in figure 73.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 74 state-of-the-art calculations considering direct photon emission from all stages of the reaction are compared with spectra and elliptic flow data from the ALICE collaboration in Pb+Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}=2.76$ TeV. Again the prompt photons from the early hard collisions dominate the high p_T region. At lower transverse photon momenta the thermal emission from the hydrodynamic medium dominates. The photons emitted during the pre-equilibrium stage are depicted by the full line. The magnitude of photon elliptic flow it depends on the time at which full chemical equilibrium between quarks and gluons is achieved (indicated by the different values of τ_chem).

Zoom In Zoom Out Reset image size

In general, it is still hard to explain the photon production yield and elliptic flow at the same time consistently in one theoretical calculation. The intuitive reason is as follows: processes that increase the yield need to occur early in the evolution when the plasma is hotter, while larger elliptic flow is reached in the later stages of the evolution. Therefore, elliptic flow is increased, if later sources are enhanced, for example in [369] the contribution of the non-equilibrium hadronic stage has been shown to be significant for low transverse momenta. In the future, it is going to be crucial to investigate emission of electromagnetic probes from well-calibrated models for the bulk evolution [368].

There are two different mechanism that lead to the production of light nuclei in heavy ion reactions. In the fragmentation regions, at very forward and backward rapidity, the spectator remnants can fragment and reach the detector as a multitude of smaller and larger nuclei. This can happen because the spectator remnants emerge from the collision in a highly excited state, and collisions among nucleons within these spectator remnants often result in their disintegration. We will not discuss this phenomenon further as the main interest here is on the production of light (anti-)nuclei at midrapidity within the hot and dense collision region. The production of nuclei is interesting because of its dependence on the nucleon–nucleon interaction, but also as a probe of possible differences in the properties of matter and antimatter.

Figure 18 shows the production yields of many particle species including light nuclei in Pb+Pb collisions at $\sqrt{{s}_{\mathrm{NN}}}=2.76$ TeV [182]. Due to the vanishing chemical potential at this high energy the yields of deuterons and anti-deuterons, helium-3, hypertriton and helium-4 and their anti-nuclei are pairwise identical. Generally, the predictions within the statistical hadronization model agree very well with the measurements for the same temperature as for all other hadron species. This poses immediately a question of current debate: how can particles that have small binding energies of a few MeV freeze-out chemically from a fireball of hot and dense strongly-interacting matter at a temperature that is many times higher (T_ch ∼ 150 MeV)?

Besides the thermal production of light nuclei, another production mechanism is proposed. The idea is that only individual hadrons are produced from the fireball, but later on the nucleons or antinucleons combine to form light nuclei. This coalescence picture [371, 372] involves calculating the overlap in phase-space of all nucleons and drawing conclusions about the abundance of light nuclei from there. In such a picture it is expected that the transverse momentum spectra of light nuclei and their anisotropic flow coefficients v_n scale according to the number of nucleons contained in a nucleus, in analogy to the recombination approach and the partonic quark number scaling discussed above (see section 5.8). Figure 75 depicts the expected yields at midrapidity in a thermal and a coalescence production approach. The rapid fall-off of the yields with beam energy is caused by reduced baryon stopping at higher energies, which means that fewer valence quarks that can form complex nuclei are present in the fireball. The highest yields are therefore expected in the beam energy range E_lab = 10 − 20A GeV.

Zoom In Zoom Out Reset image size

Due to recent increased interest in the topic, several models have been developed that aim to explain the mechanism behind this behavior. One model attempts to describe the early chemical freeze-out via the Saha equation in analogy to cosmology [373] or rate equations [374]. Another approach involves microscopic calculations of the non-equilibrium dynamics of light nuclei in the hadronic stage of the reaction [375]. At low beam energies in a baryon-rich environment, the main reactions are the nucleon catalysis reactions, while at high beam energies the pion catalysis is more important. Calculations within a hadronic transport approach suggest that the chemical equilibrium is maintained due to the high reaction cross sections during the rescattering phase in nuclear collisions [376, 377]. The kinetic and chemical decoupling almost coincide in such an approach.

Hypernuclei and their properties are of particular interest since one might be able to infer knowledge about the Λ − N interaction. For a while, the hypertriton (${}_{{\rm{\Lambda }}}^{3}{\rm{H}}$ 'lifetime puzzle' attracted much interest, since the first measurements of the lifetime of the hypertriton seemed to deviate from the lifetime of a free Λ hyperon. However, newer measurements by ALICE (see figure 76), STAR, and HADES, and an even more precise result recently reported by ALICE [379], point to an excellent agreement with the world data for the Λ lifetime.

Zoom In Zoom Out Reset image size

A rather recent idea that has already been applied with considerable success, is to measure the femtoscopic correlations of exotic hadrons in elementary and heavy ion reactions. These final-state correlation measurements can be connected to the hadronic interactions of those particles. In this manner it is possible to extract quantitative information about the interaction of Ξ, Ω and other exotic hadrons containing one or more strange quarks with other hadrons (see the recent reviews [380, 381]).

The initial state of a non-central heavy ion collision is characterized by a very large angular momentum in the center-of-mass frame. For example, a collision between two ²⁰⁸Pb nuclei at $\sqrt{{s}_{\mathrm{NN}}}=5.02$ TeV with impact parameter b = 7.5 fm commands an angular momentum $L=A\sqrt{{s}_{\mathrm{NN}}}b/2\approx 2\times {10}^{7}{\hslash }$. Only a tiny fraction of this angular momentum ends up in the central rapidity region where most observed particles are produced. However, even this small fraction endows the QGP in non-central collisions with a sizable vorticity.

Assuming thermal equilibrium, the degree of polarization, expressed as the probability for a certain spin orientation $\hat{S}$, of Λ-hyperons within the final-state hadronic gas is [382]:

$$\begin{eqnarray}&&{P}_{{\rm{\Lambda }}}(\hat{S})\propto \exp \left(\displaystyle \frac{(\vec{\omega }+{\mu }_{{\rm{\Lambda }}}\vec{B})\cdot \hat{S}}{T}\right),\end{eqnarray} \tag{ 61 }$$

where $\vec{\omega }$ is the vorticity vector of the matter, μ_Λ = (−0.6138 ± 0.0047)μ_N is the Λ magnetic moment in nuclear magnetons μ_N , and $\vec{B}$ is the magnetic field present at emission. Since the magnetic moments of the Λ and $\overline{{\rm{\Lambda }}}$ differ by their sign, a magnetic field would cause them to be oppositely polarized, whereas vorticity of the medium results in identical polarizations. Λ-hyperons are ideal probes of polarization, because the direction of their parity violating weak decay Λ → p + π⁻ is strongly aligned with the direction of the spin of the hyperon. The angular distribution of the decay proton momentum $\vec{p}=p\hat{n}$ in the Λ rest frame is given by [382]

$$\begin{eqnarray}&&\frac{{\rm{d}}w}{{\rm{d\Omega }}}=\frac{1+{\alpha }_{{\rm{\Lambda }}}{\vec{P}}_{{\rm{\Lambda }}}\cdot \hat{n}}{4\pi }\end{eqnarray} \tag{ 62 }$$

with α_Λ = 0.732 ± 0.014 [383].

The global polarization of Λ and $\overline{{\rm{\Lambda }}}$ have been measured in Au+Au collisions at RHIC over a wide energy range [384, 385] and in Pb+Pb collisions at LHC [386]. The polarization is found to be along the direction of the angular momentum in the collision, perpendicular to the reaction plane, and to grow with decreasing collision energy, presumably because a larger fraction of the angular momentum carried by the incident nucleons ends up at midrapidity. A recent compilation of data is shown in figure 77. The measured average polarizations can be converted into an estimate of the vorticity. From (61) one finds

$$\begin{eqnarray}&&\langle \omega \rangle =(\langle {P}_{{\rm{\Lambda }}}\rangle +\langle {P}_{\overline{{\rm{\Lambda }}}}\rangle )T,\end{eqnarray} \tag{ 63 }$$

which gives 〈ω〉 ≈ 10²² s⁻¹ for $\langle {P}_{{\rm{\Lambda }}}\rangle \approx \langle {P}_{\overline{{\rm{\Lambda }}}}\rangle \approx 0.02$. STAR has also measured the global polarizations of Ξ and Ω hyperons at $\sqrt{{s}_{\mathrm{NN}}}=200\,\mathrm{GeV}$, which are found to be of similar magnitude as the Λ-polarization supporting the interpretation of the phenomenon as a universal effect of QGP vorticity [388].

Zoom In Zoom Out Reset image size

The average polarizations 〈P_Λ〉 and $\langle {P}_{\overline{{\rm{\Lambda }}}}\rangle$ agree with each other at all collision energies within the experimental errors. The most precise values have been measured in Au+Au at $\sqrt{{s}_{\mathrm{NN}}}=200\,\mathrm{GeV}$, where $\langle {P}_{{\rm{\Lambda }}}\rangle -\langle {P}_{\overline{{\rm{\Lambda }}}}\rangle =0.037\pm 0.07$ [389]. This measurement allows to set an upper limit on the magnetic field at the emission time: ∣B∣ < 10¹² T [390]. This is less than 10⁻³ of the maximal magnetic field generated during the collision of the two nuclei.

In addition to the global polarization of Λ and $\overline{{\rm{\Lambda }}}$, the experiments have also observed local polarization of hyperons along the beam direction. The orientation of the polarization vector depends on the direction of the hyperon transverse momentum and shows a quadrupole pattern with respect to the beam axis [391]. This effect is now understood as a result of the shear caused by the anisotropy of the transverse flow (see [392] for a review).

Vector meson alignment is another phenomenon related to spin that has been experimentally observed. Both STAR [393] and ALICE [394] have found evidence for a global alignment of the spins of K^*0 and ϕ mesons with respect to the collision plane. Alignment is defined as the deviation of the $m=m^{\prime} =0$ component of the spin–1 density matrix ${\rho }_{m^{\prime} m}=\langle m^{\prime} | \rho | m\rangle$ from its equilibrium value ρ₀₀ = 1/3 when all spin orientations are equally likely. Different from the concept of polarization, which distinguishes between spin orientation up and down with respect to the collision plane, alignment makes no such distinction. The mechanisms that can cause a nonzero aligment are thus less constrained by symmetry than those that can cause polarization. Indeed, the experiments find much larger values ∣ρ₀₀ − 1/3∣ ≈ 0.1 − 0.2 (compared with 10⁻² for global polarization), but there is currently no generally accepted explanation for the origin of this effect.

Particle spin can be added to the hydrodynamical model of quark-gluon plasma expansion by introducing spin degrees of freedom to the fluid [395]. The equations for such a spinning, viscous fluid can be derived from kinetic theory in the usual way by applying a reduction to a limited number of moments of the momentum space distribution [140, 396] or by symmetry-based analysis of the allowed perturbations of the energy-momentum tensor and the spin current around equilibrium [397]. In such an approach the observed collision energy dependence of the global Λ polarization can be explained with reasonable assumptions about the initial conditions (see figure 78).

Zoom In Zoom Out Reset image size

QCD gauge fields are characterized by a topological quantum number, called winding number. Field configurations with different winding number are separated by an energy barrier. In the vacuum, transitions between field configurations with different winding number are possible by tunneling processes, which can be described semiclassically. The gauge field fluctuations involved in the tunneling process are nonperturbative (i.e. they cannot be expanded in a power series of the gauge coupling g) and commonly known as instantons.

While field configurations with a definite winding number break CP invariance, experimental data tell us that the QCD vacuum realized in nature contains a superposition of such field configurations that conserves the global CP symmetry of QCD for yet unknown reasons. (An axion field, if it exists, would explain this mystery in a natural way through the so-called Peccei–Quinn mechanism [398]). In the presence of electromagnetic fields, however, these winding number fluctuations can generate local CP violations via the chiral anomaly [399] of the axial current

$$\begin{eqnarray}&&{\partial }_{\mu }{j}_{5}^{\mu }=-\displaystyle \frac{{N}_{c}}{2{\pi }^{2}}\displaystyle \sum _{f}{Q}_{f}^{2}\vec{E}\cdot \vec{B},\end{eqnarray} \tag{ 64 }$$

where Q_f are the electric charges of the light quark flavors. Anomalous hydrodynamics [400, 401] adds (64) to the conservation laws for energy-monentum, baryon number, and electric charge as a fourth macroscopic equation. As usual, the conservation laws must be supportlemented with constitutive equations, in this case for the vector and axial vector currents:

$$\begin{eqnarray}\begin{array}{rcl}{j}^{\mu } & = & {{nu}}^{\mu }+\sigma {E}^{\mu }+{\sigma }_{B}{B}^{\mu }\\ {j}_{5}^{\mu } & = & {n}_{5}{u}^{\mu }+{\xi }_{E}{E}^{\mu }+{\xi }_{B}{B}^{\mu },\end{array}\end{eqnarray} \tag{ 65 }$$

where σ, σ_B , ξ_E , ξ_B are transport coefficients (σ is the usual electric conductivity), n, n₅ denote the usual and axial quark densities, and E^μ = F^{μ ν} u_ν , B^μ = ε^{μ ν α β} u_ν F_{α β} are Lorentz covariant expressions for the electric and magnetic field. It is convenient to express these transport coefficients through vector and axial vector chemical potentials, μ and μ₅ [401].

In the absence of explicit parity violation, the global axial chemical potential is μ₅ = 0, which implies that the magnetic conductivity σ_B vanishes. However, the presence of instantons implies that μ₅ fluctuates locally, which means that the electric current receives locally fluctuating contributions from the magnetic field. Owing to the motion of the colliding nuclei heavy ion collisions generate very strong, short-lived magnetic fields of the order ${eB}\sim {m}_{\pi }^{2}$ that point perpendicular to the collision plane. One thus expects electric current fluctuations perpendicular to the collision plane, which result in event-by-event separation of the net electric charge of particles emitted into the upper and lower hemispheres. This is called the chiral magnetic effect (CME) [399].

The charge separation can be understood as a direct kinematic consequence of the alignment of quark spins along (or against) the magnetic field and the alignment of spin and momentum encoded in the chirality of the quark. If the magnetic field aligns a quark spin in the direction of the field, then a right-handed quark (positive chirality) will move in the direction of the field, and a left-handed quark (negative chirality) will move against the field. A positive axial density implies a preponderance of quarks with positive chirality and vice versa. Since the magnetic moment of a quark depends on the sign of its electric charge, this leads to a current in the direction of the magnetic field, if the axial density n₅ > 0, and against the field if n₅ < 0, as is illustrated in figure 79.

Zoom In Zoom Out Reset image size

Since the coefficients σ_B and ξ_B in (65) are proportional to the densities n₅ and n, respectively, the anomalous hydrodynamic equations sustain a low-energy propagating mode corresponding to propagating coupled density-axial density fluctuations similar to the sound mode [402]. This excitation, which propagates with a speed proportional to the magnetic field B, is known under the name 'chiral magnetic wave' (CMW). Like the current fluctuations induced by the CME, the CMW could be seeded by the topological charge density fluctuations in the glasma during the earliest stage of the collision.

Another phenomenon worth mentioning is the 'chiral vortical effect' (CVE). A complete analysis shows that the constitutive equations (65) also receive contributions proportional to the vorticity vector ω^μ of the QGP. Their effect is similar to the chiral magnetic effect with axial charge fluctuations driving a fluctuating electric current along the direction of the vorticity vector. It can be similarly understood as the CME as a result of the alignment of quark spins along the QGP vorticity vector in thermal equilibrium, which we discussed in the context of global hyperon polarization. Thus, the illustration in figure 79 applies, except that the vorticity ω replaces the magnetic field B. In an off-central heavy ion collision both vectors, $\vec{\omega }$ and $\vec{B}$, point in the same direction perpendicular to the reaction plane. An overview of these phenomena can be found in [403].

A number of observables are specifically sensitive to such fluctuating electric charge separation phenomena. A possible search strategy in the context of known background effects is discussed in [404]. Quantitative predictions based on solutions of the anomalous hydrodynamical equations with reasonable assumptions for the initial axial density fluctuations can be found in [405]. The magnitude of the expected event-by-event fluctuations also depends strongly on the longevity of the magnetic field. For the parameters used in [405] the experimental signals for the chiral magnetic effect are in the range of 10⁻⁴, but other assumptions may lead to much smaller predicted values [406].

As theoretical predictions of the magnitude of observables for the chiral magnetic effect are beset with large uncertainties, experimental searches for it are of paramount importance. Owing to the parity conserving nature of QCD, the anomalous electric current must fluctuate event by event, thus all signals involve the measurement of event-by-event fluctuations. This means that other 'normal' sources that are sensitive to the orientation of the reaction plane can contribute, in particular, those involving charged resonance decays modulated by the elliptic flow of the final-state hadron distribution (see [404] for an in-depth discussion and references). Experimental studies at RHIC [407] and LHC [408] have concluded that at most a small fraction (less than 10% for RHIC) of the observed signals can be attributed to the CME.

The comparison of measurements of observables sensitive to the chiral magnetic effect in p+Pb and Pb+Pb collisions is another way to assess the size of background effects. Any magnetic field-driven effect, such as the CME, should be greatly suppressed in p+Pb collisions in comparison with Pb+Pb collisions (by a factor $(1/{Z}_{\mathrm{Pb}}^{2}=1.5\times {10}^{-4}$). Data from the CMS experiment constrain possible contributions of the CME to the Pb+Pb data to less than 7% [409].

A more sensitive search for CME signals requires a suppression or cancellation of such background effects. This motivated a comparative study of two collision systems involving nuclear isobars, ⁹⁶Zr and ⁹⁶Ru, which was carried out at RHIC. The magnetic field produced in Ru+Ru collisions is larger than that produced in Zr+Zr collisions under otherwise identical conditions, because a ⁹⁶Zr nucleus contains 40 protons, while a ⁹⁶Ru nucleus contains 44 protons. As the CME observables are proportional to the square of the magnetic field, one expects roughly a 15% difference between the two systems for the CME contribution to any observable. Great care was taken to ensure that the experimental conditions for collisions in the two isobar systems were identical, the data were subjected to a sophisticated blind analysis protocol [410], the first of its kind in the field of relativistic heavy ion physics, and the analysis was performed independently by several groups.

The results published by the STAR collaboration [411] showed no evidence for the presence of a CME contribution to any of the predefined observables with an experimental precision of ±4 × 10⁻³. Figure 80 shows the measured ratios S(Ru)/S(Zr) for each of the signature observables S considered in the analysis. All signals are in some way related to the difference Δγ for same-sign and opposite-sign charged pairs of particles of the quantity γ defined in (8), normalized to the elliptic flow anisotropy v₂ that drives the background effects. A contribution from the CME would cause this ratio to be larger than unity. Clearly, all measured ratios lie well below one, which means that they do not provide evidence for a CME contribution.

Zoom In Zoom Out Reset image size

The fact that all ratios related to Δγ/v₂ cluster around a value of 0.97 suggests that they have a common source that can be traced to a difference between the two isobars, which is not related to the nuclear charge. Indeed, such differences are known to exist: ⁹⁶Zr has a thicker neutron skin than ⁹⁶Ru (the ⁹⁶Zr nucleus has four additional neutrons), and the two nuclei have different quadrupole deformations. The experiment clearly revealed these differences in the centrality dependence of the final-state multiplicity distributions. This means that the shape differences of the two isobar nuclei are the source of the largest systematic uncertainty as background effects associated with elliptic flow do not cancel exactly. The resulting estimate for the background ratio of Δγ/v₂ is indicated by the shaded bars in figure 80. The measured value for the ratio, R_Ru/Zr = 1.013 ± 0.003 ± 0.005, is in excellent agreement with this baseline [412] and does not show evidence for a CME contribution.