Nuclear Science and Techniques

NUCLEAR PHYSICS AND INTERDISCIPLINARY RESEARCH

Recent development of hydrodynamic modeling in heavy-ion collisions

Chun Shen ，

Li Yan

Vol.31, No.12

Article number 122

Published in print 01 Dec 2020

Available online 04 Dec 2020

DOI：10.1007/s41365-020-00829-z

9300

We present a concise review of the recent development of relativistic hydrodynamics and its applications to heavy-ion collisions. Theoretical progress on the extended formulation of hydrodynamics toward out-of-equilibrium systems is addressed, with emphasis on the so-called attractor solution. Moreover, recent phenomenological improvements in the hydrodynamic modeling of heavy-ion collisions with respect to the ongoing Beam Energy Scan program, the quantitative characterization of transport coefficients in three-dimensionally expanding quark–gluon plasma, the fluid description of small colliding systems, and certain other interdisciplinary connections are discussed.

Heavy-ion collisionshydrodynamics

1 Introduction

Smashing heavy nuclei at high energies in large particle accelerators routinely creates extreme conditions to study the properties of many-body systems whose interactions are governed by quantum chromodynamics (QCD). Within a few yoctoseconds (10^-24 s), the collision systems are compressed to 10³⁰ atm and reach several trillion degrees Kelvin. A novel state of matter with deconfined quarks and gluons is formed under such extreme conditions, which is called quark–gluon plasma (QGP).

The QGP created in laboratories is a relativistic dynamical system, which expands and evolves like a nearly perfect liquid [1]. The size of the liquid droplet depends on the size of the colliding nuclei, which may vary from O(10) fm in gold–gold collisions at the Relativistic Heavy-Ion Collider (RHIC) at the Brookhaven National Laboratory, or lead–lead collisions at the Large Hadron Collider (LHC) at CERN, to O(1) fm in small colliding systems such as the proton–lead or even proton–proton collisions carried out at these facilities. The fluidity of QGP is one of the main subjects that has been explored in heavy-ion collisions. From experiments, it has been analyzed extensively through various types of long-range multiparticle correlations of the observed hadrons, known as the signatures of collective flow [2-4]. Theoretical model calculations using relativistic viscous hydrodynamics successfully characterize these flow observables, which makes relativistic hydrodynamics the "standard model" in heavy-ion collisions [5-10].

Phenomenological analyses within hydrodynamic frameworks provide the most efficient, robust, and effective tool to extract many-body QCD. For instance, as the macroscopic emergence of the interactions among quarks and gluons, transport coefficients in the QGP medium can be inferred from comparisons between hydrodynamic models and experimental data. To date, the specific shear viscosity, that is, the ratio between the shear viscosity and entropy density, η/s, has been constrained to values very close to a lower theoretical bound ħ/4πkB, suggesting that QGP is a strongly coupled medium [11, 12]. The specific bulk viscosity has been extracted as well, leading to a temperature dependence [13]. In addition to the transport coefficients, the properties of QCD are also hidden in the equation of state (EoS). These include the relations among the local energy density, pressure, entropy density, and speed of sound c_s. Certain ongoing attempts through hydrodynamic modeling have achieved results compatible with the solutions of lattice QCD [14, 15].

These efforts to study transport coefficients and the EoS based on hydrodynamics are essential for a quantitative characterization of QCD matter. In particular, a reliable hydrodynamic description of the system evolution is crucial for the search for a conjectured QCD critical point and its associated first-order phase transition between QGP and hadron gas at a finite baryon density [16]. Searching for the QCD critical point is the focus of the current Beam Energy Scan program (BES) at RHIC. However, extensions of the hydrodynamic model to cases involving finite baryon densities are challenging, especially considering the significant hydrodynamic fluctuations of baryon density associated with the QCD critical point. The difficulty stems not only from improving the model itself, such as the baryon charge in the EoS, but also from the requirement for fundamental progress in the theoretical formulation of hydrodynamics, such that stochastic hydrodynamical fluctuations may be taken into account systematically [17-23], and novel hydrodynamic modes due to the effect of critical slowing down can be included [24-27]. When nonzero baryon density is involved, the questions of how the collective behavior of the QGP is changed and, correspondingly, how the observed correlations of these generated hadrons are modified, must be answered in the hydrodynamic modeling. The success of hydrodynamics and its applications to heavy-ion collisions also yields many surprises. Although in large systems created in high-energy nucleus–nucleus collisions, various measurables with respect to the collective flow have been found to be consistent with hydrodynamic modeling [28], the application of hydrodynamic modeling to small systems such as those created in proton–nucleus collisions [29-31] is not straightforward, owing to the significant reduction in system size and high expansion rate [32, 33]. The "unreasonable effectiveness" [10, 28, 34, 35] in describing the collectivity in the small colliding systems has qualitatively modified the understanding of QCD system thermalization [36]. The condition of the onset hydrodynamics (hydrodynamization) is quite relaxed. The traditionally recognized hydrodynamic and nonhydrodynamic modes and propagation of these modes [37, 38] have been generalized largely beyond local thermal equilibrium. Out-of-equilibrium hydrodynamics, a novel concept associated with the discovery of attractor solutions in various dynamical systems [39], has been proposed as a theoretical candidate to generalize the applicability of hydrodynamics. In recent years, significant progress has been made in developing the theoretical formulation of out-of-equilibrium hydrodynamics.

The successful phenomenological application of relativistic fluid dynamics in heavy-ion collisions and the continuous support and challenges from the RHIC and LHC experiments have led to a vibrant program that unites research from traditionally separate disciplines such as string theory, computational physics, statistics, nuclear physics, and high-energy physics. Recently, the direct detection of gravitational waves from black holes and neutron star mergers [40-42] has added another interconnection with relativistic heavy-ion collisions at large baryon densities.

This review focuses on these recent developments in out-of-equilibrium hydrodynamics and highlight some of the current state-of-the-art phenomenological applications of hydrodynamic frameworks to describe the dynamics of relativistic heavy-ion collisions.

In Section 2, we review the theoretical formulation of out-of-equilibrium hydrodynamics at an introductory level. This is presented first from the extension of second-order viscous hydrodynamics to systems with large local gradients quantified by the Knudsen number in the Bjorken flow. An attractive solution from such a dynamical system emerges naturally, as a consequence of the existence of fixed points in both the free-streaming and hydrodynamic regimes. The relation between the attractor and the asymptotic hydrodynamic gradient expansion is also addressed in the context of the trans-series solution and resurgence properties in the theory of asymptotic series. An alternative approach from kinetic theory is discussed in terms of a set of moments of the phase-space distribution function. These moments are coupled through their equations of motion. The lowest orders of the equation reduce to the familiar hydrodynamic equation of motion when the gradients of the system tend to vanish. Out-of-equilibrium effects can be accounted for by higher-order moments, whose contribution to the system evolution out of equilibrium results in an effective correction of the transport coefficients.

Section 3 covers the state-of-the-art applications of (3+1)D hydrodynamics and hadronic transport framework for heavy-ion collisions at intermediate and high collision energies. The experimental programs at $\sqrt{s} \sim O (10)$ GeV are extremely exciting to map out the phase structure of QCD matter at finite net baryon densities. In the meantime, the 3D hydrodynamic framework also opens a new dimension to study event-by-event fluctuations along the longitudinal direction. Because it is difficult to calculate the transport properties of QGP from first principles, quantitative characterizations of QGP have been driven by phenomenological analyses. We will summarize the collective effort in constraining the specific shear and bulk viscosity over the past decades and highlight recent efforts toward accessing the baryon diffusion constant in QGP. While hydrodynamics becomes the standard theory to describe large heavy-ion systems, smaller collisions in p+A and p+p systems challenge the conventional picture of the validity region of hydrodynamics. We have started to see a connection between the leading development of an out-of-equilibrium hydrodynamics formulation and strong collectivity in these small systems. Finally, we highlight some interdisciplinary connections between heavy-ion physics and nuclear structure physics as well as statistics and machine learning applications.

2 Out-of-equilibrium hydrodynamics

We start with a brief introduction to the fundamental concepts of viscous hydrodynamics that have been applied in the study of high-energy heavy-ion collisions. As an essential ingredient of hydrodynamic modeling, it should be emphasized that a truncation at the second order in gradients is generally considered in these viscous hydrodynamics formulations. This is to be distinguished from certain recent developments in out-of-equilibrium hydrodynamics, which often involve gradients to infinite orders.

2.1 Viscous hydrodynamics

Hydrodynamics is a low-energy effective theory that describes the evolution of long-wavelength modes in a dynamical system. These are slow modes, commonly known as hydrodynamic modes, which obey a set of hydrodynamic equations of motion stemming from conservation laws. The conservation of energy and momentum, for instance, ∂μTμν = 0, plays a key role in determining the space–time evolution of hydrodynamic fields: Local energy density ≤, pressure 𝒫, and fluid four-velocity Uμ [7].1 The general form of the energy–momentum tensor Tμν is given in the corresponding constitutive relation as

T^{μ ν} = ϵ U^{μ} U^{ν} - (P + Π) Δ^{μ ν} + π^{μ ν} .

(1)

The projection tensor is defined as

Δ^{μ ν} = g^{μ ν} - U^{μ} U^{ν},

(2)

such that the spatial gradient can be formulated in a covariant form ▽μ = Δ^μν∂_ν. For convenience, we also define the co-moving time derivative D=Uμ∂_μ, with which a four-vector can be decomposed into a temporal component and a spatial component with respect to the fluid four-velocity Uμ. In particular, the normal derivative can be separated as ∂_μ = ▽_μ + UμD.

In addition to ideal hydrodynamics corresponding to a fluid system in local thermal equilibrium,

T_{ideal}^{μ ν} = ϵ U^{μ} U^{ν} - P Δ^{μ ν},

(3)

there are viscous corrections in the energy–momentum tensor Tμν to capture deviations of the fluid system from the local thermal equilibrium. In the framework of viscous hydrodynamics, these corrections are formulated in terms of an expansion over spatial gradients of the hydrodynamical fields. More precisely, this expansion is characterized by the Knudsen number, Kn, which is essentially the dimensionless ratio of the microscopic to the macroscopic length scale. For the bulk pressure ∏ and the shear stress tensor πμν, one has up to the first order in gradient, the Navier–-Stokes hydrodynamics,

\begin{array}{l} π^{μ ν} = 2 η^{〈} \nabla^{μ} U^{ν 〉} + O (\nabla^{2}), \\ Π = - ζ \nabla \cdot U + O (\nabla^{2}), \end{array}

(4)

where η and ζ are the shear and bulk viscosities, respectively. These are the transport coefficients determined by interactions among fluid constituents, reflecting the dynamic nature of the underlying theories. For Eq. (4), the Knudsen number can be approximately read as [32],

Kn \sim \frac{| η^{〈} \nabla^{μ} U^{ν 〉} |}{P} o r \frac{| ζ \nabla \cdot U |}{P} .

(5)

In Eq. (4) and in what follows, the brackets around tensor indices indicate a symmetric, transverse, and traceless projection of a tensor, i.e.,

^{〈} A^{μ ν 〉} = Δ^{μ ν α β} A_{α β},

(6)

where

Δ^{μ ν α β} \equiv \frac{1}{2} [Δ^{μ α} Δ^{ν β} + Δ^{μ β} Δ^{ν α}] - \frac{1}{3} Δ^{μ ν} Δ^{α β} .

(7)

In the spirit of the hydrodynamic gradient expansion, Eq. (4) can be systematically extended to higher orders. In particular, considering the fact that the resulted equations of motion from Navier–Stokes hydrodynamics are acausal, the extension to elevating the dissipative currents to dynamical degrees of freedom is necessary.2 For practical simulations, acausal modes can be remedied by using the Israel–Stewart formulation [44], with second-order gradient terms included. These terms relax to the Navier–-Stokes form, with the relaxation effect specified by the shear and bulk relaxation times, τπ and τ∏. With respect to conformal symmetry, the second-order shear stress tensor is completely determined by BRSSS hydrodynamics [45],

\begin{array}{l} π^{μ ν} = η^{〈} \nabla^{μ} U^{ν 〉} - τ_{π} [^{〈} D π^{μ ν 〉} + \frac{4}{3} π^{μ ν} \nabla \cdot U] \\ - \frac{λ_{1}}{η^{2}} π_{α}^{〈 μ} π^{ν 〉 α} - \frac{λ_{2}}{η} π_{α}^{〈 μ} Ω^{ν 〉 α} - λ_{3} Ω_{α}^{〈 μ} Ω^{ν 〉 α}, \end{array}

(8)

where, in addition to the shear relaxation time τπ, λ₁, λ₂, and λ₃ are independent second-order transport coefficients. For conformal fluids, these transport coefficients are known [12, 45] and can be parameterized as [46, 47]

η = C_{η} s, τ_{π} = \frac{C_{τ} C_{η}}{T}, λ_{1} = C_{λ_{1}} \frac{s}{T},

(9)

where the local entropy density s ∝ T³. For consistency and considering a weakly coupled system, in the current review, we shall take the evaluations from kinetic theory for a conformal system, which are [12, 48],3

C_{η} = \frac{1}{4 π}, C_{τ} = 5, C_{λ_{1}} = \frac{5}{7} C_{η} C_{τ} .

(10)

There exist other variant forms of second-order viscous hydrodynamics in addition to Eq. (8), when conformal symmetry is not guaranteed [49]. Note that Eq. (8) is consistent with the Müller–Israel–Stewart theory [44], which relaxes to Eq. (4) when the relaxation time τπ → 0. Note also that there are more tensor structures arranged in the second-order terms, such as the vorticity tensor Ω^μν, Ω^μν = ▽^μU^ν - ▽^νU^μ, whereas for Navier–Stokes hydrodynamics, only one term is involved.

