Hydrodynamization and resummed viscous hydrodynamics

Contents

1. Introduction

Relativistic viscous hydrodynamics (vHydro) stands as the primary theoretical framework for describing the spatiotemporal evolution of the rapidly expanding quark–gluon plasma (QGP) generated in ultrarelativistic heavy-ion collisions.¹ Despite its notable success, understanding how hydrodynamics can provide a reliable description of the rapidly expanding system created in these collisions poses a formidable challenge. Traditionally, hydrodynamics has been understood as a truncation of a gradient expansion.² Thus, its region of validity was assumed to be limited to small gradients relative to the inverse microscopic scales of the problem. The gradient expanded theory was historically regarded as a universal macroscopic limit inherent to microscopic theories, achievable at suitably late times. Recent investigations, however, have revealed that the gradient expansion may possess a zero radius of convergence for flow configurations pertinent to the QGP, both in strong coupling scenarios and within kinetic theory.^3,4,5,6 Consequently, constructing and refining a hydrodynamic theory by systematically incorporating higher-order terms in this series is fraught. As a result, the conventional notion that relativistic hydrodynamics is only applicable under conditions where gradients of macroscopic quantities are small, as derived from the gradient expansion, appears no longer well justified, if not altogether unnecessary. Ultimately, these revelations prompt a reevaluation of the very definition of vHydro to assess its range of applicability in the context of heavy-ion collisions.

While the phenomenological success of fluid-dynamical models was initially interpreted as an indication of the rapid isotropization and thermalization of the QGP,⁷ subsequent model calculations have suggested that such an interpretation may have been premature.^{8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25} This is due to the fact that systems that significantly depart from equilibrium may already exhibit hydrodynamic behavior through a phenomenon known as “hydrodynamization”. This novel feature is particularly relevant in rapidly expanding fluids like the QGP. In practice, the applicability of linearized vHydro is not boundless; it will eventually fail to be accurate as viscosity values become sufficiently large or when applied at very early times when there are large gradients. Nevertheless, even in such extreme scenarios, effective theories capable of describing the QGP exist, with anisotropic hydrodynamics (aHydro) being the most notable among them.^{26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44} In general, hydrodynamization is now expected to occur at a timescale ${\tau }_{\mathrm{\text{hydro}}}$, which is shorter than the corresponding timescale for isotropization, being driven by a novel dynamical attractor whose details vary according to the theory under consideration, e.g., kinetic theory, hydrodynamics, holography, etc.^{18,45,46,47,48} Such attractor solutions show that hydrodynamics displays a new degree of universality in far-from-equilibrium scenarios regardless of the details of the initial state of the system. In fact, the approach to the dynamical attractor effectively wipes out a subset of information about the specific initial condition used for the evolution, before the true equilibrium state and consequently, full thermalization, is reached.

In the realm of kinetic theory and conventional statistical mechanics, thermalization is characterized by the emergence of isotropic thermal one-particle distribution functions for the partons constituting the QGP. In the context of high-energy heavy-ion collisions, the substantial longitudinal expansion rate causes the QGP fireball’s center to gradually relax toward an approximately isotropic state, on a timescale of ${\tau }_{\mathrm{\text{iso}}}\gtrsim 4$fm/c.⁴⁹ Notably, the hydrodynamization of the fireball seems to occur on a shorter timescale, as reviewed in Ref. 25.^a For conformal systems, a crucial factor for determining proximity to universal attractor behavior is the dimensionless variable $w\equiv \tau T$.⁴⁵ For conformal fluids undergoing Bjorken expansion,⁵¹ this variable is proportional to the inverse of the Knudsen number Kn, with $1/T$ representing the microscopic relaxation timescale. For small gradients where $w\gg 1$, the system exhibits dynamics consistent with the universal hydrodynamical attractor. However, in the large gradient regime where $w\ll 1$, the system’s dynamics are dominated by nonhydrodynamic modes (i.e., modes in the linearized dynamics with nonzero frequency even for a spatially homogeneous system⁵²), whose evolution depends on the precise initial condition assumed. Considering a fixed proper time after the collision, this implies that, near the edge of the QGP, the system is more sensitive to the genuinely nonequilibrium dynamics linked to nonhydrodynamic modes. Consequently, certain nonuniversal aspects of the underlying theory, whether rooted in kinetic theory or holographic duality, start to influence the spatiotemporal evolution of the QGP. In such scenarios, a decision must be made regarding as to which underlying microscopic theory best captures the relevant physics. Given the system’s dilute nature close to the QGP’s edge (large mean free path), a kinetic theory approach appears preferable in this spatial region (large $|x|$, $|y|$ and/or $|\varsigma |$, with the $\varsigma$ being the spatial rapidity). In addition, at very early times after the nuclear pass through, the energy densities are sufficiently high to justify a weak-coupling approach that takes into account perturbative scattering processes, with modified power counting stemming from large amplitude gluon fields. For this reason, approaches such as aHydro, which are explicitly based on a kinetic theory approach are perhaps better suited for understanding the QGP at early times and in dilute regions near the edges of the plasma.

In Fig. 1, I show the relaxation time approximation (RTA) evolution of the typical pressure anisotropy (${\mathcal{𝒫}}_L/{\mathcal{𝒫}}_T$) as a function of the rescaled time variable $\bar{w}=\tau /{\tau }_{\mathrm{\text{eq}}}$. The black solid line indicates the attractor solution. The green dot-dashed lines represents the first-order (Navier–Stokes (NS)) limit of the hydrodynamical evolution. The gray dashed lines represent the evolution of the system with different assumed initial pressure anisotropies. The inset table provides a conversion of $\overline{w}$ to typical times in fm/c, assuming a shear viscosity to entropy density ratio, $\eta /s=0.2$, and an initial temperature of $T_0=500$MeV at ${\tau }_0=0.25$fm/c. As can be seen from this figure, irrespective of the initial condition, the particular solutions (gray dashed lines) approach the dynamical attractor rather quickly. I note, importantly, that the attractor solution is independent of the assumed initial conditions and shear viscosity to entropy density ratio when plotted as a function of $\overline{w}$. Second, I note that this figure demonstrates that at early times after the nuclear pass through, large pressure anisotropies are generated, independent of the assumed initial conditions. This fact was the historical basis for the introduction of aHydro, which is the focus of the second half of this paper.

Summarizing, in this paper, I will review recent progress in understanding the hydrodynamization of the QGP, including attractors in hydrodynamics theories themselves and different kinetic theory models including, for example, the RTA and effective kinetic theories (EKTs) based on quantum chromodynamics (QCD). Following this, I will discuss the formalism of aHydro and present a subset of the phenomenological tests of this framework.

2. Hydrodynamic Attractors

As mentioned above, in recent years, remarkable progress been made in understanding the evolution of systems subject to the far-from-equilibrium conditions pertinent to heavy-ion collisions, see e.g., Refs. 25, 53 and 54 for reviews. Substantial progress has occurred in elucidating the onset of hydrodynamic behavior in relativistic systems, both in strong coupling scenarios^{11,13,15,16,17,18,48,55,56,57} and in weak-coupling/kinetic theory.^{5,6,18,20,46,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,75,77,78,79,80} The introduction of hydrodynamic attractors⁴⁵ has led to an emerging framework for estimating typical values of ${\tau }_{\mathrm{\text{hydro}}}$, which is the timescale at which some form of hydrodynamics becomes applicable as a function of charged particle multiplicity. Notably, authors in Refs. 18 and 63 attempted to make estimates of this timescale.

To understand the main arguments underpinning these advancements, we first consider a simplified model of the QGP based on conformal kinetic theory^58,59,81 within the RTA.^82,83,84,85 This approximation has proven to be a powerful tool for gaining insights into the onset of hydrodynamic behavior in rapidly expanding systems.^{36,58,59,75,83,84,85} Despite the seemingly rudimentary nature of the relaxation time approach and the nonconformal nature of the QGP, the conclusions drawn are expected to possess a semi-universal character.⁸⁶ Fundamentally, the behavior seen reflects the interplay between free streaming and dissipative dynamics in an approximately conformal system undergoing rapid longitudinal expansion. As discussed below, simple conformal kinetics captures the scaling behavior anticipated from QCD at weak coupling⁸⁷ and the scaling properties exhibited by strongly coupled conformal theories based on the anti-de Sitter/conformal field theory (AdS/CFT) correspondence.¹¹ While the initial focus here is on conformal systems for simplicity, it is worth noting recent evidence indicating the existence of hydrodynamic attractors for certain moments of the one-particle distribution function in nonconformal systems.^75,74,88,89 I will discuss these most recent advances in detail later in this paper. To start with, I will review our understanding of attractors within dissipative hydrodynamics, focusing first on the conformal dynamics of systems subject to $0+1$-dimensional ($0+1$D) boost-invariant Bjorken flow.

2.1. Müller–Israel–Stewart-type theories

We start by assuming a $0+1$D system, which is transversally homogeneous and boost-invariant, following the framework presented in Ref. 51. Consequently, all variables depend solely on the longitudinal proper time, denoted as $\tau =\sqrt{t^2-z^2}$. The metric is considered “mostly minus” with space-time coordinates $x^{\mu }=(t,x,y,z)$, and the line element is given by $d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=d{t}^2-d{x}^2-d{y}^2-d{z}^2$, where $g^{\mu \nu }$ is the metric tensor in Minkowski space. The longitudinal space-time rapidity, denoted as $\varsigma$, is given by $\varsigma ={tanh}^{-1}(\frac{z}{t})$. Assuming conformality, the system has an equation of state stemming from $N_{\mathrm{\text{dof}}}$ massless degrees of freedom, which is Landau matched to the general nonequilibrium energy density. Within this framework, it follows that $\epsilon ={\epsilon }_0(T)=3P_0(T)$, and the temperature T is expressed as $\gamma {\epsilon }^{1/4}$, with $\gamma$ being proportional to $N_{\mathrm{\text{dof}}}^{-1/4}$. Additionally, in the context of a longitudinally boost-invariant system, the flow velocity is represented as $u^{\mu }=(cosh\varsigma ,0,0,sinh\varsigma )$ (Bjorken flow).

In this section, we will employ kinetic theory to derive the dynamical attractors for aHydro and second-order vHydro. To achieve this, we begin with the Boltzmann equation within the RTA⁹⁰

The momentum-and energy-independent relaxation time ${\tau }_{\mathrm{\text{eq}}}$ is defined as ${\tau }_{\mathrm{\text{eq}}}=5\eta /sT$,^91,92 where $\eta$ represents shear viscosity, T is the local effective temperature obtained through Landau matching and s is the entropy density. For simplicity, in this section I will assume a classical Boltzmann distribution for the underlying thermal distribution function.

