Recent development of hydrodynamic modeling in heavy-ion collisions

Recent development of hydrodynamic modeling in heavy-ion collisions 1 2 corresp first-author chunshen@wayne.edu chunshen@wayne.edu chunshen@wayne.edu https://orcid.org/0000-0002-6677-4784 3 corresp cliyan@fudan.edu.cn cliyan@fudan.edu.cn cliyan@fudan.edu.cn https://orcid.org/0000-0002-7764-0693 Department of Physics and Astronomy, Wayne State University, Detroit, MI 48201, USA Department of Physics and Astronomy, Wayne State University, Detroit, MI 48201, USA RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA Key Laboratory of Nuclear Physics and Ion-Beam Application (MOE) & Institute of Modern Physics Fudan University, 220 Handan Road, Shanghai 200433, China Key Laboratory of Nuclear Physics and Ion-Beam Application (MOE) & Institute of Modern Physics Fudan University, 220 Handan Road, Shanghai 200433, China 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 collisions hydrodynamics 1 Introduction 1 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 30 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 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 2 - 4 4 2-4 ]. Theoretical model calculations using relativistic viscous hydrodynamics successfully characterize these flow observables, which makes relativistic hydrodynamics the "standard mo del" in heavy-ion collisions [ 5 5 - 10 10 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 πk B , suggesting that QGP is a strongly coupled medium [ 11 11 , 12 12 ]. The specific bulk viscosity has been extracted as well, leading to a temperature dependence [ 13 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 14 , 15 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 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 17 - 23 23 17-23 ], and novel hydrodynamic modes due to the effect of critical slowing down can be included [ 24 24 - 27 27 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 28 ], the application of hydrodynamic modeling to small systems such as those created in proton–nucleus collisions [ 29 29 - 31 31 29-31 ] is not straightforward, owing to the significant reduction in system size and high expansion rate [ 32 32 , 33 33 ]. The "unreasonable effectiveness" [ 10 10 , 28 28 , 34 34 , 35 35 ] in describing the collectivity in the small colliding systems has qualitatively modified the understanding of QCD system thermalization [ 36 36 ]. The condition of the onset hydrodynamics ( hydrodynamization ) is quite relaxed. The traditionally recognized hydrodynamic and nonhydrodynamic modes and propagation of these modes [ 37 37 , 38 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 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 40 - 42 42 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 <math id="M1"><mrow><msqrt><mi>s</mi></msqrt><mo>∼</mo><mi mathvariant="script">O</mi><mo stretchy="false">(</mo><mn>10</mn><mo stretchy="false">)</mo></mrow></math> 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 1 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 2 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 7 ]. 1 1 The general form of the energy–momentum tensor T μν is given in the corresponding constitutive relation as <math display="block" id="M2"><mrow><msup><mi>T</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>=</mo><mi>ϵ</mi><msup><mi>U</mi><mi>μ</mi></msup><msup><mi>U</mi><mi>ν</mi></msup><mo>−</mo><mo stretchy="false">(</mo><mi mathvariant="script">P</mi><mo>+</mo><mi>Π</mi><mo stretchy="false">)</mo><msup><mi>Δ</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>+</mo><msup><mi>π</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mtext> </mtext><mn>.</mn></mrow></math> The projection tensor is defined as <math display="block" id="M3"><mrow><msup><mi>Δ</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>=</mo><msup><mi>g</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>−</mo><msup><mi>U</mi><mi>μ</mi></msup><msup><mi>U</mi><mi>ν</mi></msup><mtext> </mtext><mn>,</mn></mrow></math> 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, <math display="block" id="M4"><mrow><msubsup><mi>T</mi><mrow><mtext>ideal</mtext></mrow><mrow><mi>μ</mi><mi>ν</mi></mrow></msubsup><mo>=</mo><mi>ϵ</mi><msup><mi>U</mi><mi>μ</mi></msup><msup><mi>U</mi><mi>ν</mi></msup><mo>−</mo><mi mathvariant="script">P</mi><msup><mi>Δ</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mtext> </mtext><mn>,</mn></mrow></math> 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, <math display="block" id="M5"><mtable columnalign="left"><mtr><mtd><msup><mi>π</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>=</mo><mn>2</mn><msup><mi>η</mi><mo>〈</mo></msup><msup><mo>∇</mo><mi>μ</mi></msup><msup><mi>U</mi><mrow><mi>ν</mi><mo>〉</mo></mrow></msup><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mo>∇</mo><mn>2</mn></msup><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mtd></mtr><mtr><mtd><mi>Π</mi><mo>=</mo><mo>−</mo><mi>ζ</mi><mo>∇</mo><mo>⋅</mo><mi>U</mi><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mo>∇</mo><mn>2</mn></msup><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mtd></mtr></mtable></math> 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 32 ], <math display="block" id="M6"><mrow><mtext>Kn</mtext><mo>∼</mo><mfrac><mrow><mo>|</mo><msup><mi>η</mi><mo>〈</mo></msup><msup><mo>∇</mo><mi>μ</mi></msup><msup><mi>U</mi><mrow><mi>ν</mi><mo>〉</mo></mrow></msup><mo>|</mo></mrow><mi mathvariant="script">P</mi></mfrac><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>o</mi><mi>r</mi><mtext> </mtext><mtext> </mtext><mtext> </mtext><mfrac><mrow><mo>|</mo><mi>ζ</mi><mo>∇</mo><mo>⋅</mo><mi>U</mi><mo>|</mo></mrow><mi mathvariant="script">P</mi></mfrac><mtext> </mtext><mn>.</mn></mrow></math> In Eq. (4) and in what follows, the brackets around tensor indices indicate a symmetric, transverse, and traceless projection of a tensor, i.e. , <math display="block" id="M7"><mrow><msup><mrow/><mo>〈</mo></msup><msup><mi>A</mi><mrow><mi>μ</mi><mi>ν</mi><mo>〉</mo></mrow></msup><mo>=</mo><msup><mi>Δ</mi><mrow><mi>μ</mi><mi>ν</mi><mi>α</mi><mi>β</mi></mrow></msup><msub><mi>A</mi><mrow><mi>α</mi><mi>β</mi></mrow></msub><mtext> </mtext><mn>,</mn></mrow></math> where <math display="block" id="M8"><mrow><msup><mi>Δ</mi><mrow><mi>μ</mi><mi>ν</mi><mi>α</mi><mi>β</mi></mrow></msup><mo>≡</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>[</mo><mrow><msup><mi>Δ</mi><mrow><mi>μ</mi><mi>α</mi></mrow></msup><msup><mi>Δ</mi><mrow><mi>ν</mi><mi>β</mi></mrow></msup><mo>+</mo><msup><mi>Δ</mi><mrow><mi>μ</mi><mi>β</mi></mrow></msup><msup><mi>Δ</mi><mrow><mi>ν</mi><mi>α</mi></mrow></msup></mrow><mo>]</mo></mrow><mo>−</mo><mfrac><mn>1</mn><mn>3</mn></mfrac><msup><mi>Δ</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><msup><mi>Δ</mi><mrow><mi>α</mi><mi>β</mi></mrow></msup><mtext> </mtext><mn>.</mn></mrow></math> 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 2 For practical simulations, acausal modes can be remedied by using the Israel–Stewart formulation [ 44 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 45 ], <math display="block" id="M9"><mtable columnalign="left"><mtr><mtd><msup><mi>π</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>=</mo><mtext> </mtext><msup><mi>η</mi><mo>〈</mo></msup><msup><mo>∇</mo><mi>μ</mi></msup><msup><mi>U</mi><mrow><mi>ν</mi><mo>〉</mo></mrow></msup><mo>−</mo><msub><mi>τ</mi><mi>π</mi></msub><mrow><mo>[</mo><mrow><msup><mrow/><mo>〈</mo></msup><mi>D</mi><msup><mi>π</mi><mrow><mi>μ</mi><mi>ν</mi><mo>〉</mo></mrow></msup><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><msup><mi>π</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>∇</mo><mo>⋅</mo><mi>U</mi></mrow><mo>]</mo></mrow></mtd></mtr><mtr><mtd><mo>−</mo><mfrac><mrow><msub><mi>λ</mi><mn>1</mn></msub></mrow><mrow><msup><mi>η</mi><mn>2</mn></msup></mrow></mfrac><msubsup><mi>π</mi><mrow><mtext> </mtext><mi>α</mi></mrow><mrow><mo>〈</mo><mi>μ</mi></mrow></msubsup><msup><mi>π</mi><mrow><mi>ν</mi><mo>〉</mo><mi>α</mi></mrow></msup><mo>−</mo><mfrac><mrow><msub><mi>λ</mi><mn>2</mn></msub></mrow><mi>η</mi></mfrac><msubsup><mi>π</mi><mrow><mtext> </mtext><mi>α</mi></mrow><mrow><mo>〈</mo><mi>μ</mi></mrow></msubsup><msup><mi>Ω</mi><mrow><mi>ν</mi><mo>〉</mo><mi>α</mi></mrow></msup><mo>−</mo><msub><mi>λ</mi><mn>3</mn></msub><msubsup><mi>Ω</mi><mrow><mtext> </mtext><mi>α</mi></mrow><mrow><mo>〈</mo><mi>μ</mi></mrow></msubsup><msup><mi>Ω</mi><mrow><mi>ν</mi><mo>〉</mo><mi>α</mi></mrow></msup><mtext> </mtext><mn>,</mn></mtd></mtr></mtable></math> where, in addition to the shear relaxation time τ π , λ 1 , λ 2 , and λ 3 are independent second-order transport coefficients. For conformal fluids, these transport coefficients are known [ 12 12 , 45 45 ] and can be parameterized as [ 46 46 , 47 47 ] <math display="block" id="M10"><mrow><mi>η</mi><mo>=</mo><msub><mi>C</mi><mi>η</mi></msub><mi>s</mi><mn>,</mn><mtext> </mtext><msub><mi>τ</mi><mi>π</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>C</mi><mi>τ</mi></msub><msub><mi>C</mi><mi>η</mi></msub></mrow><mi>T</mi></mfrac><mn>,</mn><mtext> </mtext><msub><mi>λ</mi><mn>1</mn></msub><mo>=</mo><msub><mi>C</mi><mrow><msub><mi>λ</mi><mn>1</mn></msub></mrow></msub><mfrac><mi>s</mi><mi>T</mi></mfrac><mtext> </mtext><mn>,</mn></mrow></math> where the local entropy density s ∝ T 3 . 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 12 , 48 48 ], 3 3 <math display="block" id="M11"><mrow><msub><mi>C</mi><mi>η</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mn>4</mn><mi>π</mi></mrow></mfrac><mtext> </mtext><mn>,</mn><mtext> </mtext><msub><mi>C</mi><mi>τ</mi></msub><mo>=</mo><mn>5</mn><mtext> </mtext><mn>,</mn><mtext> </mtext><msub><mi>C</mi><mrow><msub><mi>λ</mi><mn>1</mn></msub></mrow></msub><mo>=</mo><mfrac><mn>5</mn><mn>7</mn></mfrac><msub><mi>C</mi><mi>η</mi></msub><msub><mi>C</mi><mi>τ</mi></msub><mtext> </mtext><mn>.</mn></mrow></math> There exist other variant forms of second-order viscous hydrodynamics in addition to Eq. (8), when conformal symmetry is not guaranteed [ 49 49 ]. Note that Eq. (8) is consistent with the Müller–Israel–Stewart theory [ 44 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 50 , 51 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 52 - 54 54 52-54 ]. 4 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, <math display="block" id="M12"><mrow><msup><mi>J</mi><mi>μ</mi></msup><mo>≡</mo><msubsup><mi>J</mi><mrow><mtext>i</mtext><mi>d</mi><mi>e</mi><mi>a</mi><mi>l</mi></mrow><mi>μ</mi></msubsup><mo>+</mo><msup><mi>I</mi><mi>μ</mi></msup><mo>=</mo><mi>n</mi><msup><mi>U</mi><mi>μ</mi></msup><mo>+</mo><msup><mi>I</mi><mi>μ</mi></msup><mtext> </mtext><mn>,</mn></mrow></math> with <math display="block" id="M13"><mrow><msup><mi>I</mi><mi>μ</mi></msup><mo>=</mo><mi>σ</mi><mi>T</mi><msup><mo>∇</mo><mi>μ</mi></msup><mrow><mo>(</mo><mrow><mfrac><mi>μ</mi><mi>T</mi></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mo>∇</mo><mn>2</mn></msup><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mrow></math> 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 57 ], <math display="block" id="M14"><mrow><msup><mi>I</mi><mi>μ</mi></msup><mo>=</mo><mi>σ</mi><mi>T</mi><msup><mo>∇</mo><mi>μ</mi></msup><msup><mrow><mo stretchy="false">(</mo><mn>1</mn><mo>+</mo><msub><mi>τ</mi><mi>Q</mi></msub><mi>D</mi><mo stretchy="false">)</mo></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup><mrow><mo>(</mo><mrow><mfrac><mi>μ</mi><mi>T</mi></mfrac></mrow><mo>)</mo></mrow><mn>.</mn></mrow></math> 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 2 Hydrodynamic modeling has been successfully applied to small colliding systems. For the high-multiplicity events of proton–lead [ 29 29 , 58 58 - 60 60 58-60 ], 3 He–gold, deuteron–gold [ 61 61 - 63 63 61-63 ], and even proton–proton collisions [ 30 30 ], the observed multiparticle correlations were found to be compatible with hydrodynamic predictions [ 64 64 , 65 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 9 , 46 46 , 66 66 - 81 81 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 3 ) ≪ O(Kn 2 ). In this way, gradient corrections of order higher than or equal to O (Kn 3 ) can be safely neglected. In heavy-ion collisions, this is by assumption satisfied at a time scale τ 0 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 τ 0 ≪ τ 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, τ 0 ∼1/ T [ 82 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 Q s such that <math id="M15"><mrow><msub><mi>τ</mi><mn>0</mn></msub><mo>∼</mo><msubsup><mi>α</mi><mi>s</mi><mi>x</mi></msubsup><msubsup><mi>Q</mi><mi>s</mi><mrow><mo>−</mo><mn>1</mn></mrow></msubsup></mrow></math> [ 83 83 ], where the exponent x is a constant negative number. However, given the actual values of the strong coupling constant and Q s in heavy-ion collisions, the expected τ 0 in realistic QGP systems could be rather large [ 83 83 ]. We concentrate on the weakly coupled system in the present discussion. Despite all the estimates of τ 0 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 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 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 85 ]. As the system expands and the density of the gluon field decreases, around τ ∼1/ Q s , 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 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/ Q s , 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 τ 0 , the second-order viscous hydrodynamics starts to dominate, and 𝒫 L /𝒫 T approaches unity, that is, close to local thermal equilibrium. The existence of such a universal evolution as illustrated in Fig. 1 Fig. 1 has been proved by various theoretical analyses. This uniquely implies the value of τ 0 , 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 33 , 39 39 , 46 46 , 47 47 , 67 67 , 68 68 , 71 71 , 74 74 , 81 81 , 86 86 - 89 89 86-89 ] and Gubser flow [ 69 69 , 79 79 , 80 80 , 90 90 , 91 91 ], and numerically in less symmetric systems [ 70 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 3 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 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 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, <math display="block" id="M16"><mrow><mi>τ</mi><mo>=</mo><msqrt><mrow><msup><mi>t</mi><mn>2</mn></msup><mo>−</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></msqrt><mtext> </mtext><mn>,</mn><mtext> </mtext><mtext> </mtext><mi>ξ</mi><mo>=</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mi>tanh</mi></mrow></mstyle><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo stretchy="false">(</mo><mi>t</mi><mo>/</mo><mi>z</mi><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mrow></math> 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 <math display="block" id="M17"><mrow><msup><mi>U</mi><mi>τ</mi></msup><mo>=</mo><mn>1</mn><mtext> </mtext><mn>,</mn><mtext> </mtext><mtext> </mtext><msup><mi>U</mi><mi>x</mi></msup><mo>=</mo><msup><mi>U</mi><mi>y</mi></msup><mo>=</mo><msup><mi>U</mi><mi>ξ</mi></msup><mo>=</mo><mn>0</mn><mtext> </mtext><mn>.</mn></mrow></math> Therefore, considering that the metric is rewritten as <math id="M18"><mrow><mi>d</mi><msup><mi>s</mi><mn>2</mn></msup><mo>=</mo><msub><mi>g</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msub><mi>d</mi><msup><mi>x</mi><mi>μ</mi></msup><mi>d</mi><msup><mi>x</mi><mi>ν</mi></msup><mo>=</mo><mi>d</mi><msup><mi>τ</mi><mn>2</mn></msup><mo>−</mo><msup><mi>τ</mi><mn>2</mn></msup><mi>d</mi><msup><mi>ξ</mi><mn>2</mn></msup><mo>−</mo><mi>d</mi><msubsup><mi>x</mi><mo>⊥</mo><mn>2</mn></msubsup></mrow></math> , gradients of the fluid fields are reduced, and all are related to the proper time τ ; for example, <math display="block" id="M19"><mrow><mo>∇</mo><mo>⋅</mo><mi>U</mi><mo>=</mo><mo>−</mo><mfrac><mn>3</mn><mn>4</mn></mfrac><msubsup><mi>σ</mi><mrow><mtext> </mtext><mtext> </mtext><mi>ξ</mi></mrow><mi>ξ</mi></msubsup><mo>=</mo><mfrac><mn>1</mn><mi>τ</mi></mfrac><mtext> </mtext><mn>.</mn></mrow></math> 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 <math display="block" id="M20"><mrow><msub><mo>∂</mo><mi>τ</mi></msub><mi>ϵ</mi><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>ϵ</mi><mi>τ</mi></mfrac><mo>+</mo><mfrac><mrow><msubsup><mi>π</mi><mrow><mtext> </mtext><mtext> </mtext><mi>ξ</mi></mrow><mi>ξ</mi></msubsup></mrow><mi>τ</mi></mfrac><mtext> </mtext><mn>.</mn></mrow></math> 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, <math display="block" id="M21"><mrow><mi>ϵ</mi><mo>∼</mo><msup><mi>τ</mi><mrow><mo>−</mo><mn>4</mn><mo>/</mo><mn>3</mn></mrow></msup><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mi>n</mi></munder></mstyle><mi>O</mi><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>τ</mi><mi>π</mi></msub></mrow><mi>τ</mi></mfrac></mrow><mo>)</mo></mrow></mrow><mi>n</mi></msup></mrow><mo>)</mo></mrow><mtext> </mtext><mn>.</mn></mrow></math> One may check that in the conformal fluid with respect to Bjorke symmetry, the only non-zero components of the shear stress tensor are <math id="M22"><mrow><msubsup><mi>π</mi><mi>ξ</mi><mi>ξ</mi></msubsup><mo>=</mo><mo>−</mo><mn>2</mn><msubsup><mi>π</mi><mi>x</mi><mi>x</mi></msubsup><mo>=</mo><mo>−</mo><mn>2</mn><msubsup><mi>π</mi><mi>y</mi><mi>y</mi></msubsup></mrow></math> . Similarly, the vorticity tensor Ω μν vanishes. By defining <math id="M23"><mrow><mi>π</mi><mo>=</mo><msubsup><mi>π</mi><mi>ξ</mi><mi>ξ</mi></msubsup></mrow></math> , and given all the information in the Bjorken flow, the constitutive relation for π follows from the BRSSS theory, Eq. (8), as <math display="block" id="M24"><mrow><mi>π</mi><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>η</mi><mi>τ</mi></mfrac><mo>−</mo><msub><mi>τ</mi><mi>π</mi></msub><mrow><mo>[</mo><mrow><mfrac><mrow><mtext>d</mtext><mi>π</mi></mrow><mrow><mtext>d</mtext><mi>τ</mi></mrow></mfrac><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>π</mi><mi>τ</mi></mfrac></mrow><mo>]</mo></mrow><mo>+</mo><mfrac><mrow><msub><mi>λ</mi><mn>1</mn></msub></mrow><mrow><mn>2</mn><msup><mi>η</mi><mn>2</mn></msup></mrow></mfrac><msup><mi>π</mi><mn>2</mn></msup><mtext> </mtext><mn>,</mn></mrow></math> 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 49 ] <math display="block" id="M25"><mrow><mi>π</mi><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>η</mi><mi>τ</mi></mfrac><mo>−</mo><msub><mi>τ</mi><mi>π</mi></msub><mrow><mo>[</mo><mrow><mfrac><mrow><mtext>d</mtext><mi>π</mi></mrow><mrow><mtext>d</mtext><mi>τ</mi></mrow></mfrac><mo>+</mo><msub><mi>β</mi><mi>π</mi></msub><mfrac><mi>π</mi><mi>τ</mi></mfrac></mrow><mo>]</mo></mrow><mo>−</mo><mfrac><mrow><mi>χ</mi><msubsup><mi>τ</mi><mi>π</mi><mn>2</mn></msubsup></mrow><mi>η</mi></mfrac><mfrac><mrow><msup><mi>π</mi><mn>2</mn></msup></mrow><mi>τ</mi></mfrac><mtext> </mtext><mn>,</mn></mrow></math> 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 88 , 92 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 <math display="block" id="M26"><mrow><mi>g</mi><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>≡</mo><mfrac><mrow><mtext>d</mtext><mtext> </mtext><mtext> </mtext><mtext>ln</mtext><mtext> </mtext><mi>ϵ</mi></mrow><mrow><mtext>d</mtext><mtext> </mtext><mtext> </mtext><mtext>ln</mtext><mtext> </mtext><mi>τ</mi></mrow></mfrac><mo>=</mo><mo>−</mo><mn>1</mn><mo>−</mo><mfrac><mrow><msub><mi>P</mi><mi>L</mi></msub></mrow><mi>ϵ</mi></mfrac><mtext> </mtext><mn>,</mn></mrow></math> as a function of the inverse Knudsen number <math display="block" id="M27"><mrow><mi>w</mi><mo>≡</mo><mi>τ</mi><mo>/</mo><msub><mi>τ</mi><mi>π</mi></msub><mtext> </mtext><mn>.</mn></mrow></math> Apparently, g ( w ) characterizes the decay rate of the energy density and is related to the pressure ratio. Although in the conformal case, d w /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, <math display="block" id="M28"><mtable columnalign="left"><mtr><mtd><mi>w</mi><mfrac><mrow><mtext>d</mtext><mi>g</mi></mrow><mrow><mtext>d</mtext><mi>w</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mi>g</mi><mn>4</mn></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><msup><mrow><mo>(</mo><mrow><mi>g</mi><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow><mn>2</mn></msup><mrow><mo>[</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><mn>3</mn><mi>w</mi></mrow><mn>8</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>λ</mi></msub></mrow><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mi>w</mi><mrow><mo>(</mo><mrow><mi>g</mi><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mrow><mn>16</mn></mrow><mn>9</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>τ</mi></msub></mrow></mfrac><mo>=</mo><mn>0</mn><mtext> </mtext><mn>,</mn></mtd></mtr></mtable></math> 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 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. As shown in Fig. 2 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 81 ]. This attractor solution in hydrodynamics was first noticed in the context of the Bjorken flow [ 39 39 , 46 46 ], and later in the Gubser flow (cf. Ref [ 92 92 ]). It is also realized in the solution of the kinetic equation for weakly coupled media and strongly coupled media using AdS/CFT [ 70 70 , 72 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, <math display="block" id="M29"><mrow><msup><mrow><mrow><mo>(</mo><mrow><mi>g</mi><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>−</mo><mfrac><mrow><mn>16</mn></mrow><mn>9</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>τ</mi></msub></mrow></mfrac><mo>=</mo><mn>0</mn><mtext> </mtext><mn>,</mn></mrow></math> with its two solutions given as ( g _ - < g + ) <math display="block" id="M30"><mrow><msub><mi>g</mi><mo>±</mo></msub><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>∓</mo><msqrt><mrow><mfrac><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>τ</mi></msub></mrow></mfrac></mrow></msqrt></mrow><mo>)</mo></mrow><mtext> </mtext><mn>.</mn></mrow></math> 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 33 , 74 74 ], <math display="block" id="M31"><mrow><mi>δ</mi><mi>g</mi><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>∼</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mi>τ</mi><mrow><msub><mi>τ</mi><mrow><mtext>i</mtext><mi>n</mi><mi>i</mi></mrow></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mo>−</mo><mo stretchy="false">(</mo><msub><mi>g</mi><mo>+</mo></msub><mo>−</mo><msub><mi>g</mi><mo>−</mo></msub><mo stretchy="false">)</mo></mrow></msup><mtext> </mtext><mn>.</mn></mrow></math> This power law decay qualitatively explains the observed pattern of g ( w ) evolution at early times in Fig. 2 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 →+∞, <math display="block" id="M32"><mrow><mi>g</mi><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>→</mo><msub><mi>g</mi><mrow><mtext>hyd</mtext></mrow></msub><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mtext> </mtext><mn>,</mn></mrow></math> 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 <math display="block" id="M33"><mtable columnalign="left"><mtr><mtd><msub><mi>β</mi><mrow><mtext>BRSSS</mtext></mrow></msub><mo stretchy="false">(</mo><mi>w</mi><mn>,</mn><mi>g</mi><mo stretchy="false">)</mo><mo>≡</mo><msup><mrow><mo>(</mo><mrow><mi>g</mi><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow><mn>2</mn></msup><mrow><mo>[</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><mn>3</mn><mi>w</mi></mrow><mn>8</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>λ</mi></msub></mrow><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mi>w</mi><mrow><mo>(</mo><mrow><mi>g</mi><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mrow><mn>16</mn></mrow><mn>9</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>τ</mi></msub></mrow></mfrac><mtext> </mtext><mn>.</mn></mtd></mtr></mtable></math> The root of Eq. (29) encodes the information of fixed points in Eq. (26) and in Eq. (28). Figure 3 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 <math id="M34"><mrow><mi>g</mi><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><mn>2</mn><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mrow><msub><mi>λ</mi><mn>1</mn></msub></mrow></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mn>2.08</mn></mrow></math> , which again, is a root of β BRSSS =0 in the limit of w →+∞. 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 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 46 , 71 71 , 93 93 , 94 94 ], in which one neglects the derivative in Eq. (24) for the lowest-order estimate: dg / dw ∼0, <math display="block" id="M35"><mtable columnalign="left"><mtr><mtd><msub><mi>g</mi><mrow><mtext>slow-roll</mtext></mrow></msub><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mo>+</mo><mfrac><mn>1</mn><mrow><mn>2</mn><mo>+</mo><mfrac><mrow><mn>3</mn><mi>w</mi></mrow><mn>4</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>λ</mi></msub></mrow><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow></mfrac></mrow></mfrac><mo stretchy="false">[</mo><mo>−</mo><mi>w</mi></mtd></mtr><mtr><mtd><mo>+</mo><msqrt><mrow><msup><mi>w</mi><mn>2</mn></msup><mo>+</mo><mfrac><mrow><mn>16</mn></mrow><mn>9</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>τ</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mn>4</mn><mo>+</mo><mfrac><mrow><mn>3</mn><mi>w</mi></mrow><mn>2</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>λ</mi></msub></mrow><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></msqrt><mo stretchy="false">]</mo></mtd></mtr></mtable></math> 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 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 , <math display="block" id="M36"><mrow><msub><mi>g</mi><mrow><mtext>hydro</mtext><mo stretchy="false">(</mo><mtext>w</mtext><mo stretchy="false">)</mo></mrow></msub><mo>≡</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>0</mn></mrow><mi>∞</mi></munderover></mstyle><mfrac><mrow><msub><mi>f</mi><mi>n</mi></msub></mrow><mrow><msup><mi>w</mi><mi>n</mi></msup></mrow></mfrac><mtext> </mtext><mn>.</mn></mrow></math> Substituting the expansion into Eq. (24), we find the recursion relation for the expansion coefficients 5 5 , <math display="block" id="M37"><mtable columnalign="left"><mtr><mtd><mo>−</mo><mrow><mo>(</mo><mrow><mi>n</mi><mo>−</mo><mfrac><mn>8</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow><msub><mi>f</mi><mi>n</mi></msub><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mi>n</mi></munderover></mstyle><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mfrac><mi>k</mi><mn>4</mn></mfrac></mrow><mo>)</mo></mrow><msub><mi>f</mi><mi>k</mi></msub><msub><mi>f</mi><mrow><mi>n</mi><mo>−</mo><mi>k</mi></mrow></msub><mn>.</mn><mo>+</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><msub><mi>C</mi><mi>λ</mi></msub></mrow><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><msub><mi>f</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msub></mtd></mtr><mtr><mtd><mo>+</mo><mfrac><mn>3</mn><mn>8</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>λ</mi></msub></mrow><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></munderover></mstyle><msub><mi>f</mi><mi>k</mi></msub><msub><mi>f</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn><mo>−</mo><mi>k</mi></mrow></msub><mo>+</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mfrac><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>τ</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><mn>16</mn></mrow><mn>9</mn></mfrac><msub><mi>δ</mi><mrow><mi>n</mi><mn>,0</mn></mrow></msub></mtd></mtr><mtr><mtd><mo>+</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><mrow><mo>(</mo><mrow><mn>2</mn><mo>+</mo><mfrac><mrow><msub><mi>C</mi><mi>λ</mi></msub></mrow><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><msub><mi>δ</mi><mrow><mi>n</mi><mn>,</mn><mo>−</mo><mn>1</mn></mrow></msub><mo>=</mo><mn>0</mn><mtext> </mtext><mn>.</mn></mtd></mtr></mtable></math> The leading-order solutions of the above equation are <math display="block" id="M38"><mrow><msub><mi>f</mi><mn>0</mn></msub><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mtext> </mtext><mtext>and</mtext><mtext> </mtext><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><mn>2</mn><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>λ</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mtext> </mtext><mn>,</mn></mrow></math> 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 0 =-4/3, using Eq. (32) iteratively, one is able to obtain the expansion coefficients to arbitrary order. For instance, one finds <math id="M39"><mrow><msub><mi>f</mi><mn>1</mn></msub><mo>=</mo><mfrac><mrow><mn>16</mn></mrow><mn>9</mn></mfrac><mfrac><mrow><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>τ</mi></msub></mrow></mfrac></mrow></math> . In fact, it can be shown that the magnitudes of these coefficients exhibit factorial growth, f n ∼ n !. This factorial growth can be recognized by noticing the ratio <math display="block" id="M40"><mrow><mfrac><mrow><msub><mi>f</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow><mrow><msub><mi>f</mi><mi>n</mi></msub></mrow></mfrac><mo>→</mo><msup><mi>S</mi><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo stretchy="false">(</mo><mi>n</mi><mo>+</mo><mi>β</mi><mo stretchy="false">)</mo><mo>+</mo><mi>O</mi><mrow><mo>(</mo><mrow><mfrac><mn>1</mn><mi>n</mi></mfrac></mrow><mo>)</mo></mrow><mtext> </mtext><mn>,</mn></mrow></math> for a large n in Eq. (32). The parameters S and β are real constants. With respect to Eq. (24), they are determined as <math display="block" id="M41"><mrow><mi>S</mi><mo>=</mo><mfrac><mn>3</mn><mn>2</mn></mfrac><mtext> </mtext><mn>,</mn><mtext> </mtext><mtext> </mtext><msup><mi>β</mi><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>=</mo><mo>−</mo><mfrac><mrow><mn>2</mn><msub><mi>C</mi><mi>η</mi></msub></mrow><mrow><msub><mi>C</mi><mi>τ</mi></msub></mrow></mfrac><mtext> </mtext><mn>.</mn></mrow></math> Therefore, the leading-order contribution to the coefficients at very large n is f n ∼Γ( 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 96 - 98 98 96-98 ] and the WKB approximation in quantum mechanics [ 99 99 , 100 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 <math display="block" id="M42"><mrow><mi>B</mi><mo stretchy="false">[</mo><msub><mi>g</mi><mrow><mtext>hydro</mtext></mrow></msub><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>≡</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>0</mn></mrow><mi>∞</mi></munderover></mstyle><mfrac><mrow><msub><mi>f</mi><mi>n</mi></msub></mrow><mrow><mi>n</mi><mo>!</mo></mrow></mfrac><msup><mi>z</mi><mi>n</mi></msup><mtext> </mtext><mn>.</mn></mrow></math> 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 <math display="block" id="M43"><mtable columnalign="left"><mtr><mtd><msub><mover accent="true"><mi>g</mi><mo>˜</mo></mover><mrow><mtext>hydro</mtext></mrow></msub><mo stretchy="false">(</mo><mtext>w</mtext><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>0</mn></mrow><mi>∞</mi></munderover></mstyle><mfrac><mrow><msub><mi>f</mi><mi>n</mi></msub></mrow><mrow><msup><mi>w</mi><mi>n</mi></msup></mrow></mfrac><mo>×</mo><mfrac><mn>1</mn><mrow><mi>n</mi><mo>!</mo></mrow></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>∞</mi></msubsup></mrow></mstyle><mi>d</mi><mi>z</mi><msup><mi>e</mi><mrow><mo>−</mo><mi>z</mi></mrow></msup><msup><mi>z</mi><mi>n</mi></msup></mtd></mtr><mtr><mtd><mo>=</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>∞</mi></msubsup></mrow></mstyle><mi>d</mi><mi>z</mi><msup><mi>e</mi><mrow><mo>−</mo><mi>z</mi></mrow></msup><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>0</mn></mrow><mi>∞</mi></munderover></mstyle><mfrac><mrow><msub><mi>f</mi><mi>n</mi></msub></mrow><mrow><mi>n</mi><mo>!</mo></mrow></mfrac><msup><mrow><mo>(</mo><mrow><mfrac><mi>z</mi><mi>w</mi></mfrac></mrow><mo>)</mo></mrow><mi>n</mi></msup></mtd></mtr><mtr><mtd><mo>=</mo><mi>w</mi><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>∞</mi></msubsup></mrow></mstyle><mi>d</mi><mi>z</mi><msup><mi>e</mi><mrow><mo>−</mo><mi>z</mi><mi>w</mi></mrow></msup><mi>B</mi><mo stretchy="false">[</mo><msub><mi>g</mi><mrow><mtext>hydro</mtext></mrow></msub><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mtext> </mtext><mn>.</mn></mtd></mtr></mtable></math> One may check that <math id="M44"><mrow><msub><mover accent="true"><mi>g</mi><mo>˜</mo></mover><mrow><mtext>hydro</mtext></mrow></msub><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo></mrow></math> 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 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 0 , as indicated by the accumulated poles in Fig. 4 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, <math id="M146"> <mrow> <mi mathvariant="script">B</mi><mo stretchy="false">[</mo><msub> <mi>g</mi> <mrow> <mtext>h</mtext><mi>y</mi><mi>d</mi><mi>r</mi><mi>o</mi></mrow> </msub> <mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>∼</mo><msup> <mrow> <mo stretchy="false">(</mo><mi>z</mi><mo>−</mo><mi>S</mi><mo stretchy="false">)</mo></mrow> <mrow> <mo>−</mo><mi>β</mi></mrow> </msup> </mrow> </math> , with S and β real constants given in Eq. (35), the complex ambiguity leads to <math display="block" id="M45"><mrow><mi>I</mi><mtext>m</mtext><mo stretchy="false">[</mo><msub><mover accent="true"><mi>g</mi><mo>˜</mo></mover><mrow><mtext>hydro</mtext></mrow></msub><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>∼</mo><mo>±</mo><mi>π</mi><msup><mi>e</mi><mrow><mo>−</mo><mi>S</mi><mo>/</mo><mi>w</mi></mrow></msup><msup><mi>w</mi><mi>β</mi></msup><mtext> </mtext><mn>.</mn></mrow></math> 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 <math display="block" id="M46"><mrow><mi>δ</mi><mi>g</mi><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>∼</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>S</mi><mi>w</mi></mrow></msup><msup><mi>w</mi><mi>β</mi></msup><mn>.</mn></mrow></math> 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, <math display="block" id="M47"><mrow><mi>g</mi><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>0</mn></mrow></munder></mstyle><msup><mrow><mo stretchy="false">(</mo><mi>σ</mi><msup><mi>e</mi><mrow><mo>−</mo><mi>S</mi><mi>w</mi></mrow></msup><msup><mi>w</mi><mi>β</mi></msup><mo stretchy="false">)</mo></mrow><mi>n</mi></msup><msup><mi>g</mi><mrow><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo></mrow></msup><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mrow></math> where the leading order gives g (0) ( 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 <math id="M48"><mrow><msup><mi>g</mi><mrow><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo></mrow></msup><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><msub><mo>∑</mo><mi>k</mi></msub></mstyle><msubsup><mi>f</mi><mi>k</mi><mrow><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo></mrow></msubsup><mo>/</mo><msup><mi>w</mi><mi>k</mi></msup></mrow></math> 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 101 , 102 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 wi th 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 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 103 ]. 2.2.2 From kinetic theory 3 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 95 ]. To formulate a gradient series that is compatible with the off-equilibrium system evolution, one has to solve the full transport equation [ 66 66 , 75 75 , 77 77 , 104 104 ]. With respect to Bjorken symmetry, the general form of the Boltzmann equation [ 105 105 ], <math display="block" id="M49"><mrow><msup><mi>p</mi><mi>μ</mi></msup><msub><mo>∂</mo><mi>μ</mi></msub><mi>f</mi><mo>+</mo><msubsup><mi>Γ</mi><mrow><mi>μ</mi><mi>ν</mi></mrow><mi>α</mi></msubsup><msup><mi>p</mi><mi>μ</mi></msup><msup><mi>p</mi><mi>ν</mi></msup><mfrac><mo>∂</mo><mrow><mo>∂</mo><msup><mi>p</mi><mi>α</mi></msup></mrow></mfrac><mi>f</mi><mo>=</mo><mover accent="true"><mi>C</mi><mo>^</mo></mover><mo stretchy="false">[</mo><mi>f</mi><mo stretchy="false">]</mo><mtext> </mtext><mn>,</mn></mrow></math> 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 <math display="block" id="M50"><mrow><msup><mi>p</mi><mi>μ</mi></msup><msub><mo>∂</mo><mi>μ</mi></msub><mi>f</mi><mo>+</mo><msubsup><mi>Γ</mi><mrow><mi>μ</mi><mi>ν</mi></mrow><mi>α</mi></msubsup><msup><mi>p</mi><mi>μ</mi></msup><msup><mi>p</mi><mi>ν</mi></msup><mfrac><mo>∂</mo><mrow><mo>∂</mo><msup><mi>p</mi><mi>α</mi></msup></mrow></mfrac><mi>f</mi><mtext> </mtext><mo>→</mo><mtext> </mtext><msup><mi>p</mi><mn>0</mn></msup><mrow><mo>(</mo><mrow><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>τ</mi></mrow></mfrac><mo>−</mo><mfrac><mrow><msub><mi>p</mi><mi>z</mi></msub></mrow><mi>τ</mi></mfrac></mrow><mo>)</mo></mrow><mi>f</mi><mtext> </mtext><mn>,</mn></mrow></math> 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, <math display="block" id="M51"><mrow><mrow><mo>(</mo><mrow><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>τ</mi></mrow></mfrac><mo>−</mo><mfrac><mrow><msub><mi>p</mi><mi>z</mi></msub></mrow><mi>τ</mi></mfrac></mrow><mo>)</mo></mrow><mi>f</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mi>p</mi><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><mfrac><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mi>p</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>f</mi><mrow><mtext>eq</mtext></mrow></msub><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mi>p</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>τ</mi><mi>R</mi></msub></mrow></mfrac><mtext> </mtext><mn>.</mn></mrow></math> There exists a formal and analytical solution to Eq. (43). For a relaxation time with arbitrary τ dependence, τ R ( τ ), the formal solution is [ 106 106 , 107 107 ], <math display="block" id="M52"><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mi>p</mi><mo stretchy="false">)</mo><mo>=</mo><mi>D</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><msub><mi>τ</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mi>f</mi><mo stretchy="false">(</mo><msub><mi>τ</mi><mn>0</mn></msub><mn>,</mn><msub><mi>p</mi><mo>⊥</mo></msub><mn>,</mn><msub><mi>p</mi><mi>z</mi></msub><mi>τ</mi><mo>/</mo><msub><mi>τ</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>τ</mi><mn>0</mn></msub></mrow><mi>τ</mi></msubsup></mrow></mstyle><mfrac><mrow><mi>d</mi><msup><mi>τ</mi><mo>′</mo></msup></mrow><mrow><msub><mi>τ</mi><mi>R</mi></msub><mo stretchy="false">(</mo><msup><mi>τ</mi><mo>′</mo></msup><mo stretchy="false">)</mo></mrow></mfrac><mi>D</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><msup><mi>τ</mi><mo>′</mo></msup><mo stretchy="false">)</mo><msub><mi>f</mi><mrow><mtext>e</mtext><mi>q</mi></mrow></msub><mo stretchy="false">(</mo><msqrt><mrow><msubsup><mi>p</mi><mo>⊥</mo><mn>2</mn></msubsup><mo>+</mo><msup><mrow><mo stretchy="false">(</mo><msub><mi>p</mi><mi>z</mi></msub><msup><mi>τ</mi><mo>′</mo></msup><mo>/</mo><msub><mi>τ</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></msqrt><mo>/</mo><mi>T</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mo>′</mo></msup><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext> </mtext><mn>.</mn></mrow></math> In this solution, f eq is the equilibrium distribution as a function of temperature, which is fixed via the Landau matching condition, ϵ = ϵ eq ∝ T 4 . 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 (τ 0 , p ⊥ , p z τ / τ 0 ), and the time evolution function <math display="block" id="M53"><mrow><mi>D</mi><mo stretchy="false">(</mo><msub><mi>τ</mi><mn>2</mn></msub><mn>,</mn><msub><mi>τ</mi><mn>1</mn></msub><mo stretchy="false">)</mo><mo>≡</mo><mi>exp</mi><mrow><mo>[</mo><mrow><mo>−</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>τ</mi><mn>1</mn></msub></mrow><mrow><msub><mi>τ</mi><mn>2</mn></msub></mrow></msubsup></mrow></mstyle><mi>d</mi><msup><mi>τ</mi><mo>′</mo></msup><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>R</mi></msub><mo stretchy="false">(</mo><msup><mi>τ</mi><mo>′</mo></msup><mo stretchy="false">)</mo></mrow></mfrac></mrow><mo>]</mo></mrow><mtext> </mtext><mn>.</mn></mrow></math> 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 <math display="block" id="M54"><mrow><msup><mi>T</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mfrac><mrow><msup><mi>d</mi><mn>3</mn></msup><mi>p</mi></mrow><mrow><msub><mi>E</mi><mi>p</mi></msub></mrow></mfrac><msup><mi>p</mi><mi>μ</mi></msup><msup><mi>p</mi><mi>ν</mi></msup><msub><mi>f</mi><mi>p</mi></msub><mtext> </mtext><mn>.</mn></mrow></math> In the case of Bjorken symmetry, the independent components of the energy–momentum tensor are the diagonal components, which include the local energy density, <math display="block" id="M55"><mrow><mi>ϵ</mi><mo>=</mo><msup><mi>T</mi><mrow><mn>00</mn></mrow></msup><mo>≡</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><msup><mi>d</mi><mn>3</mn></msup><mi>p</mi><msub><mi>E</mi><mi>p</mi></msub><mi>f</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mi>p</mi><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mrow></math> the longitudinal pressure <math display="block" id="M56"><mrow><msub><mi mathvariant="script">P</mi><mi>L</mi></msub><mo>=</mo><msup><mi>T</mi><mrow><mi>z</mi><mi>z</mi></mrow></msup><mo>≡</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mfrac><mrow><msup><mi>d</mi><mn>3</mn></msup><mi>p</mi></mrow><mrow><msub><mi>E</mi><mi>p</mi></msub></mrow></mfrac><msubsup><mi>p</mi><mi>z</mi><mn>2</mn></msubsup><mi>f</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mi>p</mi><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mrow></math> and the transverse pressure, <math display="block" id="M57"><mrow><msub><mi mathvariant="script">P</mi><mi>T</mi></msub><mo>=</mo><msup><mi>T</mi><mrow><mi>x</mi><mi>x</mi></mrow></msup><mo>=</mo><msup><mi>T</mi><mrow><mi>y</mi><mi>y</mi></mrow></msup><mo>≡</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mfrac><mrow><msup><mi>d</mi><mn>3</mn></msup><mi>p</mi></mrow><mrow><msub><mi>E</mi><mi>p</mi></msub></mrow></mfrac><msubsup><mi>p</mi><mo>⊥</mo><mn>2</mn></msubsup><mi>f</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mi>p</mi><mo stretchy="false">)</mo><mtext> </mtext><mn>.</mn></mrow></math> 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 <math display="block" id="M58"><mrow><mi>τ</mi><mfrac><mrow><mi>d</mi><mi>ϵ</mi></mrow><mrow><mi>d</mi><mi>τ</mi></mrow></mfrac><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mi>ϵ</mi><mo>+</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><mo stretchy="false">(</mo><msub><mi mathvariant="script">P</mi><mi>L</mi></msub><mo>−</mo><msub><mi mathvariant="script">P</mi><mi>T</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mn>0</mn><mtext> </mtext><mn>.</mn></mrow></math> 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 89 ] <math display="block" id="M59"><mrow><msub><mi mathvariant="script">P</mi><mi>L</mi></msub><mo>−</mo><msub><mi mathvariant="script">P</mi><mi>T</mi></msub><mo>=</mo><mo>−</mo><mfrac><mrow><mn>2</mn><mi>η</mi></mrow><mi>τ</mi></mfrac><mo>+</mo><mfrac><mn>4</mn><mrow><mn>3</mn><msup><mi>τ</mi><mn>3</mn></msup></mrow></mfrac><mo stretchy="false">(</mo><msub><mi>λ</mi><mn>1</mn></msub><mo>−</mo><mi>η</mi><msub><mi>τ</mi><mi>π</mi></msub><mo stretchy="false">)</mo><mo>+</mo><mi>O</mi><mrow><mo>(</mo><mrow><mfrac><mn>1</mn><mrow><msup><mi>τ</mi><mn>3</mn></msup></mrow></mfrac></mrow><mo>)</mo></mrow><mtext> </mtext><mn>.</mn></mrow></math> The <math id="M147"><mi mathvariant="script">L</mi></math> -moments The energy–momentum tensor T μν belongs to a specific set of moments of the phase-space distribution. We define the <math id="M148"><mi mathvariant="script">L</mi></math> moment as <math display="block" id="M60"><mrow><msub><mi>L</mi><mi>n</mi></msub><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><msup><mi>d</mi><mn>3</mn></msup><mi>p</mi><msub><mi>E</mi><mi>p</mi></msub><msub><mi>P</mi><mrow><mn>2</mn><mi>n</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>p</mi><mi>z</mi></msub><mo>/</mo><msub><mi>E</mi><mi>p</mi></msub><mo stretchy="false">)</mo><mi>f</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mi>p</mi><mo stretchy="false">)</mo><mn>,</mn></mrow></math> using the Legendre polynomial P n ( 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 E p is chosen such that the <math id="M149"><mi mathvariant="script">L</mi></math> moments are of the same dimension as the energy–momentum tensor. Indeed, it is straightforward to verify that <math display="block" id="M61"><mtable columnalign="left"><mtr><mtd><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mo>=</mo><msup><mi>T</mi><mrow><mn>00</mn></mrow></msup><mo>=</mo><mi>ϵ</mi><mo>=</mo><msub><mi mathvariant="script">P</mi><mi>L</mi></msub><mo>+</mo><mn>2</mn><msub><mi mathvariant="script">P</mi><mi>T</mi></msub><mtext> </mtext><mn>,</mn></mtd></mtr><mtr><mtd><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mo>=</mo><msup><mi>T</mi><mrow><mo>⊥</mo><mo>⊥</mo></mrow></msup><mo>−</mo><msup><mi>T</mi><mrow><mi>z</mi><mi>z</mi></mrow></msup><mo>=</mo><msub><mi mathvariant="script">P</mi><mi>L</mi></msub><mo>−</mo><msub><mi mathvariant="script">P</mi><mi>T</mi></msub><mtext> </mtext><mn>,</mn></mtd></mtr></mtable></math> and hence, the pressure anisotropy can be expressed in terms of the <math id="M150"><mi mathvariant="script">L</mi></math> -moments as <math display="block" id="M62"><mrow><mfrac><mrow><msub><mi mathvariant="script">P</mi><mi>L</mi></msub></mrow><mrow><msub><mi mathvariant="script">P</mi><mi>T</mi></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mo>+</mo><mn>2</mn><msub><mi mathvariant="script">L</mi><mn>1</mn></msub></mrow><mrow><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mo>−</mo><msub><mi mathvariant="script">L</mi><mn>1</mn></msub></mrow></mfrac><mtext> </mtext><mn>.</mn></mrow></math> 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 E p dependence. For instance, the E p dependence in f ( τ , p ) can be decomposed by large-error polynomials [ 108 108 ]. Although the generalized moments of the distribution function can be introduced (cf. Ref [ 49 49 , 75 75 , 109 109 ]), the <math id="M151"><mi mathvariant="script">L</mi></math> -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, <math id="M152"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> and <math id="M153"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> </mrow> </math> fully determine the components in the energy–momentum tensor T μν . With higher-order <math id="M154"><mi mathvariant="script">L</mi></math> -moments included, the description of system evolution can be even improved. An illustration is presented in Fig. 5 Fig. 5 . It is also interesting to note that the coupled equations for <math id="M158"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> and <math id="M159"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> </mrow> </math> are analogous to the so-called anisotropic hydrodynamics (ahydro) [ 110 110 , 111 111 ], where the pressure difference is considered as an individual field for out-of-equilibrium fluids. In a similar manner, higher-order <math id="M160"><mi mathvariant="script">L</mi></math> -moments play the role of viscous corrections, in comparison to viscous anisotropic hydrodynamics (vahydro) [ 112 112 ]. With respect to the analytical formal solution of the distribution function, the analytical solution of <math id="M161"><mi mathvariant="script">L</mi></math> -moments can be obtained: <math display="block" id="M63"><mrow><msub><mi mathvariant="script">L</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo><mo>=</mo><mi>D</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><msub><mi>τ</mi><mn>0</mn></msub><mo stretchy="false">)</mo><msubsup><mi mathvariant="script">L</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></msubsup><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>τ</mi><mn>0</mn></msub></mrow><mi>τ</mi></msubsup></mrow></mstyle><mfrac><mrow><mi>d</mi><msup><mi>τ</mi><mo>′</mo></msup></mrow><mrow><msub><mi>τ</mi><mi>R</mi></msub><mo stretchy="false">(</mo><msup><mi>τ</mi><mo>′</mo></msup><mo stretchy="false">)</mo></mrow></mfrac><mi>D</mi><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><msup><mi>τ</mi><mo>′</mo></msup><mo stretchy="false">)</mo><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mo stretchy="false">(</mo><msup><mi>τ</mi><mo>′</mo></msup><mo stretchy="false">)</mo><mo stretchy="false">(</mo><msup><mi>τ</mi><mo>′</mo></msup><mo>/</mo><mi>τ</mi><mo stretchy="false">)</mo><msub><mi mathvariant="script">F</mi><mi>n</mi></msub><mo stretchy="false">(</mo><msup><mi>τ</mi><mo>′</mo></msup><mo stretchy="false">)</mo></mrow></math> where the function <math id="M162"> <mrow> <msub> <mi mathvariant="script">F</mi> <mi>n</mi> </msub> <mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo></mrow> </math> is defined as <math display="block" id="M64"><mtable columnalign="left"><mtr><mtd><msub><mi mathvariant="script">F</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>≡</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>1</mn></msubsup></mrow></mstyle><mi>d</mi><mi>y</mi><msup><mrow><mo>[</mo><mrow><mn>1</mn><mo>−</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><msup><mi>x</mi><mn>2</mn></msup><mo stretchy="false">)</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mo>]</mo></mrow><mrow><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup></mtd></mtr><mtr><mtd><mo>×</mo><msub><mi>P</mi><mrow><mn>2</mn><mi>n</mi></mrow></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mi>x</mi><mi>y</mi></mrow><mrow><msup><mrow><mrow><mo>[</mo><mrow><mn>1</mn><mo>−</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><msup><mi>x</mi><mn>2</mn></msup><mo stretchy="false">)</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mo>]</mo></mrow></mrow><mrow><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup></mrow></mfrac></mrow><mo>)</mo></mrow><mtext> </mtext><mn>.</mn></mtd></mtr></mtable></math> Note that in the limit x →0, <math id="M163"> <mrow> <msub> <mi mathvariant="script">F</mi> <mi>n</mi> </msub> <mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>→</mo><mi>π</mi><msub> <mi>P</mi> <mrow> <mn>2</mn><mi>n</mi></mrow> </msub> <mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>/</mo><mn>2</mn></mrow> </math> , and in the limit x →1, <math id="M164"> <mrow> <msub> <mi mathvariant="script">L</mi> <mrow> <mi>n</mi><mo>≠</mo><mn>0</mn></mrow> </msub> <mo stretchy="false">(</mo><mn>1</mn><mo stretchy="false">)</mo><mo>=</mo><mn>0</mn></mrow> </math> . For the case of n =0, the integral above can be analytically evaluated, resulting in <math display="block" id="M65"><mrow><msub><mi mathvariant="script">F</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><mfrac><mrow><mi>i</mi><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mtext> </mtext><mi>cosh</mi></mrow></mstyle><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo></mrow><mrow><msqrt><mrow><mn>1</mn><mo>−</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></msqrt></mrow></mfrac></mrow><mo>)</mo></mrow><mtext> </mtext><mn>.</mn></mrow></math> The first term in Eq. (55) contains the free-streaming moment <math id="M66"><mrow><msubsup><mi mathvariant="script">L</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></msubsup><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo></mrow></math> , which is given by <math display="block" id="M67"><mrow><msubsup><mi mathvariant="script">L</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></msubsup><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>ϵ</mi><mn>0</mn></msub><mtext> </mtext><mfrac><mrow><msub><mi>τ</mi><mn>0</mn></msub></mrow><mi>τ</mi></mfrac><mtext> </mtext><msub><mi mathvariant="script">F</mi><mi>n</mi></msub><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>τ</mi><mn>0</mn></msub></mrow><mi>τ</mi></mfrac></mrow><mo>)</mo></mrow><mtext> </mtext><mn>,</mn></mrow></math> where ϵ 0 is the initial energy density. Eq. (55) should be regarded as an integral equation, for which <math id="M165"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> must be solved with respect to the initial condition ϵ 0 . Once <math id="M166"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> <mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo></mrow> </math> is provided, higher-order <math id="M167"> <mrow> <msub> <mi mathvariant="script">L</mi> <mi>n</mi> </msub> </mrow> </math> 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 <math id="M168"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> , which again leads to an asymptotic series [ 66 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 <math id="M169"><mi mathvariant="script">L</mi></math> -moments. By substituting the definition of <math id="M170"><mi mathvariant="script">L</mi></math> -moments into the kinetic equation Eq. (43), one obtains a first-order differential equation in which adjacent <math id="M171"><mi mathvariant="script">L</mi></math> -moments are coupled, <math display="block" id="M68"><mrow><msub><mo>∂</mo><mi>τ</mi></msub><msub><mi mathvariant="script">L</mi><mi>n</mi></msub><mo>+</mo><mfrac><mrow><msub><mi>a</mi><mi>n</mi></msub><msub><mi mathvariant="script">L</mi><mi>n</mi></msub><mo>+</mo><msub><mi>b</mi><mi>n</mi></msub><msub><mi mathvariant="script">L</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>c</mi><mi>n</mi></msub><msub><mi mathvariant="script">L</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow><mi>τ</mi></mfrac><mo>=</mo><mo>−</mo><mfrac><mrow><msub><mi mathvariant="script">L</mi><mi>n</mi></msub></mrow><mrow><msub><mi>τ</mi><mi>R</mi></msub></mrow></mfrac><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><msub><mi>δ</mi><mrow><mi>n</mi><mn>0</mn></mrow></msub><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mrow></math> where n =0, 1, 2, …. The constant coefficients a n , b n , and c n arise from the recursion relation of the Legendre polynomials, <math display="block" id="M69"><mtable columnalign="left"><mtr><mtd><msub><mi>a</mi><mi>n</mi></msub><mo>=</mo><mfrac><mrow><mn>2</mn><mo stretchy="false">(</mo><mn>14</mn><msup><mi>n</mi><mn>2</mn></msup><mo>+</mo><mn>7</mn><mi>n</mi><mo>−</mo><mn>2</mn><mo stretchy="false">)</mo></mrow><mrow><mo stretchy="false">(</mo><mn>4</mn><mi>n</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>4</mn><mi>n</mi><mo>+</mo><mn>3</mn><mo stretchy="false">)</mo></mrow></mfrac><mtext> </mtext><mn>,</mn><mtext> </mtext><msub><mi>b</mi><mi>n</mi></msub><mo>=</mo><mfrac><mrow><mo stretchy="false">(</mo><mn>2</mn><mi>n</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo><mn>2</mn><mi>n</mi><mo stretchy="false">(</mo><mn>2</mn><mi>n</mi><mo>+</mo><mn>2</mn><mo stretchy="false">)</mo></mrow><mrow><mo stretchy="false">(</mo><mn>4</mn><mi>n</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>4</mn><mi>n</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></mfrac><mtext> </mtext><mn>,</mn></mtd></mtr><mtr><mtd><msub><mi>c</mi><mi>n</mi></msub><mo>=</mo><mfrac><mrow><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mn>2</mn><mi>n</mi><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>2</mn><mi>n</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>2</mn><mi>n</mi><mo>+</mo><mn>2</mn><mo stretchy="false">)</mo></mrow><mrow><mo stretchy="false">(</mo><mn>4</mn><mi>n</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>4</mn><mi>n</mi><mo>+</mo><mn>3</mn><mo stretchy="false">)</mo></mrow></mfrac><mtext> </mtext><mn>,</mn></mtd></mtr></mtable></math> 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 <math display="block" id="M70"><mtable columnalign="left"><mtr><mtd><msub><mi>a</mi><mn>0</mn></msub><mo>=</mo><mn>4</mn><mo>/</mo><mn>3,</mn><mtext> </mtext><msub><mi>b</mi><mn>1</mn></msub><mo>=</mo><mn>0,</mn><mtext> </mtext><msub><mi>c</mi><mn>0</mn></msub><mo>=</mo><mn>2</mn><mo>/</mo><mn>3,</mn></mtd></mtr><mtr><mtd><msub><mi>a</mi><mn>1</mn></msub><mo>=</mo><mn>38</mn><mo>/</mo><mn>21,</mn><mtext> </mtext><msub><mi>b</mi><mn>2</mn></msub><mo>=</mo><mn>8</mn><mo>/</mo><mn>15,</mn><mtext> </mtext><msub><mi>c</mi><mn>1</mn></msub><mo>=</mo><mo>−</mo><mn>12</mn><mo>/</mo><mn>35.</mn></mtd></mtr></mtable></math> The equation for <math id="M172"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> is the conservation of the energy moment, as in Eq. (50), whereas higher-order <math id="M173"> <mrow> <msub> <mi mathvariant="script">L</mi> <mi>n</mi> </msub> </mrow> </math> ’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 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 τ R T 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 <math id="M174"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> and <math id="M175"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> </mrow> </math> (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 <math id="M176"><mi mathvariant="script">L</mi></math> -moments involved, improvements from these corrections are observed. With the truncation at n < 5 (with <math id="M177"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> ,…, <math id="M178"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>4</mn> </msub> </mrow> </math> 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: <math display="block" id="M71"><mrow><mi>ψ</mi><mo>=</mo><mo stretchy="false">(</mo><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mn>,</mn><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mn>,</mn><msub><mi mathvariant="script">L</mi><mn>2</mn></msub><mn>,</mn><mo>…</mo><mo stretchy="false">)</mo><mtext> </mtext><mn>.</mn></mrow></math> Correspondingly, the dynamics of free streaming are captured by a constant tri-diagonal matrix: <math display="block" id="M72"><mrow><mover accent="true"><mi>H</mi><mo>^</mo></mover><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>a</mi><mn>0</mn></msub></mrow></mtd><mtd><mrow><msub><mi>c</mi><mn>0</mn></msub></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mo>…</mo></mtd></mtr><mtr><mtd><mrow><msub><mi>b</mi><mn>1</mn></msub></mrow></mtd><mtd><mrow><msub><mi>a</mi><mn>1</mn></msub></mrow></mtd><mtd><mrow><msub><mi>c</mi><mn>1</mn></msub></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mo>…</mo></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mrow><msub><mi>b</mi><mn>2</mn></msub></mrow></mtd><mtd><mrow><msub><mi>a</mi><mn>2</mn></msub></mrow></mtd><mtd><mrow><msub><mi>c</mi><mn>2</mn></msub></mrow></mtd><mtd><mo>…</mo></mtd></mtr><mtr><mtd><mo>…</mo></mtd><mtd><mo>…</mo></mtd><mtd><mo>…</mo></mtd><mtd><mo>…</mo></mtd><mtd><mo>…</mo></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mtext> </mtext><mn>,</mn></mrow></math> such that the equation of free streaming becomes <math display="block" id="M73"><mrow><msub><mo>∂</mo><mi>ρ</mi></msub><mi>ψ</mi><mo>+</mo><mover accent="true"><mi>H</mi><mo>^</mo></mover><mi>ψ</mi><mo>=</mo><mn>0</mn><mtext> </mtext><mn>.</mn></mrow></math> For convenience, we introduce ρ = ln (τ/τ 0 ). The solution of moments can be found, provided the eigenvalues of the matrix are determined. We note that the eigenvalues of <math id="M74"><mover accent="true"><mi>H</mi><mo>^</mo></mover></math> are complex, with eigenvectors satisfying <math display="block" id="M75"><mrow><mover accent="true"><mi>H</mi><mo>^</mo></mover><msub><mi>ϕ</mi><mi>n</mi></msub><mo>=</mo><msub><mi>λ</mi><mi>n</mi></msub><msub><mi>ϕ</mi><mi>n</mi></msub><mtext> </mtext><mo>→</mo><mtext> </mtext><msub><mi>ϕ</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>ϕ</mi><mi>n</mi></msub><mo stretchy="false">(</mo><msub><mi>τ</mi><mn>0</mn></msub><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><msub><mi>λ</mi><mi>n</mi></msub><mi>ρ</mi></mrow></msup><mtext> </mtext><mn>.</mn></mrow></math> Let us order these eigenvalues by their real part, i.e. , <math display="block" id="M76"><mrow><mtext>Re</mtext><msub><mi>λ</mi><mn>0</mn></msub><mo><</mo><mtext>Re</mtext><msub><mi>λ</mi><mn>1</mn></msub><mo><</mo><mtext>Re</mtext><msub><mi>λ</mi><mn>2</mn></msub><mo><</mo><mo>…</mo><mo><</mo><mtext>Re</mtext><msub><mi>λ</mi><mi>∞</mi></msub><mn>.</mn></mrow></math> Then, the solution of moments is <math display="block" id="M77"><mtable columnalign="left"><mtr><mtd><mi>ψ</mi><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mi>n</mi></munder></mstyle><msub><mi>κ</mi><mi>n</mi></msub><msup><mi>e</mi><mrow><mo>−</mo><msub><mi>λ</mi><mi>n</mi></msub><mi>ρ</mi></mrow></msup><msub><mi>ϕ</mi><mi>n</mi></msub><mo stretchy="false">(</mo><msub><mi>τ</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>=</mo><msup><mi>e</mi><mrow><mo>−</mo><msub><mi>λ</mi><mn>0</mn></msub><mi>ρ</mi></mrow></msup><mrow><mo>[</mo><mrow><msub><mi>κ</mi><mn>0</mn></msub><msub><mi>ϕ</mi><mn>0</mn></msub><mo stretchy="false">(</mo><msub><mi>τ</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>n</mi><mo>≠</mo><mn>0</mn></mrow></munder></mstyle><msub><mi>κ</mi><mi>n</mi></msub><msup><mi>e</mi><mrow><mo>−</mo><mo stretchy="false">(</mo><msub><mi>λ</mi><mi>n</mi></msub><mo>−</mo><msub><mi>λ</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mi>ρ</mi></mrow></msup><msub><mi>ϕ</mi><mi>n</mi></msub><mo stretchy="false">(</mo><msub><mi>τ</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow></mtd></mtr><mtr><mtd><mo>→</mo><msup><mi>e</mi><mrow><mo>−</mo><msub><mi>λ</mi><mn>0</mn></msub><mi>ρ</mi></mrow></msup><msub><mi>κ</mi><mn>0</mn></msub><msub><mi>ϕ</mi><mn>0</mn></msub><mtext> </mtext><mn>,</mn></mtd></mtr></mtable></math> 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 - λ 0 ) ρ ∼ λρ ≫ 1, irrespective of the initial conditions [ 66 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, <math display="block" id="M78"><mrow><msub><mi>g</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>≡</mo><mfrac><mrow><mtext>d</mtext><mtext> </mtext><mtext>ln</mtext><mtext> </mtext><msub><mi mathvariant="script">L</mi><mi>n</mi></msub></mrow><mrow><mtext>d</mtext><mtext> </mtext><mtext>ln</mtext><mtext> </mtext><mi>τ</mi></mrow></mfrac><mtext> </mtext><mn>,</mn></mrow></math> the dominance of the ground state indicates the existence of a stable fixed point, g n ( w ) → - λ 0 , 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 Fig. 6 , the ground-state eigenvalue is plotted as a function of the truncation order. Although asymptotically, when n →∞, λ 0 →-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, <math display="block" id="M80"><mrow><msub><mi>ϕ</mi><mn>0</mn></msub><mo>:</mo><mtext> </mtext><msub><mi mathvariant="script">L</mi><mi>n</mi></msub><mo>=</mo><msub><mi>P</mi><mrow><mn>2</mn><mi>n</mi></mrow></msub><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mo>=</mo><mo stretchy="false">(</mo><mo>−</mo><msup><mrow><mn>1</mn><mo stretchy="false">)</mo></mrow><mi>n</mi></msup><mfrac><mrow><mo stretchy="false">(</mo><mn>2</mn><mi>n</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo><mo>!</mo><mo>!</mo></mrow><mrow><mo stretchy="false">(</mo><mn>2</mn><mi>n</mi><mo stretchy="false">)</mo><mo>!</mo><mo>!</mo></mrow></mfrac><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mtext> </mtext><mn>.</mn></mrow></math> In the original phase-space distribution, these <math id="M179"><mi mathvariant="script">L</mi></math> -moments with the specified ratios characterize a distribution spanning along p ⊥ and shrinking in p z . It may also be proven that when the truncation order goes to infinity, the unstable fixed point is g - →-2, by noting that a n + b n + c n = 2. The corresponding eigenvector is given by <math display="block" id="M81"><mrow><msub><mi>ϕ</mi><mi>∞</mi></msub><mo>:</mo><mtext> </mtext><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mo>=</mo><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mo>=</mo><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mo>=</mo><mo>…</mo></mrow></math> In the opposite limit with τ / τ R →∞, the moments admit a series expansion in powers of 1/ τ , <math display="block" id="M82"><mrow><msub><mi mathvariant="script">L</mi><mi>n</mi></msub><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mi>n</mi></mrow><mi>∞</mi></munderover></mstyle><mfrac><mrow><msubsup><mi>α</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mi>m</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow><mrow><msup><mi>τ</mi><mi>n</mi></msup></mrow></mfrac><mtext> </mtext><mn>.</mn></mrow></math> This structure follows from the Chapman–Enskog expansion of the kinetic theory [ 89 89 ]. Except for <math id="M180"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> , <math id="M83"><mrow><msubsup><mi>α</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></msubsup><mo>=</mo><mi>ϵ</mi><msub><mi>δ</mi><mrow><mi>n</mi><mn>0</mn></mrow></msub></mrow></math> , the expansion coefficients are directly related to the transport coefficients in viscous hydrodynamics. For instance, by noticing that the leading term in <math id="M181"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> </mrow> </math> satisfies <math id="M182"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> <mo>=</mo><msub> <mi mathvariant="script">P</mi> <mi>L</mi> </msub> <mo>−</mo><msub> <mi mathvariant="script">P</mi> <mi>T</mi> </msub> <mo>∝</mo><mi>η</mi><mo>/</mo><mi>τ</mi></mrow> </math> , one is able to identify the shear viscosity as <math display="block" id="M84"><mrow><mi>η</mi><mo>=</mo><mfrac><mrow><msub><mi>b</mi><mn>1</mn></msub></mrow><mn>2</mn></mfrac><msub><mi>τ</mi><mi>R</mi></msub><mi>ϵ</mi><mtext> </mtext><mn>.</mn></mrow></math> The coefficients in the expansion are dimensional, which can be further expressed in terms of dimensional variables in the original moment equations, <math id="M85"><mrow><msubsup><mi>α</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mi>m</mi><mo stretchy="false">)</mo></mrow></msubsup><mo>=</mo><msubsup><mi>B</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mi>m</mi><mo stretchy="false">)</mo></mrow></msubsup><mi>ϵ</mi><msubsup><mi>τ</mi><mi>R</mi><mi>n</mi></msubsup></mrow></math> , with <math id="M86"><mrow><msubsup><mi>B</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mi>m</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow></math> dimensionless. The leading-order and next-leading-order expansion coefficients can be solved analytically. The leading-order coefficients, <math id="M87"><mrow><msubsup><mi>α</mi><mi>n</mi><mrow><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow></math> , in particular, determine the asymptotic value of g n ( 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 <math id="M183"> <mrow> <msub> <mi mathvariant="script">L</mi> <mi>n</mi> </msub> <mo>∼</mo><msup> <mi>τ</mi> <mrow><mo>−</mo><mn>4</mn><mo>/</mo><mn>3</mn><mo>−</mo><mi>n</mi></mrow> </msup> </mrow> </math> . For a conformal system in which τT is a constant, on the contrary, one finds <math id="M184"> <mrow> <msub> <mi mathvariant="script">L</mi> <mi>n</mi> </msub> <mo>∼</mo><msup> <mi>τ</mi> <mrow> <mo>−</mo><mn>4</mn><mo>/</mo><mn>3</mn><mo>−</mo><mn>2</mn><mi>n</mi><mo>/</mo><mn>3</mn></mrow> </msup> </mrow> </math> . Therefore, <math display="block" id="M88"><mrow><msub><mi>g</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>∞</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mo>−</mo><mi>n</mi><mtext> </mtext><mn>,</mn><mtext> </mtext><mtext> </mtext><msub><mi>τ</mi><mi>R</mi></msub><mo>=</mo><mtext> </mtext><mtext>const</mtext><mn>.</mn><mtext> </mtext></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mo>−</mo><mfrac><mrow><mn>2</mn><mi>n</mi></mrow><mn>3</mn></mfrac><mtext> </mtext><mn>,</mn><mtext> </mtext><msub><mi>τ</mi><mi>R</mi></msub><mi>T</mi><mo>=</mo><mtext> </mtext><mtext>const</mtext><mn>.</mn><mtext> </mtext></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> 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 <math id="M185"> <mrow> <msub> <mi mathvariant="script">L</mi> <mi>n</mi> </msub> </mrow> </math> , 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 75 ]. In Fig. 7 Fig. 7 , the attractors of the first five orders of <math id="M186"><mi mathvariant="script">L</mi></math> -moments are plotted in terms of g n ( w ) for a conformal system. The <math id="M187"><mi mathvariant="script">L</mi></math> -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 <math id="M188"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> and <math id="M189"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> </mrow> </math> . The simplest truncation that satisfies the conservation of energy and momentum is <math display="block" id="M89"><mrow><msub><mo>∂</mo><mi>τ</mi></msub><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mo>+</mo><mfrac><mn>1</mn><mi>τ</mi></mfrac><mo stretchy="false">(</mo><msub><mi>a</mi><mn>0</mn></msub><msub><mi mathvariant="script">L</mi><mn>0</mn></msub><mo>+</mo><msub><mi>c</mi><mn>0</mn></msub><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mo stretchy="false">)</mo><mo>=</mo><mn>0</mn></mrow></math> <math display="block" id="M90"><mrow><msub><mo>∂</mo><mi>τ</mi></msub><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mo>+</mo><mfrac><mn>1</mn><mi>τ</mi></mfrac><mo stretchy="false">(</mo><msub><mi>a</mi><mn>1</mn></msub><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mo>+</mo><msub><mi>b</mi><mn>1</mn></msub><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>R</mi></msub></mrow></mfrac><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mtext> </mtext><mn>,</mn></mrow></math> Together with the traceless condition <math id="M91"><mrow><msubsup><mi>T</mi><mi>μ</mi><mi>μ</mi></msubsup><mo>=</mo><mn>0</mn></mrow></math> , 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 <math id="M190"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> <mo>=</mo><msub> <mi mathvariant="script">P</mi> <mi>L</mi> </msub> <mo>−</mo><msub> <mi mathvariant="script">P</mi> <mi>T</mi> </msub> </mrow> </math> in terms of viscous corrections. For the system with Bjorken expansion, at a later time if one further identifies <math id="M191"> <mrow> <msub> <mi>c</mi> <mn>0</mn> </msub> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> </mrow> </math> as the ξξ component of the shear stress tensor, <math id="M92"><mrow><msub><mi>c</mi><mn>0</mn></msub><msub><mi mathvariant="script">L</mi><mn>1</mn></msub><mo>=</mo><mi>π</mi><mo>=</mo><msubsup><mi>π</mi><mi>ξ</mi><mi>ξ</mi></msubsup></mrow></math> , the truncated moment equations indeed lead to the Israel–Stewart hydro equations of motion: <math display="block" id="M93"><mrow><mo>∂</mo><mi>ϵ</mi><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>ϵ</mi><mi>τ</mi></mfrac><mo>=</mo><mo>−</mo><mfrac><mi>π</mi><mi>τ</mi></mfrac></mrow></math> <math display="block" id="M94"><mrow><mi>π</mi><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>η</mi><mi>τ</mi></mfrac><mo>−</mo><msub><mi>τ</mi><mi>π</mi></msub><msub><mo>∂</mo><mi>τ</mi></msub><mi>π</mi><mo>−</mo><msub><mi>β</mi><mrow><mi>π</mi><mi>π</mi></mrow></msub><msub><mi>τ</mi><mi>π</mi></msub><mfrac><mi>π</mi><mi>τ</mi></mfrac><mtext> </mtext><mn>.</mn></mrow></math> 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 1 , which agrees with [ 87 87 ]. In Bjorken flow, ambiguity arises from interpreting 1/ τ as the ξξ component of the shear tensor <math id="M95"><mrow><msubsup><mi>σ</mi><mi>ξ</mi><mi>ξ</mi></msubsup></mrow></math> 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 <math display="block" id="M96"><mrow><msub><mi>β</mi><mrow><mi>π</mi><mi>π</mi></mrow></msub><msub><mi>τ</mi><mi>π</mi></msub><mfrac><mi>π</mi><mi>τ</mi></mfrac><mo>=</mo><msub><mi>τ</mi><mi>π</mi></msub><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>π</mi><mi>τ</mi></mfrac><mo>−</mo><mfrac><mn>3</mn><mn>4</mn></mfrac><mrow><mo>(</mo><mrow><msub><mi>β</mi><mrow><mi>π</mi><mi>π</mi></mrow></msub><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><msub><mi>τ</mi><mi>π</mi></msub></mrow><mi>η</mi></mfrac><msup><mi>π</mi><mn>2</mn></msup><mo>+</mo><mi>O</mi><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mn>1</mn><mi>τ</mi></mfrac></mrow><mo>)</mo></mrow></mrow><mn>3</mn></msup><mtext> </mtext><mn>,</mn></mrow></math> the constitutive relation gives rise to the BRSSS formulation [ 45 45 ], <math display="block" id="M97"><mrow><mi>π</mi><mo>=</mo><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>η</mi><mi>τ</mi></mfrac><mo>−</mo><msub><mi>τ</mi><mi>π</mi></msub><mrow><mo>(</mo><mrow><msub><mo>∂</mo><mi>τ</mi></msub><mi>π</mi><mo>+</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mfrac><mi>π</mi><mi>τ</mi></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mn>3</mn><mn>4</mn></mfrac><mrow><mo>(</mo><mrow><msub><mi>β</mi><mrow><mi>π</mi><mi>π</mi></mrow></msub><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><msub><mi>τ</mi><mi>π</mi></msub></mrow><mi>τ</mi></mfrac><msup><mi>π</mi><mn>2</mn></msup><mtext> </mtext><mn>.</mn></mrow></math> Correspondingly, the second-order transport coefficients λ 1 and β ππ are related, and the evaluation in a weakly coupled medium <math display="block" id="M98"><mrow><msub><mi>λ</mi><mn>1</mn></msub><mo>=</mo><mfrac><mn>3</mn><mn>4</mn></mfrac><mrow><mo>(</mo><mrow><msub><mi>β</mi><mrow><mi>π</mi><mi>π</mi></mrow></msub><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac></mrow><mo>)</mo></mrow><mi>η</mi><msub><mi>τ</mi><mi>π</mi></msub><mo>=</mo><mfrac><mn>5</mn><mn>7</mn></mfrac><mi>η</mi><msub><mi>τ</mi><mi>π</mi></msub><mtext> </mtext><mn>,</mn></mrow></math> is confirmed [ 48 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 <math id="M192"> <mrow> <mo>∝</mo><msub> <mi>c</mi> <mn>1</mn> </msub> <msub> <mi mathvariant="script">L</mi> <mn>2</mn> </msub> <mo>/</mo><mi>τ</mi></mrow> </math> in Eq. (74b) [ 95 95 ]. This term can also be absorbed into the collision term via a redefinition of the relaxation time: <math display="block" id="M99"><mrow><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>R</mi></msub></mrow></mfrac><mo>→</mo><mrow><mo>[</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><msub><mi>c</mi><mn>1</mn></msub><msub><mi>τ</mi><mi>R</mi></msub></mrow><mi>τ</mi></mfrac><mfrac><mrow><msub><mi mathvariant="script">L</mi><mn>2</mn></msub></mrow><mrow><msub><mi mathvariant="script">L</mi><mn>1</mn></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>R</mi></msub></mrow></mfrac><mo>≡</mo><mfrac><mrow><msubsup><mi>Z</mi><mrow><mi>η</mi><mo>/</mo><mi>s</mi></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msubsup></mrow><mrow><msub><mi>τ</mi><mi>R</mi></msub></mrow></mfrac><mtext> </mtext><mn>.</mn></mrow></math> It is then straightforward to see that because the ratio in the factor is related to g 2 ( w ), <math display="block" id="M100"><mrow><mfrac><mrow><msub><mi mathvariant="script">L</mi><mn>2</mn></msub></mrow><mrow><msub><mi mathvariant="script">L</mi><mn>1</mn></msub></mrow></mfrac><mo>=</mo><mo>−</mo><mfrac><mrow><msub><mi>b</mi><mn>2</mn></msub></mrow><mrow><mi>w</mi><mo>+</mo><msub><mi>g</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>w</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>a</mi><mn>2</mn></msub></mrow></mfrac><mtext> </mtext><mn>,</mn></mrow></math> one may solve the coupled moment equations to arbitrary orders, provided g 2 ( w ) is given precisely. In terms of the simplest truncation that dynamically involves only <math id="M193"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>0</mn> </msub> </mrow> </math> and <math id="M194"> <mrow> <msub> <mi mathvariant="script">L</mi> <mn>1</mn> </msub> </mrow> </math> , 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 2 ( w ) in the calculations, the renormalized τ R (or η/s ) in an out-of-equilibrium system is obtained. In Fig. 8 Fig. 8 , the factor Z η/s is obtained via the attractor solution of g 2 ( 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 71 , 113 113 - 115 115 113-115 ]. The results of a numerical test of the renormalization scheme are shown in Fig. 5 Fig. 5 , where given a renormalized η/s , even the solution of the two-moment equations achieves good agreement in comparison with the exact solution. 3 Phenomenological development 1 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 6 - 8 8 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 116 , 117 117 ], JAM [ 118 118 ], and SMASH [ 119 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 120 - 123 123 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 2 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 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 16 , 130 130 ] for a review). Current heavy-ion experiment programs, such as the RHIC BES program [ 131 131 - 134 134 131-134 ] and the NA61/SHINE experiment at the Super Proton Synchrotron (SPS) [ 135 135 , 136 136 ], as well as future experiments at the Facility for Antiproton and Ion Research (FAIR) [ 137 137 , 138 138 ], Nuclotron-based Ion Collider Facility [ 139 139 ], and JPARC-HI [ 140 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 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 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 <math id="M101"><mrow><msqrt><mi>s</mi></msqrt><mo>≳</mo><mn>15</mn></mrow></math> 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 125 , 126 126 ]. Off-diagonal susceptibility correlations offer additional insights into the chemical freeze-out conditions in relativistic heavy-ion collisions [ 141 141 , 142 142 ]. The three blobs in Fig. 8 Fig. 8 indicate the averaged fireball trajectories for typical Au+Au collisions at <math id="M102"><mrow><msqrt><mi>s</mi></msqrt><mo>=</mo><mn>200,19.6</mn></mrow></math> , 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. 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 6 , 7 7 ], for a review). Recently, increasing interest has shifted toward studying heavy-ion collisions in the intermediate energy regime. At <math id="M103"><mrow><msqrt><mi>s</mi></msqrt><mo>∼</mo><mi mathvariant="script">O</mi><mo stretchy="false">(</mo><mn>10</mn><mo stretchy="false">)</mo></mrow></math> GeV, heavy-ion collisions strongly violate longitudinal boost invariance and require full 3D modeling of their dynamics [ 143 143 , 144 144 ]. It is important to employ initial conditions with nontrivial rapidity dependence. Complex 3D collision dynamics can be approximated by parametric energy deposition [ 144 144 - 149 149 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 127 , 150 150 ]. There are 3D initial conditions based on hadronic transport simulations [ 151 151 - 153 153 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 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 155 , 156 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 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 ∼2 R /sinh( y beam ) [ 127 127 , 152 152 , 158 158 ]. Here, R is the nuclear radius and the beam rapidity <math id="M104"><mrow><msub><mi>y</mi><mrow><mtext>beam</mtext></mrow></msub><mo>=</mo><mtext>arccosh</mtext><mo stretchy="false">(</mo><msqrt><mi>s</mi></msqrt><mo>/</mo><mo stretchy="false">(</mo><mn>2</mn><msub><mi>m</mi><mi>N</mi></msub><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow></math> . 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 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, <math display="block" id="M105"><mrow><msub><mo>∂</mo><mi>μ</mi></msub><msup><mi>T</mi><mrow><mi>μ</mi><mi>ν</mi></mrow></msup><mo>=</mo><msubsup><mi>J</mi><mrow><mtext>source</mtext></mrow><mi>ν</mi></msubsup><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mover accent="true"><mi>x</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow></math> <math display="block" id="M106"><mrow><msub><mo>∂</mo><mi>μ</mi></msub><msubsup><mi>J</mi><mi>B</mi><mi>μ</mi></msubsup><mo>=</mo><msub><mi>ρ</mi><mrow><mi>B</mi><mn>,</mn><mtext>source</mtext></mrow></msub><mo stretchy="false">(</mo><mi>τ</mi><mn>,</mn><mover accent="true"><mi>x</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mn>.</mn></mrow></math> Such a dynamic initialization scheme was initially proposed in Refs. [ 159 159 , 160 160 ] and has been adopted by several groups [ 127 127 , 153 153 , 161 161 , 162 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 163 ]. However, at vanishing net baryon density or μ B = 0 GeV, higher-order susceptibilities have been computed by lattice QCD [ 124 124 ]. These susceptibility coefficients were used to construct nuclear matter EoS at finite baryon densities through a Taylor expansion [ 164 164 - 166 166 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 167 ]. At μ B = 0 GeV, matching EoS between the two phases was performed on the trace anomaly [ 168 168 , 169 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 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 127 , 164 164 , 170 170 ] was shown to reproduce the rapidity dependence of particle production as well as the collision energy dependence of the STAR v n {2} ( n = 2,3) flow measurements in Au+Au collisions from 200 to 7.7 GeV [ 144 144 , 171 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 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 144 , 152 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 148 , 173 173 - 175 175 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 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 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 16 , 130 130 ]. The rapid expansion of the collision system drives the fluctuations related to the QCD critical point out of equilibrium [ 24 24 , 25 25 , 176 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 177 ]. In heavy-ion physics, numerical simulations of stochastic hydrodynamics have been developed by several groups to study thermal and critical fluctuations [ 20 20 , 22 22 , 178 178 - 180 180 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 181 ], and it is often referred to as "hydro-kinetic" [ 21 21 , 23 23 - 25 25 23-25 , 176 176 , 182 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 26 , 27 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 23 , 25 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 183 ]. 3.2 Quantitative characterization of QGP transport properties 2 3.2.1 Charge diffusion 3 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 127 , 153 153 , 191 191 , 192 192 ] and detailed modeling of the QGP chemistry [ 164 164 , 193 193 - 195 195 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 49 , 196 196 , 197 197 ]. Recently, such formulations were generalized to include multiple diffusion currents with their cross couplings [ 198 198 ]. The functional forms of the transport coefficients for the diffusion matrix have been studied in transport models [ 199 199 , 200 200 ]. Additional coupling terms with the shear and bulk viscous tensors appear in the second order of the gradient expansion [ 49 49 , 197 197 ]. Phenomenologically, the net baryon diffusion current transports more baryons from forward rapidities to the mid-rapidity region [ 67 67 , 196 196 , 201 201 , 202 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 67 , 201 201 ]. The cross-diffusion between the net baryon and net strangeness induces an oscillating distribution of the net strangeness current at later times [ 198 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 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 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 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 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 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 57 , 179 179 , 195 195 ]. Quark–antiquark pairs can be thermally produced in local fluid cells. As a <math id="M107"><mrow><mi>q</mi><mover accent="true"><mi>q</mi><mo>¯</mo></mover></mrow></math> 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 <math id="M108"><mrow><mi>q</mi><mover accent="true"><mi>q</mi><mo>¯</mo></mover></mrow></math> pair hadronizes. In Ref. [ 195 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 <math id="M109"><mrow><mi>p</mi><mover accent="true"><mi>p</mi><mo>¯</mo></mover></mrow></math> 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 3 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 204 , 205 205 ]. Figure 11 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 184 , 188 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 p T 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 13 , 206 206 ]. By adopting the EoS from the latest lattice QCD calculations [ 207 207 , 208 208 ], the Duke-OSU group presented the first Bayesian inference on the temperature-dependent shear and bulk viscosities [ 189 189 , 190 190 ] using the flow measurements at the LHC. There was tension of the Bayesian extracted bulk viscosity ( ζ/s )( T ) [ 189 189 , 190 190 ] with the parameterization used in the IP-Glasma hybrid framework [ 185 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 28 , 186 186 ]. A lower T peak with a smaller ( ζ/s ) max is favored by hadron mean p T 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 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 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 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 28 , 209 209 ] and viscous hydrodynamics [ 210 210 , 211 211 ] using an effective QCD kinetic theory (EKT), KøMPøST [ 37 37 , 38 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 187 ] with an effective ( η/s ) eff = 0.12 and temperature-dependent bulk viscosity (see Fig. 11h Fig. 11h ). A direct comparison with the simulations without the EKT phase [ 28 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 212 ]. This additional pre-equilibrium phase thus leads to a 35% larger extracted QGP bulk viscosity to describe the identified particle mean p T . A recent work [ 213 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 2 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 214 , 215 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 217 ]. A more systematic study has recently been carried out using the state-of-the-art IP-Glasma + MUSIC + UrQMD hybrid framework [ 28 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 d N ch /d η 10 at mid-rapidity. On the contrary, the correlation between v 2 and the initial momentum anisotropy from the pre-hydrodynamic phase becomes stronger in the lower-multiplicity collisions [ 35 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 218 ]. Last but not least, a recent work [ 219 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 {2} from the IP-Glasma hybrid framework increases as the charged hadron multiplicity decreases [ 28 28 ], which is not seen in the flow measurements at the LHC. A comparison with the results from [ 65 65 , 217 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 2 {4} opposite to the experimental measurements [ 220 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 221 ] and for a general flow background [ 216 216 , 222 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: <math display="block" id="M110"><mrow><mi>r</mi><mo>=</mo><mfrac><mi>ζ</mi><mrow><msub><mi>τ</mi><mi>Π</mi></msub></mrow></mfrac><mfrac><mn>1</mn><mrow><mi>e</mi><mo>+</mo><mi>P</mi><mo>+</mo><mi>Π</mi></mrow></mfrac><mo>+</mo><msubsup><mi>c</mi><mi>s</mi><mn>2</mn></msubsup><mn>.</mn></mrow></math> The causality condition requires the ratio r < 1 [ 216 216 ]. Fig. 12a 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. 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 Fig. 12 . These bubbles may cause unstable cavitation inside the QGP [ 223 223 - 225 225 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 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 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 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 226 , 227 227 ]. 3.4 Interdisciplinary connections 2 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 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 3 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 34 , 229 229 , 230 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 231 - 234 234 231-234 ]. Studies in Ref. [ 235 235 - 237 237 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 87 , 238 238 ]. A recent study explored the effect of a possible octupole deformation of 208 Pb on the v 2 / v 3 ratio [ 239 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 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 241 - 243 243 241-243 ]. However, the theoretically proposed triggers based on a two-component Glauber model were not effective in STAR measurements [ 244 244 ]. The measurements were more consistent with the saturation-based initial condition models [ 245 245 , 246 246 ]. Recently, new types of correlations between elliptic flow and mean transverse momentum were proposed, which are sensitive to the nuclear deformation [ 247 247 ]. Similarly, at LHC energies, 129 Xe+ 129 Xe collisions were measured in addition to 208 Pb+ 208 Pb collisions. In contrast to the spherical 208 Pb, the 129 Xe nucleus has a prolate deformation. Hydrodynamic simulations [ 248 248 ] predicted that the elliptic flow coefficients would be 25% larger in 129 Xe+ 129 Xe than those in 208 Pb+ 208 Pb in 0–5% central collisions. This strong enhancement of v 2 {2} in central 129 Xe+ 129 Xe was confirmed by LHC measurements [ 249 249 - 251 251 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 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 <math id="M111"><mrow><msubsup><mrow/><mrow><mn>44</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Ru+ <math id="M112"><mrow><msubsup><mrow/><mrow><mn>44</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Ru and <math id="M113"><mrow><msubsup><mrow/><mrow><mn>40</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Zr+ <math id="M114"><mrow><msubsup><mrow/><mrow><mn>40</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Zr collisions. A recent work [ 253 253 ] pointed out that the neutron skin in <math id="M115"><mrow><msubsup><mrow/><mrow><mn>40</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> 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 <math id="M116"><mrow><msubsup><mrow/><mrow><mn>44</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Ru+ <math id="M117"><mrow><msubsup><mrow/><mrow><mn>44</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Ru with respect to <math id="M118"><mrow><msubsup><mrow/><mrow><mn>40</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Zr+ <math id="M119"><mrow><msubsup><mrow/><mrow><mn>40</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Zr collisions. Moreover, realistic nucleon configurations based on density functional theory calculations for <math id="M120"><mrow><msubsup><mrow/><mrow><mn>44</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Ru and <math id="M121"><mrow><msubsup><mrow/><mrow><mn>40</mn></mrow><mrow><mn>96</mn></mrow></msubsup></mrow></math> Zr nuclei showed sizable differences in the particle elliptic flow v 2 between the two isobar collision systems [ 254 254 ]. 3.4.2 Applications of statistics and machine-learning techniques 3 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 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 255 , 256 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 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. 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 258 - 260 260 258-260 ] and the breaking of flow factorization ratios [ 261 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 v n ( p T ), and the subleading component was responsible for factorization breaking in hydrodynamics [ 261 261 ]. A recent study [ 262 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 261 , 263 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 fr om human beings [ 264 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 204 ]. Bayesian analysis has been applied to constrain the QCD EoS [ 205 205 ], QGP transport coefficients [ 167 167 , 189 189 , 190 190 , 265 265 - 269 269 265-269 ], initial longitudinal fluctuations [ 270 270 ], QCD jet energy loss distribution [ 271 271 , 272 272 ], and heavy-quark diffusion coefficients [ 273 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 274 ]. More recently, supervised learning with deep convolutional neural networks (CNN) was used to identify the nature of QCD phase transition [ 275 275 , 276 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 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 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 279 - 282 282 279-282 ]. 4 Summary 1 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 70 ]. In a similar manner, attractors from kinetic theory solutions without conformal symmetry are reproduced [ 70 70 , 283 283 ]. In the context of kinetic theory, the <math id="M195"><mi mathvariant="script">L</mi></math> -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 115 ]. These observations suggest the universal existence of attractors in out-of-equilibrium systems [ 284 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 255 , 256 256 , 285 285 - 287 287 285-287 ]. In recent years, there have been a few formulations of relativistic causal hydrodynamics, known as Israel–Stewart [ 44 44 ], BRSSS [ 45 45 ], DNMR [ 49 49 ], relativistic third-order dissipative hydrodynamics [ 51 51 ], and anisotropic hydrodynamics [ 112 112 , 226 226 , 227 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 288 ] under high degrees of symmetry [ 112 112 , 289 289 - 291 291 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 292 , 293 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. We normalize the four-velocity as U μ U μ =1, corresponding to the mostly negative metric convention: g μν =(+,-,-,-). Causality and stability can be achieved in first-order viscous hydrodynamics as well but within a frame other than the choice by Landau and Lifshitz or Eckart [ 43 43 ]. These transport coefficients have different evaluations for a strongly coupled system. From the 𝒩 =4 super-YM field theory, they are [ 45 45 ] <math display="block" id="M122"><mrow><msub><mi>C</mi><mi>τ</mi></msub><mo>=</mo><mfrac><mrow><mn>2</mn><mo>−</mo><mi>log</mi><mn>2</mn></mrow><mrow><mn>2</mn><mi>π</mi></mrow></mfrac><mtext> </mtext><mn>,</mn><mtext> </mtext><msub><mi>C</mi><mi>λ</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mn>2</mn><mi>π</mi></mrow></mfrac><mtext> </mtext><mn>.</mn></mrow></math> The convergence of the hydro gradient expansion also depends on the detailed identification of the expansion parameter. For instance, the dispersion relation consisting of perturbations around equilibrium gives rise to a series expansion in terms of the wave number, which is convergent [cf. Ref [ 55 55 ]]. On the contrary, for series expansion in real space over spatial gradients, the convergence property may depend on the initial condition [ 56 56 ]. Note that this recursion relation differs from that in [ 47 47 ] by rescaling w → C τ w and 1 - f n /4. References 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 https://doi.org/10.1103/RevModPhys.89.035001 Strongly coupled quark-gluon plasma in heavy ion collisions J.Y. Ollitrault , Anisotropy as a signature of transverse collective flow . Phys. Rev. D 46 , 229 - 245 ( 1992 ). https://doi.org/10.1103/PhysRevD.46.229 https://doi.org/10.1103/PhysRevD.46.229 Anisotropy as a signature of transverse collective flow B. Alver , G. Roland , Collision geometry fluctuations and triangular flow in heavy-ion collisions . Phys. Rev. C 81 , 054905 ( 2010 ). [Erratum: Phys.Rev.C 82, 039903 (2010)] . https://doi.org/10.1103/PhysRevC.82.039903 https://doi.org/10.1103/PhysRevC.82.039903 Collision geometry fluctuations and triangular flow in heavy-ion collisions S.A. Voloshin , Collective phenomena in ultra-relativistic nuclear collisions: anisotropic flow and more . Prog. Part. Nucl. Phys . 67 , 541 - 546 ( 2012 ). https://doi.org/10.1016/j.ppnp.2012.01.025 https://doi.org/10.1016/j.ppnp.2012.01.025 Collective phenomena in ultra-relativistic nuclear collisions: anisotropic flow and more P. Romatschke , New Developments in Relativistic Viscous Hydrodynamics . Int. J. Mod. Phys. E 19 , 1 - 53 ( 2010 ). https://doi.org/10.1142/S0218301310014613 https://doi.org/10.1142/S0218301310014613 New Developments in Relativistic Viscous Hydrodynamics U. Heinz , R. Snellings , Collective flow and viscosity in relativistic heavy-ion collisions . Ann. Rev. Nucl. Part. Sci . 63 , 123 - 151 ( 2013 ). https://doi.org/10.1146/annurev-nucl-102212-170540 https://doi.org/10.1146/annurev-nucl-102212-170540 Collective flow and viscosity in relativistic heavy-ion collisions C. Gale , S. Jeon , B. Schenke , Hydrodynamic Modeling of Heavy-Ion Collisions . Int. J. Mod. Phys. A 28 , 1340011 ( 2013 ). https://doi.org/10.1142/S0217751X13400113 https://doi.org/10.1142/S0217751X13400113 Hydrodynamic Modeling of Heavy-Ion Collisions L. Yan , A flow paradigm in heavy-ion collisions . Chin. Phys. C 42 , 042001 ( 2018 ). https://doi.org/10.1088/1674-1137/42/4/042001 https://doi.org/10.1088/1674-1137/42/4/042001 A flow paradigm in heavy-ion collisions W. Florkowski , M.P. Heller , M. Spalinski , New theories of relativistic hydrodynamics in the LHC era . Rept. Prog. Phys . 81 , 046001 ( 2018 ). arXiv:1707.02282 , https://doi.org/10.1088/1361-6633/aaa091 https://doi.org/10.1088/1361-6633/aaa091 New theories of relativistic hydrodynamics in the LHC era P. Romatschke , U. Romatschke , Relativistic Fluid Dynamics In and Out of Equilibrium, Cambridge Monographs on Mathematical Physics , ( Cambridge University Press , 2019 ). https://doi.org/10.1017/9781108651998 https://doi.org/10.1017/9781108651998 Relativistic Fluid Dynamics In and Out of Equilibrium, Cambridge Monographs on Mathematical Physics E. Shuryak , Physics of Strongly coupled Quark-Gluon Plasma . Prog. Part. Nucl. Phys . 62 , 48 - 101 ( 2009 ). https://doi.org/10.1016/j.ppnp.2008.09.001 https://doi.org/10.1016/j.ppnp.2008.09.001 Physics of Strongly coupled Quark-Gluon Plasma P. Kovtun , D.T. Son , A.O. Starinets , Viscosity in strongly interacting quantum field theories from black hole physics . Phys. Rev. Lett . 94 , 111601 ( 2005 ). https://doi.org/10.1103/PhysRevLett.94.111601 https://doi.org/10.1103/PhysRevLett.94.111601 Viscosity in strongly interacting quantum field theories from black hole physics S. Ryu , J.F. Paquet , C. Shen , et al ., Importance of the Bulk Viscosity of QCD in Ultrarelativistic Heavy-Ion Collisions . Phys. Rev. Lett . 115 , 132301 ( 2015 ). https://doi.org/10.1103/PhysRevLett.115.132301 https://doi.org/10.1103/PhysRevLett.115.132301 Importance of the Bulk Viscosity of QCD in Ultrarelativistic Heavy-Ion Collisions F.G. Gardim , G. Giacalone , M. Luzum , et al ., Revealing QCD thermodynamics in ultrarelativistic nuclear collisions . Nature Phys . 16 , 615 - 619 ( 2020 ). https://doi.org/10.1038/s41567-020-0846-4 https://doi.org/10.1038/s41567-020-0846-4 Revealing QCD thermodynamics in ultrarelativistic nuclear collisions F.G. Gardim , G. Giacalone , J.Y. Ollitrault , Measuring the speed of sound of the quark-gluon plasma in ultracentral nucleus-nucleus collisions . arXiv:1909.11609 Measuring the speed of sound of the quark-gluon plasma in ultracentral nucleus-nucleus collisions M.A. Stephanov , QCD phase diagram and the critical point . Prog. Theor. Phys. Suppl . 153 , 139 - 156 ( 2004 ). QCD phase diagram and the critical point M.A. Stephanov , QCD phase diagram and the critical point . [ Int. J. Mod. Phys. A20 , 4387 ( 2005 )]. https://doi.org/10.1142/S0217751X05027965 https://doi.org/10.1142/S0217751X05027965 QCD phase diagram and the critical point L.D. Landau , E.M. Lifshitz , Fluid Mechanics , Second Edition: Volume 6 ( Course of Theoretical Physics ), 2nd Edition, Course of theoretical physics / by L. D. Landau and E. M. Lifshitz, Vol. 6 , ( Butterworth-Heinemann , 1987 ) Course of Theoretical Physics P. Kovtun , G.D. Moore , P. Romatschke , The stickiness of sound: An absolute lower limit on viscosity and the breakdown of second order relativistic hydrodynamics . Phys. Rev. D 84 , 025006 ( 2011 ). https://doi.org/10.1103/PhysRevD.84.025006 https://doi.org/10.1103/PhysRevD.84.025006 The stickiness of sound: An absolute lower limit on viscosity and the breakdown of second order relativistic hydrodynamics J. Kapusta , B. Muller , M. Stephanov , Relativistic Theory of Hydrodynamic Fluctuations with Applications to Heavy Ion Collisions . Phys. Rev. C 85 , 054906 ( 2012 ). https://doi.org/10.1103/PhysRevC.85.054906 https://doi.org/10.1103/PhysRevC.85.054906 Relativistic Theory of Hydrodynamic Fluctuations with Applications to Heavy Ion Collisions C. Young , J. Kapusta , C. Gale , et al ., Thermally Fluctuating Second-Order Viscous Hydrodynamics and Heavy-Ion Collisions . Phys. Rev. C 91 , 044901 ( 2015 ). https://doi.org/10.1103/PhysRevC.91.044901 https://doi.org/10.1103/PhysRevC.91.044901 Thermally Fluctuating Second-Order Viscous Hydrodynamics and Heavy-Ion Collisions Y. Akamatsu , A. Mazeliauskas , D. Teaney , A kinetic regime of hydrodynamic fluctuations and long time tails for a Bjorken expansion . Phys. Rev. C 95 , 014909 ( 2017 ). https://doi.org/10.1103/PhysRevC.95.014909 https://doi.org/10.1103/PhysRevC.95.014909 A kinetic regime of hydrodynamic fluctuations and long time tails for a Bjorken expansion M. Singh , C. Shen , S. McDonald , et al ., Hydrodynamic Fluctuations in Relativistic Heavy-Ion Collisions . Nucl. Phys. A 982 , 319 - 322 ( 2019 ). https://doi.org/10.1016/j.nuclphysa.2018.10.061 https://doi.org/10.1016/j.nuclphysa.2018.10.061 Hydrodynamic Fluctuations in Relativistic Heavy-Ion Collisions X. An , G. Basar , M. Stephanov , et al ., Relativistic Hydrodynamic Fluctuations . Phys. Rev. C 100 , 024910 ( 2019 ). https://doi.org/10.1103/PhysRevC.100.024910 https://doi.org/10.1103/PhysRevC.100.024910 Relativistic Hydrodynamic Fluctuations M. Stephanov , Y. Yin , Hydrodynamics with parametric slowing down and fluctuations near the critical point . Phys. Rev. D 98 , 036006 ( 2018 ). https://doi.org/10.1103/PhysRevD.98.036006 https://doi.org/10.1103/PhysRevD.98.036006 Hydrodynamics with parametric slowing down and fluctuations near the critical point X. An , G. Basar , M. Stephanov , et al ., Fluctuation dynamics in a relativistic fluid with a critical point . Phys. Rev. C 102 , 034901 ( 2020 ). https://doi.org/10.1103/PhysRevC.102.034901 https://doi.org/10.1103/PhysRevC.102.034901 Fluctuation dynamics in a relativistic fluid with a critical point K. Rajagopal , G. Ridgway , R. Weller , et al ., Hydro+ in Action: Understanding the Out-of-Equilibrium Dynamics Near a Critical Point in the QCD Phase Diagram . arXiv:1908.08539 Hydro+ in Action: Understanding the Out-of-Equilibrium Dynamics Near a Critical Point in the QCD Phase Diagram L. Du , U. Heinz , K. Rajagopal , et al ., Fluctuation dynamics near the QCD critical point . arXiv:2004.02719 Fluctuation dynamics near the QCD critical point B. Schenke , C. Shen , P. Tribedy , Running the gamut of high energy nuclear collisions . arXiv:2005.14682 Running the gamut of high energy nuclear collisions V. Khachatryan , et al ., Evidence for Collective Multiparticle Correlations in p-Pb Collisions . Phys. Rev. Lett . 115 , 012301 ( 2015 ). https://doi.org/10.1103/PhysRevLett.115.012301 https://doi.org/10.1103/PhysRevLett.115.012301 Evidence for Collective Multiparticle Correlations in p-Pb Collisions V. Khachatryan , et al ., Measurement of long-range near-side two-particle angular correlations in pp collisions at <math id="M123"><mrow><msqrt><mi>s</mi></msqrt><mo>=</mo></mrow></math> 13 TeV . Phys. Rev. Lett . 116 , 172302 ( 2016 ). https://doi.org/10.1103/PhysRevLett.116.172302 https://doi.org/10.1103/PhysRevLett.116.172302 Measurement of long-range near-side two-particle angular correlations in pp collisions at s=13 TeV C. Aidala , et al ., Creation of quark plasma droplets with three distinct geometries . Nature Phys . 15 , 214 - 220 ( 2019 ). https://doi.org/10.1038/s41567-018-0360-0 https://doi.org/10.1038/s41567-018-0360-0 Creation of quark plasma droplets with three distinct geometries H. Niemi , G. Denicol , How large is the Knudsen number reached in fluid dynamical simulations of ultrarelativistic heavy ion collisions? . arXiv:1404.7327 How large is the Knudsen number reached in fluid dynamical simulations of ultrarelativistic heavy ion collisions? A. Kurkela , U.A. Wiedemann , B. Wu , Flow in AA and pA as an interplay of fluid-like and non-fluid like excitations . Eur. Phys. J. C 79 , 965 ( 2019 ). https://doi.org/10.1140/epjc/s10052-019-7428-6 https://doi.org/10.1140/epjc/s10052-019-7428-6 Flow in AA and pA as an interplay of fluid-like and non-fluid like excitations H. Mäntysaari , B. Schenke , C. Shen , et al ., Imprints of fluctuating proton shapes on flow in proton-lead collisions at the LHC . Phys. Lett. B 772 , 681 - 686 ( 2017 ). https://doi.org/10.1016/j.physletb.2017.07.038 https://doi.org/10.1016/j.physletb.2017.07.038 Imprints of fluctuating proton shapes on flow in proton-lead collisions at the LHC B. Schenke , C. Shen , P. Tribedy , Hybrid Color Glass Condensate and hydrodynamic description of the Relativistic Heavy Ion Collider small system scan . Phys. Lett. B 803 , 135322 ( 2020 ). https://doi.org/10.1016/j.physletb.2020.135322 https://doi.org/10.1016/j.physletb.2020.135322 Hybrid Color Glass Condensate and hydrodynamic description of the Relativistic Heavy Ion Collider small system scan P. Romatschke , Do nuclear collisions create a locally equilibrated quark-gluon plasma? Eur. Phys. J. C 77 , 21 ( 2017 ). https://doi.org/10.1140/epjc/s10052-016-4567-x https://doi.org/10.1140/epjc/s10052-016-4567-x Do nuclear collisions create a locally equilibrated quark-gluon plasma? A. Kurkela , A. Mazeliauskas , J.F. Paquet , et al ., Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions to Hydrodynamics with QCD Kinetic Theory . Phys. Rev. Lett . 122 , 122302 ( 2019 ). https://doi.org/10.1103/PhysRevLett.122.122302 https://doi.org/10.1103/PhysRevLett.122.122302 Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions to Hydrodynamics with QCD Kinetic Theory A. Kurkela , A. Mazeliauskas , J.F. Paquet , et al ., Effective kinetic description of event-by-event pre-equilibrium dynamics in high-energy heavy-ion collisions . Phys. Rev. C 99 , 034910 ( 2019 ). https://doi.org/10.1103/PhysRevC.99.034910 https://doi.org/10.1103/PhysRevC.99.034910 Effective kinetic description of event-by-event pre-equilibrium dynamics in high-energy heavy-ion collisions M.P. Heller , R.A. Janik , P. Witaszczyk , Hydrodynamic Gradient Expansion in Gauge Theory Plasmas . Phys. Rev. Lett . 110 , 211602 ( 2013 ). https://doi.org/10.1103/PhysRevLett.110.211602 https://doi.org/10.1103/PhysRevLett.110.211602 Hydrodynamic Gradient Expansion in Gauge Theory Plasmas B. Abbott , et al ., Observation of Gravitational Waves from a Binary Black Hole Merger . Phys. Rev. Lett . 116 , 061102 ( 2016 ). https://doi.org/10.1103/PhysRevLett.116.061102 https://doi.org/10.1103/PhysRevLett.116.061102 Observation of Gravitational Waves from a Binary Black Hole Merger B.P. Abbott , et al ., GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence . Phys. Rev. Lett . 116 , 241103 ( 2016 ). https://doi.org/10.1103/PhysRevLett.116.241103 https://doi.org/10.1103/PhysRevLett.116.241103 GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence B. Abbott , et al ., GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral . Phys. Rev. Lett . 119 , 161101 ( 2017 ). https://doi.org/10.1103/PhysRevLett.119.161101 https://doi.org/10.1103/PhysRevLett.119.161101 GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral P. Kovtun , First-order relativistic hydrodynamics is stable . JHEP 10 , 034 ( 2019 ). https://doi.org/10.1007/JHEP10(2019)034 https://doi.org/10.1007/JHEP10(2019)034 First-order relativistic hydrodynamics is stable W. Israel , J. Stewart , Transient relativistic thermodynamics and kinetic theory . Annals Phys . 118 , 341 - 372 ( 1979 ). https://doi.org/10.1016/0003-4916(79)90130-1 https://doi.org/10.1016/0003-4916(79)90130-1 Transient relativistic thermodynamics and kinetic theory R. Baier , P. Romatschke , D.T. Son , et al ., Relativistic viscous hydrodynamics, conformal invariance, and holography . JHEP 04 , 100 ( 2008 ). https://doi.org/10.1088/1126-6708/2008/04/100 https://doi.org/10.1088/1126-6708/2008/04/100 Relativistic viscous hydrodynamics, conformal invariance, and holography M.P. Heller , M. Spalinski , Hydrodynamics Beyond the Gradient Expansion: Resurgence and Resummation . Phys. Rev. Lett . 115 , 072501 ( 2015 ). https://doi.org/10.1103/PhysRevLett.115.072501 https://doi.org/10.1103/PhysRevLett.115.072501 Hydrodynamics Beyond the Gradient Expansion: Resurgence and Resummation G. Basar , G.V. Dunne , Hydrodynamics, resurgence, and transasymptotics . Phys. Rev. D 92 , 125011 ( 2015 ). https://doi.org/10.1103/PhysRevD.92.125011 https://doi.org/10.1103/PhysRevD.92.125011 Hydrodynamics, resurgence, and transasymptotics K. Dusling , G.D. Moore , D. Teaney , Radiative energy loss and v(2) spectra for viscous hydrodynamics . Phys. Rev. C 81 , 034907 ( 2010 ). https://doi.org/10.1103/PhysRevC.81.034907 https://doi.org/10.1103/PhysRevC.81.034907 Radiative energy loss and v(2) spectra for viscous hydrodynamics G. Denicol , H. Niemi , E. Molnar , et al ., Derivation of transient relativistic fluid dynamics from the Boltzmann equation . Phys. Rev. D 85 , 114047 ( 2012 ). [Erratum: Phys.Rev.D 91, 039902 (2015)] . https://doi.org/10.1103/PhysRevD.85.114047 https://doi.org/10.1103/PhysRevD.85.114047 Derivation of transient relativistic fluid dynamics from the Boltzmann equation A. El , Z. Xu , C. Greiner , Third-order relativistic dissipative hydrodynamics . Phys. Rev. C 81 , 041901 ( 2010 ). https://doi.org/10.1103/PhysRevC.81.041901 https://doi.org/10.1103/PhysRevC.81.041901 Third-order relativistic dissipative hydrodynamics A. Jaiswal , Relativistic third-order dissipative fluid dynamics from kinetic theory . Phys. Rev. C 88 , 021903 ( 2013 ). https://doi.org/10.1103/PhysRevC.88.021903 https://doi.org/10.1103/PhysRevC.88.021903 Relativistic third-order dissipative fluid dynamics from kinetic theory H. Grad , Asymptotic Theory of the Boltzmann Equation . Physics of Fluids 6 , 147 - 181 ( 1963 ). https://doi.org/10.1063/1.1706716 https://doi.org/10.1063/1.1706716 Asymptotic Theory of the Boltzmann Equation M.P. Heller , R.A. Janik , P. Witaszczyk , Hydrodynamic gradient expansion in gauge theory plasmas . Phys. Rev. Lett . 110 , 211602 ( 2013 ). https://doi.org/10.1103/PhysRevLett.110.211602 https://doi.org/10.1103/PhysRevLett.110.211602 Hydrodynamic gradient expansion in gauge theory plasmas G.S. Denicol , J. Noronha , Divergence of the chapman-enskog expansion in relativistic kinetic theory . ( 2016 ). arXiv:1608.07869 Divergence of the chapman-enskog expansion in relativistic kinetic theory S.s. Grozdanov , P.K. Kovtun , A.O. Starinets , et al ., Convergence of the Gradient Expansion in Hydrodynamics . Phys. Rev. Lett . 122 , 251601 ( 2019 ). https://doi.org/10.1103/PhysRevLett.122.251601 https://doi.org/10.1103/PhysRevLett.122.251601 Convergence of the Gradient Expansion in Hydrodynamics M.P. Heller , A. Serantes , M. Spaliński , et al ., The hydrodynamic gradient expansion in linear response theory . arXiv:2007.05524 The hydrodynamic gradient expansion in linear response theory J.I. Kapusta , C. Plumberg , Causal Electric Charge Diffusion and Balance Functions in Relativistic Heavy Ion Collisions . Phys. Rev. C 97 , 014906 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.014906 https://doi.org/10.1103/PhysRevC.97.014906 Causal Electric Charge Diffusion and Balance Functions in Relativistic Heavy Ion Collisions S. Chatrchyan , et al ., Observation of Long-Range Near-Side Angular Correlations in Proton-Lead Collisions at the LHC . Phys. Lett. B 718 , 795 - 814 ( 2013 ). https://doi.org/10.1016/j.physletb.2012.11.025 https://doi.org/10.1016/j.physletb.2012.11.025 Observation of Long-Range Near-Side Angular Correlations in Proton-Lead Collisions at the LHC A.M. Sirunyan , et al ., Multiparticle correlation studies in pPb collisions at <math id="M124"><mrow><msqrt><mrow><msub><mi>s</mi><mrow><mtext>NN</mtext></mrow></msub></mrow></msqrt><mo>=</mo></mrow></math> 8.16 TeV . Phys. Rev. C 101 , 014912 ( 2020 ). https://doi.org/10.1103/PhysRevC.101.014912 https://doi.org/10.1103/PhysRevC.101.014912 Multiparticle correlation studies in pPb collisions at sNN= 8.16 TeV S. Acharya , et al ., Investigations of Anisotropic Flow Using Multiparticle Azimuthal Correlations in pp, p-Pb, Xe-Xe, and Pb-Pb Collisions at the LHC . Phys. Rev. Lett . 123 , 142301 ( 2019 ). https://doi.org/10.1103/PhysRevLett.123.142301 https://doi.org/10.1103/PhysRevLett.123.142301 Investigations of Anisotropic Flow Using Multiparticle Azimuthal Correlations in pp, p-Pb, Xe-Xe, and Pb-Pb Collisions at the LHC A. Adare , et al ., Measurements of elliptic and triangular flow in high-multiplicity <math id="M125"><mrow><msup><mrow/><mn>3</mn></msup></mrow></math> He <math id="M126"><mo>+</mo></math> Au collisions at <math id="M127"><mrow><msqrt><mrow><msub><mi>s</mi><mrow><msub><mrow/><mrow><mi>N</mi><mi>N</mi></mrow></msub></mrow></msub></mrow></msqrt><mo>=</mo><mn>200</mn></mrow></math> GeV . Phys. Rev. Lett . 115 , 142301 ( 2015 ). https://doi.org/10.1103/PhysRevLett.115.142301 https://doi.org/10.1103/PhysRevLett.115.142301 Measurements of elliptic and triangular flow in high-multiplicity 3He+Au collisions at sNN=200 GeV A. Adare , et al ., Measurements of mass-dependent azimuthal anisotropy in central p <math id="M128"><mo>+</mo></math> Au, d <math id="M129"><mo>+</mo></math> Au, and <math id="M130"><mrow><msup><mrow/><mn>3</mn></msup></mrow></math> He <math id="M131"><mo>+</mo></math> Au collisions at <math id="M132"><mrow><msqrt><mrow><msub><mi>s</mi><mrow><msub><mrow/><mrow><mi>N</mi><mi>N</mi></mrow></msub></mrow></msub></mrow></msqrt><mo>=</mo><mn>200</mn></mrow></math> GeV . Phys. Rev. C 97 , 064904 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.064904 https://doi.org/10.1103/PhysRevC.97.064904 Measurements of mass-dependent azimuthal anisotropy in central p+Au, d +Au, and 3He +Au collisions at sNN=200 GeV R.A. Lacey , in 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions, Long-range collectivity in small collision-systems with two- and four-particle correlations @ STAR . 2020 . arXiv:2002.11889 J.L. Nagle , A. Adare , S. Beckman , et al ., Exploiting Intrinsic Triangular Geometry in Relativistic He3+Au Collisions to Disentangle Medium Properties . Phys. Rev. Lett . 113 , 112301 ( 2014 ). https://doi.org/10.1103/PhysRevLett.113.112301 https://doi.org/10.1103/PhysRevLett.113.112301 Exploiting Intrinsic Triangular Geometry in Relativistic He3+Au Collisions to Disentangle Medium Properties M. Habich , G. Miller , P. Romatschke , et al ., Testing hydrodynamic descriptions of p+p collisions at <math id="M133"><mrow><msqrt><mi>s</mi></msqrt><mo>=</mo><mn>7</mn></mrow></math> TeV . Eur. Phys. J. C 76 , 408 ( 2016 ). https://doi.org/10.1140/epjc/s10052-016-4237-z https://doi.org/10.1140/epjc/s10052-016-4237-z Testing hydrodynamic descriptions of p+p collisions at s=7 TeV M.P. Heller , V. Svensson , How does relativistic kinetic theory remember about initial conditions? Phys. Rev. D 98 , 054016 ( 2018 ). https://doi.org/10.1103/PhysRevD.98.054016 https://doi.org/10.1103/PhysRevD.98.054016 How does relativistic kinetic theory remember about initial conditions? G.S. Denicol , C. Gale , S. Jeon , et al ., Net baryon diffusion in fluid dynamic simulations of relativistic heavy-ion collisions . Phys. Rev. C 98 , 034916 ( 2018 ). https://doi.org/10.1103/PhysRevC.98.034916 https://doi.org/10.1103/PhysRevC.98.034916 Net baryon diffusion in fluid dynamic simulations of relativistic heavy-ion collisions M. Strickland , J. Noronha , G. Denicol , Anisotropic nonequilibrium hydrodynamic attractor . Phys. Rev. D 97 , 036020 ( 2018 ). https://doi.org/10.1103/PhysRevD.97.036020 https://doi.org/10.1103/PhysRevD.97.036020 Anisotropic nonequilibrium hydrodynamic attractor A. Behtash , C.N. Cruz-Camacho , M. Martinez , Far-from-equilibrium attractors and nonlinear dynamical systems approach to the Gubser flow . Phys. Rev . 97 , 044041 ( 2018 ). https://doi.org/10.1103/PhysRevD.97.044041 https://doi.org/10.1103/PhysRevD.97.044041 Far-from-equilibrium attractors and nonlinear dynamical systems approach to the Gubser flow P. Romatschke , Relativistic Hydrodynamic Attractors with Broken Symmetries: Non-Conformal and Non-Homogeneous . JHEP 12 , 079 ( 2017 ). https://doi.org/10.1007/JHEP12(2017)079 https://doi.org/10.1007/JHEP12(2017)079 Relativistic Hydrodynamic Attractors with Broken Symmetries: Non-Conformal and Non-Homogeneous P. Romatschke , Relativistic Fluid Dynamics Far From Local Equilibrium . Phys. Rev. Lett . 120 , 012301 ( 2018 ). https://doi.org/10.1103/PhysRevLett.120.012301 https://doi.org/10.1103/PhysRevLett.120.012301 Relativistic Fluid Dynamics Far From Local Equilibrium A. Kurkela , W. van der Schee , U.A. Wiedemann , et al ., Early- and Late-Time Behavior of Attractors in Heavy-Ion Collisions . Phys. Rev. Lett . 124 , 102301 ( 2020 ). https://doi.org/10.1103/PhysRevLett.124.102301 https://doi.org/10.1103/PhysRevLett.124.102301 Early- and Late-Time Behavior of Attractors in Heavy-Ion Collisions G.S. Denicol , J. Noronha , Exact hydrodynamic attractor of an ultrarelativistic gas of hard spheres . Phys. Rev. Lett . 124 , 152301 ( 2020 ). https://doi.org/10.1103/PhysRevLett.124.152301 https://doi.org/10.1103/PhysRevLett.124.152301 Exact hydrodynamic attractor of an ultrarelativistic gas of hard spheres C. Chattopadhyay , U.W. Heinz , Hydrodynamics from free-streaming to thermalization and back again . Phys. Lett. B 801 , 135158 ( 2020 ). https://doi.org/10.1016/j.physletb.2019.135158 https://doi.org/10.1016/j.physletb.2019.135158 Hydrodynamics from free-streaming to thermalization and back again M. Strickland , The non-equilibrium attractor for kinetic theory in relaxation time approximation . JHEP 12 , 128 ( 2018 ). https://doi.org/10.1007/JHEP12(2018)128 https://doi.org/10.1007/JHEP12(2018)128 The non-equilibrium attractor for kinetic theory in relaxation time approximation J. Brewer , L. Yan , Y. Yin , Adiabatic hydrodynamization in rapidly-expanding quark-gluon plasma . Adiabatic hydrodynamization in rapidly-expanding quark-gluon plasma J.P. Blaizot , L. Yan , Fluid dynamics of out of equilibrium boost invariant plasmas . Phys. Lett. B 780 , 283 - 286 ( 2018 ). https://doi.org/10.1016/j.physletb.2018.02.058 https://doi.org/10.1016/j.physletb.2018.02.058 Fluid dynamics of out of equilibrium boost invariant plasmas A. Behtash , S. Kamata , M. Martinez , et al ., Dynamical systems and nonlinear transient rheology of the far-from-equilibrium Bjorken flow . Phys. Rev. D 99 , 116012 ( 2019 ). https://doi.org/10.1103/PhysRevD.99.116012 https://doi.org/10.1103/PhysRevD.99.116012 Dynamical systems and nonlinear transient rheology of the far-from-equilibrium Bjorken flow A. Dash , V. Roy , Hydrodynamic attractors for Gubser flow . Phys. Lett. B 806 , 135481 ( 2020 ). https://doi.org/10.1016/j.physletb.2020.135481 https://doi.org/10.1016/j.physletb.2020.135481 Hydrodynamic attractors for Gubser flow A. Behtash , S. Kamata , M. Martinez , et al ., Global flow structure and exact formal transseries of the Gubser flow in kinetic theory . arXiv:1911.06406 Global flow structure and exact formal transseries of the Gubser flow in kinetic theory M.P. Heller , R. Jefferson , M. Spaliński , et al ., Hydrodynamic attractors in phase space . arXiv:2003.07368 Hydrodynamic attractors in phase space P.M. Chesler , How big are the smallest drops of quark-gluon plasma? JHEP 03 , 146 ( 2016 ). https://doi.org/10.1007/JHEP03(2016)146 https://doi.org/10.1007/JHEP03(2016)146 How big are the smallest drops of quark-gluon plasma? R. Baier , A.H. Mueller , D. Schiff , et al ., ’Bottom up’ thermalization in heavy ion collisions . Phys. Lett . 502 , 51 - 58 ( 2001 ). https://doi.org/10.1016/S0370-2693(01)00191-5 https://doi.org/10.1016/S0370-2693(01)00191-5 ’Bottom up’ thermalization in heavy ion collisions F. Gelis , E. Iancu , J. Jalilian-Marian , et al ., The Color Glass Condensate . Ann. Rev. Nucl. Part. Sci . 60 , 463 - 489 ( 2010 ). https://doi.org/10.1146/annurev.nucl.010909.083629 https://doi.org/10.1146/annurev.nucl.010909.083629 The Color Glass Condensate T. Epelbaum , F. Gelis , Pressure isotropization in high energy heavy ion collisions . Phys. Rev. Lett . 111 , 232301 ( 2013 ). https://doi.org/10.1103/PhysRevLett.111.232301 https://doi.org/10.1103/PhysRevLett.111.232301 Pressure isotropization in high energy heavy ion collisions J. Bjorken , Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region . Phys. Rev. D 27 , 140 - 151 ( 1983 ). https://doi.org/10.1103/PhysRevD.27.140 https://doi.org/10.1103/PhysRevD.27.140 Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region G. Denicol , C. Gale , S. Jeon , et al ., Effect of initial-state nucleon-nucleon correlations on collective flow in ultra-central heavy-ion collisions . arXiv:1406.7792 Effect of initial-state nucleon-nucleon correlations on collective flow in ultra-central heavy-ion collisions S. Jaiswal , C. Chattopadhyay , A. Jaiswal , et al ., Exact solutions and attractors of higher-order viscous fluid dynamics for Bjorken flow . Phys. Rev. C 100 , 034901 ( 2019 ). https://doi.org/10.1103/PhysRevC.100.034901 https://doi.org/10.1103/PhysRevC.100.034901 Exact solutions and attractors of higher-order viscous fluid dynamics for Bjorken flow J.P. Blaizot , L. Yan , Onset of hydrodynamics for a quark-gluon plasma from the evolution of moments of distribution functions . JHEP 11 , 161 ( 2017 ). https://doi.org/10.1007/JHEP11(2017)161 https://doi.org/10.1007/JHEP11(2017)161 Onset of hydrodynamics for a quark-gluon plasma from the evolution of moments of distribution functions S.S. Gubser , A. Yarom , Conformal hydrodynamics in Minkowski and de Sitter spacetimes . Nucl. Phys. B 846 , 469 - 511 ( 2011 ). https://doi.org/10.1016/j.nuclphysb.2011.01.012 https://doi.org/10.1016/j.nuclphysb.2011.01.012 Conformal hydrodynamics in Minkowski and de Sitter spacetimes G.S. Denicol , J. Noronha , Hydrodynamic attractor and the fate of perturbative expansions in Gubser flow . Phys. Rev. D 99 , 116004 ( 2019 ). https://doi.org/10.1103/PhysRevD.99.116004 https://doi.org/10.1103/PhysRevD.99.116004 Hydrodynamic attractor and the fate of perturbative expansions in Gubser flow G.S. Denicol , J. Noronha , Analytical attractor and the divergence of the slow-roll expansion in relativistic hydrodynamics . Phys. Rev. D 97 , 056021 ( 2018 ). https://doi.org/10.1103/PhysRevD.97.056021 https://doi.org/10.1103/PhysRevD.97.056021 Analytical attractor and the divergence of the slow-roll expansion in relativistic hydrodynamics A.R. Liddle , P. Parsons , J.D. Barrow , Formalizing the slow roll approximation in inflation . Phys. Rev. D 50 , 7222 - 7232 ( 1994 ). https://doi.org/10.1103/PhysRevD.50.7222 https://doi.org/10.1103/PhysRevD.50.7222 Formalizing the slow roll approximation in inflation G.S. Denicol , J. Noronha , Analytical attractor and the divergence of the slow-roll expansion in relativistic hydrodynamics . Phys. Rev. D 97 , 056021 ( 2018 ). https://doi.org/10.1103/PhysRevD.97.056021 https://doi.org/10.1103/PhysRevD.97.056021 Analytical attractor and the divergence of the slow-roll expansion in relativistic hydrodynamics J.P. Blaizot , L. Yan , Emergence of hydrodynamical behavior in expanding ultra-relativistic plasmas . Annals Phys . 412 , 167993 ( 2020 ). https://doi.org/10.1016/j.aop.2019.167993 https://doi.org/10.1016/j.aop.2019.167993 Emergence of hydrodynamical behavior in expanding ultra-relativistic plasmas F.J. Dyson , Divergence of perturbation theory in quantum electrodynamics . Phys. Rev . 85 , 631 - 632 ( 1952 ). https://doi.org/10.1103/PhysRev.85.631 https://doi.org/10.1103/PhysRev.85.631 Divergence of perturbation theory in quantum electrodynamics G.V. Dunne , M. Ünsal , Continuity and resurgence: Towards a continuum definition of the <math id="M134"><mrow><mi>ℂ</mi><mi>ℙ</mi><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></math> model . Phys. Rev. D 87 , 025015 ( 2013 ). https://doi.org/10.1103/PhysRevD.87.025015 https://doi.org/10.1103/PhysRevD.87.025015 Continuity and resurgence: Towards a continuum definition of the ℂℙ(n−1) model J. Blaizot , E. Iancu , A. Rebhan , On the apparent convergence of perturbative QCD at high temperature . Phys. Rev. D 68 , 025011 ( 2003 ). https://doi.org/10.1103/PhysRevD.68.025011 https://doi.org/10.1103/PhysRevD.68.025011 On the apparent convergence of perturbative QCD at high temperature J. Zinn-Justin , U.D. Jentschura , Multi-instantons and exact results i: conjectures, wkb expansions, and instanton interactions . Annals of Physics 313 , 197 - 267 ( 2004 ). https://doi.org/10.1016/j.aop.2004.04.004 https://doi.org/10.1016/j.aop.2004.04.004 Multi-instantons and exact results i: conjectures, wkb expansions, and instanton interactions J. Zinn-Justin , U.D. Jentschura , Multi-instantons and exact results ii: specific cases, higher-order effects, and numerical calculations . Annals of Physics 313 , 269 - 325 ( 2004 ). https://doi.org/10.1016/j.aop.2004.04.003 https://doi.org/10.1016/j.aop.2004.04.003 Multi-instantons and exact results ii: specific cases, higher-order effects, and numerical calculations I. Aniceto , G. Basar , R. Schiappa , A Primer on Resurgent Transseries and Their Asymptotics . Phys. Rept . 809 , 1 - 135 ( 2019 ). arXiv:1802.10441 , https://doi.org/10.1016/j.physrep.2019.02.003 https://doi.org/10.1016/j.physrep.2019.02.003 A Primer on Resurgent Transseries and Their Asymptotics I. Aniceto , R. Schiappa , Nonperturbative Ambiguities and the Reality of Resurgent Transseries . Commun. Math. Phys . 335 , 183 - 245 ( 2015 ). https://doi.org/10.1007/s00220-014-2165-z https://doi.org/10.1007/s00220-014-2165-z Nonperturbative Ambiguities and the Reality of Resurgent Transseries J.P. Blaizot , L. Yan , Analytical attractor for Bjorken expansion . arXiv:2006.08815 Analytical attractor for Bjorken expansion M.P. Heller , A. Kurkela , M. Spaliński , et al ., Hydrodynamization in kinetic theory: Transient modes and the gradient expansion . Phys. Rev. D 97 , 091503 ( 2018 ). https://doi.org/10.1103/PhysRevD.97.091503 https://doi.org/10.1103/PhysRevD.97.091503 Hydrodynamization in kinetic theory: Transient modes and the gradient expansion S. Groot , W. Leeuwen , C. van Weert , Relativistic kinetic theory: principles and applications , ( North-Holland Pub. Co. , 1980 ) Relativistic kinetic theory: principles and applications G. Baym , Thermal equilibration in ultra-relativistic heavy-ion collisions . Physics Letters B 138 , 18 - 22 ( 1984 ). https://doi.org/10.1016/0370-2693(84)91863-X https://doi.org/10.1016/0370-2693(84)91863-X Thermal equilibration in ultra-relativistic heavy-ion collisions W. Florkowski , R. Ryblewski , M. Strickland , Testing viscous and anisotropic hydrodynamics in an exactly solvable case . Phys. Rev. C 88 , 024903 ( 2013 ). https://doi.org/10.1103/PhysRevC.88.024903 https://doi.org/10.1103/PhysRevC.88.024903 Testing viscous and anisotropic hydrodynamics in an exactly solvable case D. Bazow , G. Denicol , U. Heinz , et al ., Nonlinear dynamics from the relativistic Boltzmann equation in the Friedmann-Lemaître-Robertson-Walker spacetime . Phys. Rev. D 94 , 125006 ( 2016 ). https://doi.org/10.1103/PhysRevD.94.125006 https://doi.org/10.1103/PhysRevD.94.125006 Nonlinear dynamics from the relativistic Boltzmann equation in the Friedmann-Lemaître-Robertson-Walker spacetime L. Tinti , G. Vujanovic , J. Noronha , et al ., Resummed hydrodynamic expansion for a plasma of particles interacting with fields . Phys. Rev. D 99 , 016009 ( 2019 ). https://doi.org/10.1103/PhysRevD.99.016009 https://doi.org/10.1103/PhysRevD.99.016009 Resummed hydrodynamic expansion for a plasma of particles interacting with fields M. Martinez , M. Strickland , Dissipative Dynamics of Highly Anisotropic Systems . Nucl. Phys . 848 , 183 - 197 ( 2010 ). https://doi.org/10.1016/j.nuclphysa.2010.08.011 https://doi.org/10.1016/j.nuclphysa.2010.08.011 Dissipative Dynamics of Highly Anisotropic Systems W. Florkowski , R. Ryblewski , Highly-anisotropic and strongly-dissipative hydrodynamics for early stages of relativistic heavy-ion collisions . Phys. Rev . 83 , 034907 ( 2011 ). https://doi.org/10.1103/PhysRevC.83.034907 https://doi.org/10.1103/PhysRevC.83.034907 Highly-anisotropic and strongly-dissipative hydrodynamics for early stages of relativistic heavy-ion collisions D. Bazow , U.W. Heinz , M. Strickland , Second-order (2+1)-dimensional anisotropic hydrodynamics . Phys. Rev . 90 , 054910 ( 2014 ). https://doi.org/10.1103/PhysRevC.90.054910 https://doi.org/10.1103/PhysRevC.90.054910 Second-order (2+1)-dimensional anisotropic hydrodynamics M. Lublinsky , E. Shuryak , How much entropy is produced in strongly coupled Quark-Gluon Plasma (sQGP) by dissipative effects? Phys. Rev. C 76 , 021901 ( 2007 ). https://doi.org/10.1103/PhysRevC.76.021901 https://doi.org/10.1103/PhysRevC.76.021901 How much entropy is produced in strongly coupled Quark-Gluon Plasma (sQGP) by dissipative effects? A. Behtash , C. Cruz-Camacho , S. Kamata , et al ., Non-perturbative rheological behavior of a far-from-equilibrium expanding plasma . Phys. Lett. B 797 , 134914 ( 2019 ). https://doi.org/10.1016/j.physletb.2019.134914 https://doi.org/10.1016/j.physletb.2019.134914 Non-perturbative rheological behavior of a far-from-equilibrium expanding plasma S. Kamata , M. Martinez , P. Plaschke , et al ., Hydrodynamization and non-equilibrium Green’s functions in kinetic theory . arXiv:2004.06751 Hydrodynamization and non-equilibrium Green’s functions in kinetic theory S.A. Bass , et al ., Microscopic models for ultrarelativistic heavy ion collisions . Prog. Part. Nucl. Phys . 41 , 255 - 369 ( 1998 ). https://doi.org/10.1016/S0146-6410(98)00058-1 https://doi.org/10.1016/S0146-6410(98)00058-1 Microscopic models for ultrarelativistic heavy ion collisions M. Bleicher , et al ., Relativistic hadron hadron collisions in the ultrarelativistic quantum molecular dynamics model . J. Phys . 25 , 1859 - 1896 ( 1999 ). arXiv:hep-ph/9909407 , https://doi.org/10.1088/0954-3899/25/9/308 https://doi.org/10.1088/0954-3899/25/9/308 Relativistic hadron hadron collisions in the ultrarelativistic quantum molecular dynamics model Y. Nara , N. Otuka , A. Ohnishi , et al ., Study of relativistic nuclear collisions at AGS energies from p + Be to Au + Au with hadronic cascade model . Phys. Rev . 61 , 024901 ( 2000 ). https://doi.org/10.1103/PhysRevC.61.024901 https://doi.org/10.1103/PhysRevC.61.024901 Study of relativistic nuclear collisions at AGS energies from p + Be to Au + Au with hadronic cascade model J. Weil , et al ., Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions . Phys. Rev . 94 , 054905 ( 2016 ). https://doi.org/10.1103/PhysRevC.94.054905 https://doi.org/10.1103/PhysRevC.94.054905 Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions H. Song , S.A. Bass , U. Heinz , et al ., 200 A GeV Au+Au collisions serve a nearly perfect quark-gluon liquid . Phys. Rev. Lett . 106 , 192301 ( 2011 ). [Erratum: Phys.Rev.Lett. 109, 139904 (2012)] . https://doi.org/10.1103/PhysRevLett.106.192301 https://doi.org/10.1103/PhysRevLett.106.192301 200 A GeV Au+Au collisions serve a nearly perfect quark-gluon liquid H. Song , S.A. Bass , U. Heinz , et al ., Hadron spectra and elliptic flow for 200 A GeV Au+Au collisions from viscous hydrodynamics coupled to a Boltzmann cascade . Phys. Rev. C 83 , 054910 ( 2011 ). [Erratum: Phys.Rev.C 86, 059903 (2012)] . https://doi.org/10.1103/PhysRevC.83.054910 https://doi.org/10.1103/PhysRevC.83.054910 Hadron spectra and elliptic flow for 200 A GeV Au+Au collisions from viscous hydrodynamics coupled to a Boltzmann cascade C. Shen , U. Heinz , P. Huovinen , et al ., Radial and elliptic flow in Pb+Pb collisions at the Large Hadron Collider from viscous hydrodynamic . Phys. Rev. C 84 , 044903 ( 2011 ). https://doi.org/10.1103/PhysRevC.84.044903 https://doi.org/10.1103/PhysRevC.84.044903 Radial and elliptic flow in Pb+Pb collisions at the Large Hadron Collider from viscous hydrodynamic U. Heinz , C. Shen , H. Song , The viscosity of quark-gluon plasma at RHIC and the LHC . AIP Conf. Proc . 1441 , 766 - 770 ( 2012 ). https://doi.org/10.1063/1.3700674 https://doi.org/10.1063/1.3700674 The viscosity of quark-gluon plasma at RHIC and the LHC A. Bazavov , et al ., Chiral crossover in QCD at zero and non-zero chemical potentials . Phys. Lett . 795 , 15 - 21 ( 2019 ). https://doi.org/10.1016/j.physletb.2019.05.013 https://doi.org/10.1016/j.physletb.2019.05.013 Chiral crossover in QCD at zero and non-zero chemical potentials A. Andronic , P. Braun-Munzinger , K. Redlich , et al ., Decoding the phase structure of QCD via particle production at high energy . Nature 561 , 321 - 330 ( 2018 ). https://doi.org/10.1038/s41586-018-0491-6 https://doi.org/10.1038/s41586-018-0491-6 Decoding the phase structure of QCD via particle production at high energy L. Adamczyk , et al ., Bulk Properties of the Medium Produced in Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program . Phys. Rev . 96 , 044904 ( 2017 ). https://doi.org/10.1103/PhysRevC.96.044904 https://doi.org/10.1103/PhysRevC.96.044904 Bulk Properties of the Medium Produced in Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program C. Shen , B. Schenke , Dynamical initial state model for relativistic heavy-ion collisions . Phys. Rev . 97 , 024907 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.024907 https://doi.org/10.1103/PhysRevC.97.024907 Dynamical initial state model for relativistic heavy-ion collisions C. Shen , in 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions, Studying QGP with flow: A theory overview . 2020 . arXiv:2001.11858 Y. Aoki , G. Endrodi , Z. Fodor , et al ., The Order of the quantum chromodynamics transition predicted by the standard model of particle physics . Nature 443 , 675 - 678 ( 2006 ). https://doi.org/10.1038/nature05120 https://doi.org/10.1038/nature05120 The Order of the quantum chromodynamics transition predicted by the standard model of particle physics A. Bzdak , S. Esumi , V. Koch , et al ., Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan . arXiv:1906.00936 Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan H. Caines , in Proceedings, 44th Rencontres de Moriond on QCD and High Energy Interactions: La Thuile, Italy, March 14-21, 2009 , The RHIC Beam Energy Scan: STAR’S Perspective . 2009 , pp. 375 - 378 . arXiv:0906.0305 The RHIC Beam Energy Scan: STAR’S Perspective B. Mohanty , STAR experiment results from the beam energy scan program at RHIC . J. Phys . 38 , 124023 ( 2011 ). arXiv:1106.5902 , https://doi.org/10.1088/0954-3899/38/12/124023 https://doi.org/10.1088/0954-3899/38/12/124023 STAR experiment results from the beam energy scan program at RHIC J.T. Mitchell , The RHIC Beam Energy Scan Program: Results from the PHENIX Experiment . Nucl. Phys . 904 - 905 , 903c - 906c ( 2013 ). https://doi.org/10.1016/j.nuclphysa.2013.02.161 https://doi.org/10.1016/j.nuclphysa.2013.02.161 The RHIC Beam Energy Scan Program: Results from the PHENIX Experiment G. Odyniec , Future of the beam energy scan program at RHIC . EPJ Web Conf . 95 , 03027 ( 2015 ). https://doi.org/10.1051/epjconf/20149503027 https://doi.org/10.1051/epjconf/20149503027 , https://doi.org/10.1051/epjconf/20159503027 https://doi.org/10.1051/epjconf/20159503027 Future of the beam energy scan program at RHIC M. Gazdzicki , Ion Program of Na61/Shine at the CERN SPS . J. Phys . 36 , 064039 ( 2009 ). https://doi.org/10.1088/0954-3899/36/6/064039 https://doi.org/10.1088/0954-3899/36/6/064039 Ion Program of Na61/Shine at the CERN SPS N. Abgrall , et al ., NA61/SHINE facility at the CERN SPS: beams and detector system . JINST 9 , P06005 ( 2014 ). https://doi.org/10.1088/1748-0221/9/06/P06005 https://doi.org/10.1088/1748-0221/9/06/P06005 NA61/SHINE facility at the CERN SPS: beams and detector system P. Spiller , G. Franchetti , The FAIR accelerator project at GSI . Nucl. Instrum. Meth . 561 , 305 - 309 ( 2006 ). https://doi.org/10.1016/j.nima.2006.01.043 https://doi.org/10.1016/j.nima.2006.01.043 The FAIR accelerator project at GSI T. Ablyazimov , et al ., Challenges in QCD matter physics -The scientific programme of the Compressed Baryonic Matter experiment at FAIR . Eur. Phys. J . 53 , 60 ( 2017 ). https://doi.org/10.1140/epja/i2017-12248-y https://doi.org/10.1140/epja/i2017-12248-y Challenges in QCD matter physics -The scientific programme of the Compressed Baryonic Matter experiment at FAIR A.N. Sissakian , A.S. Sorin , The nuclotron-based ion collider facility (NICA) at JINR: New prospects for heavy ion collisions and spin physics . J. Phys . 36 , 064069 ( 2009 ). https://doi.org/10.1088/0954-3899/36/6/064069 https://doi.org/10.1088/0954-3899/36/6/064069 The nuclotron-based ion collider facility (NICA) at JINR: New prospects for heavy ion collisions and spin physics T. Sakaguchi , High density matter physics at J-PARC-HI . PoS 2018 , 189 ( 2019 ). https://doi.org/10.22323/1.347.0189 https://doi.org/10.22323/1.347.0189 High density matter physics at J-PARC-HI R. Bellwied , J. Noronha-Hostler , P. Parotto , et al ., Freeze-out temperature from net-kaon fluctuations at energies available at the BNL Relativistic Heavy Ion Collider . Phys. Rev. C 99 , 034912 ( 2019 ). https://doi.org/10.1103/PhysRevC.99.034912 https://doi.org/10.1103/PhysRevC.99.034912 Freeze-out temperature from net-kaon fluctuations at energies available at the BNL Relativistic Heavy Ion Collider R. Bellwied , S. Borsanyi , Z. Fodor , et al ., Off-diagonal correlators of conserved charges from lattice QCD and how to relate them to experiment . Phys. Rev. D 101 , 034506 ( 2020 ). https://doi.org/10.1103/PhysRevD.101.034506 https://doi.org/10.1103/PhysRevD.101.034506 Off-diagonal correlators of conserved charges from lattice QCD and how to relate them to experiment J. Noronha , Collective effects in nuclear collisions: theory overview . Nucl. Phys. A 982 , 78 - 84 ( 2019 ). https://doi.org/10.1016/j.nuclphysa.2018.11.017 https://doi.org/10.1016/j.nuclphysa.2018.11.017 Collective effects in nuclear collisions: theory overview C. Shen , S. Alzhrani , A collision geometry-based 3D initial condition for relativistic heavy-ion collisions . arXiv:2003.05852 A collision geometry-based 3D initial condition for relativistic heavy-ion collisions T. Hirano , U.W. Heinz , D. Kharzeev , et al ., Hadronic dissipative effects on elliptic flow in ultrarelativistic heavy-ion collisions . Phys. Lett. B 636 , 299 - 304 ( 2006 ). https://doi.org/10.1016/j.physletb.2006.03.060 https://doi.org/10.1016/j.physletb.2006.03.060 Hadronic dissipative effects on elliptic flow in ultrarelativistic heavy-ion collisions P. Bozek , W. Broniowski , J. Moreira , Torqued fireballs in relativistic heavy-ion collisions . Phys. Rev. C 83 , 034911 ( 2011 ). arXiv:1011.3354 , https://doi.org/10.1103/PhysRevC.83.034911 https://doi.org/10.1103/PhysRevC.83.034911 Torqued fireballs in relativistic heavy-ion collisions P. Bozek , W. Broniowski , The torque effect and fluctuations of entropy deposition in rapidity in ultra-relativistic nuclear collisions . Phys. Lett. B 752 , 206 - 211 ( 2016 ). https://doi.org/10.1016/j.physletb.2015.11.054 https://doi.org/10.1016/j.physletb.2015.11.054 The torque effect and fluctuations of entropy deposition in rapidity in ultra-relativistic nuclear collisions P. Bozek , W. Broniowski , Longitudinal decorrelation measures of flow magnitude and event-plane angles in ultrarelativistic nuclear collisions . Phys. Rev. C 97 , 034913 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.034913 https://doi.org/10.1103/PhysRevC.97.034913 Longitudinal decorrelation measures of flow magnitude and event-plane angles in ultrarelativistic nuclear collisions A. Sakai , K. Murase , T. Hirano , Rapidity decorrelation of anisotropic flow caused by hydrodynamic fluctuations . arXiv:2003.13496 Rapidity decorrelation of anisotropic flow caused by hydrodynamic fluctuations A. Bialas , A. Bzdak , V. Koch , Stopped nucleons in configuration space . Acta Phys. Polon. B 49 , 103 ( 2018 ). https://doi.org/10.5506/APhysPolB.49.103 https://doi.org/10.5506/APhysPolB.49.103 Stopped nucleons in configuration space L. Pang , Q. Wang , X.N. Wang , Effects of initial flow velocity fluctuation in event-by-event (3+1)D hydrodynamics . Phys. Rev. C 86 , 024911 ( 2012 ). https://doi.org/10.1103/PhysRevC.86.024911 https://doi.org/10.1103/PhysRevC.86.024911 Effects of initial flow velocity fluctuation in event-by-event (3+1)D hydrodynamics I.A. Karpenko , P. Huovinen , H. Petersen , et al ., Estimation of the shear viscosity at finite net-baryon density from <math id="M135"><mrow><mi>A</mi><mo>+</mo><mi>A</mi></mrow></math> collision data at <math id="M136"><mrow><msqrt><mrow><msub><mi>s</mi><mrow><mtext>NN</mtext></mrow></msub></mrow></msqrt><mo>=</mo><mn>7.7</mn><mo>−</mo><mn>200</mn></mrow></math> GeV . Phys. Rev . 91 , 064901 ( 2015 ). https://doi.org/10.1103/PhysRevC.91.064901 https://doi.org/10.1103/PhysRevC.91.064901 Estimation of the shear viscosity at finite net-baryon density from A+A collision data at sNN=7.7−200 GeV L. Du , U. Heinz , G. Vujanovic , Hybrid model with dynamical sources for heavy-ion collisions at BES energies . Nucl. Phys. A 982 , 407 - 410 ( 2019 ). https://doi.org/10.1016/j.nuclphysa.2018.09.015 https://doi.org/10.1016/j.nuclphysa.2018.09.015 Hybrid model with dynamical sources for heavy-ion collisions at BES energies R. Anishetty , P. Koehler , L.D. McLerran , Central Collisions Between Heavy Nuclei at Extremely High-Energies: The Fragmentation Region . Phys. Rev. D 22 , 2793 ( 1980 ). https://doi.org/10.1103/PhysRevD.22.2793 https://doi.org/10.1103/PhysRevD.22.2793 Central Collisions Between Heavy Nuclei at Extremely High-Energies: The Fragmentation Region M. Li , J.I. Kapusta , Large Baryon Densities Achievable in High Energy Heavy Ion Collisions Outside the Central Rapidity Region . Phys. Rev. C 99 , 014906 ( 2019 ). https://doi.org/10.1103/PhysRevC.99.014906 https://doi.org/10.1103/PhysRevC.99.014906 Large Baryon Densities Achievable in High Energy Heavy Ion Collisions Outside the Central Rapidity Region L.D. McLerran , S. Schlichting , S. Sen , Spacetime picture of baryon stopping in the color-glass condensate . Phys. Rev. D 99 , 074009 ( 2019 ). https://doi.org/10.1103/PhysRevD.99.074009 https://doi.org/10.1103/PhysRevD.99.074009 Spacetime picture of baryon stopping in the color-glass condensate M. Attems , Y. Bea , J. Casalderrey-Solana , et al ., Holographic Collisions across a Phase Transition . Phys. Rev. Lett . 121 , 261601 ( 2018 ). https://doi.org/10.1103/PhysRevLett.121.261601 https://doi.org/10.1103/PhysRevLett.121.261601 Holographic Collisions across a Phase Transition C. Shen , B. Schenke , Initial state and hydrodynamic modeling of heavy-ion collisions at RHIC BES energies . PoS 2017 , 006 ( 2018 ). https://doi.org/10.22323/1.311.0006 https://doi.org/10.22323/1.311.0006 Initial state and hydrodynamic modeling of heavy-ion collisions at RHIC BES energies M. Okai , K. Kawaguchi , Y. Tachibana , et al ., New approach to initializing hydrodynamic fields and mini-jet propagation in quark-gluon fluids . Phys. Rev. C 95 , 054914 ( 2017 ). https://doi.org/10.1103/PhysRevC.95.054914 https://doi.org/10.1103/PhysRevC.95.054914 New approach to initializing hydrodynamic fields and mini-jet propagation in quark-gluon fluids C. Shen , G. Denicol , C. Gale , et al ., A hybrid approach to relativistic heavy-ion collisions at the RHIC BES energies . Nucl. Phys. A 967 , 796 - 799 ( 2017 ). https://doi.org/10.1016/j.nuclphysa.2017.06.008 https://doi.org/10.1016/j.nuclphysa.2017.06.008 A hybrid approach to relativistic heavy-ion collisions at the RHIC BES energies Y. Akamatsu , M. Asakawa , T. Hirano , et al ., Dynamically integrated transport approach for heavy-ion collisions at high baryon density . Phys. Rev. C 98 , 024909 ( 2018 ). https://doi.org/10.1103/PhysRevC.98.024909 https://doi.org/10.1103/PhysRevC.98.024909 Dynamically integrated transport approach for heavy-ion collisions at high baryon density Y. Kanakubo , Y. Tachibana , T. Hirano , Unified description of hadron chemistry from dynamical core-corona initialization . Phys. Rev. C 101 , 024912 ( 2020 ). https://doi.org/10.1103/PhysRevC.101.024912 https://doi.org/10.1103/PhysRevC.101.024912 Unified description of hadron chemistry from dynamical core-corona initialization C. Ratti , Lattice QCD and heavy ion collisions: a review of recent progress . Rept. Prog. Phys . 81 , 084301 ( 2018 ). https://doi.org/10.1088/1361-6633/aabb97 https://doi.org/10.1088/1361-6633/aabb97 Lattice QCD and heavy ion collisions: a review of recent progress A. Monnai , B. Schenke , C. Shen , Equation of state at finite densities for QCD matter in nuclear collisions . Phys. Rev. C 100 , 024907 ( 2019 ). https://doi.org/10.1103/PhysRevC.100.024907 https://doi.org/10.1103/PhysRevC.100.024907 Equation of state at finite densities for QCD matter in nuclear collisions J. Noronha-Hostler , P. Parotto , C. Ratti , et al ., Lattice-based equation of state at finite baryon number, electric charge and strangeness chemical potentials . Phys. Rev. C 100 , 064910 ( 2019 ). https://doi.org/10.1103/PhysRevC.100.064910 https://doi.org/10.1103/PhysRevC.100.064910 Lattice-based equation of state at finite baryon number, electric charge and strangeness chemical potentials P. Parotto , M. Bluhm , D. Mroczek , et al ., QCD equation of state matched to lattice data and exhibiting a critical point singularity . Phys. Rev. C 101 , 034901 ( 2020 ). https://doi.org/10.1103/PhysRevC.101.034901 https://doi.org/10.1103/PhysRevC.101.034901 QCD equation of state matched to lattice data and exhibiting a critical point singularity J.F. Paquet , et al ., in 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions, Revisiting Bayesian constraints on the transport coefficients of QCD . 2020 . arXiv:2002.05337 P. Huovinen , P. Petreczky , QCD Equation of State and Hadron Resonance Gas . Nucl. Phys. A 837 , 26 - 53 ( 2010 ). https://doi.org/10.1016/j.nuclphysa.2010.02.015 https://doi.org/10.1016/j.nuclphysa.2010.02.015 QCD Equation of State and Hadron Resonance Gas J.S. Moreland , R.A. Soltz , Hydrodynamic simulations of relativistic heavy-ion collisions with different lattice quantum chromodynamics calculations of the equation of state . Phys. Rev. C 93 , 044913 ( 2016 ). https://doi.org/10.1103/PhysRevC.93.044913 https://doi.org/10.1103/PhysRevC.93.044913 Hydrodynamic simulations of relativistic heavy-ion collisions with different lattice quantum chromodynamics calculations of the equation of state C. Shen , B. Schenke , Dynamical initialization and hydrodynamic modeling of relativistic heavy-ion collisions . Nucl. Phys. A 982 , 411 - 414 ( 2019 ). https://doi.org/10.1016/j.nuclphysa.2018.08.007 https://doi.org/10.1016/j.nuclphysa.2018.08.007 Dynamical initialization and hydrodynamic modeling of relativistic heavy-ion collisions L. Adamczyk , et al ., Harmonic decomposition of three-particle azimuthal correlations at energies available at the BNL Relativistic Heavy Ion Collider . Phys. Rev. C 98 , 034918 ( 2018 ). https://doi.org/10.1103/PhysRevC.98.034918 https://doi.org/10.1103/PhysRevC.98.034918 Harmonic decomposition of three-particle azimuthal correlations at energies available at the BNL Relativistic Heavy Ion Collider L. Adamczyk , et al ., Beam Energy Dependence of the Third Harmonic of Azimuthal Correlations in Au+Au Collisions at RHIC . Phys. Rev. Lett . 116 , 112302 ( 2016 ). https://doi.org/10.1103/PhysRevLett.116.112302 https://doi.org/10.1103/PhysRevLett.116.112302 Beam Energy Dependence of the Third Harmonic of Azimuthal Correlations in Au+Au Collisions at RHIC H. Li , L. Yan , Pseudorapidity dependent hydrodynamic response in heavy-ion collisions . Phys. Lett. B 802 , 135248 ( 2020 ). https://doi.org/10.1016/j.physletb.2020.135248 https://doi.org/10.1016/j.physletb.2020.135248 Pseudorapidity dependent hydrodynamic response in heavy-ion collisions R. Franco , M. Luzum , Rapidity-dependent eccentricity scaling in relativistic heavy-ion collisions . arXiv:1910.14598 Rapidity-dependent eccentricity scaling in relativistic heavy-ion collisions S. McDonald , S. Jeon , C. Gale , in 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions, Exploring Longitudinal Observables with 3+1D IP-Glasma . 2020 . arXiv:2001.08636 Y. Akamatsu , D. Teaney , F. Yan , et al ., Transits of the QCD critical point . Phys. Rev. C 100 , 044901 ( 2019 ). https://doi.org/10.1103/PhysRevC.100.044901 https://doi.org/10.1103/PhysRevC.100.044901 Transits of the QCD critical point E. Lifshitz , L. Pitaevskii , Statistical Physics: Theory of the Condensed State, no. pt. 2 in Course of Theoretical Physics , ( Elsevier Science , 2013 ) Statistical Physics: Theory of the Condensed State, no. pt. 2 in Course of Theoretical Physics M. Bluhm , et al ., Dynamics of critical fluctuations: Theory - phenomenology - heavy-ion collisions . arXiv:2001.08831 Dynamics of critical fluctuations: Theory - phenomenology - heavy-ion collisions A. De , C. Plumberg , J.I. Kapusta , Calculating Fluctuations and Self-Correlations Numerically for Causal Charge Diffusion in Relativistic Heavy-Ion Collisions . arXiv:2003.04878 Calculating Fluctuations and Self-Correlations Numerically for Causal Charge Diffusion in Relativistic Heavy-Ion Collisions M. Nahrgang , M. Bluhm , Modeling the diffusive dynamics of critical fluctuations near the QCD critical point . arXiv:2007.10371 Modeling the diffusive dynamics of critical fluctuations near the QCD critical point S. Pratt , Calculating n-point charge correlations in evolving systems . Phys. Rev. C 101 , 014914 ( 2020 ). https://doi.org/10.1103/PhysRevC.101.014914 https://doi.org/10.1103/PhysRevC.101.014914 . Calculating n-point charge correlations in evolving systems M. Martinez , T. Schäfer , Stochastic hydrodynamics and long time tails of an expanding conformal charged fluid . Phys. Rev. C 99 , 054902 ( 2019 ). https://doi.org/10.1103/PhysRevC.99.054902 https://doi.org/10.1103/PhysRevC.99.054902 Stochastic hydrodynamics and long time tails of an expanding conformal charged fluid S. Pratt , Calculating <math id="M137"><mi>n</mi></math> -Point Charge Correlations in Evolving Systems . Phys. Rev. C 101 , 014914 ( 2020 ). https://doi.org/10.1103/PhysRevC.101.014914 https://doi.org/10.1103/PhysRevC.101.014914 Calculating n-Point Charge Correlations in Evolving Systems C. Gale , S. Jeon , B. Schenke , et al ., Event-by-event anisotropic flow in heavy-ion collisions from combined Yang-Mills and viscous fluid dynamics . Phys. Rev. Lett . 110 , 012302 ( 2013 ). https://doi.org/10.1103/PhysRevLett.110.012302 https://doi.org/10.1103/PhysRevLett.110.012302 Event-by-event anisotropic flow in heavy-ion collisions from combined Yang-Mills and viscous fluid dynamics B. Schenke , C. Shen , P. Tribedy , Features of the IP-Glasma . Nucl. Phys. A 982 , 435 - 438 ( 2019 ). https://doi.org/10.1016/j.nuclphysa.2018.08.015 https://doi.org/10.1016/j.nuclphysa.2018.08.015 Features of the IP-Glasma B. Schenke , C. Shen , P. Tribedy , in 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions, Bulk properties and multi-particle correlations in large and small systems . 2020 . arXiv:2001.09949 C. Gale , J.F. Paquet , B. Schenke , et al ., in 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions, Probing Early-Time Dynamics and Quark-Gluon Plasma Transport Properties with Photons and Hadrons . 2020 . arXiv:2002.05191 H. Niemi , K. Eskola , R. Paatelainen , Event-by-event fluctuations in a perturbative QCD + saturation + hydrodynamics model: Determining QCD matter shear viscosity in ultrarelativistic heavy-ion collisions . Phys. Rev. C 93 , 024907 ( 2016 ). https://doi.org/10.1103/PhysRevC.93.024907 https://doi.org/10.1103/PhysRevC.93.024907 Event-by-event fluctuations in a perturbative QCD + saturation + hydrodynamics model: Determining QCD matter shear viscosity in ultrarelativistic heavy-ion collisions J.E. Bernhard , J.S. Moreland , S.A. Bass , et al ., Applying Bayesian parameter estimation to relativistic heavy-ion collisions: simultaneous characterization of the initial state and quark-gluon plasma medium . Phys. Rev. C 94 , 024907 ( 2016 ). https://doi.org/10.1103/PhysRevC.94.024907 https://doi.org/10.1103/PhysRevC.94.024907 Applying Bayesian parameter estimation to relativistic heavy-ion collisions: simultaneous characterization of the initial state and quark-gluon plasma medium J.E. Bernhard , J.S. Moreland , S.A. Bass , Bayesian estimation of the specific shear and bulk viscosity of quark-gluon plasma . Nature Phys . 15 , 1113 - 1117 ( 2019 ). https://doi.org/10.1038/s41567-019-0611-8 https://doi.org/10.1038/s41567-019-0611-8 Bayesian estimation of the specific shear and bulk viscosity of quark-gluon plasma M. Martinez , M.D. Sievert , D.E. Wertepny , et al ., Initial state fluctuations of QCD conserved charges in heavy-ion collisions . arXiv:1911.10272 Initial state fluctuations of QCD conserved charges in heavy-ion collisions M. Martinez , M.D. Sievert , D.E. Wertepny , et al ., Toward Initial Conditions of Conserved Charges Part II: The ICCING Monte Carlo Algorithm . arXiv:1911.12454 Toward Initial Conditions of Conserved Charges Part II: The ICCING Monte Carlo Algorithm S. Pratt , Identifying the Charge Carriers of the Quark-Gluon Plasma . Phys. Rev. Lett . 108 , 212301 ( 2012 ). https://doi.org/10.1103/PhysRevLett.108.212301 https://doi.org/10.1103/PhysRevLett.108.212301 Identifying the Charge Carriers of the Quark-Gluon Plasma S. Pratt , C. Plumberg , Evolving Charge Correlations in a Hybrid Model with both Hydrodynamics and Hadronic Boltzmann Descriptions . Phys. Rev . 99 , 044916 ( 2019 ). https://doi.org/10.1103/PhysRevC.99.044916 https://doi.org/10.1103/PhysRevC.99.044916 Evolving Charge Correlations in a Hybrid Model with both Hydrodynamics and Hadronic Boltzmann Descriptions S. Pratt , C. Plumberg , Determining the Diffusivity for Light Quarks from Experiment . arXiv:1904.11459 Determining the Diffusivity for Light Quarks from Experiment A. Monnai , Dissipative Hydrodynamic Effects on Baryon Stopping . Phys. Rev. C 86 , 014908 ( 2012 ). https://doi.org/10.1103/PhysRevC.86.014908 https://doi.org/10.1103/PhysRevC.86.014908 Dissipative Hydrodynamic Effects on Baryon Stopping A. Jaiswal , B. Friman , K. Redlich , Relativistic second-order dissipative hydrodynamics at finite chemical potential . Phys. Lett. B 751 , 548 - 552 ( 2015 ). https://doi.org/10.1016/j.physletb.2015.11.018 https://doi.org/10.1016/j.physletb.2015.11.018 Relativistic second-order dissipative hydrodynamics at finite chemical potential J.A. Fotakis , M. Greif , C. Greiner , et al ., Diffusion processes involving multiple conserved charges: A study from kinetic theory and implications to the fluid-dynamical modeling of heavy ion collisions . Phys. Rev. D 101 , 076007 ( 2020 ). https://doi.org/10.1103/PhysRevD.101.076007 https://doi.org/10.1103/PhysRevD.101.076007 Diffusion processes involving multiple conserved charges: A study from kinetic theory and implications to the fluid-dynamical modeling of heavy ion collisions M. Greif , J.A. Fotakis , G.S. Denicol , et al ., Diffusion of conserved charges in relativistic heavy ion collisions . Phys. Rev. Lett . 120 , 242301 ( 2018 ). https://doi.org/10.1103/PhysRevLett.120.242301 https://doi.org/10.1103/PhysRevLett.120.242301 Diffusion of conserved charges in relativistic heavy ion collisions J.B. Rose , M. Greif , J. Hammelmann , et al ., Cross-conductivity: novel transport coefficients to constrain the hadronic degrees of freedom of nuclear matter . arXiv:2001.10606 Cross-conductivity: novel transport coefficients to constrain the hadronic degrees of freedom of nuclear matter M. Li , C. Shen , Longitudinal Dynamics of High Baryon Density Matter in High Energy Heavy-Ion Collisions . Phys. Rev. C 98 , 064908 ( 2018 ). https://doi.org/10.1103/PhysRevC.98.064908 https://doi.org/10.1103/PhysRevC.98.064908 Longitudinal Dynamics of High Baryon Density Matter in High Energy Heavy-Ion Collisions L. Du , U. Heinz , (3+1)-dimensional dissipative relativistic fluid dynamics at non-zero net baryon density . Comput. Phys. Commun . 251 , 107090 ( 2020 ). https://doi.org/10.1016/j.cpc.2019.107090 https://doi.org/10.1016/j.cpc.2019.107090 (3+1)-dimensional dissipative relativistic fluid dynamics at non-zero net baryon density D. Kharzeev , Can gluons trace baryon number? Phys. Lett. B 378 , 238 - 246 ( 1996 ). https://doi.org/10.1016/0370-2693(96)00435-2 https://doi.org/10.1016/0370-2693(96)00435-2 Can gluons trace baryon number? J. Novak , K. Novak , S. Pratt , et al ., Determining Fundamental Properties of Matter Created in Ultrarelativistic Heavy-Ion Collisions . Phys. Rev. C 89 , 034917 ( 2014 ). https://doi.org/10.1103/PhysRevC.89.034917 https://doi.org/10.1103/PhysRevC.89.034917 Determining Fundamental Properties of Matter Created in Ultrarelativistic Heavy-Ion Collisions S. Pratt , E. Sangaline , P. Sorensen , et al ., Constraining the Eq. of State of Super-Hadronic Matter from Heavy-Ion Collisions . Phys. Rev. Lett . 114 , 202301 ( 2015 ). https://doi.org/10.1103/PhysRevLett.114.202301 https://doi.org/10.1103/PhysRevLett.114.202301 Constraining the Eq. of State of Super-Hadronic Matter from Heavy-Ion Collisions S. Ryu , J.F. Paquet , C. Shen , et al ., Effects of bulk viscosity and hadronic rescattering in heavy ion collisions at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider . Phys. Rev. C 97 , 034910 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.034910 https://doi.org/10.1103/PhysRevC.97.034910 Effects of bulk viscosity and hadronic rescattering in heavy ion collisions at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider S. Borsanyi , Z. Fodor , C. Hoelbling , et al ., Full result for the QCD equation of state with 2+1 flavors . Phys. Lett. B 730 , 99 - 104 ( 2014 ). https://doi.org/10.1016/j.physletb.2014.01.007 https://doi.org/10.1016/j.physletb.2014.01.007 Full result for the QCD equation of state with 2+1 flavors A. Bazavov , et al ., Equation of state in (2+1)-flavor QCD . Phys. Rev. D 90 , 094503 ( 2014 ). https://doi.org/10.1103/PhysRevD.90.094503 https://doi.org/10.1103/PhysRevD.90.094503 Equation of state in (2+1)-flavor QCD B. Schenke , P. Tribedy , R. Venugopalan , Fluctuating Glasma initial conditions and flow in heavy ion collisions . Phys. Rev. Lett . 108 , 252301 ( 2012 ). https://doi.org/10.1103/PhysRevLett.108.252301 https://doi.org/10.1103/PhysRevLett.108.252301 Fluctuating Glasma initial conditions and flow in heavy ion collisions B. Schenke , S. Jeon , C. Gale , Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics . Phys. Rev. Lett . 106 , 042301 ( 2011 ). https://doi.org/10.1103/PhysRevLett.106.042301 https://doi.org/10.1103/PhysRevLett.106.042301 Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics J.F. Paquet , C. Shen , G.S. Denicol , et al ., Production of photons in relativistic heavy-ion collisions . Phys. Rev. C 93 , 044906 ( 2016 ). https://doi.org/10.1103/PhysRevC.93.044906 https://doi.org/10.1103/PhysRevC.93.044906 Production of photons in relativistic heavy-ion collisions J. Liu , C. Shen , U. Heinz , Pre-equilibrium evolution effects on heavy-ion collision observables . Phys. Rev. C 91 , 064906 ( 2015 ). [Erratum: Phys.Rev.C 92, 049904 (2015)] . https://doi.org/10.1103/PhysRevC.91.064906 https://doi.org/10.1103/PhysRevC.91.064906 Pre-equilibrium evolution effects on heavy-ion collision observables T. Nunes da Silva , D. Chinellato , M. Hippert , et al ., Pre-hydrodynamic evolution and its signatures in final-state heavy-ion observables . arXiv:2006.02324 Pre-hydrodynamic evolution and its signatures in final-state heavy-ion observables P. Romatschke , Light-Heavy Ion Collisions: A window into pre-equilibrium QCD dynamics? Eur. Phys. J. C 75 , 305 ( 2015 ). https://doi.org/10.1140/epjc/s10052-015-3509-3 https://doi.org/10.1140/epjc/s10052-015-3509-3 Light-Heavy Ion Collisions: A window into pre-equilibrium QCD dynamics? G. Giacalone , B. Schenke , C. Shen , Observable signatures of initial state momentum anisotropies in nuclear collisions . Phys. Rev. Lett . 125 , 192301 ( 2020 ). https://doi.org/10.1103/PhysRevLett.125.192301 https://doi.org/10.1103/PhysRevLett.125.192301 Observable signatures of initial state momentum anisotropies in nuclear collisions F.S. Bemfica , M.M. Disconzi , J. Noronha , Causality of the Einstein-Israel-Stewart Theory with Bulk Viscosity . Phys. Rev. Lett . 122 , 221602 ( 2019 ). https://doi.org/10.1103/PhysRevLett.122.221602 https://doi.org/10.1103/PhysRevLett.122.221602 Causality of the Einstein-Israel-Stewart Theory with Bulk Viscosity R.D. Weller , P. Romatschke , One fluid to rule them all: viscous hydrodynamic description of event-by-event central p+p, p+Pb and Pb+Pb collisions at <math id="M138"><mrow><msqrt><mi>s</mi></msqrt><mo>=</mo><mn>5.02</mn></mrow></math> TeV . Phys. Lett. B 774 , 351 - 356 ( 2017 ). https://doi.org/10.1016/j.physletb.2017.09.077 https://doi.org/10.1016/j.physletb.2017.09.077 One fluid to rule them all: viscous hydrodynamic description of event-by-event central p+p, p+Pb and Pb+Pb collisions at s=5.02 TeV J. Sousa , J. Noronha , M. Luzum , in 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions, System response to the initial energy-momentum tensor in relativistic heavy-ion collisions . 2020 . arXiv:2002.12735 W. Zhao , C.M. Ko , Y.X. Liu , et al ., Probing the partonic degrees of freedom in high-multiplicity <math id="M139"><mrow><mi>p</mi><mo>−</mo><mtext>Pb</mtext></mrow></math> collisions at <math id="M140"><mrow><msqrt><mrow><msub><mi>s</mi><mrow><mi>N</mi><mi>N</mi></mrow></msub></mrow></msqrt><mo>=</mo><mn>5.02</mn><mtext>TeV</mtext></mrow></math> . Phys. Rev. Lett . 125 , 072301 ( 2020 ). https://doi.org/10.1103/PhysRevLett.125.072301 https://doi.org/10.1103/PhysRevLett.125.072301 Probing the partonic degrees of freedom in high-multiplicity p−Pb collisions at sNN=5.02TeV W. Zhao , Y. Zhou , K. Murase , et al ., Searching for small droplets of hydrodynamic fluid in proton-proton collisions at the LHC . arXiv:2001.06742 Searching for small droplets of hydrodynamic fluid in proton-proton collisions at the LHC S. Floerchinger , E. Grossi , Causality of fluid dynamics for high-energy nuclear collisions . JHEP 08 , 186 ( 2018 ). https://doi.org/10.1007/JHEP08(2018)186 https://doi.org/10.1007/JHEP08(2018)186 Causality of fluid dynamics for high-energy nuclear collisions F.S. Bemfica , M.M. Disconzi , V. Hoang , et al ., Nonlinear Constraints on Relativistic Fluids Far From Equilibrium . arXiv:2005.11632 Nonlinear Constraints on Relativistic Fluids Far From Equilibrium K. Rajagopal , N. Tripuraneni , Bulk Viscosity and Cavitation in Boost-Invariant Hydrodynamic Expansion . JHEP 03 , 018 ( 2010 ). https://doi.org/10.1007/JHEP03(2010)018 https://doi.org/10.1007/JHEP03(2010)018 Bulk Viscosity and Cavitation in Boost-Invariant Hydrodynamic Expansion M. Habich , P. Romatschke , Onset of cavitation in the quark-gluon plasma . JHEP 12 , 054 ( 2014 ). https://doi.org/10.1007/JHEP12(2014)054 https://doi.org/10.1007/JHEP12(2014)054 Onset of cavitation in the quark-gluon plasma M. Byres , S. Lim , C. McGinn , et al ., The Skinny on Bulk Viscosity and Cavitation in Heavy Ion Collisions . Phys. Rev. C 101 , 044902 ( 2020 ). https://doi.org/10.1103/PhysRevC.101.044902 https://doi.org/10.1103/PhysRevC.101.044902 The Skinny on Bulk Viscosity and Cavitation in Heavy Ion Collisions M. Alqahtani , M. Nopoush , R. Ryblewski , et al ., Anisotropic hydrodynamic modeling of 2.76 TeV Pb-Pb collisions . Phys. Rev. C 96 , 044910 ( 2017 ). https://doi.org/10.1103/PhysRevC.96.044910 https://doi.org/10.1103/PhysRevC.96.044910 Anisotropic hydrodynamic modeling of 2.76 TeV Pb-Pb collisions M. McNelis , D. Bazow , U. Heinz , (3+1)-dimensional anisotropic fluid dynamics with a lattice QCD equation of state . Phys. Rev. C 97 , 054912 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.054912 https://doi.org/10.1103/PhysRevC.97.054912 (3+1)-dimensional anisotropic fluid dynamics with a lattice QCD equation of state J. Berges , M.P. Heller , A. Mazeliauskas , et al ., Thermalization in QCD: theoretical approaches, phenomenological applications, and interdisciplinary connections . arXiv:2005.12299 Thermalization in QCD: theoretical approaches, phenomenological applications, and interdisciplinary connections B. Schenke , R. Venugopalan , Eccentric protons? Sensitivity of flow to system size and shape in p+p, p+Pb and Pb+Pb collisions . Phys. Rev. Lett . 113 , 102301 ( 2014 ). https://doi.org/10.1103/PhysRevLett.113.102301 https://doi.org/10.1103/PhysRevLett.113.102301 Eccentric protons? Sensitivity of flow to system size and shape in p+p, p+Pb and Pb+Pb collisions H. Mäntysaari , Review of proton and nuclear shape fluctuations at high energy . arXiv:2001.10705 Review of proton and nuclear shape fluctuations at high energy R. Shneor , et al ., Investigation of proton-proton short-range correlations via the C-12(e, e-prime pp) reaction . Phys. Rev. Lett . 99 , 072501 ( 2007 ). https://doi.org/10.1103/PhysRevLett.99.072501 https://doi.org/10.1103/PhysRevLett.99.072501 Investigation of proton-proton short-range correlations via the C-12(e, e-prime pp) reaction O. Hen , G. Miller , E. Piasetzky , et al ., Nucleon-Nucleon Correlations, Short-lived Excitations, and the Quarks Within . Rev. Mod. Phys . 89 , 045002 ( 2017 ). https://doi.org/10.1103/RevModPhys.89.045002 https://doi.org/10.1103/RevModPhys.89.045002 Nucleon-Nucleon Correlations, Short-lived Excitations, and the Quarks Within M. Duer , et al ., Direct Observation of Proton-Neutron Short-Range Correlation Dominance in Heavy Nuclei . Phys. Rev. Lett . 122 , 172502 ( 2019 ). https://doi.org/10.1103/PhysRevLett.122.172502 https://doi.org/10.1103/PhysRevLett.122.172502 Direct Observation of Proton-Neutron Short-Range Correlation Dominance in Heavy Nuclei A. Schmidt , et al ., Probing the core of the strong nuclear interaction . Nature 578 , 540 - 544 ( 2020 ). https://doi.org/10.1038/s41586-020-2021-6 https://doi.org/10.1038/s41586-020-2021-6 Probing the core of the strong nuclear interaction M. Alvioli , H.J. Drescher , M. Strikman , A Monte Carlo generator of nucleon configurations in complex nuclei including Nucleon-Nucleon correlations . Phys. Lett. B 680 , 225 - 230 ( 2009 ). https://doi.org/10.1016/j.physletb.2009.08.067 https://doi.org/10.1016/j.physletb.2009.08.067 A Monte Carlo generator of nucleon configurations in complex nuclei including Nucleon-Nucleon correlations W. Broniowski , M. Rybczynski , Two-body nucleon-nucleon correlations in Glauber models of relativistic heavy-ion collisions . Phys. Rev. C 81 , 064909 ( 2010 ). https://doi.org/10.1103/PhysRevC.81.064909 https://doi.org/10.1103/PhysRevC.81.064909 Two-body nucleon-nucleon correlations in Glauber models of relativistic heavy-ion collisions M. Alvioli , H. Holopainen , K. Eskola , et al ., Initial state anisotropies and their uncertainties in ultrarelativistic heavy-ion collisions from the Monte Carlo Glauber model . Phys. Rev. C 85 , 034902 ( 2012 ). https://doi.org/10.1103/PhysRevC.85.034902 https://doi.org/10.1103/PhysRevC.85.034902 Initial state anisotropies and their uncertainties in ultrarelativistic heavy-ion collisions from the Monte Carlo Glauber model C. Shen , Z. Qiu , U. Heinz , Shape and flow fluctuations in ultracentral Pb + Pb collisions at the energies available at the CERN Large Hadron Collider . Phys. Rev. C 92 , 014901 ( 2015 ). https://doi.org/10.1103/PhysRevC.92.014901 https://doi.org/10.1103/PhysRevC.92.014901 Shape and flow fluctuations in ultracentral Pb + Pb collisions at the energies available at the CERN Large Hadron Collider P. Carzon , S. Rao , M. Luzum , et al ., Possible octupole deformation of 208Pb and the ultracentral v2 to v3 puzzle . arXiv:2007.00780 Possible octupole deformation of 208Pb and the ultracentral v2 to v3 puzzle U.W. Heinz , A. Kuhlman , Anisotropic flow and jet quenching in ultrarelativistic U + U collisions . Phys. Rev. Lett . 94 , 132301 ( 2005 ). https://doi.org/10.1103/PhysRevLett.94.132301 https://doi.org/10.1103/PhysRevLett.94.132301 Anisotropic flow and jet quenching in ultrarelativistic U + U collisions S.A. Voloshin , Testing the Chiral Magnetic Effect with Central U+U collisions . Phys. Rev. Lett . 105 , 172301 ( 2010 ). https://doi.org/10.1103/PhysRevLett.105.172301 https://doi.org/10.1103/PhysRevLett.105.172301 Testing the Chiral Magnetic Effect with Central U+U collisions M. Rybczynski , W. Broniowski , G. Stefanek, Influence of initial fluctuations on geometry measures in relativistic U+U and Cu+Au collisions . Phys. Rev. C 87 , 044908 ( 2013 ). https://doi.org/10.1103/PhysRevC.87.044908 https://doi.org/10.1103/PhysRevC.87.044908 G. Stefanek, Influence of initial fluctuations on geometry measures in relativistic U+U and Cu+Au collisions A. Goldschmidt , Z. Qiu , C. Shen , et al ., Collision geometry and flow in uranium + uranium collisions . Phys. Rev. C 92 , 044903 ( 2015 ). https://doi.org/10.1103/PhysRevC.92.044903 https://doi.org/10.1103/PhysRevC.92.044903 Collision geometry and flow in uranium + uranium collisions L. Adamczyk , et al ., Azimuthal anisotropy in U+U and Au+Au collisions at RHIC . Phys. Rev. Lett . 115 , 222301 ( 2015 ). https://doi.org/10.1103/PhysRevLett.115.222301 https://doi.org/10.1103/PhysRevLett.115.222301 Azimuthal anisotropy in U+U and Au+Au collisions at RHIC J.S. Moreland , J.E. Bernhard , S.A. Bass , Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions . Phys. Rev. C 92 , 011901 ( 2015 ). https://doi.org/10.1103/PhysRevC.92.011901 https://doi.org/10.1103/PhysRevC.92.011901 Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions B. Schenke , P. Tribedy , R. Venugopalan , Initial-state geometry and fluctuations in Au + Au, Cu + Au, and U + U collisions at energies available at the BNL Relativistic Heavy Ion Collider . Phys. Rev. C 89 , 064908 ( 2014 ). https://doi.org/10.1103/PhysRevC.89.064908 https://doi.org/10.1103/PhysRevC.89.064908 Initial-state geometry and fluctuations in Au + Au, Cu + Au, and U + U collisions at energies available at the BNL Relativistic Heavy Ion Collider G. Giacalone , Observing the deformation of nuclei with relativistic nuclear collisions . arXiv:1910.04673 Observing the deformation of nuclei with relativistic nuclear collisions G. Giacalone , J. Noronha-Hostler , M. Luzum , et al ., Hydrodynamic predictions for 5.44 TeV Xe+Xe collisions . Phys. Rev. C 97 , 034904 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.034904 https://doi.org/10.1103/PhysRevC.97.034904 Hydrodynamic predictions for 5.44 TeV Xe+Xe collisions S. Acharya , et al ., Anisotropic flow in Xe-Xe collisions at <math id="M141"><mrow><msqrt><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>s</mi></mstyle><mrow><mtext>NN</mtext></mrow></msub></mrow></msqrt><mstyle mathvariant="bold" mathsize="normal"><mo>=</mo><mn>5</mn></mstyle><mstyle mathvariant="bold" mathsize="normal"><mo>.</mo><mn>4</mn><mn>4</mn></mstyle></mrow></math> TeV . Phys. Lett. B 784 , 82 - 95 ( 2018 ). https://doi.org/10.1016/j.physletb.2018.06.059 https://doi.org/10.1016/j.physletb.2018.06.059 Anisotropic flow in Xe-Xe collisions at sNN=5.44 TeV A.M. Sirunyan , et al ., Charged-particle angular correlations in XeXe collisions at <math id="M142"><mrow><msqrt><mrow><msub><mi>s</mi><mrow><msub><mrow/><mrow><mtext>NN</mtext></mrow></msub></mrow></msub></mrow></msqrt><mo>=</mo></mrow></math> 5.44 TeV . Phys. Rev. C 100 , 044902 ( 2019 ). https://doi.org/10.1103/PhysRevC.100.044902 https://doi.org/10.1103/PhysRevC.100.044902 Charged-particle angular correlations in XeXe collisions at sNN= 5.44 TeV G. Aad , et al ., Measurement of the azimuthal anisotropy of charged-particle production in <math id="M143"><mrow><mi>X</mi><mi>e</mi><mo>+</mo><mi>X</mi><mi>e</mi></mrow></math> collisions at <math id="M144"><mrow><msqrt><mrow><msub><mi>s</mi><mrow><mtext>NN</mtext></mrow></msub></mrow></msqrt><mo>=</mo><mn>5.44</mn></mrow></math> TeV with the ATLAS detector . Phys. Rev. C 101 , 024906 ( 2020 ). https://doi.org/10.1103/PhysRevC.101.024906 https://doi.org/10.1103/PhysRevC.101.024906 Measurement of the azimuthal anisotropy of charged-particle production in Xe+Xe collisions at sNN=5.44 TeV with the ATLAS detector B.A. Li , L.W. Chen , C.M. Ko , Recent Progress and New Challenges in Isospin Physics with Heavy-Ion Reactions . Phys. Rept . 464 , 113 - 281 ( 2008 ). https://doi.org/10.1016/j.physrep.2008.04.005 https://doi.org/10.1016/j.physrep.2008.04.005 Recent Progress and New Challenges in Isospin Physics with Heavy-Ion Reactions J. Hammelmann , A. Soto-Ontoso , M. Alvioli , et al ., Influence of the neutron-skin effect on nuclear isobar collisions at RHIC . arXiv:1908.10231 Influence of the neutron-skin effect on nuclear isobar collisions at RHIC H.J. Xu , X. Wang , H. Li , et al ., Importance of isobar density distributions on the chiral magnetic effect search . Phys. Rev. Lett . 121 , 022301 ( 2018 ). https://doi.org/10.1103/PhysRevLett.121.022301 https://doi.org/10.1103/PhysRevLett.121.022301 Importance of isobar density distributions on the chiral magnetic effect search D. Bazow , U.W. Heinz , M. Strickland , Massively parallel simulations of relativistic fluid dynamics on graphics processing units with CUDA . Comput. Phys. Commun . 225 , 92 - 113 ( 2018 ). https://doi.org/10.1016/j.cpc.2017.01.015 https://doi.org/10.1016/j.cpc.2017.01.015 Massively parallel simulations of relativistic fluid dynamics on graphics processing units with CUDA L.G. Pang , H. Petersen , X.N. Wang , Pseudorapidity distribution and decorrelation of anisotropic flow within the open-computing-language implementation CLVisc hydrodynamics . Phys. Rev. C 97 , 064918 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.064918 https://doi.org/10.1103/PhysRevC.97.064918 Pseudorapidity distribution and decorrelation of anisotropic flow within the open-computing-language implementation CLVisc hydrodynamics https://indico.cern.ch/event/433345/contributions/2321600/ https://indico.cern.ch/event/433345/contributions/2321600/ R.S. Bhalerao , J.Y. Ollitrault , S. Pal , et al ., Principal component analysis of event-by-event fluctuations . Phys. Rev. Lett . 114 , 152301 ( 2015 ). https://doi.org/10.1103/PhysRevLett.114.152301 https://doi.org/10.1103/PhysRevLett.114.152301 Principal component analysis of event-by-event fluctuations A. Mazeliauskas , D. Teaney , Subleading harmonic flows in hydrodynamic simulations of heavy ion collisions . Phys. Rev. C 91 , 044902 ( 2015 ). https://doi.org/10.1103/PhysRevC.91.044902 https://doi.org/10.1103/PhysRevC.91.044902 Subleading harmonic flows in hydrodynamic simulations of heavy ion collisions A. Sirunyan , et al ., Principal-component analysis of two-particle azimuthal correlations in PbPb and pPb collisions at CMS . Phys. Rev. C 96 , 064902 ( 2017 ). https://doi.org/10.1103/PhysRevC.96.064902 https://doi.org/10.1103/PhysRevC.96.064902 Principal-component analysis of two-particle azimuthal correlations in PbPb and pPb collisions at CMS A. Mazeliauskas , D. Teaney , Fluctuations of harmonic and radial flow in heavy ion collisions with principal components . Phys. Rev. C 93 , 024913 ( 2016 ). https://doi.org/10.1103/PhysRevC.93.024913 https://doi.org/10.1103/PhysRevC.93.024913 Fluctuations of harmonic and radial flow in heavy ion collisions with principal components M. Hippert , J.G.P. Barbon , D. Dobrigkeit Chinellato , et al ., Probing the structure of the initial state of heavy-ion collisions with pT-dependent flow fluctuations . arXiv:2006.13358 Probing the structure of the initial state of heavy-ion collisions with pT-dependent flow fluctuations F.G. Gardim , F. Grassi , P. Ishida , et al ., pT-dependent particle number fluctuations from principal-component analyses in hydrodynamic simulations of heavy-ion collisions . Phys. Rev. C 100 , 054905 ( 2019 ). https://doi.org/10.1103/PhysRevC.100.054905 https://doi.org/10.1103/PhysRevC.100.054905 pT-dependent particle number fluctuations from principal-component analyses in hydrodynamic simulations of heavy-ion collisions Z. Liu , W. Zhao , H. Song , Principal Component Analysis of collective flow in Relativistic Heavy-Ion Collisions . Eur. Phys. J. C 79 , 870 ( 2019 ). https://doi.org/10.1140/epjc/s10052-019-7379-y https://doi.org/10.1140/epjc/s10052-019-7379-y Principal Component Analysis of collective flow in Relativistic Heavy-Ion Collisions E. Sangaline , S. Pratt , Toward a deeper understanding of how experiments constrain the underlying physics of heavy-ion collisions . Phys. Rev. C 93 , 024908 ( 2016 ). https://doi.org/10.1103/PhysRevC.93.024908 https://doi.org/10.1103/PhysRevC.93.024908 Toward a deeper understanding of how experiments constrain the underlying physics of heavy-ion collisions J.E. Bernhard , P.W. Marcy , C.E. Coleman-Smith , et al ., Quantifying properties of hot and dense QCD matter through systematic model-to-data comparison . Phys. Rev. C 91 , 054910 ( 2015 ). https://doi.org/10.1103/PhysRevC.91.054910 https://doi.org/10.1103/PhysRevC.91.054910 Quantifying properties of hot and dense QCD matter through systematic model-to-data comparison J.F. Paquet , C. Shen , G. Denicol , et al ., Phenomenological constraints on the bulk viscosity of QCD . Nucl. Phys. A 967 , 429 - 432 ( 2017 ). https://doi.org/10.1016/j.nuclphysa.2017.06.024 https://doi.org/10.1016/j.nuclphysa.2017.06.024 Phenomenological constraints on the bulk viscosity of QCD J.S. Moreland , J.E. Bernhard , S.A. Bass , Bayesian calibration of a hybrid nuclear collision model using p-Pb and Pb-Pb data at energies available at the CERN Large Hadron Collider . Phys. Rev. C 101 , 024911 ( 2020 ). https://doi.org/10.1103/PhysRevC.101.024911 https://doi.org/10.1103/PhysRevC.101.024911 Bayesian calibration of a hybrid nuclear collision model using p-Pb and Pb-Pb data at energies available at the CERN Large Hadron Collider J. Auvinen , K.J. Eskola , P. Huovinen , et al ., Temperature dependence of η/s of strongly interacting matter: effects of the equation of state and the parametric form of (η/s)(T) . arXiv:2006.12499 Temperature dependence of η/s of strongly interacting matter: effects of the equation of state and the parametric form of (η/s)(T) W. Ke , J.S. Moreland , J.E. Bernhard , et al ., Constraints on rapidity-dependent initial conditions from charged particle pseudorapidity densities and two-particle correlations . Phys. Rev. C 96 , 044912 ( 2017 ). https://doi.org/10.1103/PhysRevC.96.044912 https://doi.org/10.1103/PhysRevC.96.044912 Constraints on rapidity-dependent initial conditions from charged particle pseudorapidity densities and two-particle correlations Y. He , L.G. Pang , X.N. Wang , Bayesian extraction of jet energy loss distributions in heavy-ion collisions . Phys. Rev. Lett . 122 , 252302 ( 2019 ). https://doi.org/10.1103/PhysRevLett.122.252302 https://doi.org/10.1103/PhysRevLett.122.252302 Bayesian extraction of jet energy loss distributions in heavy-ion collisions R. Soltz , Bayesian extraction of <math id="M145"><mover accent="true"><mi>q</mi><mo>^</mo></mover></math> with multi-stage jet evolution approach . PoS 2018 , 048 ( 2019 ). https://doi.org/10.22323/1.345.0048 https://doi.org/10.22323/1.345.0048 Bayesian extraction of q^ with multi-stage jet evolution approach Y. Xu , J.E. Bernhard , S.A. Bass , et al ., Data-driven analysis for the temperature and momentum dependence of the heavy-quark diffusion coefficient in relativistic heavy-ion collisions . Phys. Rev. C 97 , 014907 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.014907 https://doi.org/10.1103/PhysRevC.97.014907 Data-driven analysis for the temperature and momentum dependence of the heavy-quark diffusion coefficient in relativistic heavy-ion collisions S. Bass , A. Bischoff , J. Maruhn , et al ., Neural networks for impact parameter determination . Phys. Rev. C 53 , 2358 - 2363 ( 1996 ). https://doi.org/10.1103/PhysRevC.53.2358 https://doi.org/10.1103/PhysRevC.53.2358 Neural networks for impact parameter determination L.G. Pang , K. Zhou , N. Su , et al ., An equation-of-state-meter of quantum chromodynamics transition from deep learning . Nature Commun . 9 , 210 ( 2018 ). https://doi.org/10.1038/s41467-017-02726-3 https://doi.org/10.1038/s41467-017-02726-3 An equation-of-state-meter of quantum chromodynamics transition from deep learning Y.L. Du , K. Zhou , J. Steinheimer , et al ., Identifying the nature of the QCD transition in relativistic collision of heavy nuclei with deep learning . arXiv:1910.11530 Identifying the nature of the QCD transition in relativistic collision of heavy nuclei with deep learning J. Steinheimer , L. Pang , K. Zhou , et al ., A machine learning study to identify spinodal clumping in high energy nuclear collisions . JHEP 12 , 122 ( 2019 ). https://doi.org/10.1007/JHEP12(2019)122 https://doi.org/10.1007/JHEP12(2019)122 A machine learning study to identify spinodal clumping in high energy nuclear collisions H. Huang , B. Xiao , H. Xiong , et al ., Applications of deep learning to relativistic hydrodynamics . arXiv:1801.03334 Applications of deep learning to relativistic hydrodynamics Y.T. Chien , R. Kunnawalkam Elayavalli , Probing heavy ion collisions using quark and gluon jet substructure . arXiv:1803.03589 Probing heavy ion collisions using quark and gluon jet substructure Y.S. Lai , Automated Discovery of Jet Substructure Analyses . arXiv:1810.00835 Automated Discovery of Jet Substructure Analyses P.T. Komiske , E.M. Metodiev , J. Thaler , Energy Flow Networks: Deep Sets for Particle Jets . JHEP 01 , 121 ( 2019 ). https://doi.org/10.1007/JHEP01(2019)121 https://doi.org/10.1007/JHEP01(2019)121 Energy Flow Networks: Deep Sets for Particle Jets L.G. Pang , K. Zhou , X.N. Wang , Interpretable deep learning for nuclear deformation in heavy ion collisions . arXiv:1906.06429 Interpretable deep learning for nuclear deformation in heavy ion collisions W. Florkowski , E. Maksymiuk , R. Ryblewski , Coupled kinetic equations for fermions and bosons in the relaxation-time approximation . Phys. Rev. C 97 , 024915 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.024915 https://doi.org/10.1103/PhysRevC.97.024915 Coupled kinetic equations for fermions and bosons in the relaxation-time approximation G. Giacalone , A. Mazeliauskas , S. Schlichting , Hydrodynamic attractors, initial state energy and particle production in relativistic nuclear collisions . Phys. Rev. Lett . 123 , 262301 ( 2019 ). https://doi.org/10.1103/PhysRevLett.123.262301 https://doi.org/10.1103/PhysRevLett.123.262301 Hydrodynamic attractors, initial state energy and particle production in relativistic nuclear collisions H. Marrochio , J. Noronha , G.S. Denicol , et al ., Solutions of Conformal Israel-Stewart Relativistic Viscous Fluid Dynamics . Phys. Rev. C 91 , 014903 ( 2015 ). https://doi.org/10.1103/PhysRevC.91.014903 https://doi.org/10.1103/PhysRevC.91.014903 Solutions of Conformal Israel-Stewart Relativistic Viscous Fluid Dynamics C. Shen , Z. Qiu , H. Song , et al ., The iEBE-VISHNU code package for relativistic heavy-ion collisions . Comput. Phys. Commun . 199 , 61 - 85 ( 2016 ). https://doi.org/10.1016/j.cpc.2015.08.039 https://doi.org/10.1016/j.cpc.2015.08.039 The iEBE-VISHNU code package for relativistic heavy-ion collisions J. Noronha-Hostler , J. Noronha , F. Grassi , Bulk viscosity-driven suppression of shear viscosity effects on the flow harmonics at energies available at the BNL Relativistic Heavy Ion Collider . Phys. Rev. C 90 , 034907 ( 2014 ). https://doi.org/10.1103/PhysRevC.90.034907 https://doi.org/10.1103/PhysRevC.90.034907 Bulk viscosity-driven suppression of shear viscosity effects on the flow harmonics at energies available at the BNL Relativistic Heavy Ion Collider W. Florkowski , R. Ryblewski , M. Strickland , Anisotropic Hydrodynamics for Rapidly Expanding Systems . Nucl. Phys . 916 , 249 - 259 ( 2013 ). https://doi.org/10.1016/j.nuclphysa.2013.08.004 https://doi.org/10.1016/j.nuclphysa.2013.08.004 Anisotropic Hydrodynamics for Rapidly Expanding Systems W. Florkowski , E. Maksymiuk , R. Ryblewski , et al ., Exact solution of the (0+1)-dimensional Boltzmann equation for a massive gas . Phys. Rev. C 89 , 054908 ( 2014 ). https://doi.org/10.1103/PhysRevC.89.054908 https://doi.org/10.1103/PhysRevC.89.054908 Exact solution of the (0+1)-dimensional Boltzmann equation for a massive gas G.S. Denicol , W. Florkowski , R. Ryblewski , et al ., Shear-bulk coupling in nonconformal hydrodynamics . Phys. Rev. C 90 , 044905 ( 2014 ). https://doi.org/10.1103/PhysRevC.90.044905 https://doi.org/10.1103/PhysRevC.90.044905 Shear-bulk coupling in nonconformal hydrodynamics W. Florkowski , R. Ryblewski , M. Strickland , et al ., Leading-order anisotropic hydrodynamics for systems with massive particles . Phys. Rev. C 89 , 054909 ( 2014 ). https://doi.org/10.1103/PhysRevC.89.054909 https://doi.org/10.1103/PhysRevC.89.054909 Leading-order anisotropic hydrodynamics for systems with massive particles Y.X. Zhang , et al ., Comparison of heavy-ion transport simulations: Collision integral in a box . Phys. Rev. C 97 , 034625 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.034625 https://doi.org/10.1103/PhysRevC.97.034625 Comparison of heavy-ion transport simulations: Collision integral in a box A. Ono , et al ., Comparison of heavy-ion transport simulations: Collision integral with pions and Δ resonances in a box . Phys. Rev. C 100 , 044617 ( 2019 ). https://doi.org/10.1103/PhysRevC.100.044617 https://doi.org/10.1103/PhysRevC.100.044617 Comparison of heavy-ion transport simulations: Collision integral with pions and Δ resonances in a box We thank C. Gale, U. Heinz, J. P. Blaizot, J. F. Paquet, and B. Schenke for fruitful discussions. We thank the JETSCAPE Collaboration for providing the preliminary results shown in Fig. 11 Fig. 11 . 10.1007/s41365-020-00829-z This work was supported in part by the U.S. Department of Energy (DOE) (No. DE-SC0013460), the National Science Foundation (NSF) (No. PHY-2012922), the National Natural Science Foundation of China (No. 11975079),and the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration. 2020-06-30 2020-10-19 2020-10-22 2020-12-04 2020-12-01 2021-04-10 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2020 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-020-00829-z http://dx.doi.org/10.1007/s41365-020-00829-z . 2020 Nuclear Science and Techniques 12 31