Extension to even higher orders has been considered in the literature (cf. [50, 51]), with more tensor structures introduced with correspondingly new transport coefficients. As a consequence of the increasing number of tensor structures, it is expected at the n-th order, the number of new transport coefficients, or the number of new tensor structures scales as n!. This factorial increase essentially affects the convergence properties of the hydrodynamic gradient expansion, such that it is rather asymptotic than convergent [52-54].4 In addition to the shear channel, the asymptotic property of the hydro gradient expansion exists in the bulk and diffusion channels as well. In principle, the applicability of the classical framework of hydrodynamics relies on the analysis of gradient expansion.

In a similar way, the charge conservation gives, ∂μJμ = 0, where the conserved current of hydrodynamics undergoes dissipative corrections as well,

J^{μ} \equiv J_{i d e a l}^{μ} + I^{μ} = n U^{μ} + I^{μ},

(12)

with

I^{μ} = σ T \nabla^{μ} (\frac{μ}{T}) + O (\nabla^{2}),

(13)

and σ is the corresponding conductivity of the conserved charge. For heavy-ion collisions, the net baryon number, which is related to the QCD critical behavior, is commonly considered in hydrodynamic analysis. As in the shear channel, acausal modes can be avoided by extending the constitutive relation in Eq. (13) to the Cattaneo equation with a finite relaxation time τQ [57],

I^{μ} = σ T \nabla^{μ} {(1 + τ_{Q} D)}^{- 1} (\frac{μ}{T}) .

(14)

For most theoretical analyses carried out with respect to QGP in high-energy nucleus–nucleus collisions, with up to second-order viscous corrections, the aforementioned equations provide the essential ingredient of a successful phenomenological model that captures the system evolution. Together with the EoS provided by lattice QCD simulations, e.g., 𝒫 = 𝒫(e), numerical solutions to the hydrodynamic modeling give rise to the space–time evolution of the hydrodynamic fields, which eventually reaches freeze-out and yields the observed particles in experiments. More details on phenomenological modeling will be provided later in Section 3.

2.2 Hydrodynamization and out-of-equilibrium fluid dynamics

Hydrodynamic modeling has been successfully applied to small colliding systems. For the high-multiplicity events of proton–lead [29, 58-60], ³He–gold, deuteron–gold [61-63], and even proton–proton collisions [30], the observed multiparticle correlations were found to be compatible with hydrodynamic predictions [64, 65]. Compared to nucleus–nucleus collisions, in these systems, the created QGP fireball is expected to be small and short-lived. Even though the application of hydrodynamics suffers from small system size and finite evolution time, which in turn gives large spatial gradients, the question of why hydrodynamics is "unreasonably" successful in small colliding systems is one of the recent focuses in the heavy-ion community, which has motivated extensive theoretical development of out-of-equilibrium hydrodynamics [9, 46, 66-81].

When applied to realistic simulations of heavy-ion collisions, the framework of viscous hydrodynamics assumes a valid truncation of the gradient expansion at second order in the gradients. The validity of the truncation requires that the QGP system is locally close to thermal equilibrium, such that in the equation of motion O(Kn³) ≪ O(Kn²). In this way, gradient corrections of order higher than or equal to O(Kn³) can be safely neglected. In heavy-ion collisions, this is by assumption satisfied at a time scale τ₀ when the QGP created from heavy-ion collisions approaches local thermal equilibrium (thermalization), or when the system evolution starts to be captured by second-order fluid dynamics (hydrodynamization). It should be emphasized that hydrodynamization is a more relaxed condition than thermalization, which does not require isotropization between the longitudinal and transverse pressures, and the finite pressure difference 𝒫L-𝒫T is accounted for by dissipative effects.

Compared to the lifetime τf of QGP in heavy-ion collisions, it is obvious that only when τ₀≪τf can the fluid dynamics description of the system evolution be reliable. For a strongly coupled QGP medium, a theoretical analysis based on AdS/CFT estimates that hydrodynamization scales inversely to temperature, τ₀∼1/T [82]. On the contrary, if the QGP system is weakly coupled, and the fundamental interactions are scattering events among quarks and gluons characterized by perturbative QCD, the onset of hydrodynamics is related to the strong coupling constant αs and saturation scale Qs such that $τ_{0} \sim α_{s}^{x} Q_{s}^{- 1}$ [83], where the exponent x is a constant negative number. However, given the actual values of the strong coupling constant and Qs in heavy-ion collisions, the expected τ₀ in realistic QGP systems could be rather large [83]. We concentrate on the weakly coupled system in the present discussion.

Despite all the estimates of τ₀ in various theories, how the systems created in heavy-ion collisions evolve toward fluids is an outstanding question. A schematic illustration of the early stages of system evolution in heavy-ion collisions, considering more realistic situations that are compatible with QCD, is shown in Fig. 1. The created medium is believed to experience at first the stage of the classical gluon field. This is a color glass condensate (CGC) picture [84], in which the system evolution is dominated by a saturated gluon field. The longitudinal and transverse pressures in the gluon field are highly anisotropic [85]. As the system expands and the density of the gluon field decreases, around τ∼1/Qs, a kinetic theory description becomes available for gluons. In the kinetic theory description, isotropization is eventually achieved via scattering among quarks and gluons, against the effect of longitudinal expansion. In Fig. 1, isotropization is characterized in terms of the ratio of longitudinal to transverse pressure (blue lines), 𝒫L/𝒫T. The blue dashed lines and the solid line correspond to various initializations due to quantum fluctuations around 1/Qs, which lead to different pressure anisotropies 𝒫L/𝒫T when the kinetic theory description starts. At late times, regardless of the arbitrary initial conditions, the evolution of pressure anisotropy becomes universal. At around τ₀, the second-order viscous hydrodynamics starts to dominate, and 𝒫L/𝒫T approaches unity, that is, close to local thermal equilibrium.

Fig. 1.

Schematic illustration of the early-time stages of the system evolution in high-energy heavy-ion collisions, describable respectively starting from initial time of heavy-ion collisions by classical gauge theory for the gluon field, then around τ∼1/Qs by kinetic theory for quarks and gluons, and by fluid dynamics for hydrodynamic variables such as pressure and energy density after hydrodynamization. The blue lines indicate the longitudinal to transverse pressure ratio, 𝒫L/𝒫T, characterizing how far the system is from the local thermal equilibrium, i.e., 𝒫L/𝒫T=1.

The existence of such a universal evolution as illustrated in Fig. 1 has been proved by various theoretical analyses. This uniquely implies the value of τ₀, irrespective of the initial conditions. Moreover, it provides a novel and extended description of the system evolution that applies to out-of-equilibrium system evolution. This is the fundamental idea of out-of-equilibrium hydrodynamics. This universal evolution, which is dubbed an "attractor," can be shown beyond the characterization of hydrodynamics, including infinite orders in the gradient expansion. To date, studies based on attractor solutions have been found to be theoretically feasible in certain highly symmetric expanding systems, such as the Bjorken flow [33, 39, 46, 47, 67, 68, 71, 74, 81, 86-89] and Gubser flow [69, 79, 80, 90, 91], and numerically in less symmetric systems [70]. It has been solved with respect to the equation of motion from hydrodynamics as well as kinetic theory. From both perspectives, we shall present a pedestrian introduction to the derivation of the attractor solutions.

2.2.1 From fluid dynamics

As an example, we first present the analysis with respect to conformal fluids and Bjorken symmetry. Bjorken symmetry is a good approximation for the system created in high-energy heavy-ion collisions in its very early stages [86]. It describes the longitudinal boost-invariant expansion of the medium along the beam axis (which we identify as z), while transverse expansions along the x and y axes are ignored. Boost-invariant symmetry is motivated by observations in high-energy collisions in the mid-rapidity region [86], and considering the fact that longitudinal expansion dominates at the stage shortly after collisions allows one to ignore transverse expansions.

Instead of using the usual Minkowski space–time coordinates t and z, with the new Milne coordinates,

τ = \sqrt{t^{2} - z^{2}}, ξ = \tanh^{- 1} (t / z),

(15)

the Bjorken flow is simplified such that hydrodynamic fields can be written independently of the space–time rapidity ξ. In the Milne coordinate system, the fluid four-velocity under Bjorken symmetry is fully determined as

U^{τ} = 1, U^{x} = U^{y} = U^{ξ} = 0 .

(16)

Therefore, considering that the metric is rewritten as $d s^{2} = g_{μ ν} d x^{μ} d x^{ν} = d τ^{2} - τ^{2} d ξ^{2} - d x_{⊥}^{2}$ , gradients of the fluid fields are reduced, and all are related to the proper time τ; for example,

\nabla \cdot U = - \frac{3}{4} σ_{ξ}^{ξ} = \frac{1}{τ} .

(17)

This relation implies a unique expression of the Knudsen number for the Bjorken symmetry, that is, Kn^-1∼τ/τπ.

The equation of motion of hydrodynamics, ∂μTμν=0, becomes

\partial_{τ} ϵ = - \frac{4}{3} \frac{ϵ}{τ} + \frac{π_{ξ}^{ξ}}{τ} .

(18)

Without dissipation, the above equation has an analytic solution for the energy density, ideal hydrodynamic evolution, ϵ∼τ^-4/3. It should be noted that the exponent -4/3 is the characteristic decay rate of the energy density in the ideal and conformal fluid experiencing Bjorken expansion. With dissipative corrections, the expected full solution of the energy density consists of gradient corrections,

ϵ \sim τ^{- 4 / 3} (1 + \sum_{n} O {(\frac{τ_{π}}{τ})}^{n}) .

(19)

One may check that in the conformal fluid with respect to Bjorke symmetry, the only non-zero components of the shear stress tensor are $π_{ξ}^{ξ} = - 2 π_{x}^{x} = - 2 π_{y}^{y}$ . Similarly, the vorticity tensor Ω^μν vanishes. By defining $π = π_{ξ}^{ξ}$ , and given all the information in the Bjorken flow, the constitutive relation for π follows from the BRSSS theory, Eq. (8), as

π = - \frac{4}{3} \frac{η}{τ} - τ_{π} [\frac{d π}{d τ} + \frac{4}{3} \frac{π}{τ}] + \frac{λ_{1}}{2 η^{2}} π^{2},

(20)

which is a nonlinear first-order differential equation. Without the constraint of conformal symmetry, the equation of motion is not unique. For instance, it can also be formulated according to the DNMR approach as [49]

π = - \frac{4}{3} \frac{η}{τ} - τ_{π} [\frac{d π}{d τ} + β_{π} \frac{π}{τ}] - \frac{χ τ_{π}^{2}}{η} \frac{π^{2}}{τ},

(21)

where βπ and χ are the transport coefficients in this alternate formulation. It should be noted that βπ and χ are related to the second-order transport coefficients appearing in BRSSS hydrodynamics (cf. Eq. (78)). It is worth mentioning that, for some evaluations of the transport coefficients, such as constant τπ, the coupled equations of motion in Eq. (20) and Eq. (21) can be solved analytically [88, 92].

To proceed, one needs to solve the coupled equations, Eqs. (18) and (20) and to construct the hydro gradient expansion accordingly. A convenient way to do so is to introduce

g (w) \equiv \frac{d ln ϵ}{d ln τ} = - 1 - \frac{P_{L}}{ϵ},

(22)

as a function of the inverse Knudsen number

w \equiv τ / τ_{π} .

(23)

Apparently, g(w) characterizes the decay rate of the energy density and is related to the pressure ratio. Although in the conformal case, dw/dτ>0 is satisfied, w is more significant than purely a time scale. A small w could indicate an early time of the system evolution, but it can also be interpreted as the system being far away from local equilibrium. On the contrary, a large w could indicate a late time, but it also corresponds to systems close to the local equilibrium.

After some algebra, rewriting the coupled equations in terms of g(w), one realizes a nonlinear first-order differential equation,

\begin{array}{l} w \frac{d g}{d w} (1 + \frac{g}{4}) + {(g + \frac{4}{3})}^{2} [1 + \frac{3 w}{8} \frac{C_{λ}}{C_{η}}] \\ + w (g + \frac{4}{3}) - \frac{16}{9} \frac{C_{η}}{C_{τ}} = 0, \end{array}

(24)

where the transport coefficients have been substituted by the parameterization constants.

With the initial conditions of the function g(w) given, we can solve Eq. (24). Numerical solutions are presented in Fig. 2, together with the results obtained from the viscous hydrodynamic solution with viscous corrections up to the first, second, and 50th order (solid blue, orange, and red lines, respectively). In comparison with the full solutions, it is apparent that hydrodynamics with truncated viscous corrections are valid only when w≫1, as anticipated. Moreover, including more viscous corrections does not help to improve the solution out of equilibrium. This is a direct consequence of the fact that hydro gradient expansion is not convergent, as we shall detail later.

Fig. 2.

(Color online) Numerical solution of g(w) to the conformal second-order fluid dynamics with respect to various initial conditions (dashed lines). The solid black line corresponds to the attractor solution, starting from the free-streaming stable fixed point g₊ of the evolving system when w→0⁺. Solutions from random initial conditions (dashed lines) converge toward the attractor solution when w∼1. The leading order slow-roll approximation of the attractor solution is indicated by the solid purple line Eq. (30). Results with corrections from first order, second order, and up to 50th order in viscous hydrodynamics are plotted as solid blue, orange, and red lines, respectively.

As shown in Fig. 2, full numerical solutions starting from various initial conditions merge to a single curve, indicated by the solid black line, around w≲1. This solid black line is known as the attractor solution, which in a dynamical system is often referred to as a solution collected in the phase space irrespective of initial conditions [81]. This attractor solution in hydrodynamics was first noticed in the context of the Bjorken flow [39, 46], and later in the Gubser flow (cf. Ref [92]). It is also realized in the solution of the kinetic equation for weakly coupled media and strongly coupled media using AdS/CFT [70, 72]. The attractor solution offers a valid and universal description of the system evolution, even if w≲1. Therefore, it extends the applicability of hydrodynamics to out-of-equilibrium systems.