In kinetic theory the covariantly conserved energy–momentum tensor is given by

with $\int dP=\int d^{3}p/{(2\pi )}^{3}E$ being the appropriate Lorentz-invariant integration measure.⁹⁰ The local energy density is obtained via $\epsilon =u_{\mu }{u}_{\nu }{T}^{\mu \nu }$ and the shear stress tensor is obtained by projecting with a transverse projector where ${\mathrm{\Delta }}_{\alpha \beta }^{\mu \nu }=({\mathrm{\Delta }}_{\alpha }^{\mu }{\mathrm{\Delta }}_{\beta }^{\nu }+{\mathrm{\Delta }}_{\beta }^{\mu }{\mathrm{\Delta }}_{\alpha }^{\nu })/2-{\mathrm{\Delta }}^{\mu \nu }{\mathrm{\Delta }}_{\alpha \beta }/3$ is the tensor projector orthogonal to the flow constructed using ${\mathrm{\Delta }}^{\mu \nu }=g^{\mu \nu }-u^{\mu }{u}^{\nu }$.

Bjorken symmetry and conformal invariance allow for the reduction of the energy–momentum conservation laws derived from the first moment of the Boltzmann equation to a single dynamical equation

involving the energy density and $\pi \equiv {\pi }_{\varsigma }^{\varsigma }$. In second-order hydrodynamic frameworks such as the Müller–Israel–Stewart (MIS)^93,94,95 and Denicol–Niemi–Molnár–Rischke (DNMR)^96,97 approaches, the 14-moment approximation for the single-particle distribution function is employed to derive a differential equation for the evolution of $\pi$. For a $0+1$D system undergoing Bjorken expansion, the dynamical for $\pi$ equation takes the following form : where $\dot{}=d/d\tau$. In RTA, ${\beta }_{\pi \pi }=38/21$ and ${\tau }_{\pi }={\tau }_{\mathrm{\text{eq}}}$ in the complete second-order calculation of DNMR,^{91,96,97,98,99} while in MIS ${\beta }_{\pi \pi }=4/3$ and ${\tau }_{\pi }=6{\tau }_{\mathrm{\text{eq}}}/5$.^100,b Solving Eqs. (4) and (5) allows one to determine the dynamical evolution of a viscous fluid described by second-order hydrodynamics. From these simple equations we can understand the emergence of hydrodynamic attractor behavior, as demonstrated in Ref. 45.

2.2. $0+1$D conformal aHydro

Later in this paper we will discuss the general $3+1$D nonconformal aHydro formalism. Here, we only present the formalism for a $0+1$D conformal system. In the $0+1$D case, aHydro requires only one anisotropy direction and parameter, $\widehat{n}$ and $\xi$, respectively. This leads to a distribution function of the form^101,102

where $\mathrm{\Lambda }$ can be interpreted as the local “transverse temperature” field and $-1<\xi <\infty$ is the anisotropy field. For a $0+1$D system boost-invariant system, both of these fields only depend on the longitudinal proper time. Using this form, the integrals defining the conformal energy density, transverse pressure and longitudinal pressure factorize, resulting in with^27,103

$$\mathcal{ℛ}(\xi )=\frac{1}{2}\left[\frac{1}{1+\xi }+\frac{arctan\sqrt{\xi }}{\sqrt{\xi }}\right],$$(7)

$${\mathcal{ℛ}}_T(\xi )=\frac{3}{2\xi }\left[\frac{1+({\xi }^2-1)\mathcal{ℛ}(\xi )}{\xi +1}\right],$$(8)

$${\mathcal{ℛ}}_L(\xi )=\frac{3}{\xi }\left[\frac{(\xi +1)\mathcal{ℛ}(\xi )-1}{\xi +1}\right],$$(9)

which satisfy $3\mathcal{ℛ}=2{\mathcal{ℛ}}_T+{\mathcal{ℛ}}_L$ (the conformal isotropic pressure is $P_0=\epsilon /3$). In these expressions, L and T represent the directions parallel and perpendicular to $\widehat{n}$, respectively. Conventionally, the anisotropy direction is aligned along the beam line direction in heavy-ion applications ($\widehat{n}=\widehat{z}$). Employing Landau matching, one obtains $\epsilon ={\epsilon }_0(T)$. For a conformal system with a momentum-and energy-independent relaxation time, this leads to To dynamically evolve the system, we need an equation of motion for $\xi$, since $\mathrm{\Lambda }$ is already connected to the temperature (energy density) using Eq. (10).

To do this, we will employ the following moment of the Boltzmann distribution³³ :

which will be important for the aHydro approach. Using the Boltzmann equation in RTA (1), the equation of motion for this moment is

We note that $I^{\mu \nu \lambda }$ is symmetric with respect to interchanges of $\mu$, $\nu$ and $\lambda$ and traceless in any pair of indices (massless particles/conformal invariance). In an isotropic system, one finds $I_{xxx}=I_{yyy}=I_{zzz}=I_0$ with

Using the aHydro form for the one-particle distribution function, one obtains with, e.g., $I_{uuu}\equiv u_{\mu }{u}_{\nu }{u}_{\lambda }{I}^{\mu \nu \lambda }$, etc. Taking the $zz$ projection of Eq. (12) minus one-third of the sum of its $xx$, $yy$ and $zz$ projections gives our second equation of motion³² which can be used to obtain the proper-time evolution of the anisotropy parameter.

2.2.1. Connection with shear stress tensor and the inverse Reynolds number

To enable a more transparent comparison between the equations of motion for aHydro and those of standard vHydro, one can make a change of variables in Eq. (15) to express it in terms of the shear stress tensor component $\pi$. This can be achieved by using $\pi =P_0-{\mathcal{𝒫}}_L$, yielding

In the left panel of Fig. 2, we present the dependence of $\overline{\pi }$ on $\xi$, determined by Eq. (16). In the right panel, we present the dependence of $\xi$ on $\overline{\pi }$ as obtained through the numerical inversion of Eq. (16). It is important to note that in aHydro, $\overline{\pi }$ is constrained within $-2/3<\overline{\pi }<1/3$, which is a restriction associated with the naturally occurring positivity of longitudinal and transverse pressures in aHydro. Next we note that, for a 0+1D boost-invariant system, the magnitude of $\overline{\pi }$ is proportional to the inverse Reynolds number

As a consequence, a series in $\overline{\pi }$ can be understood as an expansion in $R_{\pi }^{-1}$.

We will also need the relation between the time derivatives of $\pi$ and $\xi$. This relation can be obtained from Eq. (16)

which upon using Eqs. (16) and (4) gives where ${\overline{\pi }}^{\prime }\equiv d\overline{\pi }/d\xi$.

Plugging (19) into (15), one obtains

with and the understanding that $\xi =\xi (\overline{\pi })$ with $\xi (\overline{\pi })$ being the inverse function of $\overline{\pi }(\xi )$ (shown in the right panel of Fig. 2). Written in this form, we can see explicitly that the aHydro second-moment equation sums an infinite number of terms in the expansion in the inverse Reynolds number (17). In the next section we will expand this equation in powers of the inverse Reynolds number through second order in order to compare it to standard vHydro.

2.2.2. Small $\xi$ expansion

To establish the final connection to standard vHydro, one can expand Eq. (20) in terms of $\xi$ around $\xi =0$. To accomplish this, we need the $\xi$ expansions of the various functions in order to construct an explicit inversion and rewrite the equations exclusively in terms of $\overline{\pi }$. At the second order in $\xi$, the expression is

Inverting the relationship between $\overline{\pi }$ and $\xi$ to second order in $\overline{\pi }$ gives which results in

Applying this to the equation of motion (20) and keeping only terms through ${\pi }^2$ gives

where on the left-hand side, we employed the fact that the energy density can be eliminated by expressing it in terms of the transport coefficients as and we have relabeled ${\tau }_{\mathrm{\text{eq}}}\to {\tau }_{\pi }$ in order to cast the equations in second-order vHydro form. Note that, to linear order in $\pi$, Eq. (25) agrees with the previously obtained RTA second-order vHydro results.^{91,96,97,98,99} However, one can go beyond this truncation by extending the power series expansions used above or simply solving the nonlinear equation of motion (20) numerically without Taylor expansion.⁶¹

2.3. Comparison of attractors in vHydro and aHydro models

In this section, we will explore the hydrodynamic attractor behavior of aHydro and contrast it with the corresponding outcomes in MIS and DNMR theories. In each of these cases, the dynamics of the system is governed by solving the differential equations for $\epsilon$ and $\pi$. However, to establish a connection with the literature, we adopt the methodology introduced in Ref. 45 and introduce the dimensionless “time” variable

with which one may define the amplitude which is related to the shear correction $\pi$ as follows : As a consequence, a solution for the proper-time evolution of the energy density uniquely determines the dependence of the amplitude $\phi$ on w. Additionally, it is important to note that the positive energy condition imposes a bound on $\phi$ within the region $0\le \phi \le 1$.¹⁰⁴

2.3.1. The hydrodynamic attractor in the DNMR framework

The change of variables from $\{\epsilon ,\pi \}\to \{w,\phi \}$ is convenient because it allows one to express the coupled set of first-order ordinary differential equations (ODEs) for $\{\epsilon ,\pi \}$ in terms of a single first-order ODE for $\phi (w)$.⁴⁵ In the case of MIS and DNMR, this procedure gives

where ${\phi }^{\prime }=d\phi (w)/dw$, $c_{\pi }\equiv {\tau }_{\pi }T$ and $c_{\eta }=\eta /s$ (with $c_{\pi }=5c_{\eta }$ in the cases considered here). After defining the rescaled variable $\overline{w}=w/c_{\pi }$ one obtains The form above makes it clear that the solution only depends on the ratio $c_{\eta /\pi }\equiv c_{\eta }/c_{\pi }=(\eta /s)/({\tau }_{\pi }T)$ and the value chosen for ${\beta }_{\pi \pi }$. To connect these equations with the RTA Boltzmann one must set $c_{\eta /\pi }=1/5$.

Using the MIS value ${\beta }_{\pi \pi }=4/3$ one obtains

which agrees with Eq. (9) of Ref. 45; however, for RTA this value for ${\beta }_{\pi \pi }$ is incorrect. Using the correct value for ${\beta }_{\pi \pi }=38/21$ obtained using the complete second-order formalism of DNMR and again neglecting quadratic terms in $\pi$ one obtains^c