Numerically, the attractor solution can be solved with respect to a special initial condition, corresponding to the free-streaming stable fixed point. Analytically, the emergence of an attractor solution in the nonlinear differential equation can be understood in terms of either the evolution of (pseudo) fixed points or the Borel resummation of the hydrodynamic gradients.

Fixed point analysis. Eq. (24) has singularities at w=0 and infinity, between which the solution is expected to be analytic. Hence, it is worth examining the two extremes: the far-from-equilibrium extreme with w→0⁺ and the close-to-local-equilibrium extreme with w→+∞. In particular, one should concentrate on the stable fixed points in these extremes that govern the system evolution.

Because w→0⁺ is equivalent to setting τπ→∞, which corresponds to infinitely weak interactions among fluid constituents, the far-from-equilibrium extreme can also be interpreted as system evolution determined entirely by expansion. This is known as free streaming. In the limit of small w, the nonlinear differential equation is reduced to an algebraic equation,

{(g + \frac{4}{3})}^{2} - \frac{16}{9} \frac{C_{η}}{C_{τ}} = 0,

(25)

with its two solutions given as (g_{_}-<g₊)

g_{\pm} = - \frac{4}{3} (1 \mp \sqrt{\frac{C_{η}}{C_{τ}}}) .

(26)

These are the two fixed points of free streaming, and g₊ corresponds to a stable point while g_- leads to an unstable fixed point. That is, if the system is initialized at g₊, it will stay at g₊ for a purely expanding system. In pure free-streaming systems, perturbations around the stable fixed point decay with time, and the decay rate scales as a power law [33, 74],

δ g (w) \sim {(\frac{τ}{τ_{i n i}})}^{- (g_{+} - g_{-})} .

(27)

This power law decay qualitatively explains the observed pattern of g(w) evolution at early times in Fig. 2.

In the opposite limit that w→+∞, the Knudsen number is small enough that the system approaches the hydrodynamic regime. Given Eq. (19) and the definition of g(w), it is not difficult to notice that in the limit of w→+∞,

g (w) \to g_{hyd} = - \frac{4}{3},

(28)

which is the energy density decay rate of ideal hydrodynamics. This is the hydrodynamic fixed point, to be reached as long as hydrodynamization is realized in a system.

In fact, a convenient way to reveal the properties of these fixed points and to understand how these fixed points emerge from Eq. (24) in both extremes, is to define effectively a beta function. By considering the root of dg/dw = 0, one effectively defines

\begin{array}{l} β_{BRSSS} (w, g) \equiv {(g + \frac{4}{3})}^{2} [1 + \frac{3 w}{8} \frac{C_{λ}}{C_{η}}] \\ + w (g + \frac{4}{3}) - \frac{16}{9} \frac{C_{η}}{C_{τ}} . \end{array}

(29)

The root of Eq. (29) encodes the information of fixed points in Eq. (26) and in Eq. (28). Figure 3 shows an illustration of the beta function for various values of w, from a small value (w=0.01) in the vicinity of free streaming, to a large value (w=10) in the hydrodynamic regime. The root of the beta function appears as the intersection of the solved line with the x-axis, and the stable fixed point corresponds to the crossing with a negative slope. Indeed, as w increases from w=0.01, where the crossing gives rise to a stable fixed point at approximately g₊ -1.165 (and an unstable fixed point at approximately g_- -1.502, at which dg/dw>0), the stable fixed point moves smoothly toward the hydro fixed point at g_hyd=-4/3. It should be noted that the hydro fixed point is highly stable because it is related to the crossing with an infinite negative slope. Accordingly, the unstable fixed point evolves from the free-streaming system to a super-unstable fixed point in the hydro regime, around $g = - \frac{4}{3} (1 + \frac{2 C_{η}}{C_{λ_{1}}}) = - 2.08$ , which again, is a root of β_BRSSS=0 in the limit of w→+∞.

Fig. 3.

Numerical evaluation of the beta function Eq. (29) for different values of w. Each intersection with the x-axis represents the (pseudo) fixed point at the corresponding w. The blue dot indicates the location of the hydrodynamic fixed point, corresponding to w→ ∞. A similar figure can be found in [95].

One may consider the accumulation of all these solved stable fixed points from the beta function to form an analytical representation of the attractor. This is shown in Fig. 1 as the solid purple line. Indeed, this approximation correctly captures the system evolution in both extremes, but deviates when w∼1, even though this deviation is relatively small. This approximation procedure coincides with the leading-order approximation using the slow-roll expansion [46, 71, 93, 94], in which one neglects the derivative in Eq. (24) for the lowest-order estimate: dg/dw∼0,

\begin{array}{l} g_{slow-roll} (w) = - \frac{4}{3} + \frac{1}{2 + \frac{3 w}{4} \frac{C_{λ}}{C_{η}}} [- w \\ + \sqrt{w^{2} + \frac{16}{9} \frac{C_{η}}{C_{τ}} (4 + \frac{3 w}{2} \frac{C_{λ}}{C_{η}})}] \end{array}

(30)

This is also compatible with the adiabatic evolution of a ground-state eigenmode (slowest mode), as determined by the linear system of the original coupled hydrodynamic equation of motion [76].

Hydro gradient expansion, trans-series, and resurgence. The hydrodynamic gradient expansion is a power series in terms of the Knudsen number, that is, 1/w,

g_{hydro (w)} \equiv \sum_{n = 0}^{\infty} \frac{f_{n}}{w^{n}} .

(31)

Substituting the expansion into Eq. (24), we find the recursion relation for the expansion coefficients5 ,

\begin{array}{l} - (n - \frac{8}{3}) f_{n} + \sum_{k = 0}^{n} (1 - \frac{k}{4}) f_{k} f_{n - k} . + (1 + \frac{C_{λ}}{C_{η}}) f_{n + 1} \\ + \frac{3}{8} \frac{C_{λ}}{C_{η}} \sum_{k = 0}^{n + 1} f_{k} f_{n + 1 - k} + (1 - \frac{C_{η}}{C_{τ}}) \frac{16}{9} δ_{n,0} \\ + \frac{2}{3} (2 + \frac{C_{λ}}{C_{η}}) δ_{n, - 1} = 0 . \end{array}

(32)

The leading-order solutions of the above equation are

f_{0} = - \frac{4}{3} and - \frac{4}{3} (1 + \frac{2 C_{η}}{C_{λ}}),

(33)

in agreement with those fixed points in the hydrodynamic regime found earlier from the root of the beta function, as anticipated. The super-unstable fixed point is nonphysical and is not expected in realistic solutions. Starting from the hydrodynamic fixed point f₀=-4/3, using Eq. (32) iteratively, one is able to obtain the expansion coefficients to arbitrary order. For instance, one finds $f_{1} = \frac{16}{9} \frac{C_{η}}{C_{τ}}$ . In fact, it can be shown that the magnitudes of these coefficients exhibit factorial growth, fn∼n!. This factorial growth can be recognized by noticing the ratio

\frac{f_{n + 1}}{f_{n}} \to S^{- 1} (n + β) + O (\frac{1}{n}),

(34)

for a large n in Eq. (32). The parameters S and β are real constants. With respect to Eq. (24), they are determined as

S = \frac{3}{2}, β^{- 1} = - \frac{2 C_{η}}{C_{τ}} .

(35)

Therefore, the leading-order contribution to the coefficients at very large n is fn∼Γ(n+β)/S^n+β∼n!. The factorial growth of the expansion coefficients results in the well-known statement that the hydrodynamic gradient expansion is rather asymptotic than convergent, which has a vanishing radius of convergence. Asymptotic series are commonly seen in physics, such as the perturbative expansion in quantum field theory [96-98] and the WKB approximation in quantum mechanics [99, 100].

One way to reveal hidden information in the gradient expansion, especially the emergence of nonhydrodynamic contributions from the hydrodynamic equation of motion, is to apply the Borel resummation technique. With respect to the hydrodynamic gradient expansion Eq. (31), the Borel transform defines a new series as

B [g_{hydro} (z)] \equiv \sum_{n = 0}^{\infty} \frac{f_{n}}{n!} z^{n} .

(36)

With an extra factor of 1/n! introduced, this series now has a finite radius of convergence. It can be shown that this new convergent series is related to the original hydrodynamic gradient expansion via an inverse Laplace transform, such that the Borel resummation of the hydrodynamic gradient expansion is obtained as

\begin{array}{l} {\tilde{g}}_{hydro} (w) = \sum_{n = 0}^{\infty} \frac{f_{n}}{w^{n}} \times \frac{1}{n!} \int_{0}^{\infty} d z e^{- z} z^{n} \\ = \int_{0}^{\infty} d z e^{- z} \sum_{n = 0}^{\infty} \frac{f_{n}}{n!} {(\frac{z}{w})}^{n} \\ = w \int_{0}^{\infty} d z e^{- z w} B [g_{hydro} (z)] . \end{array}

(37)

One may check that ${\tilde{g}}_{hydro} (w)$ is a solution to the nonlinear differential equation Eq. (24). Without singularities, the integration can be performed in a straightforward manner. For an asymptotic series, such as the hydrodynamic gradient expansion, we are considering, however, singularities in the Borel transform are expected on the real axis, which gives rise to additional contributions.

In Fig. 4, the singularity structure of the Borel transform is shown on the complex Borel plane, where a series of poles on the real axis are observed. This structure is estimated numerically by a symmetric Padé approximation of the Borel transform up to truncation order n=300. Note that the leading pole (the pole closest to the origin) lies at z=S, which is not sensitive to truncation orders. Ideally, the Borel transform of the original asymptotic series should result in a branch cut on the real axis starting from z₀, as indicated by the accumulated poles in Fig. 4. To avoid the branch cut on the real axis, the integral contour connecting zero and infinity in Eq. (37) must be analytically continued to the complex plane. Upon the integration contour considered above or below the real axis, there exists a complex ambiguity. This ambiguity results in the hydrodynamic gradient expansion in a complex term. With respect to the leading pole, $B [g_{h y d r o} (z)] \sim {(z - S)}^{- β}$ , with S and β real constants given in Eq. (35), the complex ambiguity leads to

Fig. 4.

Singularity structure of the Borel transform Eq. (36) from the Padé approximant. All poles are on the real axis, implying a branch cut from the one closest to the origin located at 3/2 to +∞. A similar figure can be found in [47].

I m [{\tilde{g}}_{hydro} (w)] \sim \pm π e^{- S / w} w^{β} .

(38)

Because the solution g(w) is real and definite, this complex ambiguity in g_hydro must be cancelled by a term with the same exponentially suppressed factor. That is, the singularity of the Borel resummation and the reality condition implies the existence of an extra contribution to the hydro gradient expansion, which is exponentially suppressed.

The existence of such extra terms with exponential factors can also be proved through a small perturbation around the hydro gradient expansion. In Eq. (24), assuming g(w) = g_hydro(w) + δg(w), we find an equation for δg(w) to be solved asymptotically by

δ g (w) \sim e^{- S w} w^{β} .

(39)

In fact, the complex ambiguity arises not only from the leading singularity, and the cancellation of all the complex ambiguities requires the complete form of the solution to be a trans-series, rather than a simple power series expansion,

g (w) = \sum_{n = 0} {(σ e^{- S w} w^{β})}^{n} g^{(n)} (w),

(40)

where the leading order gives g⁽⁰⁾(w)=g_hydro(w). In the trans-series, σ is a complex parameter denoting the order of trans-series expansion, whose real part is related to the initial condition. At each order, the function $g^{(n)} (w) = \sum_{k} f_{k}^{(n)} / w^{k}$ itself is an asymptotic series, which with respect to the Borel regression gives rise to a complex ambiguity to be cancelled by a factor similar to the next-order term in the trans-series. The imaginary part of σ is determined as a consequence of the cancellation. This is a typical property of the resurgent theory [101, 102], where different orders in the trans-series are related via the cancellation of complex ambiguities. Mathematically, it is not surprising to expect the trans-series solution with an exponentially suppressed factor for a nonlinear first-order differential equation.

For an asymptotically large w, the higher-order transient contributions in the trans-series are suppressed. However, in the small-w regime, that is, the nonperturbative regime of the hydrodynamic gradient expansion, these higher-order contributions become important. This can be shown, for instance, through the reconstruction of the attractor solution from the Borel resummation of the trans-series solution. A key step in the procedure is the identification of the real part of the σ parameter, corresponding to the initial condition that determines the attractor solution, which, as we discussed before, corresponds to g_hydro(0⁺) → g₊. Given Reσ, one must resum the trans-series order by order, according to the detailed cancellation rules provided by the resurgence relations. For a conformal fluid, it has been shown numerically that an attractor does emerge, provided higher orders in the trans-series are taken into account [46]. For certain fluid systems with analytical solutions, the procedure can be proved explicitly by the Borel resummation of the trans-series to infinite order [103].

2.2.2 From kinetic theory

The discussion in the previous section relies on an equation of motion provided by hydrodynamics, where viscous corrections are introduced up to the second order. Although a series expansion can be generated from the equation to arbitrary order, based on the equation with second-order viscous (gradient) corrections, this series expansion does not consistently and completely capture the information of the off-equilibrium physics. For instance, it would not be surprising to realize that in Eq. (31), the expansion coefficients with n≥3 are modified once in the original hydrodynamic equation of motion; third-order and higher-order viscous corrections are taken into account [95]. To formulate a gradient series that is compatible with the off-equilibrium system evolution, one has to solve the full transport equation [66, 75, 77, 104].

With respect to Bjorken symmetry, the general form of the Boltzmann equation [105],

p^{μ} \partial_{μ} f + Γ_{μ ν}^{α} p^{μ} p^{ν} \frac{\partial}{\partial p^{α}} f = \hat{C} [f],

(41)

is simplified. In particular, in the Milne coordinates, in the slide of z=0 (or ξ=0), the left-hand side of the kinetic equation reduces to

p^{μ} \partial_{μ} f + Γ_{μ ν}^{α} p^{μ} p^{ν} \frac{\partial}{\partial p^{α}} f \to p^{0} (\frac{\partial}{\partial τ} - \frac{p_{z}}{τ}) f,

(42)

where the phase-space distribution becomes a function of (τ, p). As expected, this corresponds to the same kinematic domain as in the fluid dynamics discussed in the previous section. We now consider a relaxation time approximation for the collision kernel to further simplify the kinetic equation,

(\frac{\partial}{\partial τ} - \frac{p_{z}}{τ}) f (τ, p) = - \frac{f (τ, p) - f_{eq} (τ, p)}{τ_{R}} .

(43)

There exists a formal and analytical solution to Eq. (43). For a relaxation time with arbitrary τ dependence, τR(τ), the formal solution is [106, 107],

f (τ, p) = D (τ, τ_{0}) f (τ_{0}, p_{⊥}, p_{z} τ / τ_{0}) + \int_{τ_{0}}^{τ} \frac{d τ^{'}}{τ_{R} (τ^{'})} D (τ, τ^{'}) f_{e q} (\sqrt{p_{⊥}^{2} + {(p_{z} τ^{'} / τ_{0})}^{2}} / T (τ^{'})) .

(44)

In this solution, f_eq is the equilibrium distribution as a function of temperature, which is fixed via the Landau matching condition, ϵ=ϵ_eq∝ T⁴. Apparently, as a consequence of conformal symmetry, f_eq does not depend on the chemical potential. The function in the first term is the free-streaming solution (when the collision kernel vanishes) of the kinetic equation f_FS(τ, p) = f(τ₀, p_⊥, pzτ/τ₀), and the time evolution function

D (τ_{2}, τ_{1}) \equiv \exp [- \int_{τ_{1}}^{τ_{2}} d τ^{'} \frac{1}{τ_{R} (τ^{'})}] .

(45)

The conservation of energy and momentum is implied in the kinetic equation, ∂μTμν = 0, where the energy–momentum tensor is defined in kinetic theory as

T^{μ ν} = \int \frac{d^{3} p}{E_{p}} p^{μ} p^{ν} f_{p} .

(46)

In the case of Bjorken symmetry, the independent components of the energy–momentum tensor are the diagonal components, which include the local energy density,

ϵ = T^{00} \equiv \int d^{3} p E_{p} f (τ, p),

(47)

the longitudinal pressure

P_{L} = T^{z z} \equiv \int \frac{d^{3} p}{E_{p}} p_{z}^{2} f (τ, p),

(48)

and the transverse pressure,

P_{T} = T^{x x} = T^{y y} \equiv \frac{1}{2} \int \frac{d^{3} p}{E_{p}} p_{⊥}^{2} f (τ, p) .

(49)

It should be noted that with respect to the conformal symmetry, є=𝒫L+2𝒫T. In terms of these components, the conservation of energy and momentum gives

τ \frac{d ϵ}{d τ} + \frac{4}{3} ϵ + \frac{2}{3} (P_{L} - P_{T}) = 0 .

(50)

From the above equation, one notices that the ideal hydrodynamic equation of motion is recovered when the pressures are isotropized, 𝒫L=𝒫T. In the case of viscous hydrodynamics, the pressure anisotropy corresponds to small viscous corrections. With respect to the BRSSS form of viscous hydrodynamics, [89]

P_{L} - P_{T} = - \frac{2 η}{τ} + \frac{4}{3 τ^{3}} (λ_{1} - η τ_{π}) + O (\frac{1}{τ^{3}}) .

(51)

The L-moments The energy–momentum tensor Tμν belongs to a specific set of moments of the phase-space distribution. We define the $L$ moment as

L_{n} = \int d^{3} p E_{p} P_{2 n} (p_{z} / E_{p}) f (τ, p),

(52)

using the Legendre polynomial Pn(x). As a result of the parity symmetry in the Bjorken expansion, moments associated with odd orders of the Legendre polynomials vanish. In addition to the Legendre polynomial that specifies asymmetry in the phase space, the weight Ep is chosen such that the $L$ moments are of the same dimension as the energy–momentum tensor. Indeed, it is straightforward to verify that

\begin{array}{l} L_{0} = T^{00} = ϵ = P_{L} + 2 P_{T}, \\ L_{1} = T^{⊥ ⊥} - T^{z z} = P_{L} - P_{T}, \end{array}

(53)

and hence, the pressure anisotropy can be expressed in terms of the $L$ -moments as

\frac{P_{L}}{P_{T}} = \frac{L_{0} + 2 L_{1}}{L_{0} - L_{1}} .

(54)

Legendre polynomials provide a complete set of decomposition in the angular dependence with respect to Bjorken symmetry, but the reconstruction of the phase-space distribution also requires a complete mode decomposition for the Ep dependence. For instance, the Ep dependence in f(τ, p) can be decomposed by large-error polynomials [108]. Although the generalized moments of the distribution function can be introduced (cf. Ref [49, 75, 109]), the $L$ -moments are sufficient for the description of system hydrodynamization, especially concerning the evolution of pressures and energy density. For instance, with respect to conformal symmetry, $L_{0}$ and $L_{1}$ fully determine the components in the energy–momentum tensor Tμν. With higher-order $L$ -moments included, the description of system evolution can be even improved. An illustration is presented in Fig. 5.

Fig. 5.

(Color online) Pressure anisotropy solved from the truncated moment equations with a truncation including two, three, four, and five

L

-moments, in comparison with the exact solution of the kinetic equation (open symbols). The solution with respect to viscous hydrodynamics, or the truncated moment equation at

L_{0}

and

L_{1}

, but using a renormalized η/s is shown as the gray dashed line.

It is also interesting to note that the coupled equations for $L_{0}$ and $L_{1}$ are analogous to the so-called anisotropic hydrodynamics (ahydro) [110, 111], where the pressure difference is considered as an individual field for out-of-equilibrium fluids. In a similar manner, higher-order $L$ -moments play the role of viscous corrections, in comparison to viscous anisotropic hydrodynamics (vahydro) [112].

With respect to the analytical formal solution of the distribution function, the analytical solution of $L$ -moments can be obtained:

L_{n} (τ) = D (τ, τ_{0}) L_{n}^{(0)} (τ) + \int_{τ_{0}}^{τ} \frac{d τ^{'}}{τ_{R} (τ^{'})} D (τ, τ^{'}) L_{0} (τ^{'}) (τ^{'} / τ) F_{n} (τ^{'})

(55)

where the function $F_{n} (x)$ is defined as

\begin{array}{l} F_{n} (x) \equiv \frac{1}{2} \int_{- 1}^{1} d y {[1 - (1 - x^{2}) y^{2}]}^{1 / 2} \\ \times P_{2 n} (\frac{x y}{{[1 - (1 - x^{2}) y^{2}]}^{1 / 2}}) . \end{array}

(56)

Note that in the limit x→0, $F_{n} (0) \to π P_{2 n} (0) / 2$ , and in the limit x→1, $L_{n \neq 0} (1) = 0$ . For the case of n=0, the integral above can be analytically evaluated, resulting in

F_{0} (x) = \frac{1}{2} (x - \frac{i \cosh^{- 1} (x)}{\sqrt{1 - x^{2}}}) .

(57)

The first term in Eq. (55) contains the free-streaming moment $L_{n}^{(0)} (τ)$ , which is given by

L_{n}^{(0)} (τ) = ϵ_{0} \frac{τ_{0}}{τ} F_{n} (\frac{τ_{0}}{τ}),

(58)

where ϵ₀ is the initial energy density. Eq. (55) should be regarded as an integral equation, for which $L_{0}$ must be solved with respect to the initial condition ϵ₀. Once $L_{0} (τ)$ is provided, higher-order $L_{n}$ can be determined accordingly.

Eq. (55) allows one to study the solution of energy density in powers of 1/w=τR/τ, i.e., the gradient expansion at late time. For instance, gradient expansion can be generated using integration by parts in the integral equation of $L_{0}$ , which again leads to an asymptotic series [66]. With the asymptotic series given, by analogy with the hydrodynamic gradient expansion, the emergence of a trans-series solution, and hence the resurgence phenomenon, etc., can be obtained following the standard procedure of Borel resummation.

Alternatively, one may start from a set of coupled equations of $L$ -moments. By substituting the definition of $L$ -moments into the kinetic equation Eq. (43), one obtains a first-order differential equation in which adjacent $L$ -moments are coupled,

\partial_{τ} L_{n} + \frac{a_{n} L_{n} + b_{n} L_{n - 1} + c_{n} L_{n + 1}}{τ} = - \frac{L_{n}}{τ_{R}} (1 - δ_{n 0}),

(59)

where n=0, 1, 2, …. The constant coefficients an, bn, and cn arise from the recursion relation of the Legendre polynomials,

\begin{array}{l} a_{n} = \frac{2 (14 n^{2} + 7 n - 2)}{(4 n - 1) (4 n + 3)}, b_{n} = \frac{(2 n - 1) 2 n (2 n + 2)}{(4 n - 1) (4 n + 1)}, \\ c_{n} = \frac{(1 - 2 n) (2 n + 1) (2 n + 2)}{(4 n + 1) (4 n + 3)}, \end{array}

(60)

reflecting the geometric nature of the Bjorken expansion. These can also be understood as the limiting case of the Clebsch–-Gordan coefficients without coupling between transverse and odd-parity modes. The first several constants are

\begin{array}{l} a_{0} = 4 / 3, b_{1} = 0, c_{0} = 2 / 3, \\ a_{1} = 38 / 21, b_{2} = 8 / 15, c_{1} = - 12 / 35. \end{array}

(61)

The equation for $L_{0}$ is the conservation of the energy moment, as in Eq. (50), whereas higher-order $L_{n}$ ’s bring in corrections. In the vicinity of the local equilibrium, these are viscous corrections.

To solve the time evolution of these moments, in comparison with the exact solution from the kinetic equation, one would expect truncating the coupled equations at a finite order as a good approximation. In Fig. 5, the numerically solved pressure anisotropy 𝒫L/𝒫T is plotted as a function of τ/τR, for the case of a conformal medium in which τ_RT is a constant with respect to a specified initial condition, 𝒫L/𝒫T≈0.49. The exact solution to the kinetic equation is shown as open symbols, compared to which the simplest truncation of the moment equations, at order n<2 involving moments $L_{0}$ and $L_{1}$ (dotted line), already captures the bulk property of time evolution. Deviations from the lowest-order approximation are significant only when τ/τ_R<1 and negligible in the hydrodynamic regime when τ/τ_R≫1. With higher-order $L$ -moments involved, improvements from these corrections are observed. With the truncation at n<5 (with $L_{0}$ ,…, $L_{4}$ taken into consideration), a satisfactory description of the pressure anisotropy over the entire time evolution is obtained.

In fact, the effectiveness of the truncation of the moment equations is guaranteed by the existence of fixed points and can be studied analytically in the limiting cases. If one considers free streaming, i.e., focuses on the very early time limit τ/τR→0 in which collisions are effectively suppressed by expansion, the solutions of moments are analytical. This is obvious in the formal solution, Eq. (55). One may also recast the set of moment equations into a matrix form with respect to the following vector:

ψ = (L_{0}, L_{1}, L_{2}, \dots) .

(62)

Correspondingly, the dynamics of free streaming are captured by a constant tri-diagonal matrix:

\hat{H} = (\begin{matrix} a_{0} & c_{0} & 0 & 0 & \dots \\ b_{1} & a_{1} & c_{1} & 0 & \dots \\ 0 & b_{2} & a_{2} & c_{2} & \dots \\ \dots & \dots & \dots & \dots & \dots \end{matrix}),

(63)

such that the equation of free streaming becomes

\partial_{ρ} ψ + \hat{H} ψ = 0 .

(64)

For convenience, we introduce ρ = ln (τ/τ₀). The solution of moments can be found, provided the eigenvalues of the matrix are determined. We note that the eigenvalues of $\hat{H}$ are complex, with eigenvectors satisfying

\hat{H} ϕ_{n} = λ_{n} ϕ_{n} \to ϕ_{n} (τ) = ϕ_{n} (τ_{0}) e^{- λ_{n} ρ} .

(65)

Let us order these eigenvalues by their real part, i.e.,

Re λ_{0} < Re λ_{1} < Re λ_{2} < \dots < Re λ_{\infty} .

(66)

Then, the solution of moments is

\begin{array}{l} ψ (τ) = \sum_{n} κ_{n} e^{- λ_{n} ρ} ϕ_{n} (τ_{0}) \\ = e^{- λ_{0} ρ} [κ_{0} ϕ_{0} (τ_{0}) + \sum_{n \neq 0} κ_{n} e^{- (λ_{n} - λ_{0}) ρ} ϕ_{n} (τ_{0})] \\ \to e^{- λ_{0} ρ} κ_{0} ϕ_{0}, \end{array}

(67)

where κn are constant coefficients fixed by the initial condition. Finally, the evolution of the moment is dominated by the ground-state mode on a timescale (λn-λ₀)ρ∼λρ ≫ 1, irrespective of the initial conditions [66]. The gap between the eigenvalues is of the order of unity, Δλ∼1.

If we generalize the definition of the function g(w) in Eq. (22), for all the moments,

g_{n} (w) \equiv \frac{d ln L_{n}}{d ln τ},