Following Ref. 45, the underlying dynamical attractor can be deduced from Eq. (30) through a procedure similar to the slow-roll expansion in cosmology.¹⁰⁵ This process can be described as follows: a small parameter $\delta$ is formally introduced as a pre-factor in the term $\overline{w}\phi {\phi }^{\prime }$ in (31). It is assumed that the solution of the differential equation $\phi (\overline{w};\delta )$ can be expressed as a power series expansion $\phi (\overline{w};\delta )={\phi }_0(\overline{w})+{\phi }_1(\overline{w})\delta +\mathcal{𝒪}({\delta }^2)$. After considering all orders, one may then take the limit $\delta \to 1$. The 0th-order truncation is obtained by solving the simple quadratic equation

Out of the two possible solutions, only one is stable and remains finite in the large $\overline{w}$ limit and is

It is possible to compute higher-order corrections in the slow-roll expansion, however, it is better to define a numerically determined attractor solution using the boundary condition ${lim}_{\overline{w}\to 0}\overline{w}\phi {\phi }^{\prime }=0$, which then implies that⁴⁵

Numerical solution of Eq. (31) subject to this boundary condition allows one to construct the attractor without having to resort to the slow-roll expansion. In the next section, we generalize this analysis to determine the dynamical attractor in aHydro.

2.4. The hydrodynamic attractor in the aHydro framework

I now discuss how to obtain the hydrodynamic attractor in aHydro. This can be achieved by combining the two first-order aHydro differential equations into a single second-order differential equation written in terms of $\phi$ and w. To obtain the necessary aHydro dynamical equation, one can use the following identity :

and (20). We first express Eq. (20) in terms of $\phi$ and w and then using $\tau {\partial }_{\tau }log\epsilon =4(\phi -1)=-4/3+\pi /\epsilon$, one finds Inserting this into Eq. (37) and converting to $\overline{w}$ gives our final result for the aHydro attractor equation Note that above $\xi =\xi (\overline{\pi })=\xi (4\phi -8/3)$ and likewise for ${\overline{\pi }}^{\prime }$. It is worth noting that the aHydro equation above does not contain $c_{\eta /\pi }$ since it cancels explicitly. The aHydro attractor solution remains universal when plotted as a function of $\overline{w}$. This universality holds true for second-order vHydro approximations as well.

Equation (39), when expanded in gradients (powers of $1/\overline{w}$ here), has zero radius of convergence.¹⁰⁶ Thus, the numerical solution of the differential equation (39) may be considered an all-orders resummation of the gradient series, as in MIS theory.⁴⁵ However, it is important to highlight that the right-hand side of Eq. (39) involves a sum of an infinite number of terms in the inverse Reynolds number. This conceptual difference sets it apart from DNMR, which derived their equations of motion assuming a perturbative series in $R_{\pi }^{-1}$.

In the case of aHydro, solving even the 0th-order approximation in the slow-roll expansion requires a numerical approach. Hence, we directly proceed to solve the differential equation numerically. Once again, for this purpose, the attractor solution is determined by imposing the same boundary condition as before at $\overline{w}=0$. Utilizing the numerical solution of the slow-roll equation, one finds that in aHydro the attractor boundary condition is

With this boundary condition, we simply numerically solve Eq. (39). Note that the boundary condition above guarantees the positivity of the longitudinal pressure of the attractor solution as $\overline{w}\to 0$.

In the left panel of Fig. 3 we compare the aHydro, MIS and DNMR attractors to the corresponding quantity obtained from the exact solution to the 0+1D RTA Boltzmann equation.^83,84 We also include a curve showing the NS result⁴⁵

As the left panel of Fig. 3 demonstrates, the aHydro attractor solution closely resembles the exact RTA attractor.^46,63,83,84 Given that aHydro involves resummation not only in the Knudsen number but also in the inverse Reynolds number, the remarkable agreement with the exact kinetic theory result suggests that the inverse Reynolds number resummation is critical. This observation could serve as a guide for developing alternative approaches to far-from-equilibrium hydrodynamics that do not depend on perturbative expansions in both the Knudsen and the inverse Reynolds number. Finally, in the right panel of Fig. 3, I present the attractor solution for the pressure ratio which can be obtained using In this figure, we have also included a result from Ref. 61 for the aHydro equations truncated at third order in the inverse Reynolds number. As can be seen from this figure, only the resummed aHydro attractor gives positive longitudinal pressure for all $\overline{w}$.

In Fig. 4, I show the DNMR and aHydro attractors expressed in terms of ${\mathcal{𝒫}}_L/{\mathcal{𝒫}}_T$ along with specific numerical solutions corresponding to different assumed values of the initial plasma anisotropy ratio. As can be seen from this figure, all numerical solutions collapse to their respective attractors for $\overline{w}\gtrsim 2$. This indicates that information about the precise initial conditions for ${\mathcal{𝒫}}_L/{\mathcal{𝒫}}_T$ quickly becomes ignorable. In this figure, we also see the feature that truncations in the inverse Reynolds number lead to attractors with negative longitudinal pressure at early times.

3. Attractors in Nonequilibrium Kinetic Theory

Having established the existence of hydrodynamic attractors in conformal $0+1$D dissipative hydrodynamics, I now turn to a discussion of the existence of dynamical attractors in nonequilibrium kinetic theory. I review the solutions obtained for both conformal (massless)⁶³ and nonconformal (massive) systems.^75,88,89 In kinetic theory, one finds that attractors exist at very early times for a large set of moments of the Boltzmann equation implying the existence of an attractor for the entire one-particle distribution function.^63,d

3.1. Nonconformal Boltzmann equation in RTA

We continue with the Boltzmann equation in RTA

where f is the one-particle distribution function, $p^{\mu }$ is the particle four-momentum and C is the collision kernel with $u^{\mu }$ being the four-velocity of the local rest frame (LRF) and $a\cdot b\equiv a^{\mu }{b}_{\mu }$. The quantity ${\tau }_{\mathrm{\text{eq}}}$ appearing above is the relaxation time, which will be precisely specified below. For the equilibrium distribution, we will follow Ref. 85 and assume a Boltzmann distribution^e Herein, we will assume Bjorken flow, in which case in Milne coordinates one has $u^{\tau }=1$ and $u^{x,y,\varsigma }=0$, where $\tau$ is the longitudinal proper time, $\tau =\sqrt{t^2-z^2}$ and $\varsigma$ is the spatial rapidity, $\varsigma ={tanh}^{-1}(z/t)$.

3.1.1. Thermodynamic variables

For a single-component massive gas obeying Boltzmann statistics, the equilibrium thermodynamic quantities are

with $\widehat{m}\equiv m/T$ and $K_n$ being modified Bessel functions of the second kind. Above $𝔫$ is the number density, s is the entropy density, $\epsilon$ is the energy density and P is the pressure. These satisfy $\epsilon +P=Ts$ and, from the above relations, one can determine the speed of sound squared

3.1.2. Relaxation time for a massive gas

For a massive system, the shear viscosity $\eta$ can be expressed as^37,90,109

with and $L_n(x)$ being modified Struve functions. For fixed specific shear viscosity, $\bar{\eta }\equiv \eta /s$, using $\epsilon +P=Ts$ one obtains with Note that, in the massless limit, $m\to 0$, one has $\kappa (\widehat{m})\to 12$, $\epsilon \to 3P$ and $\gamma \to 1$, giving the usual conformal RTA relaxation time For small $\widehat{m}$, one has and in the large $\widehat{m}$ limit, one has

In Fig. 5, I plot $\gamma (\widehat{m})$. As can be seen from this figure, $\gamma (\widehat{m})$ goes to unity in the massless limit and grows linearly at large $m/T$, which corresponds either to fixed temperature and large mass or fixed mass and small temperature. The fact that $\gamma (\widehat{m})\ge 1$ implies that a massive gas always relaxes more slowly to equilibrium than a massless one in physical units, however, it is unclear a priori how things will change as a function of the rescaled time $\overline{\tau }\equiv \tau /{\tau }_{\mathrm{\text{eq}}}$. We note that the strong enhancement of the relaxation time at low temperatures modifies the asymptotic approach to equilibrium.

3.1.3. Exact solution for the distribution function and its solution

Here, I review the derivation of the exact solution presented in Ref. 85 and derive the integral equation obeyed by all moments of the distribution function. I will also generalize the results contained in that reference to the full set of integral moments. For details concerning the derivation of the exact solution in the nonconformal case, I refer the reader to Ref. 75.

In the case of one-dimensional boost-invariant expansion ($0+1$D), all scalar quantities depend only on the longitudinal proper time $\tau$. To describe boost-invariant $0+1$D dynamics, one can introduce a spacelike vector $z^{\mu }$, which is orthogonal to the fluid four-velocity $u^{\mu }$ in all frames and corresponds to the z-direction in the LRF of the matter.^28,29 The requirement of boost invariance implies that $f(x,p)$ may depend only on three variables, $\tau$, w and $p_T$,^110,111 with the boost-invariant variable w defined by^f

Using w and $p_T$ one can define Using these variables, one can write the energy and the longitudinal momentum as

$$E=\frac{vt+wz}{{\tau }^2},$$(57)

$$p_L=\frac{wt+vz}{{\tau }^2},$$(58)

and the Lorentz-invariant momentum integration measure becomes

When written in terms of these variables, the $0+1$D RTA Boltzmann equation takes a particularly simple form^83,84,85

with ${\tau }_{\mathrm{\text{eq}}}$ specified in Eq. (50) and the equilibrium distribution function given by The exact solution to Eq. (60) is^{82,83,84,85,112,113,114} where $f_0(w,p_T)$ is the initial distribution function specified at $\tau ={\tau }_0$ and the damping function D is defined as

Here, I will assume that the initial distribution function $f_0$ can be expressed in spheroidally deformed form^101,102

where ${\xi }_0$ is the initial anisotropy parameter and ${\mathrm{\Lambda }}_0$ is the initial transverse momentum scale. For $-1<{\xi }_0<0$, this corresponds to an initially prolate distribution in the LRF and, conversely, for ${\xi }_0>0$ this corresponds to an initially oblate distribution function. For ${\xi }_0=0$, one obtains an isotropic Boltzmann distribution function as the initial condition.

3.1.4. The integral equation obeyed by all moments

To understand the emergence of the kinetic attractor, I introduce the following moments of the one-particle distribution function^63,67

In principle, powers of $p_T^2$ could also appear in a general moment, however, such moments can be expressed as a linear combination of the two-index moment appearing above using $p^2=m^2$ to write $p_T^2={(p\cdot u)}^2-{(p\cdot z)}^2-m^2$.

Some specific cases of ${\mathcal{ℳ}}^{nl}$ map to familiar quantities, e.g., $n=1$ and $l=0$ maps to the number density $𝔫={\mathcal{ℳ}}^{10}$, $n=2$ and $l=0$ maps to the energy density and $n=0$ and $l=1$ maps to the longitudinal pressure, $P_L$. The transverse pressure, $P_T$, can be obtained by using $p_T^2={(p\cdot u)}^2-{(p\cdot z)}^2-m^2$ to obtain $P_T={\mathcal{ℳ}}^{20}-{\mathcal{ℳ}}^{01}-m^2{\mathcal{ℳ}}^{00}$.

For a Boltzmann equilibrium distribution function, these moments reduce to

I will present results for these general moments scaled by their equilibrium values, i.e.,

In the late-time limit ($\tau \to \infty$), if the system approaches equilibrium, then ${\overline{\mathcal{ℳ}}}^{nl}\to 1$.

In the general case, using the boost-invariant variables introduced earlier, one finds that the general moments can be expressed as

Taking a general moment of Eq. (62), one obtains Evaluating the integrals necessary results in where with where ${}_2{F}_1$ is a hypergeometric function.

Finally, for a momentum-and energy-independent relaxation time, specializing to the case $n=2$ and $l=0$ and requiring conservation of energy $\epsilon (\tau )={\epsilon }_{\mathrm{\text{eq}}}(T)$, also known as Landau matching, we obtain the following integral equation :

This is the integral equation obtained originally in Ref. 85 with the understanding that ${\tilde{H}}^{20}={\widetilde{\mathcal{\mathscr{H}}}}_2$ defined therein. This equation can be solved iteratively for $T(\tau )$ and, once converged to the desired numerical accuracy, the solution can be used in Eq. (69) to compute all moments.

In Fig. 6, I present the resulting solutions for the scaled moments ${\overline{M}}_{nl}$ as a function of rescaled time in the conformal limit ($m=0$). The black solid line indicates the attractor solution for each moment and the nonsolid lines indicate particular solutions with a range of anisotropic initial conditions. As can be seen from this figure, there exists an attractor for all moments shown. This implies that there is an attractor for the entire one-particle distribution function. I remind the reader that $M^{20}=\epsilon$, $M^{10}=𝔫$ and $M^{01}=P_L$. All other moments represent different modes of the distribution function not associated with basic quantities. I note that in the case of ${\overline{M}}^{20}$, it should be equal to 1 at all times due to energy conservation.

In Fig. 7, I present the resulting solutions for the scaled moments ${\overline{M}}_{nl}$ as a function of rescaled time for a nonconformal system ($m=1$GeV). In this figure, I varied the initialization time while holding the initial momentum-space anisotropy and energy density fixed. This is done in order to search for the existence of an early-time attractor or pull-back attractor. As before, the black solid line indicates the attractor solution for each moment and the nonsolid lines indicate particular solutions with different initialization times. As can be seen from this figure, there exists an attractor for all moments with $l>0$. For moments with $l=0$ there are indications that the early-time behavior is only semi-universal. This semi-universality does not detract from the usefulness of the phenomenon, however, it does result in a small inherent uncertainties if one wants to use the attractor in order to extrapolate to early times.

3.2. Attractors in QCD kinetic theory

Having discussed results obtained within the RTA, I present some results obtained in QCD using EKT with an emphasis on what this may imply for freeze-out prescriptions in heavy-ion phenomenology. This section is based on the findings reported in Ref. 71. The use of fluid dynamics is instrumental in understanding ultrarelativistic nuclear collisions.^53,115,116 Within this framework, only select degrees of freedom evolve dynamically, namely those stemming from the energy–momentum tensor $T^{\mu \nu }$. These encompass local temperatures, velocities, and, within vHydro, the shear and bulk viscous tensors. Despite this, experimental measurements do not directly capture fluid-dynamic variables; instead, they measure distributions of particles that have undergone“freeze-out” and stream unimpeded to the detectors. The angular and momentum distributions of these particles provide crucial details about the fluid’s characteristics.¹¹⁷

Extensive research has focused on assessing the ability of different formulations of viscous fluid dynamics to accurately portray the time evolution of energy–momentum tensor components. Notably, it has been found that hydrodynamic constitutive equations, which correlate the stress tensor with flow field gradients, are well fulfilled in systems that are significantly far from equilibrium, particularly in cases characterized by highly symmetrical flow. The freeze-out procedure’s validity far from equilibrium has received considerably less scrutiny. One challenge lies in the fact that models with simplistic momentum dependence, like kinetic theory in RTA, offer limited insights into freeze-out procedure validity. Conversely, strongly coupled models do not possess quasiparticle structure and hence lack underlying particle distributions altogether.

In this section, I discuss the possibility of reconstructing particle distributions from the energy–momentum tensor within EKT for weak-coupling QCD.¹¹⁸ The detailed momentum-dependent structure of the EKT collision kernel allows for a meaningful examination of the freeze-out procedure in this theoretically clean case. Focusing again on $0+1$D far-from-equilibrium Bjorken flow, I compare various moments of the distribution function with predictions from hydrodynamic freeze-out prescriptions. Interestingly, one finds qualitative similarities between EKT and RTA kinetic theory, as well as Israel–Stewart (IS)-type hydrodynamics,^45,46,66 both in early (early-time or pull-back attractor) and late times (late-time or hydrodynamic attractor).

As discussed in the previous section, this attractor behavior extends beyond the energy–momentum tensor components to other integral moments of the one-particle distribution function. While commonly used freeze-out prescriptions well reproduce low-order moments of the distribution at late times, they may not be accurate at early times or when addressing moments that are sensitive to significantly higher momenta than the temperature scale. To study this one can use numerical implementations of the EKT introduced in Refs. 20, 118 and 119. In parametrically isotropic systems, EKT gives a leading-order accurate description (in ${\alpha }_s$) of the QCD time evolution of the one-particle distribution function and allows for a numerical realization of the so-called bottom-up thermalization scenario.⁸⁷ In practice, we solved the EKT Boltzmann equation for a gluonic plasma undergoing one-dimensional Bjorken expansion with transverse translational symmetry such that the effective Boltzmann equation reads¹²⁰

where $f(p)$ is the gluonic one-particle distribution function (per degree of freedom). The elastic scattering term ${\mathcal{𝒞}}_{2\leftrightarrow 2}$ and the effective inelastic term ${\mathcal{𝒞}}_{1\leftrightarrow 2}$ include physics of dynamical screening and Landau–Pomeranchuk–Migdal suppression and, in order to find the form of the collision kernels, self-energy and ladder resummations are required.^20,118,119

To numerically solve Eq. (73), one can discretize $n(p)=p^{2}f(p)$ on an optimized momentum-space grid and employ Monte Carlo sampling to compute the integrals within the elastic and inelastic collisional kernels. The methodology used is based on the methods outlined in Refs. 20 and 119. The method ensures exact conservation of energy, while accurately addressing particle number violation originating from inelastic contributions to the collisional kernel. Leveraging the azimuthal symmetry inherent in Bjorken flow, one can reduce the momentum-space discretization to a two-dimensional grid. In Ref. 71, we utilized 250 points in the p-direction and 2000 points in the $x=cos\theta$ direction. The momenta p were distributed logarithmically within the ranges $[0.02,45]\mathrm{\Lambda }$, where $\mathrm{\Lambda }$ denotes the typical energy scale of the initial condition. In the results presented in this section we used a ’t Hooft coupling $\lambda =N_c{g}^2=10$, corresponding to a specific shear viscosity of $\bar{\eta }=\eta /s\approx 0.62$.^18,121

We track the temporal evolution of the same set of integral moments considered in RTA⁶³

where $p=|p|$. Notably, the energy density is expressed as $\epsilon =N_{\mathrm{\text{dof}}}{\mathcal{ℳ}}^{20}$, longitudinal pressure as $P_L=N_{\mathrm{\text{dof}}}{\mathcal{ℳ}}^{01}$ and number density as $n=N_{\mathrm{\text{dof}}}{\mathcal{ℳ}}^{10}$. Other moments lack interpretation in terms of conventional hydrodynamic moments often discussed in literature.

We scaled these moments by their corresponding equilibrium values with ${\overline{\mathcal{ℳ}}}^{nl}(\tau )\equiv {\mathcal{ℳ}}^{nl}(\tau )/{\mathcal{ℳ}}_{\mathrm{\text{eq}}}^{nl}(\tau )$. Using a Bose distribution, one obtains

The temperature T here corresponds to the temperature of an equilibrium system with the same energy density, given by $T={(30\epsilon /N_{\mathrm{\text{dof}}}{\pi }^2)}^{1/4}$.

Our simulations were initialized with two types of initial conditions: either spheroidally deformed thermal initial conditions, referred to as “RS” initial conditions,¹⁰¹ or nonthermal color-glass-condensate (CGC)-inspired initial conditions.²⁰ In the former scenario, the initial gluonic one-particle distribution function takes the following form :

where $-1<{\xi }_0<\infty$ encodes the initial momentum anisotropy and ${\mathrm{\Lambda }}_0$ is a temperature-like scale which sets the magnitude of the initial average transverse momentum. In the second case, we take for the form of the initial gluonic one-particle distribution This latter form has been utilized in several prior studies (see, for instance, Refs. 20, 65 and 122–124) and is motivated by the saturation framework. Here, the initial average transverse momentum scale ${\tilde{\mathrm{\Lambda }}}_0$ is connected to the saturation scale, which is ${\tilde{\mathrm{\Lambda }}}_0={〈p_T〉}_0\approx 1.8Q_s$.^125,126,127 The overall constant A is determined by adjusting the initial energy density to match an expected value ${\tau }_0{\epsilon }_0=0.358\nu Q_s^3/\lambda$ from a classical Yang–Mills simulation.¹²⁷

In Fig. 8, I present the evolution of the scaled moments ${\overline{\mathcal{ℳ}}}^{nl}$ with $n\in \{0,1,2\}$ and $m\in \{0,1,2,3\}$. As shown in Fig. 8, evidence of a nonequilibrium attractor emerges both at early and late times across all moments. The integral moments are plotted against a rescaled time variable $\bar{w}=\tau /{\tau }_R(\tau )$, which represents the system’s age in units of the instantaneous interaction time ${\tau }_R(\tau )$. This interaction timescale ${\tau }_R(\tau )=4\pi \bar{\eta }/T(\tau )$ varies in time. The chosen time scaling ensures that, as long as the system is closely described by hydrodynamics near thermal equilibrium, ${\overline{\mathcal{ℳ}}}^{01}$ will eventually be described by the first-order gradient expansion, ${\overline{\mathcal{ℳ}}}^{01}=1-(120\zeta (5)/{\pi }^5){\tau }_R/\tau$.^5,63 This behavior remains unaffected by microscopic details or specific values of macroscopic hydrodynamic parameters. A similar convergence was also observed in RTA and it occurs before the system is effectively characterized by the hydrodynamic gradient expansion.