(68)

the dominance of the ground state indicates the existence of a stable fixed point, gn(w) → -λ₀, regardless of the order n. Similarly, an unstable fixed point is also expected, corresponding to λ_∞. This is very similar to the observations from the fluid dynamics analysis, although both the stable free-streaming fixed point and the unstable fixed point depend weakly on the truncation order. For truncation at n<2, the stable fixed point is g₊ ≈ -0.929, whereas the unstable fixed point is g_- ≈ -2.213, which is analogous to the fixed-point analysis from fluid dynamics. This is not accidental, but rather a direct consequence of the fact that the simplest truncation of moment equations leads to second-order viscous hydrodynamics. In Fig. 6, the ground-state eigenvalue is plotted as a function of the truncation order. Although asymptotically, when n→∞, λ₀→-1, truncating at finite orders only gives a small correction. This observation guarantees the validity of moment truncation in the free-streaming limit. Corresponding to the stable fixed point, or the ground-state eigenvalue, the ground-state eigenvector is determined to fix the ratios between moments,

Fig. 6.

Ground-state eigenvalue of the matrix

\hat{H}

, with respect to truncation orders. The asymptotic value is 1, while λ₀≈-0.929 when truncating at n<2. Figure taken from [95].

ϕ_{0} : L_{n} = P_{2 n} (0) L_{0} = (- {1)}^{n} \frac{(2 n - 1)!!}{(2 n)!!} L_{0} .

(69)

In the original phase-space distribution, these $L$ -moments with the specified ratios characterize a distribution spanning along p_⊥ and shrinking in pz. It may also be proven that when the truncation order goes to infinity, the unstable fixed point is g_-→-2, by noting that an+bn+cn = 2. The corresponding eigenvector is given by

ϕ_{\infty} : L_{0} = L_{1} = L_{1} = \dots

(70)

In the opposite limit with τ/τ_R→∞, the moments admit a series expansion in powers of 1/τ,

L_{n} = \sum_{m = n}^{\infty} \frac{α_{n}^{(m)}}{τ^{n}} .

(71)

This structure follows from the Chapman–Enskog expansion of the kinetic theory [89]. Except for $L_{0}$ , $α_{n}^{(0)} = ϵ δ_{n 0}$ , the expansion coefficients are directly related to the transport coefficients in viscous hydrodynamics. For instance, by noticing that the leading term in $L_{1}$ satisfies $L_{1} = P_{L} - P_{T} \propto η / τ$ , one is able to identify the shear viscosity as

η = \frac{b_{1}}{2} τ_{R} ϵ .

(72)

The coefficients in the expansion are dimensional, which can be further expressed in terms of dimensional variables in the original moment equations, $α_{n}^{(m)} = B_{n}^{(m)} ϵ τ_{R}^{n}$ , with $B_{n}^{(m)}$ dimensionless. The leading-order and next-leading-order expansion coefficients can be solved analytically. The leading-order coefficients, $α_{n}^{(n)}$ , in particular, determine the asymptotic value of gn(w). For instance, for a constant relaxation time, taking into account the time dependence of the energy density in an ideal fluid, ϵ∼τ^-4/3, one finds $L_{n} \sim τ^{- 4 / 3 - n}$ . For a conformal system in which τT is a constant, on the contrary, one finds $L_{n} \sim τ^{- 4 / 3 - 2 n / 3}$ . Therefore,

g_{n} (\infty) = {\begin{matrix} - \frac{4}{3} - n, τ_{R} = const . \\ - \frac{4}{3} - \frac{2 n}{3}, τ_{R} T = const . \end{matrix}

(73)

These asymptotics represent the hydrodynamic fixed points of the moments of different orders, which the moments would eventually approach, irrespective of the initial conditions. Correspondingly, attractors are smooth solutions that connect from the free-streaming fixed point and the hydrodynamic fixed point. Because the hydrodynamic fixed points differ in order for different $L_{n}$ , there are infinitely many attractors from the coupled moment equations, or the original kinetic equations. This is also observed in other forms of moments [75]. In Fig. 7, the attractors of the first five orders of $L$ -moments are plotted in terms of gn(w) for a conformal system.

Fig. 7.

(Color online) Attractor solution of gn for n=0,1,…,4, calculated from the numerical solution of the coupled moment equations with a truncation order of n=20. A conformal system wherein τRT is a constant is considered, for which the asymptotic values of gn(w) in the hydrodynamic regime depend on the order n, according to Eq. (73). Figure taken from Ref. [95].

The L-moments and variants of hydrodynamics. As we have discussed, the conservation of energy and momentum is only a subset of the coupled moment equations and is the leading one. Note that conservation of energy and momentum ∂μTμν=0 involves the first two moments $L_{0}$ and $L_{1}$ . The simplest truncation that satisfies the conservation of energy and momentum is

\partial_{τ} L_{0} + \frac{1}{τ} (a_{0} L_{0} + c_{0} L_{1}) = 0

(74a)

\partial_{τ} L_{1} + \frac{1}{τ} (a_{1} L_{1} + b_{1} L_{1}) = - \frac{1}{τ_{R}} L_{1},

(74b)

Together with the traceless condition $T_{μ}^{μ} = 0$ , all components in Tμν can thus be determined.

In the hydrodynamic regime, with τ/τR→∞, one would expect Eqs. (74) to be identified as the hydrodynamic equations of motion. In fact, Eq. (74a) is simply the conservation of energy and momentum, ∂μTμν=0, for a Bjorken expanding system. Eq. (74b), however, generalizes the constitutive relation that interprets the pressure anisotropy $L_{1} = P_{L} - P_{T}$ in terms of viscous corrections. For the system with Bjorken expansion, at a later time if one further identifies $c_{0} L_{1}$ as the ξξ component of the shear stress tensor, $c_{0} L_{1} = π = π_{ξ}^{ξ}$ , the truncated moment equations indeed lead to the Israel–Stewart hydro equations of motion:

\partial ϵ + \frac{4}{3} \frac{ϵ}{τ} = - \frac{π}{τ}

(75a)

π = - \frac{4}{3} \frac{η}{τ} - τ_{π} \partial_{τ} π - β_{π π} τ_{π} \frac{π}{τ} .

(75b)

In obtaining the above equation, we consider conformal EoS є=3𝒫 and use the shear viscosity from Eq. (72). By doing so, we find the second-order transport coefficient τπ=τR=5η/sT. The transport coefficient βππ is precisely a₁, which agrees with [87]. In Bjorken flow, ambiguity arises from interpreting 1/τ as the ξξ component of the shear tensor $σ_{ξ}^{ξ}$ or the expansion rate ▽·U, as indicated in Eq. (17); therefore, the hydrodynamic constitutive relation from Eq. (74b) is not unique. For instance, if one splits the last term in Eq. (75b) as

β_{π π} τ_{π} \frac{π}{τ} = τ_{π} \frac{4}{3} \frac{π}{τ} - \frac{3}{4} (β_{π π} - \frac{4}{3}) \frac{τ_{π}}{η} π^{2} + O {(\frac{1}{τ})}^{3},

(76)

the constitutive relation gives rise to the BRSSS formulation [45],

π = - \frac{4}{3} \frac{η}{τ} - τ_{π} (\partial_{τ} π + \frac{4}{3} \frac{π}{τ}) + \frac{3}{4} (β_{π π} - \frac{4}{3}) \frac{τ_{π}}{τ} π^{2} .

(77)

Correspondingly, the second-order transport coefficients λ₁ and βππ are related, and the evaluation in a weakly coupled medium

λ_{1} = \frac{3}{4} (β_{π π} - \frac{4}{3}) η τ_{π} = \frac{5}{7} η τ_{π},

(78)

is confirmed [48].

The inclusion of higher-order moments in the coupled equation will improve the quantitative characterization of the system evolution, as has already been noted. These corrections due to higher-order moments can be explicitly incorporated in the equations, even for the simplest truncation, by one additional term $\propto c_{1} L_{2} / τ$ in Eq. (74b) [95]. This term can also be absorbed into the collision term via a redefinition of the relaxation time:

\frac{1}{τ_{R}} \to [1 + \frac{c_{1} τ_{R}}{τ} \frac{L_{2}}{L_{1}}] \frac{1}{τ_{R}} \equiv \frac{Z_{η / s}^{- 1}}{τ_{R}} .

(79)

It is then straightforward to see that because the ratio in the factor is related to g₂(w),

\frac{L_{2}}{L_{1}} = - \frac{b_{2}}{w + g_{2} (w) + a_{2}},

(80)

one may solve the coupled moment equations to arbitrary orders, provided g₂(w) is given precisely. In terms of the simplest truncation that dynamically involves only $L_{0}$ and $L_{1}$ , as in the case of viscous hydrodynamics, Eq. (79) effectively gives a renormalized relaxation time τR. Because η/s ∝ τR, this procedure equivalently renormalizes the ratio of shear viscosity to the entropy density, η/s, by a multiplicative renormalization factor Zη/s. If one further considers the attractor as a generic generalization to an out-of-equilibrium hydrodynamic mode, and hence substitutes the attractor solution of g₂(w) in the calculations, the renormalized τR (or η/s) in an out-of-equilibrium system is obtained. In Fig. 8, the factor Zη/s is obtained via the attractor solution of g₂(w). Unless the system is close to equilibrium, τ/τR ≫ 1, the out-of-equilibrium effects will reduce the value of η/s, thus making the out-of-equilibrium system closer to an ideal fluid [71, 113-115]. The results of a numerical test of the renormalization scheme are shown in Fig. 5, where given a renormalized η/s, even the solution of the two-moment equations achieves good agreement in comparison with the exact solution.

Fig. 8.

Renormalization of η/s by out-of-equilibrium effects. Figure taken from Ref. [77].

3 Phenomenological development

Relativistic hydrodynamics, incorporated with a lattice QCD-based EoS, viscosity, and initial state fluctuations, has been used as a precision tool to understand the dynamics of strongly coupled QGP and experimental flow observables (see reviews [6-8]). Fluid dynamics serves as a universal long-wavelength description of the system’s macroscopic degrees of freedom from the QGP to the hadronic gas phase. This strongly coupled description naturally breaks down as the system becomes increasingly dilute within its hadronic phase at low temperatures. One must then transit to a microscopic transport description. The numerical realizations of hadronic transport models are UrQMD [116, 117], JAM [118], and SMASH [119]. Such a hydrodynamics–hadronic-transport hybrid theoretical framework has successfully described and even predicted various types of flow correlation measurements with remarkable precision [120-123].

In this section, we review the recent phenomenological developments in modeling the full 3D dynamics at intermediate collision energies, current state-of-the-art constraints on the QGP transport properties, understanding collective behavior in small systems, and interdisciplinary cross-talk with other fields of science.

3.1 Hydrodynamic perspectives on Beam Energy Scan and longitudinal dynamics

Quantifying the phase structure of QCD matter is one of the primary questions in relativistic heavy-ion physics. First-principles lattice QCD calculations have established that hadron resonance gas (HRG) transitions to the QGP phase as a smooth crossover at vanishing net baryon density [129]. Meanwhile, many model calculations conjectured the presence of a first-order transition accompanied by a critical point at some finite net baryon density in the QCD phase diagram (see e.g., [16, 130] for a review). Current heavy-ion experiment programs, such as the RHIC BES program [131-134] and the NA61/SHINE experiment at the Super Proton Synchrotron (SPS) [135, 136], as well as future experiments at the Facility for Antiproton and Ion Research (FAIR) [137, 138], Nuclotron-based Ion Collider Facility [139], and JPARC-HI [140], routinely produce hot and dense QCD matter to probe an extensive temperature and baryon chemical potential region in the phase diagram. Measurements from such a beam energy scan of heavy-ion collisions provide a unique opportunity to quantitatively study the nature of the QCD phase transition from hadron gas to QGP at different net baryon densities.

Figure 9 presents our current somewhat limited knowledge of the nuclear matter phase diagram. A recent lattice QCD calculation of higher-order susceptibilities at μB = 0 allows for a Taylor series extrapolation of the thermodynamic quantities to moderate finite μB [124]. This work showed that the phase transition from HRG to QGP remains as a smooth crossover to μB∼250–300 MeV. The phase diagram region with μB < 300 MeV corresponds to mid-rapidity heavy-ion collisions with a collision energy $\sqrt{s} ≳ 15$ GeV. The estimated crossover line agrees with the chemical freeze-out temperatures and net baryon chemical potentials extracted from the STAR BES hadron yield measurements [125, 126]. Off-diagonal susceptibility correlations offer additional insights into the chemical freeze-out conditions in relativistic heavy-ion collisions [141, 142]. The three blobs in Fig. 8 indicate the averaged fireball trajectories for typical Au+Au collisions at $\sqrt{s} = 200,19.6$ , and 5 GeV. A fireball created in the lower collision energy can probe a larger-μB but lower-T region of the QCD phase diagram.

Fig. 9.

(Color online) Sketch of QCD phase diagram together with current and future heavy-ion experimental programs. As relativistic heavy-ion collisions evolve from QGP to the hadron gas phase, they explore the phase diagram of hot and dense QCD matter. Lattice QCD calculations identified a smooth crossover between QGP and hadron gas for μB < 250 MeV [124]. The chemical freeze-out points extracted from thermal fits at RHIC [125, 126] are shown. The three blobs represent the fireball trajectories of Au+Au collisions at RHIC BES energies mapped onto the QCD phase diagram event by event using the dynamical framework [127]. Their brightness is proportional to the fireball space–time volume weighted by T⁴. A crossover phase transition is assumed in the simulations. This figure was adapted from [128].

To establish definitive links between observables and structures in the phase diagram, detailed dynamical modeling of all stages of heavy-ion collisions is required. Precise flow measurements of the hadronic final state, together with phenomenological studies, can elucidate the collective aspects of the baryon-rich QGP and extract the QGP transport properties, such as its viscosity and charge diffusion coefficients. Because relativistic heavy-ion collisions have complex dynamics in multiple stages, a fully integrated theoretical framework is required to provide reliable estimates of the dynamical evolution of the collisions and all relevant sources of fluctuations.

Over the past decade, extensive phenomenological studies have focused on relativistic heavy-ion collisions at the LHC and the top RHIC energies (see e.g., [6, 7], for a review). Recently, increasing interest has shifted toward studying heavy-ion collisions in the intermediate energy regime. At $\sqrt{s} \sim O (10)$ GeV, heavy-ion collisions strongly violate longitudinal boost invariance and require full 3D modeling of their dynamics [143, 144]. It is important to employ initial conditions with nontrivial rapidity dependence. Complex 3D collision dynamics can be approximated by parametric energy deposition [144-149]. Nontrivial dynamics could be included to model the energy loss in individual nucleon–nucleon (NN) collision impacts. Initial state models have been built based on classical string deceleration [127, 150]. There are 3D initial conditions based on hadronic transport simulations [151-153]. These provide interesting correlations between the longitudinal energy distribution and flow velocity. The earlier work of Anishetty, Koehler, and McLerran (AKM) in 1980 [154], which found that nuclei were significantly compressed and excited when they collided at extreme relativistic energies, was pursued only sporadically. This AKM picture was generalized recently to understand early-stage baryon stopping from the CGC-based approaches in the fragmentation region [155, 156]. The incoming nucleons or valence quarks are decelerated and compressed by the strong shock wave of small-x gluons. Compared to the previous phenomenological approaches, this CGC-based approach consists of ab initio calculations of baryon stopping at high energies. This approach becomes less applicable at intermediate collision energies. Last but not least, there are recent theoretical developments from a holographic approach at intermediate couplings to understand the initial energy density and baryon charge distributions [157].

At a lower collision energy, the longitudinal Lorentz contraction becomes weaker on the colliding nuclei. The overlapping time for the two nuclei to pass through each other is significant, τ_overlap∼2R/sinh(y_beam) [127, 152, 158]. Here, R is the nuclear radius and the beam rapidity $y_{beam} = arccosh (\sqrt{s} / (2 m_{N}))$ . Therefore, it is important to understand the pre-equilibrium dynamics during this period. Dynamical initialization schemes have been developed to model this extended interaction region in heavy-ion collisions (see Fig. 10). They interweave the initial collision stage with hydrodynamics on a local basis, as the two nuclei pass through each other. The initial state energy–momentum and conserved charge density currents are treated as source inputs to the hydrodynamic fields at different times,

Fig. 10.

(Color online) Illustration of relativistic heavy-ion collisions at an intermediate collision energy. When the two colliding nuclei overlap with each other early, dynamical initialization connects the initial state collision impact with hydrodynamics.

\partial_{μ} T^{μ ν} = J_{source}^{ν} (τ, \vec{x})

(81)

\partial_{μ} J_{B}^{μ} = ρ_{B, source} (τ, \vec{x}) .

(82)

Such a dynamic initialization scheme was initially proposed in Refs. [159, 160] and has been adopted by several groups [127, 153, 161, 162].

Solving the equations of motion of hydrodynamics at low energies requires an EoS, which describes the thermodynamic relations of nuclear matter at finite baryon density. The current lattice QCD techniques cannot directly compute such an EoS, because of the sign problem [163]. However, at vanishing net baryon density or μB = 0 GeV, higher-order susceptibilities have been computed by lattice QCD [124]. These susceptibility coefficients were used to construct nuclear matter EoS at finite baryon densities through a Taylor expansion [164-166]. These estimated EoS are reliable within the region where μB/T≲2 in the phase diagram. For the region where the temperature is below ∼ 150 MeV, the lattice QCD EoS is glued with an EoS for HRG. To ensure energy and momentum conservation in the hydrodynamics plus hadronic transport approaches, the particle species in the HRG EoS need to be the same as those in the transport model. A mismatch in the particle content of the HRG EoS could lead to 5–10% variations in the particle yields and flow observables [167]. At μB = 0 GeV, matching EoS between the two phases was performed on the trace anomaly [168, 169]. At a finite baryon density, susceptibility coefficients are matched individually, and then the thermal pressure is constructed. Full-fledged hydrodynamics plus hadronic transport simulations with an EoS at finite baryon density, NEoS, have been applied to heavy-ion collisions at intermediate collision energies [164]. That study found that the enforcement of strangeness neutrality improved the description of relative particle yields for multistrange particles measured in Pb+Pb collisions at the top SPS energy.

The dynamical initialization and EoS at finite baryon density are two essential components to enable hybrid simulations for heavy-ion collisions at intermediate collision energies. A fully integrated framework [127, 164, 170] was shown to reproduce the rapidity dependence of particle production as well as the collision energy dependence of the STAR vn{2} (n = 2,3) flow measurements in Au+Au collisions from 200 to 7.7 GeV [144, 171]. Remarkably, this preliminary calculation can produce a similar nonmonotonic collision energy dependence present in the experimental triangular flow data measured at the RHIC BES phase I [172], without the need for a critical point in the phase diagram. Therefore, it is essential to understand the interplay among the duration of dynamical initialization, the variation in the speed of sound, and the T and μB dependence of the specific shear viscosity. This work demonstrated a critical role of theoretical modeling in elucidating the origin of the nonmonotonic behavior seen in the RHIC measurements. Phenomenological studies of precise anisotropic flow measurements from the upcoming analysis of RHIC BES II will be able to further constrain the μB dependence of the QGP shear and bulk viscosities [144, 152].

Modeling relativistic heavy-ion collisions beyond the Bjorken boost invariance offers us a new dimension to study the fluctuations of flow anisotropy as a function of particle rapidity. The flow correlations between different rapidity regions reflect the longitudinal fluctuations of local energy density profiles in the full 3D dynamics [148, 173-175]. The fluctuations of the anisotropic flow coefficients along the longitudinal direction elucidate how the shape of the fireball varies as a function of rapidity. The longitudinal fluctuations will cause the initial eccentricity vectors to fluctuate from one space–time rapidity to another. Consequently, the anisotropic flow coefficients decorrelate as a function of particle rapidity. A recent systematic study [149] showed that initial state fluctuations and thermal fluctuations in hydrodynamics were equally important for understanding the centrality dependence of flow decorrelation measurements at the LHC. The RHIC BES program will systematically study the collision energy dependence of the rapidity flow fluctuations from 200 GeV down to 7.7 GeV. The program offers a unique opportunity to study the interplay between thermal production and collision transport (stopping mechanisms). Such measurements combined with phenomenological modeling could lead to strong constraining power for our understanding of the longitudinal dynamics in heavy-ion collisions [127]. Future measurements with identified particles have the potential to shed further light on the initial distributions of conserved charges (net baryons, strangeness, and electric charges) at different collision energies.

The realistic dynamical simulations of heavy-ion collisions at the RHIC BES energies lay the foundation to quantitatively understand the out-of-equilibrium stochastic fluctuations when the collision systems evolve close to a conjectured QCD critical point in the phase diagram. Near the QCD critical point, the relaxation times of the critical fluctuations grow large and rapidly become comparable with the system size [16, 130]. The rapid expansion of the collision system drives the fluctuations related to the QCD critical point out of equilibrium [24, 25, 176]. Therefore, these "critical slowing down" dynamics require realistic event-by-event simulations to address quantitatively the extent to which they are out of equilibrium.

There are two primary approaches to this problem in the literature. One is to explicitly evolve these fluctuations in a stochastic hydrodynamics framework [177]. In heavy-ion physics, numerical simulations of stochastic hydrodynamics have been developed by several groups to study thermal and critical fluctuations [20, 22, 178-180]. However, because these stochastic fluctuations are local δ-functions in the coordinate space, one needs to take care of the numerical cut-off dependence to ensure that the simulation results are physical. The other approach studies the deterministic evolution of the two-point correlation function of the fluctuations. This approach was pioneered by Andreev in the 1970s in the nonrelativistic case [181], and it is often referred to as "hydro-kinetic" [21, 23-25, 176, 182].

The hydrokinetic approach should be consistent with the stochastic hydrodynamics approach for the two-point correlation function. A side-by-side comparison between these two approaches would be extremely valuable in improving both theories. On the one hand, the renormalization of the hydrodynamic equations is well controlled in the hydrokinetic framework. It can guide the stochastic hydrodynamics on how to introduce a UV cut-off to regulate multiplicative noise in simulations. On the other hand, stochastic hydrodynamics can estimate the effect of nonlinearity in the fluid dynamics on the evolution of the two-point function, which is neglected in the hydrokinetic formalism. Meanwhile, it is straightforward to access higher-order correlation functions in stochastic hydrodynamics.

Recent studies investigated the deterministic "hydrokinetic" formalism on simplified (1+1)D hydrodynamic backgrounds [26, 27]. They found that the feedback contributions to the thermodynamic quantities from the out-of-equilibrium fluctuations are on the order of 10^-4, which can be safely neglected in the simulations.

To deliver quantitative predictions for experimental signals of the critical fluctuations at the RHIC BES II, we need to further develop the theoretical frameworks in the following directions. In the "hydrokinetic" approach, the effects of flow gradients on the critical fluctuations must be addressed quantitatively with realistic 3D event-by-event hydrodynamic simulations. Although the works [23, 25] derived the equations of motion for two-point correlation functions under a general flow background, there are substantial challenges to implementing these equations in the state-of-the-art hydrodynamic framework. In addition, a theoretical formulation is needed to map these two-point correlation functions from the coordinate space to momentum correlations among particles, which is an essential step to providing theoretical predictions for measurements. Finally, generalization to n-point correlation functions (n > 2) could enhance the sensitivity of the experimental signals, but this requires substantial effort on the theoretical side [183].

3.2 Quantitative characterization of QGP transport properties

3.2.1 Charge diffusion

The fluid dynamics of the bulk QGP medium are driven by local pressure gradients. The RHIC BES program allows us to study the evolution of nonvanishing conserved charge currents, which is also dependent on the gradients of chemical potentials. Therefore, the interplay among different gradient forces in the dynamics of conserved charge currents can elucidate novel transport processes inside the medium, namely the charge diffusion constants and heat conductivity. These QGP transport coefficients are to date poorly constrained but as important as the specific shear and bulk viscosities. Because individual quarks carry multiple quantum charges, the diffusion currents of conserved charges are coupled with each other. A quantitative understanding of how multiple conserved charge currents diffuse in and out of the expanding fluid cells has led to the development of the next generation of dynamical frameworks. Such a framework can unravel the detailed chemical aspects of the QGP. This topic has been stimulating interest in developing realistic initial conditions for conserved charge distributions [127, 153, 191, 192] and detailed modeling of the QGP chemistry [164, 193-195].

The net baryon diffusion is driven by the local gradients of the net baryon chemical potential with respect to temperature in the hydrodynamic evolution. Causal Israel–Stewart-like equations of motion for the net baryon diffusion current were derived based on the Grad’s 14-moment and Chapman–Enskog methods [49, 196, 197]. Recently, such formulations were generalized to include multiple diffusion currents with their cross couplings [198]. The functional forms of the transport coefficients for the diffusion matrix have been studied in transport models [199, 200]. Additional coupling terms with the shear and bulk viscous tensors appear in the second order of the gradient expansion [49, 197].

Phenomenologically, the net baryon diffusion current transports more baryons from forward rapidities to the mid-rapidity region [67, 196, 201, 202]. The gradients of μB/T act against local pressure gradients and decelerate baryon charges with respect to the bulk fluid cells along the longitudinal direction. Therefore, the shape of the rapidity distribution of the net protons shows a strong sensitivity to baryon diffusion [67, 201]. The cross-diffusion between the net baryon and net strangeness induces an oscillating distribution of the net strangeness current at later times [198]. It will be interesting to see how this pattern is mapped to the final state hadron correlations, such as Kaons and Λ. Measurements of identified particle rapidity distributions will play an important role in unraveling the charge diffusion processes in heavy-ion collisions.

The net baryon diffusion process in the hydrodynamic phase can only transport the net baryon charges by ∼1 unit in rapidity [67]. As the bulk fluid rapidly explodes along the z direction, transporting net baryon charge back to the mid-rapidity region during hydrodynamic evolution is very challenging at high energies. It turned out to be difficult to reproduce the small but nonzero net proton rapidity distribution at 200 GeV measured by the BRAMHS Collaboration [201] through baryon diffusion only. The measurement suggests that there is a large baryon stopping at the early stage of heavy-ion collisions. Allowing the baryon charge to fluctuate to the string junction [203] in the initial state, we can reproduce the net baryon rapidity distributions at 200 GeV. In fact, this model can consistently reproduce the net proton distribution measured by the STAR Collaboration down to 7.7 GeV at mid-rapidity [144]. The rapidity distribution in RHIC BES phase II will further help to constrain the initial state baryon stopping in this phenomenological model. Because the net proton rapidity distribution is sensitive to both the initial state stopping and baryon diffusion [201], independent experimental observables are needed to disentangle these two effects.

Recently, the charge balance functions of identified hadrons were proposed as independent observables to constrain the charge diffusion constants of QGP [57, 179, 195]. Quark–antiquark pairs can be thermally produced in local fluid cells. As a $q \bar{q}$ pair randomly walks through the medium, it develops a finite correlation length that is controlled by the local fluid expansion and charge diffusion constants. This spatial correlation eventually maps to momentum correlations among hadrons into which the $q \bar{q}$ pair hadronizes. In Ref. [195], the authors found that a larger diffusion constant in the QGP medium leads to wider azimuthal and rapidity distributions for the K⁺K^- and $p \bar{p}$ correlations. Meanwhile, the underlying hydrodynamic flow acts as a boost to these correlation functions, which results in a focus on the azimuthal distributions of balance functions. Therefore, measurements of balance functions can provide an independent constraint on the QGP charge diffusion constants.

3.2.2 Specific shear and bulk viscosities

The shear and bulk viscosities, normalized by entropy density, characterize the dissipation of energy––momentum currents in a medium. Extracting the specific QGP viscosity has been the main theme of phenomenological studies in relativistic heavy-ion collisions over the past 20 years. Because of the strongly coupled nature of QGP, it is difficult to calculate these transport coefficients from first principles. Meanwhile, phenomenological constraints drawn from comparisons with precise anisotropic flow measurements have led our field to understand the properties of QGP. Adopting Bayesian statistical analysis has become a standard approach to systematically constrain the QGP transport properties [204, 205].

Figure 11 summarizes the collective effort of quantifying the QGP transport properties over the past eight years. As theoretical tools are being developed rapidly to include more realistic physics, the extraction of the QGP transport properties becomes more systematic. Chronologically, before the implementations of bulk viscous effects in the dynamical models, the saturation-based IP-Glasma and EKRT initial conditions prefer an effective shear viscosity of 0.12–0.20 in the hydrodynamic phase to achieve a simultaneous description of all orders of harmonic flow coefficients [184, 188]. Owing to large pressure gradients and finite initial radial flow in the IP-Glasma initial conditions, a temperature-dependent bulk viscosity is essential to balance the consequent strong flow and reproduce the mean pT measurements in the heavy-ion collisions. At the same time, the introduction of the bulk viscosity reduced the extracted QGP shear viscosity by almost 50% [13, 206]. By adopting the EoS from the latest lattice QCD calculations [207, 208], the Duke-OSU group presented the first Bayesian inference on the temperature-dependent shear and bulk viscosities [189, 190] using the flow measurements at the LHC. There was tension of the Bayesian extracted bulk viscosity (ζ/s)(T) [189, 190] with the parameterization used in the IP-Glasma hybrid framework [185] in 2018. This difference was greatly reduced in the 2019 updated simulations. On the one hand, the significant changes in the IP-Glasma hybrid framework stem from allowing the peak temperature of bulk viscosity to drop from T_peak = 180 MeV to 160 MeV [28, 186]. A lower T_peak with a smaller (ζ/s)_max is favored by hadron mean pT measurements in peripheral Pb+Pb and Au+Au collisions, in which the maximum temperatures at the starting time of hydrodynamic simulations are close to, or even below, 180 MeV. On the other hand, a more flexible prior parameterization of (ζ/s)(T) in the JETSCAPE preliminary Bayesian analysis allows a large QGP bulk viscosity in the posterior distribution for the model to reproduce flow measurements [167].

Fig. 11.

(Color online) Summary of phenomenological constraints on the temperature-dependent QGP specific shear and bulk viscosities from 2012 to 2019. This figure is taken from [128]. The QGP viscosities are constrained by the collective flow measurements at the top RHIC and LHC collision energies. The results from open-source dynamical frameworks using three different initialization models are compared. The IP-Glasma initial conditions with the hydrodynamics plus hadronic transport approach continuously improve the constraints on the QGP viscosity [13, 28, 184-187]. The blue hatched area in (h) indicates the variation of (ζ/s)(T) with (upper) and without (lower) a pre-equilibrium KøMPøST effective kinetic theory (EKT) phase [187]. The result obtained using the EKRT initial condition is from Ref. [188]. The Bayesian extracted QGP (η/s)(T) and (ζ/s)(T) [189, 190] were constrained by flow measurements in p+Pb and Pb+Pb collisions at the LHC. The orange bands indicate a 90% confidence level. The preliminary Bayesian analysis calibrated with combined RHIC and LHC flow measurements was presented by the JETSCAPE Collaboration [167].

The preliminary Bayesian analysis from the JETSCAPE Collaboration demonstrated the first simultaneous calibration using flow observables at the top RHIC and LHC energies, which differ by an order of magnitude. Such a combined analysis showed strong constraints on the temperature dependence of the QGP shear and bulk viscosities around the crossover temperature region [167]. Heavy-ion collisions at the RHIC and LHC offer a wide dynamical range for the theoretical framework to explore the parameter space. At the same time, different types of off-equilibrium corrections at the particlization stage introduce sizable theoretical uncertainties to the Bayesian extraction, demanding more theoretical work in this direction. That work also emphasized that performing closure tests in Bayesian analysis was an essential step before extracting any physical information from experimental measurements [167].

The next-generation Bayesian analysis requires improving the state-of-the-art theoretical framework to access new physics in dynamical simulations. A recent study bridges the event-by-event IP-Glasma initial state [28, 209] and viscous hydrodynamics [210, 211] using an effective QCD kinetic theory (EKT), KøMPøST [37, 38]. The EKT in the pre-equilibrium stage can drive the collision system sufficiently near the local thermal equilibrium and smoothly matches the full energy–momentum tensor to viscous hydrodynamics. Remarkably, this framework can quantitatively reproduce a variety of flow measurements in heavy-ion collisions from 200 to 5020 GeV [187] with an effective (η/s)_eff = 0.12 and temperature-dependent bulk viscosity (see Fig. 11h). A direct comparison with the simulations without the EKT phase [28] showed that the conformal EKT generates a faster expansion than viscous hydrodynamics at this early stage. A similar finding was observed earlier for free-streaming dynamics [212]. This additional pre-equilibrium phase thus leads to a 35% larger extracted QGP bulk viscosity to describe the identified particle mean pT. A recent work [213] pointed out that the breaking of conformality in the pre-equilibrium stage could be important in the extraction of bulk viscosity. These studies demonstrated the significant phenomenological impact of a realistic model of the early stage of heavy-ion collisions on constraining the QCD bulk viscosity.

3.3 Challenges and opportunities in small systems

The RHIC and LHC collide a variety of nuclei, which offers us measurements to study the collectivity as a function of the collision system size. As the size of a QGP droplet shrinks, the lifetime for the strongly coupled hydrodynamic evolution becomes shorter. Therefore, the final state particles’ momentum distributions and correlations can reveal more information about the early-stage dynamics [214, 215]. On the theoretical side, small systems are phenomenological grounds for understanding how heavy-ion collisions achieve macroscopic hydrodynamic behavior, chemical, and kinetic equilibration from states far out of equilibrium. Nonhydrodynamic modes can potentially play an important role in pushing the hydrodynamic framework to its limits. At the same time, the increasing roles of fluctuations and nonflow correlations stress the ability of experiments to unambiguously identify flow signatures.

The collective flow and its hydrodynamic description are most robust in the central Pb+Pb and Au+Au collisions. Therefore, by first constraining the model parameters with flow measurements in heavy-ion collisions and then extrapolating to small systems, we can provide a stringent parameter-free test for hydrodynamic models. This approach was pioneered in Ref. [217]. A more systematic study has recently been carried out using the state-of-the-art IP-Glasma + MUSIC + UrQMD hybrid framework [28]. This work achieved remarkable success in describing the system size dependence of the flow measurements over more than two orders of magnitude in particle multiplicity. With a single set of model parameters, this theoretical framework can quantitatively describe particle production, radial and anisotropic flow observables, and multiparticle correlations from Pb+Pb to p+Pb collisions at the LHC and from Au+Au to p+Au collisions at the top RHIC energy. This work demonstrated that the universal hydrodynamic response to collision geometry dominated the flow production in collisions with dN^ch/dη 10 at mid-rapidity. On the contrary, the correlation between v₂ and the initial momentum anisotropy from the pre-hydrodynamic phase becomes stronger in the lower-multiplicity collisions [35]. This work provided phenomenological evidence that flow in low-multiplicity collisions can elucidate the early dynamics of the collisions. A similar finding was shown in a study of the correlation between the system’s elliptic flow coefficient and the initial energy–momentum tensor [218]. Last but not least, a recent work [219] demonstrated that the quark coalescence is essential to reproduce the measured the momentum anisotropy of produced hadrons with transverse momenta 3 < p_T < 6 GeV in small systems. That work indicated the collective behavior was originated from partonic degrees of freedom in high multiplicity p-Pb collisions at the LHC.

The anisotropic flow in high-energy p+p collisions still challenges our understanding of the underlying dynamics in these small systems. The computed two-particle cumulant v₂{2} from the IP-Glasma hybrid framework increases as the charged hadron multiplicity decreases [28], which is not seen in the flow measurements at the LHC. A comparison with the results from [65, 217] suggests that the momentum anisotropy from the pre-equilibrium Glasma phase might be too strong in p+p collisions. Moreover, recent work pointed out that pure hydrodynamic evolution introduces an excessively strong nonlinear cubic response of elliptic flow to the initial eccentricity in p+p collisions, which results in a positive four-particle cumulant C₂{4} opposite to the experimental measurements [220]. Elucidating the real dynamics in p+p collisions requires more theoretical progress on quantifying the contributions from pre-equilibrium dynamics, thermal and longitudinal fluctuations, nontrivial correlations from local conservation laws, and nonhydrodynamic modes.

Finally, it is instructive to quantify the extent to which small collision systems have been pushing the hydrodynamic framework to its limits. Conditions for the nonlinear causality bounds of second-order hydrodynamics were derived for radially expanding systems [221] and for a general flow background [216, 222]. These nonlinear causality bounds set strong constraints on the maximal allowed viscous pressure in dynamical simulations. One can check the nonlinear causality condition, including the bulk viscous pressure, by computing the following ratio:

r = \frac{ζ}{τ_{Π}} \frac{1}{e + P + Π} + c_{s}^{2} .

(83)

The causality condition requires the ratio r < 1 [216]. Fig. 12a compares the distribution of this ratio from individual fluid cell in a typical 0–5% p+Au collision with that from a 30–40% Au+Au collisions. Because of the strong expansion rate and the consequent large negative bulk viscous pressure, the causality conditions in a typical p+A collision are ∼ 20% closer to the bound than those in an A+A collision.

Fig. 12.

Upper panel: Color contour plots for the initial event-by-event energy densities of U+U, Au+Au, Ru+Ru, O+O, ³He+Au, d+Au, p+Au, and p+p collisions at 200 GeV from IP-Glasma initial conditions [28]. Panel (a): Check the nonlinear causality condition with bulk viscous pressure [216] in typical 0–5% p+Au and 30–40% Au+Au collisions at 200 GeV. Panel (b): Distribution of the ratio of bulk viscous pressure over thermal pressure, ∏/P, in fluid cells for 0–5% p+Au and 30–40% Au+Au collisions. Panel (d): Effect of regulating negative total pressure on the elliptic flow coefficients in 0–5% p+Au and 30–40% Au+Au collisions at 200 GeV [35].

The large pressure gradients and the consequent violent expansion in small systems can result in negative total (thermal + bulk viscous) pressure in a significant fraction of fluid cells, as shown in Fig. 12. These bubbles may cause unstable cavitation inside the QGP [223-225]. It would be preferable to switch the fluid dynamic description to a dilute transport approach before the total pressure becomes negative. In the absence of such an advanced theoretical framework, we can roughly estimate the phenomenological impact of these negative-pressure regions on final flow observables by numerically regulating the size of the bulk viscous pressure to be less than the thermal pressure. This modification leads to a sizable variation in the elliptic flow in 0–5% p+Au collisions in Fig. 12c. For p_T < 1.5 GeV, the associated theoretical uncertainty is comparable to the effects resulting from varying the second-order transport coefficients [35]. The situation in semi-peripheral A+A collisions is much better, with a negligible effect on the final flow observables. Whereas this ad hoc numerical regulation can only provide us with a rough estimate, Fig. 12c demonstrates that the phenomenological descriptions of small systems with the standard second-order viscous hydrodynamics approach the limits of the model. Anisotropic hydrodynamics, which reproduce the free-streaming limit at large viscosity, are a good theoretical tool for providing more robust guidance on this issue in small systems [226, 227].

3.4 Interdisciplinary connections

Relativistic heavy-ion collisions embrace a richness of physics that expands multiple energy scales. At the early stages of the collisions, the dynamics of approaching hydrodynamic behavior from far out of equilibrium have strong connections with the reheating of the early universe following inflation and over-occupied cold atomic gases (see a recent review [228]). The collision geometry in heavy-ion collisions can be used as a powerful microscope to image structures inside nuclei and nucleons. Meanwhile, it has become an active topic to apply cutting-edge machine learning techniques to study the complex dynamics of heavy-ion collisions.

3.4.1 Nuclear structure physics

Although the bulk evolution in relativistic heavy-ion collisions is above 100 MeV in temperature, nuclear structure physics at a much lower energy scale can still play an important role in precision studies of flow observables. The Monte-Carlo Glauber model is a key element in computing the initial states for hydrodynamical modeling of ultrarelativistic heavy-ion collisions. The spatial configurations of nucleons inside the colliding nuclei are inputs based on nuclear structure physics for each collision event. Hydrodynamics can then efficiently transform the shape of the initial energy density profile to the momentum anisotropy of the final state particles. Therefore, the measurements of anisotropic flow coefficients and their fluctuations provide a tool to image the event-by-event shape fluctuations of the colliding nuclei. Furthermore, small collision systems, such as proton–lead collisions, can image the subnucleon fluctuations inside protons in a similar fashion [34, 229, 230].

The structure of nuclei has been approximated for a long time by independent particle models in which the nucleons inside a nucleus are treated as a collection of free and point-like fermions. The colliding nuclei in the relativistic heavy-ion collisions are highly Lorentz-boosted along the beam direction. In the lab frame, because of the relativistic time dilation effect, individual nucleons are frozen in their spatial positions as the two nuclei collide with each other. The spatial positions of the nucleons inside the nucleus are usually sampled independently from parametric Woods–Saxon distributions. More realistically, correlations among nucleons at short distances are not negligible and are known as short-range correlations (SRC). These NN correlations were unambiguously observed in a series of dedicated experiments [231-234]. Studies in Ref. [235-237] showed that including realistic NN correlations had sizable effects on the generated initial eccentricity of the energy profile. The NN correlations have sizable effects on the spectrum of anisotropic flow coefficients in central heavy-ion collisions, in which the geometric distortion of the overlapping region from the impact parameter is minimized. The ratios of the elliptic flow to triangular flow were found to be sensitive to the spatial configurations of the colliding nuclei with and without SRC [87, 238]. A recent study explored the effect of a possible octupole deformation of ²⁰⁸Pb on the v₂/v₃ ratio [239].

Furthermore, collisions with deformed nuclei are particularly interesting as one can use the nucleus’ intrinsic deformation as an additional control to study the hydrodynamic conversion from spatial eccentricity to momentum anisotropy and inform the transport properties of the QGP. At the top RHIC energy, the deformed U+U collisions were studied with Au+Au collisions. Full-overlap U+U collisions have the potential to study the hydrodynamic behavior of elliptic flow in large and dense collision systems, in addition to the nonlinear path length dependence of radiative parton energy loss [240]. Experimentally, such a study requires to select "tip–tip" collisions, defined to occur when the major axes of the uranium nuclei lie parallel with the beam direction, from the "body–body" events, where the major axis of each nucleus is perpendicular to the beam direction [241-243]. However, the theoretically proposed triggers based on a two-component Glauber model were not effective in STAR measurements [244]. The measurements were more consistent with the saturation-based initial condition models [245, 246]. Recently, new types of correlations between elliptic flow and mean transverse momentum were proposed, which are sensitive to the nuclear deformation [247]. Similarly, at LHC energies, ¹²⁹Xe+¹²⁹Xe collisions were measured in addition to ²⁰⁸Pb+²⁰⁸Pb collisions. In contrast to the spherical ²⁰⁸Pb, the ¹²⁹Xe nucleus has a prolate deformation. Hydrodynamic simulations [248] predicted that the elliptic flow coefficients would be 25% larger in ¹²⁹Xe+¹²⁹Xe than those in ²⁰⁸Pb+²⁰⁸Pb in 0–5% central collisions. This strong enhancement of v₂{2} in central ¹²⁹Xe+¹²⁹Xe was confirmed by LHC measurements [249-251].

Stable heavy nuclei are also neutron rich. The neutron density profiles usually have larger RMS radii than those of proton densities, which is referred to as the neutron skin of the nucleus. Measurements of neutron skins in heavy nuclei are critical for studying the EoS of neutron-rich nuclear matter [252]. In relativistic heavy-ion collisions, the neutron skin has a negligible effect on regular flow observables. However, the RHIC proposed the use of isobar collisions to minimize the background flow signals in search of chiral magnetic effects (CME). The success of this experimental program relies on a precision-level understanding of the flow background in $_{44}^{96}$ Ru+ $_{44}^{96}$ Ru and $_{40}^{96}$ Zr+ $_{40}^{96}$ Zr collisions. A recent work [253] pointed out that the neutron skin in $_{40}^{96}$ Zr can result in a stronger magnetic field in peripheral collisions and consequently lead to reductions by a factor of 2 in the difference in CME signals from $_{44}^{96}$ Ru+ $_{44}^{96}$ Ru with respect to $_{40}^{96}$ Zr+ $_{40}^{96}$ Zr collisions. Moreover, realistic nucleon configurations based on density functional theory calculations for $_{44}^{96}$ Ru and $_{40}^{96}$ Zr nuclei showed sizable differences in the particle elliptic flow v₂ between the two isobar collision systems [254].

3.4.2 Applications of statistics and machine-learning techniques

Extracting quantitative physics from the complex dynamics in relativistic heavy-ion collisions requires the production and analysis of a high volume of data, both from experiments and numerical simulations. The data acquisition rates at the RHIC and LHC experiments grow factorially with the detectors’ updates. Precise measurements, such as six- and eight-particle azimuthal correlations [29], are pushing large-scale and high-performance numerical simulations. Parallelization using graphics processing units (GPUs) is employed to speed up the event-by-event simulations [255, 256]. Meanwhile, we need to adopt advanced statistical tools, such as unsupervised learning, Bayesian inference, and deep learning (DL) techniques, to systematically analyze the high volume of simulation data to extract the physical properties of the QGP. As shown in Fig. 13, the dynamical modeling of heavy-ion collisions involves multiple model parameters, each of which can influence several experimental measurements. Therefore, it is difficult to identify a single experimental observable to constrain one physical parameter. Machine learning is a collection of generic algorithms that allow computers to find nontrivial correlations and patterns in data samples with minimum bias.

Fig. 13.

(Color online) Illustration of the workflow of statistical analysis, such as Bayesian inference, in relativistic heavy-ion collisions. The individual physical parameters that pass through the physical model can influence multiple experimental observables. Statistical analysis provides a systematic way to constrain multiple model parameters using a collection of measurements. Trained model emulators, together with Monte-Carlo Markov Chain processes, are often applied in statistical analysis to explore the high-dimensional parameter space efficiently. This figure is modified based on Bass’s talk at Quark Matter 2017 [257].

Principle component analysis (PCA) has been widely used to study the roles of fluctuations and correlations in relativistic heavy-ion collisions. PCA is a statistical technique for extracting the dominant components in fluctuating data by transforming a set of correlated variables into independent variables via orthogonal transformations. This method was first introduced to analyze the event-by-event fluctuations of anisotropic flow coefficients [258-260] and the breaking of flow factorization ratios [261]. The PCA procedure in that study found two dominant contributions to the two-particle correlation function. The leading component was identified with the event plane anisotropic flow vn(p_T), and the subleading component was responsible for factorization breaking in hydrodynamics [261]. A recent study [262] showed that the subleading principal components of anisotropic flow can reveal details of the hydrodynamic response to small-scale structures in the initial density profiles. Similar studies have also been performed to understand the event-by-event fluctuations in particle multiplicity and radial flow [261, 263]. PCA has recently been applied in unsupervised learning to test whether a machine can directly discover anisotropic flow coefficients from high-volume simulation data without explicit instructions from human beings [264]. Because PCA can disentangle and extract the dominant components in data, it is a standard technique to perform dimensional reduction in statistical analysis as well as model emulation.

Bayesian inference or Bayesian analysis is based on the Bayes’ theorem to derive posterior distributions of model parameters by constraining prior information with experimental data. It relies on a solid theoretical framework and precise experimental measurements. As discussed in Sec. 3.2.2, it is a systematic way to constrain high-dimensional model parameters using multiple experimental measurements. Over the past few years, Bayesian analysis has gradually become a standard tool to perform phenomenological extractions of the QGP properties from experimental measurements. The MADAI Collaboration initiated a community effort of applying Bayesian analysis to heavy-ion physics [204]. Bayesian analysis has been applied to constrain the QCD EoS [205], QGP transport coefficients [167, 189, 190, 265-269], initial longitudinal fluctuations [270], QCD jet energy loss distribution [271, 272], and heavy-quark diffusion coefficients [273]. These works have paved the way for precision physics in relativistic heavy-ion collisions. Bayesian analysis is a powerful tool to systematically extract information from experimental data using a well-established model. However, it has a limitation, as the produced posterior distributions are influenced by the defined subjective priors. There is no unique and unbiased method for choosing a prior in Bayesian analysis.

As an outlook, the next generation of Bayesian analysis in heavy-ion physics is evolving toward drawing global constraints with multiple subfields, such as combining bulk flow observables with electromagnetic and QCD jet probes. Using unified theoretical frameworks and observables that probe multiple scales in heavy-ion collisions, more reliable and systematic information can be extracted from Bayesian analysis.

DL is a branch of machine learning methods that are based on artificial neural networks to capture highly correlated patterns/features from big data. It has achieved tremendous success in science and technology. The great advantage of DL over conventional methods is its ability to extract hidden features from highly dynamical and complex nonlinear systems, such as relativistic heavy-ion collisions.

The application of neural networks was pioneered in heavy-ion physics more than 20 years ago to determine the impact parameter of every heavy-ion collision based on final particle momentum distributions [274]. More recently, supervised learning with deep convolutional neural networks (CNN) was used to identify the nature of QCD phase transition [275, 276]. The complex dynamics of heavy-ion collisions hide the experimental signals of a first-order phase transition and potentially enhanced fluctuations near a critical point. These works demonstrated that the deep CNN could provide a powerful and efficient "decoder" to extract information about the QCD EoS from final particle momentum distributions. If the QGP fluid transitions to HRG through a first-order phase transition, the conserved net baryon density can clump together in space because of the spinodal decomposition. This phenomenon leaves characteristic imprints on the spatial net density distribution in every collision event, which can be detected by DL techniques [277]. Meanwhile, the impact of spinodal decomposition on the measurable particle momentum space information is still challenging for DL methods to recognize.

DL techniques have also been applied to learn and mimic the nonlinear dynamics in relativistic heavy-ion collisions. A deep neural network was designed to learn and capture the main features of relativistic hydrodynamics [278]. By treating the initial energy density and flow velocity as inputs, the trained neural network can quantitatively reproduce the realistic event-by-event hydrodynamical evolution. This work demonstrated that DL could speed up event-by-event simulations of heavy-ion collisions by orders of magnitude by replacing real hydrodynamic simulations with neural network predictions.

DL is gaining popularity in a variety of aspects of heavy-ion physics to hunt for hidden experimental signals of important physics [279-282].

4 Summary

Relativistic heavy-ion collisions interconnect nuclear and high-energy physics. Experiments at the RHIC and LHC push the field to evolve rapidly and introduce many surprises in the process. The relativistic expanding QGP is a unique fluid to study emergent many-body physics of strong QCD interactions. Meanwhile, the QGP shares many universal collective features with other strongly coupled systems in condensed matter and cosmology. Therefore, understanding the properties of the QGP will not only advance our knowledge about the many-body aspects of strong interactions, but also enrich the crosstalk with other fields in physics.

To date, theoretical studies on the out-of-equilibrium fluid dynamics in heavy-ion collisions are mostly limited to highly symmetric systems. These are QGP media experiencing either Bjorken expansion (0+1)D or Gubser expansion (1+1)D and are further simplified with respect to the conformal EoS. However, in realistic heavy-ion collisions, neither of these symmetry conditions is rigorously satisfied. For instance, Bjorken boost invariance is apparently broken in proton–lead collisions. It is a necessary step to examine the robustness of these attractors in realistic QGP in heavy-ion collisions, such that out-of-equilibrium hydrodynamics can provide a theoretical foundation for interpreting the collective behavior in small colliding systems.

Some generalized discussions have been carried out on various aspects. By numerically solving the second-order viscous hydrodynamics without Bjorken symmetry and with a nonconformal EoS, attractors are identified in multiple channels [70]. In a similar manner, attractors from kinetic theory solutions without conformal symmetry are reproduced [70, 283]. In the context of kinetic theory, the $L$ -moments can be generalized in 3+1D expanding systems as well, by replacing the weight of Legendre polynomials with spherical harmonic functions. Attractors are found in the induced coupled moment equations [115]. These observations suggest the universal existence of attractors in out-of-equilibrium systems [284].

Apart from the out-of-equilibrium extensions, macroscopic fluid dynamics is an effective, robust, and efficient description of the bulk dynamics of heavy-ion collisions. It is a bridge that connects fundamental QCD theory and experiments. Phenomenological studies combined with high-precision flow measurements at the RHIC and LHC have been driving the field toward an era of precision. The adoption of modern statistics and machine learning techniques, such as Bayesian analysis and DL, have become a popular and standard approach to extract the QGP transport properties systematically. The next challenge lies in how to quantify the model uncertainty in the theoretical framework. We are expecting increasing input from first-principles calculations to reduce the theoretical uncertainties. Analytical solutions based on Bjorken and Gubser flow profiles have been widely adopted to validate hydrodynamic frameworks [255, 256, 285-287]. In recent years, there have been a few formulations of relativistic causal hydrodynamics, known as Israel–Stewart [44], BRSSS [45], DNMR [49], relativistic third-order dissipative hydrodynamics [51], and anisotropic hydrodynamics [112, 226, 227]. These theories differ from each other based on the number of velocity gradient terms included in the equations of motion for the dissipative tensors. Systematic comparisons among different hydrodynamic theories have been carried out by comparing the results to the exact solution of the Boltzmann equation [288] under high degrees of symmetry [112, 289-291]. Future extensions of such comparisons to full (3+1)D will help to standardize the fluid dynamic model for relativistic heavy-ion collisions. A similar community-wide effort has been carried out for different transport frameworks in the field of low-energy nuclear physics [292, 293]. Flow observables from the RHIC BES program and future FAIR/NICA experiments bring us to a new era of full 3D dynamics beyond the boost-invariant approximation. With the development of dynamic initialization schemes that interweave the 3D collision dynamics with fluid simulations, we are starting to quantify the initial baryon stopping and study the collectivity of QGP in a baryon-rich environment. This framework provides a reliable baseline to search for critical point signals in the upcoming RHIC BES II measurements. Flow measurements in small systems offer a window to study the early-stage dynamics of QGP. Understanding the collective origin in small systems has led to state-of-the-art theory development of rapid hydrodynamization as well as in- and out-of-equilibrium relativistic hydrodynamics.

References

[1]

E. Shuryak,

Strongly coupled quark-gluon plasma in heavy ion collisions

. Rev. Mod. Phys. 89, 035001 (2017). https://doi.org/10.1103/RevModPhys.89.035001