While the late-time attractor behavior for the longitudinal pressure, ${\overline{\mathcal{ℳ}}}^{01}$, has previously been noted in simplified kinetic theories, e.g., RTA, QCD EKT evolution allows one to examine the extent to which this attractive behavior governs the overall shape of the distribution function. The timescale at which various solutions to Eq. (73) converge to the attractor is approximately $\tau /{\tau }_R\sim 0.5$. While all theories ultimately converge onto a single curve, the timescale for individual solutions to collapse to the attractor depends on specifics of the model. In Ref. 66, two distinct patterns were identified. In RTA kinetic theory and IS hydrodynamics, the transition to the attractor occurs via a power law, with the scale determined by the initial time ${\tau }_0$. Consequently, a unique attractor emerges at arbitrarily early times, discernible through initialization with decreasing ${\tau }_0$. Conversely, in AdS/CFT, the transition to the attractor occurs solely on the timescale dictated by the decay of the quasinormal modes.

4. Two Phenomenological Implications of Attractors

We now highlight two phenomenological implications of attractors. The first is based on the ability to use attractors to extrapolate backwards in time and the second stems from the detailed QCD EKT simulations presented in the previous section.

4.1. Entropy generation constraints

As mentioned above, for a system undergoing Bjorken expansion, conformality, along with energy–momentum conservation ${\partial }_{\mu }{T}^{\mu \nu }=0$, leads to the following evolution equation for the energy density :

where $P_L$ is the longitudinal pressure of the system and $\tau$ is the longitudinal proper time. In first-order NS vHydro, the longitudinal pressure is given by with second-order viscous corrections being of order ${(\bar{\eta }/T\tau )}^2$,⁵³ where $\eta$ is the shear viscosity, s is the entropy density and $\bar{\eta }\equiv \eta /s$ is the specific shear viscosity. The hydrodynamic attractor concept⁴⁵ envisions an all-order resummation of these terms into a constitutive relation of the form $P_L/e=f(w)$, where $w(\tau )\propto \tau T_{\mathrm{\text{eff}}}/\bar{\eta }$.²⁵ The effective temperature $T_{\mathrm{\text{eff}}}$ is defined through the nonequilibrium energy density and the equation of state, $sT=\frac{4}{3}\epsilon =b_{\mathrm{\text{qgp}}}{T}^4$, where $b_{\mathrm{\text{qgp}}}\simeq 17.6$ is estimated from the lattice equation of state at high temperatures.^128,129

In the conventional view, hydrodynamics is applicable only when the Knudsen number, defined as the ratio of microscopic scales like $\ell$ (such as the mean free path of a gas) to macroscopic scales L associated with spatial gradients of conserved quantities, becomes significantly smaller than one. This Knudsen number is given by $\mathrm{\text{Kn}}=\ell /L$. In this context, we have $w\sim 1/\mathrm{\text{Kn}}$, meaning that the hydrodynamic regime becomes valid at later times when $w\gg 1$. Corrections are then organized in powers of $w^{-1}$ (or equivalently, in powers of Kn). If the appropriate attractor constitutive relation is in play, the energy density can be expressed as

where ${(\epsilon {\tau }^{4/3})}_{\infty }\propto {(\tau s)}_{\infty }^{4/3}$ normalizes the entropy in the system and $\mathcal{ℰ}\to 1$ at late times.

Different types of simulations, encompassing the RTA,^{46,61,63,67,74,75,130} QCD kinetic theory^{18,20,64,65,66,71} and the AdS/CFT correspondence,^{11,13,15,16,17,18,48,56,57} suggest that the simplified interpolating expression in Eq. (80) effectively captures the overall dynamics. The behavior of $\mathcal{ℰ}(w)$ is highly constrained by its characteristics at both early and late times. Specifically, at late times, the entropy per rapidity ${(\tau s)}_{\infty }$ remains constant, whereas at early times when $\tau \to 0$, the energy per rapidity ${(\tau \epsilon (\tau ))}_0$ remains constant, dictating the behavior as $\mathcal{ℰ}\simeq C_{\infty }^{-1}{w}^{4/9}$. The constant $C_{\infty }$ exhibits a mild dependence on the underlying theory and is approximately unity. Simulations within QCD kinetic theory suggest $C_{\infty }\simeq 0.87$, while the AdS/CFT correspondence yields $C_{\infty }=1.06$.^11,45

The hydrodynamic attractor, therefore, becomes very useful, allowing one to determine the energy density with an accuracy of 20% during the initial stages, using the measured charged particle multiplicity, nearly independent of the underlying theoretical framework.¹³¹ This precision level serves to significantly constrain models for the initial state, including those based on the CGC.¹³² The second noteworthy significance of the attractor solution is its role as a criterion for the initiation of hydrodynamics. This criterion has undergone validation through detailed studies, particularly highlighting that the hydrodynamic gradient expansion becomes applicable at a time ${\tau }_{\mathrm{\text{hydro}}}$ when $w\simeq 1$ or larger. To translate this theoretical criterion into an experimental context, it is important to recognize that the total entropy per rapidity in hydrodynamics is directly linked to the hydrodynamic multiplicity.⁶⁵

where $b_{\mathrm{\text{hrg}}}\simeq 8.3$ is based on the particle composition of the hadron resonance gas model and the lattice equation of state, with a 15% uncertainty.⁶⁵ The temperature during later times is ascertained through the constant entropy or multiplicity, where $d{N}_{\mathrm{\text{ch}}}/d\eta \propto \tau T^3$, resulting in where we have defined the opacity of the system We note that the multiplicity factor $N_0$ is primarily determined by the properties of the equation of state, which is known from lattice QCD. Consequently the uncertainties in $N_0$ are only of order 20%. The system will have a strong hydrodynamic response if ${\tau }_{\mathrm{\text{hydro}}}/{\tau }_R\ll 1$⁶⁴ which amounts to a requirement that opacity is large $\chi \gg 1$.

4.2. Freeze-out in heavy-ion collisions

Transforming hydrodynamic fields into particle distributions requires a “freeze-out” procedure. While the energy–momentum tensor relies solely on the first momentum-integral moments of the distribution function, particle distributions span all moments. Therefore, translating results from hydrodynamic simulations into particle distributions requires the incorporation of additional information via theoretical assumptions. The standard approach involves assuming a distribution function with a near-equilibrium configuration, wherein deviations from equilibrium stem from formally small corrections. These corrections’ structure is dictated by the linearized collision kernel’s response, within some assumed kinetic theory, to infinitesimal strain.^133,134 Since the freeze-out procedure strongly affects the phenomenological analysis and conclusions about the matter created in ultrarelativistic heavy-ion collisions, it is of great interest to assess how well justified the theoretical assumptions about the shape of the nonequilibrium distribution functions are. The need for such an assessment becomes increasingly important in the case of small systems, e.g., peripheral nucleus–nucleus collisions, proton–nucleus and high-multiplicity proton–proton collisions, where hydrodynamical descriptions are being applied to situations that remain far from equilibrium throughout their dynamical evolution.

Although hydrodynamics does not explicitly describe higher moments of the distribution functions, it is customary to deduce the entire distribution’s shape solely from the shear components of the energy–momentum tensor. For a given $T^{\mu \nu }$, the linearized viscous correction to the one-particle distribution function, $\delta f$, can be locally determined based on some assumptions about the collision kernel. Here, I will present results obtained with two forms for $\delta f$. First, the (i) quadratic ansatz

This form arises from a broad array of models, including RTA with momentum-independent relaxation time, the momentum diffusion approximation, scalar field theory and the EKT in the leading-log approximation.¹³⁴ Here, $\bar{\mathrm{\Pi }}=\mathrm{\Pi }/\epsilon =1/3-T^{zz}/\epsilon$ represents the shear viscous correction to the longitudinal pressure, normalized by the energy density. However, QCD EKT exhibits additional structure; for large momenta $p\gg T$, it reduces to a power law form of the (ii) LPM ansatz This $p^{1.5}$ power law is numerically close to $\propto p^{1.38}$ that is found to describe the full EKT in the relevant momentum range in Ref. 134.

One can also explore a straightforward (iii) aHydro freeze-out ansatz, which does not rely on linearization around equilibrium. Instead, in this approach, the nonequilibrium distribution function is approximated by a spheroidally deformed Bose distribution, $f(p)=f_{\mathrm{\text{Bose}}}(\sqrt{p^2+\xi p_z^2}/\mathrm{\Lambda })$.^26,27,101 To assess this method, we determine $\xi (\tau )$ such that the energy–momentum tensor of the ansatz reproduces the full simulation, enabling predictions for higher-order moments.

In Fig. 9, I compare the different moments obtained from the prescriptions mentioned above with the QCD EKT attractor solution. At late times ($\tau \gtrsim 5{\tau }_R$), all prescriptions describe the low-order moments within a few percent. However, some discrepancy persists even at $\tau \sim 20{\tau }_R$ between the quadratic ansatz (i) and the QCD EKT results. This disparity gradually worsens at earlier times, particularly around $\tau \sim {\tau }_R$, where corrections to the longitudinal pressure become substantial ($P_L/P_L^{\mathrm{\text{eq}}}\sim 65\%$). For instance, ${\mathcal{ℳ}}^{11}$ exhibits an approximately 20% discrepancy between EKT and both linearized ansatz. This disagreement increases for higher moments and at earlier times. Conversely, we observe a remarkably good agreement between the aHydro ansatz and the QCD EKT results at all times.

In the phenomenological analysis of nuclear collisions, a critical step is the freeze-out procedure, wherein hydrodynamic fields are transformed into particle distributions. Presently, the quadratic ansatz (i) is the most commonly used form. However, this approach assumes linear deviations from thermal equilibrium, a stark contrast to the far-from-equilibrium conditions prevalent in current phenomenological applications, especially in modeling small systems (see, e.g., Refs. 40 and 135,136,137,138,139). To investigate whether these linearized procedures remain quantitatively predictive in far-from-equilibrium scenarios, I have compared them with far-from-equilibrium simulations of QCD EKT. The findings shown in Fig. 9 indicate that the nonlinear aHydro freeze-out ansatz performs better in reconstructing moments of the distribution function compared to linearized ansatz in far-from-equilibrium systems.

5. 3+1D aHydro

The results in the first part of this paper demonstrated that the aHydro formulation of dissipative hydrodynamics resums an infinite number of terms in the expansion in the inverse Reynolds number. As a result, aHydro best reproduces the evolution of all moments of the one-particle distribution function and also provides a superior formulation when it comes to freeze-out. I now describe phenomenological application of 3+1D aHydro, allowing for a nonconformal realistic equation of state and full dynamics in both the transverse plane and longitudinal (spatial rapidity) direction.

To model heavy-ion collisions one must obtain the evolution equations for the resulting 3+1D configurations and include the nonconformality of QCD consistent with a realistic lattice-based equation of state. To do this one can use aHydro,^26,27,140 specifically quasiparticle aHydro (aHydroQP) in which one assumes that the QGP consists of massive relativistic quasiparticles with temperature-dependent masses $m(T)$.¹⁴⁰ In aHydroQP the system is assumed to obey a relativistic Boltzmann equation with $m(T)$ determined from lattice QCD computations of QCD thermodynamics. Since the masses are temperature dependent, the Boltzmann equation contains an additional force term on the left-hand side related to gradients in the temperature

With being the collisional kernel in RTA. Above, $u^{\mu }$ is the four-velocity associated with the LRF of the matter and Latin indices such as $i\in \{x,y,z\}$ are spatial indices.

For a gas of massive quasiparticles, the relaxation time is given by¹⁴⁰

where ${\widehat{m}}_{\mathrm{\text{eq}}}=m/T$, $\bar{\eta }=\eta /s$ is the specific shear viscosity, $\epsilon$ is the energy density, P is the pressure, which is fixed by the equation of state. The definition of the special function $I_{3,2}({\widehat{m}}_{\mathrm{\text{eq}}})$ can be found in Refs. 140 and 141. The effective temperature $T(\tau )$ is determined by requiring energy–momentum conservation. For a momentum-and energy-independent relaxation time this requires that the nonequilibrium kinetic energy densities calculated from f are equal to the equilibrium kinetic energy density calculated from the equilibrium distribution, $f_{\mathrm{\text{eq}}}(T,m)$.

In all leading-order aHydro phenomenological works to date, one assumes that the distribution is parametrized by a diagonal anisotropy tensor in the LRF^g

As indicated above, in the LRF, the argument of the distribution function can be expressed in terms of three independent momentum-anisotropy parameters ${\alpha }_i$ and a scale parameter $\lambda$, which are space-time-dependent fields. Herein, we will assume that $f_{\mathrm{\text{eq}}}$ is given by a Boltzmann distribution.

In order to determine the space-time evolution of the fields $u$, $\alpha$ and $\lambda$ one must obtain seven dynamical equations. The first aHydroQP equation of motion is obtained from the first moment of the left-hand side of the quasiparticle Boltzmann equation (88), which reduces to ${\partial }_{\mu }{T}^{\mu \nu }$. In the RTA, however, the first moment of the collisional kernel on the right-hand side results in a constraint that must be satisfied in order to conserve energy and momentum, i.e., $\int dP{p}^{\mu }C[f]=0$. This constraint can be enforced by expressing the effective temperature in terms of the microscopic parameters $\lambda$ and $\alpha$. As a consequence, computing the first moment of the Boltzmann equation gives the partial differential equation resulting from energy–momentum conservation

where For the second equation of motion, one can perform a similar procedure using the second moment of the quasiparticle Boltzmann equation with $I^{\mu \nu \lambda }=\int dP{p}^{\mu }{p}^{\nu }{p}^{\lambda }f$ and $J^{\mu }=\int dP{p}^{\mu }f$.

5.1. The equation of state for aHydroQP

For a system of massive particles obeying Boltzmann statistics, the equilibrium energy density, pressure and entropy density are given by

where ${\widehat{m}}_{\mathrm{\text{eq}}}=m/T$, $K_i$ are modified Bessel functions of the second kind and $\widehat{N}=N_{\mathrm{\text{dof}}}/2{\pi }^2$, with $N_{\mathrm{\text{dof}}}$ being the number of degrees of freedom present in the theory under consideration.

In the quasiparticle approach, one assumes the mass to be temperature dependent. This results in a change in the bulk variables in Eqs. (95). Due to the temperature dependence of the quasiparticle mass, one cannot simply insert $m(T)$ into the bulk variables since this will not be thermodynamically consistent. This is due to the fact that the entropy density may be obtained in two ways: $s_{\mathrm{\text{eq}}}=({\epsilon }_{\mathrm{\text{eq}}}+P_{\mathrm{\text{eq}}})/T$ and $s_{\mathrm{\text{eq}}}=\partial P_{\mathrm{\text{eq}}}/\partial T$. In order to guarantee that both methods result in the same expression for the entropy density, the energy–momentum tensor definition must include a background field $B_{\mathrm{\text{eq}}}$, i.e.,

where $B_{\mathrm{\text{eq}}}$ is the additional background contribution. As a result, in an equilibrium gas of massive quasiparticles, the bulk thermodynamic variables for the gas become

$${\epsilon }_{\mathrm{\text{eq}}}(T,m)={\epsilon }_{\mathrm{\text{kinetic}}}+B_{\mathrm{\text{eq}}},$$(97)

$$P_{\mathrm{\text{eq}}}(T,m)=P_{\mathrm{\text{kinetic}}}-B_{\mathrm{\text{eq}}},$$(98)

$$s_{\mathrm{\text{eq}}}(T,m)=s_{\mathrm{\text{kinetic}}}.$$(99)

To determine the temperature dependence of $B_{\mathrm{\text{eq}}}$ one requires thermodynamic consistency

To fix the equation of state, one can determine $m(T)$ using This can be solved based on the equilibrium energy density and pressure determined from lattice QCD calculations.³⁷ Once $m(T)$ is fixed, the background field $B(T)$ can be determined using the requirement of thermodynamic consistency. The resulting effective mass scaled by T extracted from continuum extrapolated Wuppertal–Budapest lattice data¹²⁸ can be found in Refs. 37 and 143. At high temperatures ($T\sim 0.6$GeV) the scaled mass is $\propto T$ in agreement with the expected high-temperature behavior of QCD.¹⁴³ Once $m(T)$ is fixed, one can self-consistently determine all transport coefficients.

In Fig. 10, I show the result for the bulk viscosity to entropy density ratio. I note that the peak value of $\zeta /s$ is smaller than obtained from some other analyses, but consistent with Bayesian extractions and AdS/CFT-inspired calculations of this ratio. As can be seen from Fig. 10, the peak value obtained using aHydroQP is $\sim 0.05$, whereas in some other fits to data,¹⁴⁴ the peak value can be as large as $\zeta /s\sim 0.5$. Finally, I point out that, in practice, one does not have to encode this in the dynamical equations since it is automatically incorporated. In aHydroQP one does not have this as an independent function to fit to data. In fact, not only this transport coefficient, but an infinite number of nonconformal transport coefficients are self-consistently included by the aHydroQP resummation, without the need for arbitrary parametrizations.

5.2. Leading-order aHydroQP evolution equations

The evolution equations for $u^{\mu }$, $\lambda$ and ${\alpha }_i$ are obtained from moments of the quasiparticle Boltzmann equation. These can be expressed compactly by introducing a timelike vector $u^{\mu }$ which is normalized as $u^{\mu }{u}_{\mu }=1$ and three spacelike vectors $X_i^{\mu }$ which are individually normalized as $X_i^{\mu }{X}_{\mu ,i}=-1$.^28,29 These vectors are mutually orthogonal and obey $u_{\mu }{X}_i^{\mu }=0$ and $X_{\mu ,i}{X}_j^{\mu }=0$ for $i\ne j$. The four equations resulting from the first moment of the Boltzmann equation are

where $D_u\equiv u^{\mu }{\partial }_{\mu }$ and $D_i\equiv X_i^{\mu }{\partial }_{\mu }$. The expansion scalars are ${\theta }_u={\partial }_{\mu }{u}^{\mu }$ and ${\theta }_i={\partial }_{\mu }{X}_i^{\mu }$. Expressions for the basis vectors, derivative operators and expansion scalars appearing above can be found in Refs. 33, 37, 40 and 44. The quantities $\epsilon$ and $P_i$ are the energy density and pressures obtained using the aHydro ansatz for the one-particle distribution function including the background contribution $B(T)$ necessary to enforce thermodynamic consistency

$$\epsilon ={\epsilon }_{\mathrm{\text{kinetic}}}(\lambda ,\alpha ,m)+B(\lambda ,\alpha ),$$(104)

$$P_i=P_{i,\mathrm{\text{kinetic}}}(\lambda ,\alpha ,m)-B(\lambda ,\alpha ).$$(105)

The three equations resulting from the second moment of the Boltzmann equation are

with³³ where $\widehat{m}=m/\lambda$ and $\alpha ={\alpha }_x{\alpha }_y{\alpha }_z$.

Equations (102), (103) and (106) provide seven partial differential equations for $u$, $\alpha$ and $\lambda$, which can be solved numerically. To determine the local effective temperature, one make uses of Landau matching; requiring the equilibrium and nonequilibrium energy densities in the LRF to be equal and solving for T. This system of partial differential equations are evolved until the effective temperature in the entire simulation volume falls below a given freeze-out temperature, $T_{\mathrm{\text{FO}}}$.

Initial transverse and longitudinal profiles

For the results reviewed herein, the initial energy density distribution in the transverse plane is computed from a “tilted” profile.¹⁴⁷ The distribution used is a linear combination of smooth Glauber wounded-nucleon and binary-collision density profiles, with a binary-collision mixing factor of $\chi =0.15$. In the longitudinal direction, we used a profile with a central plateau and Gaussian “tails”, resulting in a longitudinal profile function of the form

We fit $\mathrm{\Delta }\varsigma$ and ${\sigma }_{\varsigma }$ to the pseudorapidity distribution of charged hadron production, with the results being $\mathrm{\Delta }\varsigma =2.3$ and ${\sigma }_{\varsigma }=1.6$^43,44,148 at Large Hadron Collider (LHC) energies and $\mathrm{\Delta }\varsigma =1.4$ and ${\sigma }_{\varsigma }=1.4$ in $\sqrt{s_{\mathrm{\text{NN}}}}=200$GeV collisions.¹⁴⁹

The resulting initial energy density at a given transverse position $x_T$ and spatial rapidity $\varsigma$ was computed using

where $W_{A,B}(x_T)$ are the wounded nucleon densities for nuclei A and B, $C(x_T)$ is the binary-collision density, and $g(\varsigma )$ is where $y_N=log(2\sqrt{s_{\mathrm{\text{NN}}}}/(m_p+m_n))$ is the nucleon momentum rapidity.¹⁴⁷

5.3. Freeze-out and hadronic decays

From the $3+1$D aHydroQP solution one extracts a three-dimensional freeze-out hypersurface at a fixed energy density (temperature). For this purpose, we assume the same fluid anisotropy and scale parameters for all hadronic species. We also assume that all hadrons are created in chemical equilibrium. Employing an extended Cooper–Frye prescription,¹⁴⁰ one translates the underlying hydrodynamic evolution parameters, including flow velocity, anisotropy parameters and scale, into explicit “primordial” hadronic distribution functions on this hypersurface. The aHydroQP parameters on the freeze-out hypersurface are then processed by a modified version of T herminator 2, which generates hadronic configurations using Monte Carlo sampling.¹⁵⁰ Following the sampling of primordial hadrons, hadronic decays are accounted for using standard routines in T herminator 2. The source code for aHydroQP and the custom version of T herminator 2 utilized are publicly accessible.¹⁵¹

5.4. Phenomenological applications of leading-order aHydroQP

Next, I discuss the phenomenological use of leading-order aHydroQP with a diagonal anisotropy tensor. This includes application at $\sqrt{s_{\mathrm{\text{NN}}}}=2.76$TeV, 5.02TeV and 200GeV. I present comparisons with data from both LHC and Relativistic Heavy-Ion Collider (RHIC).

5.4.1. 2.76TeV Pb–Pb collisions

In Ref. 44, 2.76TeV Pb–Pb collisions were considered. I summarize the results found here. In Fig. 11, I show the spectra of pions, kaons and protons as a function of the transverse momentum $p_T$ in four centrality classes 0–5%, 5–10%, 10–20% and 20–30%. As can be seen from these comparisons, the aHydroQP model shows good agreement with the experimental data over the entire $p_T$ range shown, with some discrepancies at high $p_T$ in relatively higher centrality classes as shown in panel (d). From the spectra, one can compute the average transverse momentum for identified hadrons ${〈p_T〉}_i$ and the total charged particle multiplicity $dN/d\eta$. In the left panel of Fig. 12, I show the multiplicity as a function of pseudorapidity. In the right panel, I show the average transverse momentum as a function of centrality for pions, kaons and protons. On the left, I show the multiplicity in five centrality classes (0–5%, 5–10%, 10–20%, 20–30% and 30–40%, from top to bottom), which demonstrates that aHydroQP describes the multiplicity quite well compared to available experimental data. In the right panel, I show ${〈p_T〉}_i$ as a function of centrality. Again, we see that aHydroQP agrees with the data quite well.

In Fig. 13, the resulting identified elliptic flow is displayed as a function of $p_T$ in two centrality classes, 20–30% and 30–40%. As observed in this figure, aHydroQP describes the data well up to $p_T\sim 2$GeV. However, at higher $p_T$, aHydroQP predictions deviate from the data, which is a common occurrence in hydrodynamic calculations since at high transverse momentum hard physics starts to dominate. Next, comparisons of Hanbury–Brown–Twiss (HBT) radii ratios $R_{\mathrm{\text{out}}}/R_{\mathrm{\text{side}}}$, $R_{\mathrm{\text{out}}}/R_{\mathrm{\text{long}}}$ and $R_{\mathrm{\text{side}}}/R_{\mathrm{\text{long}}}$ are presented as a function of the pair relative momentum, $k_T$. The top row of Fig. 14 shows the HBT ratios in the 0–5% centrality class, while the bottom row displays the ratios in the 20–30% centrality class. In both centrality classes, aHydroQP predictions agree well with experimental data. In Ref. 44, we provide HBT ratios in more centrality classes, along with results for $R_{\mathrm{\text{side}}}$, $R_{\mathrm{\text{long}}}$ and $R_{\mathrm{\text{out}}}$, and more comparisons to experimental data.

5.4.2. 5.02TeV Pb–Pb collisions

In this section, I present comparisons between aHydroQP predictions and 5.02TeV Pb–Pb collision data gathered by the ALICE Collaboration. This work appeared originally in Ref. 148. I focus on two key free parameters: the initial central temperature $T_0$ and the specific shear viscosity $\eta /s$. These parameters are determined by fitting to the transverse momentum spectra of pions, kaons and protons in the 0–5% and 30–40% centrality classes. The obtained parameters from the spectra fit were: $T_0=630$MeV and $\eta /s=0.159$. I note that the initial temperature derived is only slightly higher than that found at 2.76TeV, which was $T_0^{2.76\mathrm{\text{TeV}}}=600$MeV.⁴⁴ The best-fit value for $\eta /s$ matches that found at 2.76TeV.⁴⁴

I begin by presenting aHydroQP’s predictions for the transverse momentum distribution of identified hadrons in 5.02TeV Pb–Pb collisions. These predictions are compared with experimental data from the ALICE Collaboration.^157,158 In Fig. 15, I present the combined spectra of pions, kaons and protons as a function of transverse momentum across six distinct centrality classes. At low centralities, aHydroQP demonstrates remarkable consistency with the data, as depicted in Fig. 15(a). Conversely, for more peripheral collisions, satisfactory agreement is observed only for $p_T\lesssim 1$GeV.

In Fig. 16, the $K/\pi$ ratio (left column) and $p/\pi$ ratio (right column) are presented as a function of $p_T$ in three distinct centrality classes and compared to experimental data. Notably, the agreement between aHydroQP and the data at $p_T\lesssim 1$GeV is very good across all centrality bins. Specifically, in the 0–5% centrality class, the agreement between aHydroQP and the data for the $K/\pi$ ratio persists up to $p_T\sim 2.5$GeV, while for $p/\pi$, it extends up to $p_T\sim 1.5$GeV. The aHydroQP predictions for the integrated $K/\pi$ and $p/\pi$ ratios as a function of centrality are shown in the left and right panels of Fig. 17, respectively. In both panels, we see that aHydroQP describes the centrality dependence of the integrated multiplicity ratios quite well across a wide range of centralities.

5.4.3. 200GeV Au–Au collisions

Next, I review results obtained using aHydroQP in 200GeV Au–Au collisions.¹⁴⁹ Drawing from our prior investigation of 2.76 TeV collisions at LHC,⁴⁴ we adopted a fixed switching (freeze-out) temperature of $T_{\mathrm{\text{FO}}}=130$MeV. This choice leaves only the shear viscosity to entropy density ratio $\eta /s$ and the initial central temperature $T_0$ (the center of the system for a $b=0$ collision) as independent parameters. Similarly to before, we assume that $\bar{\eta }$ is independent of the temperature. To determine $T_0$ and $\bar{\eta }$, we conducted comparisons between model predictions and observed pion, proton and kaon spectra in the 0–5% and 30–40% centrality classes. These comparisons yielded $T_0=455$MeV at ${\tau }_0=0.25$fm/c and $\bar{\eta }=0.179$. The resulting fits to the pion, kaon and proton spectra are depicted in Fig. 18 and compared with experimental data from the PHENIX Collaboration.¹⁵⁹ As demonstrated by this figure, the model effectively describes the identified particle spectra with this parameter set. However, in high centrality classes, we observe that the model underestimates hadron production at large transverse momenta, $p_T\gtrsim 1.5$GeV.

In Fig. 19, I show our findings for the identified particle multiplicities as a function of centrality. From top to bottom, the particles depicted are ${\pi }^{+}$, $K^{+}$, p, $\varphi$ and ${\mathrm{\Omega }}^{+}+{\mathrm{\Omega }}^{-}$. The data for ${\pi }^{+}$, $K^{+}$ and p are from the PHENIX Collaboration,¹⁵⁹ while the data for the $\varphi$ meson are from the PHENIX Collaboration.¹⁶⁰ The ${\mathrm{\Omega }}^{+}+{\mathrm{\Omega }}^{-}$ data are from the STAR Collaboration.¹⁶¹ The aHydroQP model results are binned according to the centrality bins used by the PHENIX Collaboration for ${\pi }^{+}$, $K^{+}$ and p. As illustrated by this figure, aHydroQP, coupled with our customized version of Therminator 2, effectively replicates the centrality dependence of observed identified particle multiplicities. This is particularly noteworthy given that we utilized a single iso-thermal switching (freeze-out) temperature, which is relatively low ($T_{\mathrm{\text{FO}}}=130$MeV), yet are able to reasonably reproduce observed identified particle multiplicities not only for central collisions but across many centrality classes.

Finally, in Fig. 20, I present the aHydroQP model results for the charged particle multiplicity as a function of pseudorapidity in 200GeV Au–Au collisions. In this figure, I compare to the data reported by the PHOBOS Collaboration.¹⁶² In the figure, panel (a) shows centrality classes in the range 0–25% and panel (b) shows centrality classes in the range 25–50%. As can be seen from this figure, once the initial central temperature and shear viscosity to entropy density ratio are fit to the identified spectra, aHydroQP is able to describe the centrality and pseudorapidity dependence of the charged particle multiplicity quite well.

5.4.4. 5.02 and 8.16TeV p-Pb collisions

Next, I discuss recent results obtained for p-Pb collisions at 5.02TeV and 8.16TeV in Ref. 167, which used the same setup as the RHIC and LHC results presented above. Here, I focus on 5.02TeV and results for 8.16TeV can be found in Ref. 167. In this case we, again, use a tilted profile and fit the central width $\mathrm{\Delta }\varsigma$ to the pseudorapidity distribution of charged hadron production. This gives $\mathrm{\Delta }\varsigma =1.8$ for $\sqrt{s_{\mathrm{\text{NN}}}}=5.02$TeV collisions and $\mathrm{\Delta }\varsigma =2.1$ for $\sqrt{s_{\mathrm{\text{NN}}}}=8.16$TeV collisions. The width of the Gaussian tails, ${\sigma }_{\varsigma }$, is largely unconstrained based on available p-Pb data. As a result, we have used ${\sigma }_{\varsigma }=1.6$, which was used previously in Pb–Pb collisions at both $\sqrt{s_{\mathrm{\text{NN}}}}=5.02$TeV and 2.76TeV.^43,44,148

The resulting initial energy density at a given transverse position $x_{\perp }$ and spatial rapidity $\varsigma$ was computed using¹³⁵

where $W_A(x_{\perp })$ is the wounded-nucleon density for nucleus A,¹⁶⁸$C(x_{\perp })$ is the binary-collision density,¹⁶⁸ and $g(\varsigma )$ is the tilt function introduced earlier. To compute $W_p$, we parametrize the proton overlap function as^169,170 with $n=1.85$ and $r_p=0.975$fm.

To evaluate the tuning of aHydroQP in p-Pb collisions, Figs. 21 and 22 present comparisons of our aHydroQP p-Pb results with standard soft hadron observables measured at LHC energies. Figure 21 presents our aHydroQP results for the charged particle multiplicity as a function of pseudorapidity, $\eta \equiv \frac{1}{2}ln((p+p_z)/(p-p_z))$, at $\sqrt{s_{\mathrm{\text{NN}}}}=$ 5.02 and $\sqrt{s_{\mathrm{\text{NN}}}}=8.16$TeV. In this figure, comparisons are made with data from the ALICE^163,165 and CMS¹⁶⁶ Collaborations. The top panel displays the model results alongside ALICE and CMS data, while the bottom panel shows the relative error (data/model) of the aHydroQP results. As evident from Fig. 21, the aHydroQP model can reproduce the charged particle multiplicity measurements of the ALICE and CMS Collaborations to within approximately 5% in the reported pseudorapidity range. This level of agreement is comparable to or better than all model comparisons presented in Refs. 163, 165 and 166.

In Fig. 22, the computed transverse momentum spectra in p-Pb collisions at $\sqrt{s_{\mathrm{\text{NN}}}}=5.02$TeV for pions, kaons and protons are compared to data from the ALICE¹⁶⁵ and CMS Collaborations.¹⁶⁶ In the left panel, comparisons of the identified spectra obtained using aHydroQP with experimental data are shown, while the three right panels show the relative error (data/model). As shown in this figure, aHydroQP is capable of reproducing the experimental observations for $p_T\lesssim 1.5$GeV reasonably well. The remaining discrepancies are comparable to those encountered with other models presented in Refs. 165 and 166, thereby instilling confidence in the ability of aHydroQP to provide a reasonable description of min-bias p-Pb collisions. For the details concerning the p-Pb results presented in this section, see Ref. 167.

6. Phenomenological Results from Viscous aHydro

Finally, I would like to highlight recent phenomenological results that make use of viscous aHydro (vaHydro or VAH). VAH is a model that goes beyond leading-order aHydro by including nonspheroidal corrections in a linearized manner akin to second-order vHydro.^{31,35,171,172,173,174} In a recent paper¹⁷⁴ a comprehensive Bayesian analysis using the JETSCAPE framework was performed, replacing the free streaming and hydrodynamical model components by the VAH model. Prior to showing the main results of this study, I quickly sketch the VAH framework. For details, I refer the reader to Refs. 31, 35 and 171–174.

To setup VAH, one starts by decomposing the spatial projector ${\mathrm{\Delta }}^{\mu \nu }$ into projectors along the beam direction ($z^{\mu }$) and a transverse projector ${\mathrm{\Xi }}^{\mu \nu }$ such that ${\mathrm{\Delta }}^{\mu \nu }={\mathrm{\Xi }}^{\mu \nu }-z^{\mu }{z}^{\nu }$. Multiplying these projectors by distinct longitudinal and transverse pressures, the energy–momentum tensor can then decomposed as³⁸

This decomposition allows VAH to evolve the longitudinal and transverse pressures $P_L$ and $P_T$ separately, treating them at the same level as the thermal pressure in standard vHydro. Importantly, it is not by assumed that their differences from the thermal pressure and from one another are small.

The evolution equations for the energy density and the flow velocity are obtained from the conservation laws for energy and momentum

The equilibrium pressure $P(\epsilon )$ is taken from lattice QCD calculations conducted by the HotQCD Collaboration.¹⁷⁵ The dynamical evolution equations for the nonequilibrium flows $P_L$, $P_T$, $W_{Tz}^{(\mu }{z}^{\nu )}$ and ${\pi }_T^{\mu \nu }$ are derived under the assumption that the fluid’s microscopic physics can be described by the relativistic Boltzmann equation with a medium-dependent mass,^40,172,176 as was the case with aHydroQP. The relaxation times for $P_L$ and $P_T$ are expressed in terms of those for the bulk and shear viscous pressures,¹⁷² which are parametrized as

The functions ${\beta }_{\pi }$, ${\beta }_{\mathrm{\Pi }}$, along with all necessary anisotropic transport coefficients, are computed within the quasiparticle kinetic theory model in Ref. 172.

For the temperature-dependent specific shear and bulk viscosities, $(\eta /s)(T)$ and $(\zeta /s)(T)$, the authors of Ref. 174 used the same parametrizations as the JETSCAPE Collaboration,¹⁷⁷ namely

with and with

In addition to the eight parameters related to the viscosities that are included in the above parametrizations, there was one additional parameter inferred from the experimental data, namely the initial pressure ratio $R={(P_L/P_T)}_0$ at the time ${\tau }_0$ that VAH was initialized

Above $P_0\equiv P({\epsilon }_0)$ is the equilibrium pressure at the initial proper time. The authors assumed that the initial bulk viscous pressure was zero. The initial flow profile was assumed to be given by Bjorken flow in the longitudinal direction, with zero initial transverse velocities, and the initial residual shear stresses $W_{Tz}^{\mu }$ and ${\pi }_T^{\mu \nu }$ were taken to be zero. At freeze-out, as done in aHydroQP, the authors sampled from an anisotropic momentum-space distribution that is nonnegative for all momentum.

The VAH model replaces the free-streaming and relativistic vHydro stages by a single VAH module, coupled to its own anisotropic freeze-out/particlization routine. This eliminates the free-streaming time parameter (and an associated discontinuity in the equation of state¹⁷⁷) and extends the sensitivity of the model to the QGP viscosities into the far-from-equilibrium stage. This allows VAH to be applied at much higher temperatures than in other available hydrodynamical evolution models.

In Fig. 23, I show results from Refs. 174 and 180 for the $p_T$-spectra of identified hadrons, $\pi$ (red), K (green), p (blue), predicted by the best-fit parameters of the calibrated JETSCAPE SIMS model with Chapman–Enskog (CE) particlization^178,177 (top row) and the calibrated VAH model¹⁷⁴ (bottom row), for 2.76TeV Pb–Pb collisions at three collision centralities (from left to right: 0–5%, 5–10%, 10–20%¹⁷⁹). As can be seen from this figure the VAH best-fit results agree better with the experimental data than the calibrated JETSCAPE SIMS result. This is quantified by the model/data ratios which are closer to unity when using the calibrated VAH model than when using the calibrated JETSCAPE SIMS result.

In Fig. 24, I show the VAH results^174,180 for the $p_T$-differential elliptic $v_2\{\mathrm{\text{SP}}\}$ (left), triangular $v_3\{\mathrm{\text{SP}}\}$ (middle) and quadrangular $v_4\{\mathrm{\text{SP}}\}$ (right) flows for identified pions, kaons and protons. Results predicted by the best-fit parameters of the calibrated JETSCAPE SIMS model with CE particlization^177,178 are indicated by dotted lines and the calibrated VAH model¹⁷⁴ by solid lines. Results are for 2.76TeV Pb–Pb collisions at 20–30% centrality. The experimental data are from the ALICE Collaboration.^181,182 As this figure demonstrates, the VAH predictions exhibit better agreement with the data compared to the JETSCAPE SIMS predictions, across all three particlization models examined in Refs. 177 and 178 and at all available collision centralities.

Finally, from the Bayesian analysis performed in Ref. 174 it is possible to extract confidence intervals on the extracted transport coefficients. In Fig. 25, I show the VAH results^174,180 for the prior (gray) and posterior (colored) credibility intervals for the temperature-dependent specific bulk (left) and shear (right) viscosities inferred from the VAH model¹⁷⁴ (top) and from the JETSCAPE SIMS model¹⁷⁸ (bottom). In the subpanels, the authors plotted the Kullback–Leibler divergence $D_{\mathrm{\text{KL}}}$, which quantifies the information gain provided by the experimental data. From this figure, we see that $D_{\mathrm{\text{KL}}}$ is nonzero up to temperatures of around 700MeV in the VAH model and, at low temperatures, the VAH model provides tighter constraints than JETSCAPE SIMS. I note, importantly, that the temperature dependence of the bulk viscosity to entropy density ratio extracted using VAH (see the 60% confidence interval) is consistent with the result obtained using leading-order nonlinearized aHydroQP, in which it was found that the bulk viscosity to entropy density ratio was below 5% at all temperatures (see Fig. 10).

7. Conclusions and Outlook

In this paper, I have made a connection between our understanding of nonequilibrium attractors and the phenomenological predictions of aHydro. The motivation for aHydro comes from the understanding that a high degree of momentum-space anisotropy exists at early times in the QGP’s lifetime and also near the transverse edges. This can be seen easily from the attractor since small $w=\tau T$ corresponds either to small $\tau$ (early times) or low temperature regions (transverse and longitudinal edges).

I demonstrated that aHydro provides the most accurate reproduction of exact solutions to the Boltzmann equation for 0+1D systems and the associated attractor for the one-particle distribution function. Although not discussed in detail here, this also holds true for nontrivial conformal flows like Gubser flow^183,184 and nonconformal systems.³³ I also discussed tests of different freeze-out formulations using QCD EKT simulations in which it was found that, like in RTA, an attractor for the full one-particle distribution function exists and its form is better described by the anisotropic aHydro ansatz that linearized vHydro ansatz. This is particularly important for peripheral collisions and small systems, where the lifetime of the QGP may be short and the system may not have time enough to approach a near-equilibrium configuration for which linearized approaches are applicable.

Taking this together, I then reviewed various phenomenological applications of aHydroQP, including to Pb–Pb collisions at 2.76TeV and 5.02TeV, Au–Au collisions at 200GeV and finally p-Pb collisions at 5.02TeV and 8.16TeV. In all cases, one finds good agreement between the aHydroQP results and experimental data for spectra, multiplicities, elliptic flow, HBT radii, etc. I closed by highlighting some very recent results using the VAH framework, which relies on single anisotropy parameter (spheroidal deformation) plus linearized corrections in the spirit of second-order vHydro. The results shown indicate that the calibrated VAH did a better job at reproducing experimental observations than the calibrated JETSCAPE SIMS framework. In addition, the VAH model allows one to dispense with the free streaming phase and, therefore, access early times when the temperature is large. This allows VAH to, in principle, constrain QGP transport coefficients at much higher temperatures that standard simulation chains.

Acknowledgments

I thank my collaborators for their contributions to the work reported herein. This work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics (Nuclear Theory) under contract number DE-SC0013470.

Notes

^a It is worth mentioning that studies utilizing either the 2PI formalism or holography suggest that, at the highest temperatures explored in heavy-ion collisions, an equation of state may be established well before pressure isotropization.^22,24,50

^b I emphasize that the DNMR results are the correct values resulting from a complete treatment, while the MIS equations discard certain terms. Despite this, MIS is frequently used as a standard reference point for discussions.

^c As demonstrated in Eq. (25), aHydro naturally reproduces this equation when truncated at leading order in $\xi$ (linear order in the inverse Reynolds number).

^d For a recent analytic understanding of this, see Ref. 107.

^e It is possible to investigate the emergence of attractors with underlying Fermi–Dirac or Bose statistics using the exact solution presented in Ref. 108.

^f In Eq. (55), z is the spatial coordinate, which is not to be confused with the basis vector $z^{\mu }$.

^g In Ref. 142, a method for nonperturbatively including off-diagonal components of the anisotropy tensor was presented, however, to date it has not been applied to phenomenology.