Hydrodynamic attractors in heavy ion collisions: a review

A review of the recent progress of relativistic hydrodynamic attractors is presented, with a focus on applications in heavy ion collisions and the quark gluon plasma. Pedagogical introductions to the effective descriptions relevant for attractors in high energy physics, namely hydrodynamics, holography and kinetic theory, are followed by highlights of some recent advances.


Introduction
Relativistic hydrodynamics [1][2][3] is an incredibly useful tool in effective theories, capturing the underlying symmetries of a system in the long wavelength limit. It has been successfully applied at a vast swathe of scales: from the smallest in understanding quantum chromodynamics (QCD) at high energies [4], to condensed matter systems [5] and to the largest in cosmological applications [6,7].
The application of hydrodynamics is not always readily apparent. In fact, in the high energy context one of the most unexpected surprises to come from evaluating the data at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) was that the quark-gluon plasma (QGP), the extreme state of matter composed from liberated constituents of neutrons and protons following collisions of heavy ions, behaves like the "most perfect fluid" and as such is amenable to a description via relativistic hydrodynamics. Abundant evidence points to the QGP exhibiting collective flow, such as the observation of elliptic flow [8,9] and the particularly low value of the ratio of the shear viscosity to entropy density [10], close to the value predicted from holography [11][12][13]. For further mechanisms and features of the thermalization of the QGP, see [14] and references therein.
As such, the hydrodynamic description of the QGP had great phenomenological success. Active research deals with the hydrodynamization, i.e. the approach of the post-collision a e-mail: alexander.soloviev9810@gmail.com (corresponding author) state to a fluid description, see [15] for a review. It is worth reiterating how surprising this is since hydrodynamics is an effective theory based on the standard assumption that the system is near thermal equilibrium, although what exactly sets the regime of applicability of hydrodynamics was not clear. In typical heavy ion collisions, the system is certainly not in or even near equilibrium, which naturally beckons the question: why is hydrodynamics so unreasonably effective?
Recently, there has been progress in understanding how systems with little hope of being in equilibrium become nonetheless well-approximated by hydrodynamics. One proposed solution is due to the emergence of attractor phenomena. In a word, the study of attractors stems from the observation that solutions of evolution equations with different initial conditions evolve to the same late time behavior, which is called the (hydrodynamic) attractor. Quite generally, nonequilibrium solutions of different microscopic theories decay to attractors, losing knowledge of their initial conditions well before reaching thermal equilibrium. Specifically, heavy ion collisions have a mix of hydrodynamic and nonhydrodynamic modes, the latter of which decay rapidly in the expanding medium, leaving behind just the long wavelength modes, well described by hydrodynamics. For the heavy ion community, attractors came to the forefront with the remarkable discovery of the hydrodynamic attractor in an expanding boost-invariant background by Heller and Spalinski [16], building on the earlier results of Borel resummed hydrodynamic gradient expansion in holographic Bjorken flow [17,18].
Spurred on by the characterization of the hydrodynamic attractors, attention was drawn to attractors in different descriptions relevant for heavy ion collisions. For instance, similar attractors were found in kinetic theory and in holography, even in non-conformal and non-homogeneous settings [19,20], which indicates that attractor solutions should exist for a wide range of theories independent of the symmetries of those theories. In hindsight, this should not be so surprising, since the hydrodynamic limit is found in so many diverse physical settings. For instance, hydrodynamic attractors have also been explored in holographic models of QGP [21] and in Gauss-Bonnet gravity [22], which could have well been anticipated due to the fluid-gravity correspondence [23].
A key feature of the attractor is that when it is present, the system "forgets" about its initial conditions. Therefore, attractors can be used phenomenologically to distinguish between observables which are dependent on initial conditions and those which are not. This can lead to the identification of universal behavior in heavy ion collisions. Recently, the potential of hydrodynamic attractors as a tool for phenomenological studies has begun to be appreciated. For instance, they have been used to relate initial state energies to the produced particle multiplicity [24,25]. Moreover, recent studies examined the hydrodynamic attractor at finite chemical potential [26,27], tracing the hydrodynamic evolution as trajectories through the QCD phase diagram. One of the more exciting ideas relevant to heavy ion physics to come to fruition at least in part due to attractors is the Kurkela-Mazeliauskas-Paquet-Schlichtling-Teaney (KøMPøST) framework [28,29], essentially linking the early time dynamics, captured by a kinetic description, to the late time hydrodynamic behavior. Another practical application of hydrodynamic attractors is as a viable arena to study resurgence [30,31]. A comprehensive discussion of resurgence is beyond the scope of the present review.
Of course, it would be remiss to not mention that attractors have a rich, storied history in dynamical systems. There is considerable value in trying to understand the attractor studied in heavy ion collisions in the more mathematical language of dynamical systems. Progress in this direction has already been made, see [32][33][34].
A related story to the hydrodynamic attractor is the study of non-thermal fixed points [35][36][37] in heavy ion physics (as well as in a variety of other contexts [38]), which are essentially an observation that at early enough times, the distribution function exhibits a special scaling behavior. This early time regime is currently the subject of numerical study, such as in a simplified model of the glasma in 2 + 1 dimensions [39]. This should be considered as a distinct feature of the weakly coupled degrees of freedom at early times in the QGP, different from the hydrodynamic attractors that we focus on here.
Since hydrodynamic attractors have been found in a wide range of (effective) theories pertinent to heavy ion collisions, many works have focussed on comparing features of attractor behavior in relativistic hydrodynamics, kinetic theory and holography [20,40]. These theories vary not only in their frameworks, but also in their description of weak and strong coupling regimes. The benefits of this approach are clear, giving ample evidence of the importance of the attractors. This review will present in pedagogical detail the attractor story for the trio of theories mentioned previously, first by describing an illustrative example, followed by reviewing current progress in the field.
As such, the organization of the paper is as follows. In Sect. 2, we will discuss the prototypical example of the hydrodynamic attractor in the simplest setting of (0 + 1)dimensional Bjorken expansion, before reviewing further progress. Next, in Sect. 3, we will present the attractor as it pertains to holographic systems. Then, we will turn our attention to attractors in kinetic theory in Sect. 4, especially focussing on the relaxed time approximation. Finally, in Sect. 5, we turn to the recent applications of attractors to heavy ion collisions. A table of useful abbreviations is included in Appendix.

Hydrodynamic attractors
In this section, we will review the advances in studying the attractor in expanding relativistic hydrodynamics. As is usually the case, it is best to start simple and then discuss more complicated examples. We will refresh the reader with a discussion of relativistic hydrodynamics and various models of curing the acausality, before reaching the simplest and quintessential example of an attractor in 0 + 1d Müller-Israel-Stewart (MIS) Bjorken flow. We will finish by going into some detail of some advances in the field.

Primer on relativistic hydrodynamics
Hydrodynamics represents the course-grained description of microscopic dynamics. The macroscopic objects of interest are the energy momentum tensor, T μν , and the conserved currents, J μ . The dynamical equations are generated by considering the conservation of the energy momentum tensor, as well as the conservation of the currents, schematically ∇ μ J μ = 0. Generally, the fluid energy momentum tensor can be decomposed into ideal and dissipative contributions The ideal energy-momentum tensor is where ε is the energy density, P is the pressure and Δ μν = u μ u ν + g μν is the spatial projection operator. Note that the four-velocity, u μ , is timelike normalized, u μ u μ = −1. The dissipative tensor, Π μν , generically contains higher deriva-tive corrections of u μ and the temperature, T . This organization is known as the gradient expansion. We will ultimately need to specify a reference frame for the fluid. Here, we choose to work in the Landau frame, where the energy density is the eigenvalue associated with u μ via T μν u ν = −εu μ . This leads to the requirement that the dissipation tensor is orthogonal to the flow, i.e.
At this point, the standard approach would be to specify the dissipation tensor order by order in a gradient expansion with arbitrary transport coefficients. These transport coefficients are part of the required input for a complete hydrodynamic description, which can be determined from the underlying microscopic theory. To first order in derivatives, we have two distinct tensor structures, which we write as follows where the two transport coefficients are the shear viscosity, η, and the bulk viscosity, ζ . Higher orders lead to more possible tensor structures with additional transport coefficients. The shear tensor is transverse, traceless, and is explicitly given by The causal hydrodynamic evolution of a system can be described via (1) . There is one complication: the four equations in (1) contain five variables: the energy density, ε, the pressure, P, and the four velocity, u μ . As mentioned before, the last needs to satisfy the timelike normalization condition u μ u μ = −1, meaning that u μ has three independent components. Thus, we are shy one equation.
To fill this gap, we will need to specify an equation of state, relating the pressure to the energy density, P = P(ε). In general, the equation of state depends on the system that one is interested in modelling. One of the simplest choices (which we will largely adopt here) is the conformal equation of state ε = 3P, which means that the energy momentum tensor is traceless, In the conformal case, the energy density can only depend on temperature, which means that ε ∼ T 4 .

Truncation schemes
The issue with the choice of (5) for the dissipation tensor is that resulting equations become acausal [41,42]. A straightforward way to see this, e.g. [2], is by considering metric perturbations of the form, g μν = η μν + h μν , and seeing how the stress tensor (2) changes due to the metric perturbation in linear response. In essence, we are computing the retarded correlator, G μν,αβ R , via where the dots obscure terms higher order in perturbations. Focussing simply on the shear channel, we see that the retarded correlator is [1] G 0x,0x From the denominator, we identify the diffusive pole, ω = −i η ε+P k 2 . The group velocity, |dω/dk|, is given by which grows without bound as k gets large.
There are a few ways to cure such acausality. A simple and popular approach is to implement the MIS equations [43][44][45] (for more recent developments, see Baier-Romatschke-Son-Starinets-Stephanov (BRSSS) [46]), where the dissipation tensor satisfies a relaxation equation. Here, as we will be interested in the conformal case of non-zero shear viscosity, but vanishing bulk viscosity (as well as higher order transport coefficients), the MIS equation in this case is then where τ π is the relaxation time. The previous case (5) corresponds to vanishing relaxation time τ π (as well as ζ = 0). How does this cure the acausality? For small only timedependent perturbations, we see that which can be compared directly to (5) in the limit that ω, ζ → 0. This simply amounts to replacing η → η iτ π ω−1 . Revisiting the computation for the retarded correlator, we see that the diffusive mode is found when Solving for ω and computing the group velocity, we find one of the modes has dispersion The group velocity now saturates to a maximum value, namely lim k→∞ dω dk = η τ π (ε + P) (17) and as such, there is no issue with causality, so long as the right hand side of (17) is smaller than 1. A more recent improvement to the MIS equations is known as Denicol-Niemi-Molnar-Rischke (DNMR) effective theory [47], which is a second order hydrodynamic theory using Grad's 14-moment approximation for the single particle distribution function. To first order in gradients, DNMR is similar to MIS (13), in that (taking u μ = (1, 0, 0, 0)) we see that the general form the equations take to this order is [48] where for DNMR we have β ππ = 38/21 and τ π = τ eq [47], while for MIS this is β ππ = 4/3 and τ π = 6τ eq /5 [49]. Another truncation scheme goes under the name of anisotropic hydrodynamics (aHydro) [50][51][52][53]. Essentially, the idea is to formulate hydrodynamics in momentum space with a view towards taking into account the large anisotropies found in the QGP. One takes a distribution function with an explicit anisotropy parameter, ξ (see Sect. 4 for a primer on kinetic theory and (163) for the form of the anisotropic distribution function). As outlined in Sect. 4, the Landau matching condition determines the temperature. The evolution of ξ is determined by considering a higher moment of the Boltzmann equation. In this way, aHydro resums an infinite number of terms in an expansion in inverse Reynolds number [48] whereΠ = Π/ε. A comparison to the other truncation schemes can be made by rewriting the evolution equation for ξ as a function of Π and expanding for small anisotropy. This results in agreement with (18) to linear order in Π . However, note that at higher order Π 2 , the term λ 1 Π 2 /(2τ π η 2 ) has a numerical disagreement of the transport coefficient, λ aHydro 1 = ητ π /7, with λ DNMR 1 = 5λ aHydro 1 [48,54]. The differences between the two approaches arises from the fact that aHydro includes higher order terms compared to DNMR at the same level of truncation.
Which scheme to use? In terms of conceptual simplicity, the MIS framework tends to get more mileage since it is rather straightforward to include additional higher order transport coefficients. MIS and DNMR differ mainly in their far-from-equilibrium transient effects, as they have the same late time Navier-Stokes limit [27]. It was found in [48] that DNMR is more accurate than MIS when compared to the exact solution in the relaxed time approximation (RTA), as it captures systematically second order hydrodynamic contributions, while aHydro approximates the exact attractor the best. Likewise, in the context of Gubser flow where the exact attractor is known, aHydro did a better job in capturing the attractor compared to MIS and DNMR [32,55].
As an aside, it is worthwhile to underline that the MIS equations represent a particularly easy model to work with. However, if one were inclined to include non-hydrodynamic modes to their model building, such as the oscillatory behavior of N = 4 super Yang-Mills theory via a hydrodynamic model, the Heller-Janik-Spalinski-Witaszczyk (HJSW) model could be considered a natural generalization to the MIS equations [56][57][58]. In this case, the dissipation tensor, Π μν total = Π μν MIS + Π μν HSJW , can be decomposed into a piece evolving with the MIS equations (13) and another piece satisfying 1 where Ω = Ω R + iΩ I is the complex frequency, meant to mimic the quasinormal mode frequency. The values of the real and imaginary part can be fixed by comparing to N = 4 Super Yang-Mills (SYM) theories [57,59]. The presence of these non-hydrodynamic modes is guaranteed by the appearance of the second derivative in the above. In principle, one can further generalize by considering higher derivatives on the left hand side of (20), although how to specify the extra initial conditions is not clear.

0 + 1d Bjorken flow
Now let's get our hands dirty with an example computation. We will review the hydrodynamic attractor in Bjorken flow [16], a (0 + 1)-dimensional expanding transversally homogeneous system with longitudinal boost-invariance. Bjorken flow [60] describes a simple model for heavy ion collisions, where the collision axis is along the z-axis and the system is otherwise homogeneous. Essentially, the nuclei are modelled as infinite flat sheets in the transverse (x, y) directions. Furthermore, the matter produced in the forward light cone is boost invariant. The background metric is given in Milne coordinates We can relate Milne coordinates (τ, x, y, ξ) to Minkowski coordinates (t, x, y, z) by observing that the proper time and the spacetime rapidity are respectively. Note that the metric (21) describes a space which is Ricci flat. For future convenience, we also note that the nonzero Christoffel symbols are We will assume all variables depend only on the proper time, τ , which means that (1) and (13) take a particularly simple form: where φ ≡ −Π ξ ξ . Thus, the conservation equation, (25), and the MIS equation, (26), represent the set of equations describing the relativistic evolution of a viscous fluid in a Bjorken expanding background. We proceed by expressing the shear viscosity and the relaxation time in terms of the energy density: where C η and C τ Π are dimensionless constants. As we will see shortly, the relaxation time is the hydrodynamic transport coefficient which indicates the decay time of the nonhydrodynamic mode. Taking insight from holography [46], we will take these constants to be We can then eliminate φ to arrive to arrive at an expression solely for the energy density We can recover the expressions found in [16] by remembering that in a conformal theory the energy density is related to the temperature via ε ∼ T 4 .
Following [16,17], we can now introduce dimensionless quantities which leads to the following equation for f = f (ω) We can now numerically solve the above equation for a variety of initial conditions, which we plot in Fig. 1. Irrespective of the choice of initial conditions, the large ω behavior converges to a single curve, which is known as the hydrodynamic attractor. In this concrete sense, the attractor is said to "forget" about its initial conditions. The approach to the attractor can be found if we solve the linearized form of (31) for δ f . We find that From here it clear to see that the perturbation decays exponentially to the attractor. Note that another often used dimensionless quantity (see e.g. [15]), known as the pressure anisotropy, is defined as  [16] which is normalized with respect to the equilibrium pressure, P = E/3. The two dimensionless quantities are related simply via We can now examine the hydrodynamic gradient expansion as an asymptotic expansion in ω of (31), given by It is straightforward to see that the resulting series has terms growing factorially [18], indicating that the series is indeed divergent, see Fig. 2.
It is worthwhile to point out that the choice of variables T (τ ) ↔ f (ω) has nothing to do with the appearance of the attractor, which is inherent to the system. For the case considered above, the attractor can alternatively be captured in a study of the phase space [58], by studying for instance the evolution of (τ 0 T (τ ), τ 2 0Ṫ (τ )) as a function of τ . As food for thought, the evolution of a set of initial conditions in phase space can be viewed as dynamical dimensional reduction of the phase space regions. This might be a helpful way to think about the attractor in more realistic scenarios, e.g. in the case of relaxed symmetries.
Another method to characterize the attractor is known as the slow-roll approximation, originally studied in the cosmological setting [61]. Formally, it amounts to expanding (31) in derivatives, by adding a small parameter, δ, to the f f term [48]. One then expands (31) in powers of δ This is a useful computational technique, as at a given order, it contains derivative terms to all orders, although it should be noted that the slow-roll approximation diverges [62]. The Bjorken case is one of the more simple examples due to the assumed symmetries and low dimensionality. The discussion presented lends itself naturally to generalizations, which we will review here after a brief aside.

Quick introduction to trans-series
Trans-series are an interesting mathematical development with real physics application. In some sense, they represent a generalization to Taylor series, taking into account not just the polynomial powers, but also the non-perturbative structure, such as exponentials and logarithms. A related topic, resurgence, is the curing of a perturbatively divergent series with non-perturbative terms [63]. Resurgence can be used to deduce nonperturbative effects via analytic continuation from a perturbative series. For an introduction to trans-series with a more mathematical flavor, see [64], and for more on resurgence, see for example [31].
Here, we will follow the discussion of trans-series relevant to hydrodynamics [30], noting that the transseries of N = 4 SYM in Bjorken flow have been explored in [65].
As we see in (32), there is the non-hydrodynamic mode, exponentially suppressed at late times, analogous to quasinormal modes in gravitational theory [57] and to instantons in quantum field theories [66] . It becomes clear that instead of the naive expansion (35), one should seek a trans-series ansatz of the form where σ is the trans-series expansion parameter, which counts the "non-perturbative" order of the trans-series and is generally complex and is the so-called "instanton fugacity". Plugging this ansatz into (31) leads to equations at each order in σ. For n = 0, we have the previous case discussed, namely the naive Taylor series. There is only freedom to fix one parameter, f (1) 0 = 1, and all coefficients are then determined by their recursion relation, which is a generic feature of trans-series [67]. The resulting non-hydrodynamic terms in the trans-series correspond to exponentially damped modes times a gradient expansion. It is also worthwhile to point out that the overall constant is known as the Stokes constant, which when fixed resolves the ambiguity in the choice of contour when doing the inverse Borel transform [22] (note that in the example we consider in Sect. 2.5, we would have two Stokes coefficients, while they are infinite for holography).
Do trans-series tell us anything about early times, i.e. when ω is small? The exponentially suppressed terms at asymptotically late times can become large compared to the powers of ω for early times, with such terms jockeying for position in size. Trans-asymptotic rearrangement rewrites the transseries to better study the early time behavior [68]. For the example considered in [30], we should rewrite the expansion (38) as where Inserting the reordered ansatz into (31), one can solve order by order, finding that for k > 1 each equation contains just a linear term in F k . The imaginary part of the parameter σ is fixed in terms of the Stokes parameter by the requirement that f (ω) is real. The real part of σ is arbitrary [69], thus parameterizing the different possible initial conditions at early times.

Gubser flow
A straightforward generalization of Bjorken flow was introduced by Gubser in [70], which has been studied at length in various settings. Like Bjorken flow, Gubser flow is boostinvariant with the additional aspect that there is azimuthally symmetric radial expansion. In fact, a radially expanding Bjorken flow can be mapped to Gubser flow via an appropriate Weyl rescaling. It can be thought of as a comoving flow in conformal theories. The Gubser flow attractor has been widely studied as it is closer to the set up relevant in heavy ion collisions [32,34,[71][72][73][74].
In terms of the de Sitter coordinates (ρ, θ, φ, η), the Gubser flow metric is which relates to the Milne coordinates via where q −1 is the characteristic length scale of the transverse direction of the system. The time variable is −∞ < ρ < ∞.
That the system is expanding is straightforward to see by computing One then finds that the conservation of the energy momentum tensor, ∇ μ T μν = 0, satisfies a single equation, namely [32,75] where the normalized shear stress is π ≡ −π ηη /(ε + P) (note that this parameterization is different from φ introduced below (26)). As mentioned above, the evolution of π can be specified in a number of ways, the simplest of which is arguably MIS (13), which in Gubser coordinates reads [32] Using conformal symmetry, we can write τ π = c/T with c = 5η/s. As in the Bjorken case, the two equations can be combined into one equation by solving for π in (45). Similar to the discussion in the Sect. 2.2, a compact way to write the above equations is by introducing the dimensionless time ω = tanh ρ/T and the anisotropy parameter [32] which is related via the conservation equations to the normalized shear stress as π = 3A + 2. We thus have As a useful aside, note that DNMR in Gubser flow yield similar evolution equation for A(ω) [32] 3ω( The qualitative differences between the previously discussed Bjorken flow and Gubser flow were established in a number of studies, although it is worth pointing out that the gradient expansion of Gubser flow was found to diverge factorially similarly to the Bjorken case [72]. One of the major distinctions between Gubser flow and the previously discussed case of Bjorken flow is the rich range of behaviors Gubser flow experiences. During its evolution, a system with Gubser dynamics undergoes early time free-streaming, progressing to intermediate time hydrodynamization and finally at late times free streaming [74]. The final stage is absent from 1 + 1d Bjorken flow, as a system with Gubser flow expands more rapidly than in a Bjorken flow and interactions cannot keep up. In examining the local entropy current and entropy density for the two flows [73,76], it was found that Bjorken entropy eventually saturates, while for Gubser flow, the entropy continues to grow for late times, marking the ever-growing departure from thermal equilibrium of the system.
Seeing how Gubser flow has more features closer to real heavy ion collisions compared to Bjorken flow, such as transverse expansion, one can speculate that with further relaxed symmetries and the inclusion of transverse dynamics, the attractor story will more resemble Gubser flow.

Phase transitions
Since hydrodynamic attractors seem to play a role in simpler models of heavy ion collisions, it would be prudent to study the interplay between the attractor and other known features of the QGP, such as spontaneous symmetry breaking. A step in this direction was recently explored in [77] by studying superfluid Bjorken flow, where a conformal fluid in a MIS framework was coupled to a U (1) scalar field undergoing a phase transition via a temperature dependent mass. In this case, the action reads where Σ = ρe iψ is the complex scalar field, ρ is the modulus, ψ is the phase, the gauge covariant derivative is and the metric is the usual Bjorken metric (21). The phase transition is characterized by parameterizing the temperature-dependent mass via m 2 ∼ (T − T c ). The complete set of equations arise from the energy momentum tensor conservation, the conservation of the U (1) current, the scalar equation of motion and the MIS equations. Note that the energy momentum tensor and the current can be decomposed into a sum of normal and superfluid components, e.g.
As a momentary aside, it is worthwhile to show that the effective fluid action [78] indeed generates the hydrodynamic equations via diffeomorphism invariance. We define the temperature as β μ = 1 T u μ and the four-velocity, u μ , is timelike normalized, u μ u μ = −1. If we vary (53), using the thermodynamic identity involving the entropy, s = ∂ p ∂ T , and ε + p = sT , we find where we recognize the ideal fluid energy momentum tensor in the last line Under diffeomorphisms x μ → x μ + ξ μ , we have that which directly leads to the conservation of the stress tensor after partial integration, ∇ μ T μν = 0. Since we recover the usual ideal hydrodynamic equations from the variation of (53), we will refer to (53) as the fluid action. Dissipation is added to the energy momentum tensor, T μν = T μν ideal + π μν , and the current j μ = j μ ideal + q μ , with relaxation akin to the MIS formulation, where τ π and τ q are the relaxation times. Furthermore, the equation for the condensate, ρ, and the phase, ψ, also have added noise We now determine the form of the dissipation by requiring positive entropy production. The entropy is given by s = (ε + P − μN )/T , where μ is the chemical potential. The divergence of the entropy current reads Enforcing positive entropy production leads to the following identification of dissipative structures where η, ζ, κ, κ 1 and κ 2 are positive functions of T and μ. For a conformal system, p ∼ T 4 , the bulk viscosity vanishes, ζ = 0, and the transport coefficients can be parameterized via Defining the pressure anisotropy χ := −3π η η /(4T 4 ) and setting ψ = 0 as well as μ = 0, for simplicity, the equations of motion arë τṪ The relatively simple system described by (50) and the above equations turns out to have a rich set of behaviors. Starting in the unbroken phase with T > T c , the system approaches the usual hydrodynamic attractor (T ∝ τ −1/3 ), as can be seen in Fig. 3. During this time, the condensate, ρ ≡ |Σ|, value is (close to) zero. Then, some time later, the system reheats and exponentially quickly develops a non-zero value for the condensate. This corresponds to a symmetry-breaking fixed point of the system. The fixed points at (±ρ , t ≡ T /T c ) can be computed directly from the equations of motion by recognizing that there should be no dissipation at the fixed point, i.e. χ = 0, andρ = 0, which means the equation of motion for ρ leads to Plugging this into the equation for the temperature results in which in the range 1 ≥ t ≥ 0 always has a real solution, irrespective of the choice of constants.
Of course, to make contact to phase transitions found in the evolution of the QGP, like the chiral phase transition, one would need to include the full SU (2) L ×SU (2) R → SU (2) V symmetry breaking pattern, following the set-up studied in [79]. As such, this study serves as a playground for exploring further interplay between hydrodynamic attractors and other features relevant in the QGP.

Attractor of hybrid models
The QGP has a rich structure of strongly and weakly coupled degrees of freedom, interacting in a range of energies. To study the interplay between these two sectors, as a first step it would be worthwhile to study a hybrid fluid model in a Bjorken background [80], where the two fluids represent the hydrodynamic modes of the hard and soft sectors. The central idea is to see the effect of coupling two viscous fluids which naturally exhibit attractor behavior to one another. A priori, it is not clear how the full system will behave or whether a new attractor will emerge. The discussion here will deviate from [80], but arrive at the same coupling equations studied in the Minkowski background of [81].
We can study the coupling of two fluids via the fluid action discussed previously. The action of the first fluid is given by Fig. 3 Evolution of the system for a variety of initial conditions, which are shown in green. The inset shows the approach of the system to the usual hydrodynamic attractor for certain initial conditions, before veering away and approaching one of two possible fixed points at the end of the system's evolution, shown in red and cyan. Figure from [77] (53), while the second fluid is What is the combined pressure of these two subsystems, (53) and (73), if they are brought in contact? Taking inspiration from Dalton's law [82], which states that the total pressure is given by the sum of partial pressures for fixed volume and temperature, i.e.
as well as Agamat's law [83], which states that the total volume is given by the sum of partial volumes with fixed pressure and temperature, i.e.
we make the ansatz that the total action is where the combined action is given by the total pressure Note that if the effective metrics, g μν andg μν , are the same as the physical background metric, G μν , we then have that the total pressure as just the sum of the individual pressures, which is just a statement of Dalton's law (74). If the pressures are constant and equal, then the integral of the volume elements can be identified as the partial volumes (i.e. d d x −g =Ṽ ), whose sum is the total volume of the system, which is Agamat's law (75). As a perhaps useful aside, there seem to be too many metrics. As will be shown later, the effective metrics, g μν and g μν , encode the interaction between the two subsystems. The physical background metric, G μν , is the metric that could be made dynamical à la [84] via Here, we will consider the background metric G μν to not have any dynamics, ultimately choosing the Minkowski metric as the background metric, G μν = η μν . We will assume that each subsystem is in local thermal equilibrium with a temperature (Killing) vectorβ μ = 1 T 2ũ μ [45], which means the temperature is given by T 2 = (−β μg μνβ ν ) −1/2 . This leads to the following identification when we are in hydrostatic equilibrium [85] T We will refer to T as the physical temperature, while T 1 and T 2 will be referred to as the effective temperatures of each subsystem.
We can use the thermodynamic identity dP = SdT (and also dP = s 1 dT 1 and dP = s 2 dT 2 ) to compute the total entropy density. We thus have which agrees with the result of the entropy in [81]. Thus, the entropy is essentially additive (upto interactions hidden in the effective metric determinants).
Finally, we point out that the two fluids are coupled only through the effective metrics, i.e. varying with respect to the effective metrics, g μν andg μν , the microscopic detail is given by where the covariant derivative is defined with respect to each subsystem, i.e.
and similarly for the tilded equations. Thus, one fluid only feels the other through the effective metric. In other words, the interactions are encoded by the effective metrics.
We can now compute the full energy-momentum tensor of (76), where we find We see that the change in the effective metrics, g μν andg μν , encompasses the interactions between the two fluids. We introduce interactions by deforming the effective metrics in such a way that the total energy momentum tensor, T μν , is conserved in the physical background, G μν , i.e. ∇ (B) μ T μν = 0. First, we mention that it is straightforward to verify that (84) can be rewritten as and similarly for the other subsystem. Then, using the previous identity, we see that the two subsystems are coupled via the lowest order coupling rules where we use the shorthand t = t μν G μν andt =t μν G μν . Note that the tensor coupling, γ , and trace coupling, γ , have mass dimension −4. Higher order couplings, like O(γ 2 , γ γ , (γ ) 2 ), can be included in a straightforward manner, see [81] for a general discussion. Using (86) and (87), we can read off the full energy momentum tensor which is conserved in the background metric: where t ·t = t μν G ναt αβ G βμ . Note that stress tensors have their indices raised/lowered by their effective metric, e.g. t μ ν = t μα g αν . The coupling rules (87) are chosen in this way to ensure that the full system has a conserved energymomentum tensor in the physical background, G μν = η μν , namely To solve the coupling equations (87), we make the ansatz that the effective metrics will take the diagonal form, moti- Fig. 4 Anisotropy as a function of proper time. Note that χ = A/6, where A is the anisotropy parameter defined in (33). The thick lines denote the attractor, while the thin lines denote solutions with different initial conditions. The more viscous system is in red, while the blue lines are the less viscous system. Finally the orange line denotes the emergent system attractor. Figure from [80] vated by the diagonal form of the background metric, G μν = η μν : Likewise, the two stress tensors are assumed to take the form and similarly for the other, tilded subsystem. Π μν is taken to evolve using the MIS equation (13) in each sector with each system having their own shear viscosities, η,η. To distinguish the fluids and in mind to describe the weak/strong coupling sectors found in the QGP, we set the tilded subsystem to be more strongly coupled to have Cη = 1/4π to have a smaller shear viscosity than the more weakly coupled untilded sector Thus the untilded weakly coupled sector can be thought of as governed by kinetic theory, while the tilded strongly coupled sector is governed by a holography.
It is worth noting that in the decoupling limit, γ, γ → 0, each subsystem contains an attractor. As can be seen in Fig. 4, with interactions turned on, the combined system surprisingly has an emergent attractor as well. This emergent attractor is parameterized by the energy densities of each subsystem, γ ε and γε. Interestingly, the scenario modelled here is reminiscent of bottom-up thermalization [86], in that at early times the strongly coupled subsystem is dominant, while at later times it is the weakly coupled subsystem which dominates. Intriguingly, though the individual subsystems were taken to be conformal (t μν g μν = 0 =t μνg μν ), the full system is in general not, T μν G μν = 0, becoming conformal at late times.
The hydrodynamization time for each subsystem can be defined (see [87,88] for similar definitions) when the pressure is sufficiently close to first order hydrodynamics, namely and similarly for the other sector. It was found in a variety of cases that the hydrodynamization time for the harder sector is longer than for the strongly interacting sector. One can compare the ratio of the hard to soft hydrodynamization times via the ratio R hd = τ hd /τ hd , where it is clear that depending on the initial conditions, R hd can vary by an order of magnitude, perhaps indicating that the hydrodynamization scenario in small-system collisions [89][90][91] may be quite different to their heavier ion counterpart. The hybrid coupling described in this section lends itself naturally as a framework to study more nuanced examples, closer to the physics of the QGP. For instance, instead of fluids, one can consider the coupling between an ultraviolet sector described by kinetic theory, while the low energy degrees of freedom are described holographically. Such a framework is known as semiholography, a term coined in [92] and introduced to heavy ion collisions in [93]. It would be interesting to see the behavior of the attractor with the democratic metric coupling in addition to the scalar coupling between sectors [94][95][96], which drives the energy from the hard to the soft sector.

Holographic attractor
The connection between hydrodynamic and holographic theories is firmly established, e.g. see [97] for a discussion bridging BRSSS hydrodynamics to a holographic theory. Moreover, boost-invariant holography had been studied by numerous groups over the years [17,98,99], going back to the linearized analysis of [100]. The holographic setting was also where the Borel resummed hydrodynamic gradient expansion was first considered [17,18]. As such, after finding the hydrodynamic attractor [16], the stage was set to explore the attractor story in the strong coupling context. An important early work in this direction was due to Romatschke [101], whose example we will reproduce and expand on in Sect. 3.1. We will then turn our attention to some of the novel directions explored in holographic attactors.

An illustrative example
Here we will follow [101]. In the holographic context, one aims to solve a classical gravity problem which is dual to a strongly coupled quantum field theory. We will want this dual theory to be comparable to the hydrodynamic theory. The dual theory is encoded in the boundary conditions of the gravity equations, which we will outline shortly. We begin with the vacuum Einstein equations where R μν is the Ricci tensor and R ≡ R μν g μν is the Ricci scalar. The cosmological constant in AdS is given by Λ = −(d − 1)(d − 2)/2L 2 = −6, since we work in d = 5 dimensions and with the AdS radius, L, set to unity.
The metric in Eddington-Finkelstein coordinates is given by where ξ is the spacetime rapidity. We take the functions A, Σ and B to be functions of the proper time, τ , and the holographic radius, r . With this metric ansatz, the Einstein equations (96) can be helpfully written as a set of nested ordinary differential equations [102] where the prime denotes the radial derivative, h ≡ ∂ r h, and the dot denotes the radial derivative along the null geodesiċ Note that (101) and (102) are constraint equations, while the other equations are to be solved numerically.
Near the boundary at r = ∞, the asymptotic expansion of the metric functions is given by Σ(τ, r ) = r ∞ n=0 a n (τ )r −n .
We want the metric at the boundary to take the form of the Bjorken metric, (21), in order for the dual theory to exhibit boost invariance. This implies that the boundary conditions at r → ∞ are This determines the leading coefficients in the expansion of the metric functions, a 0 = 1 and s 0 = 1. Furthermore, it is easy to see that the metric (97) remains invariant under the diffeomorphism, r → r + f (τ ). This residual gauge freedom can be fixed by setting a 1 = 0.
The boundary stress tensor is related to the metric coefficient a 4 via holographic renormalization [103,104] by The Einstein equations are then solved using numerical techniques due to Chesler and Yaffe (namely, the spectral method) [102,105,106]. There is also publicly available code in [107,108].
Various initial conditions at τ = τ 0 were generated in [101] by taking Σ to be where Then choosing the profile ε ∼ τ α , one can generate initial data by changing α andc.
Sorting through the initial data, there is a clear late-time attractor as can be seen in Fig. 5. However, at earlier times the story is not so clear due to the well-known oscillatory nature of N = 4 SYM, which arises primarily from the existence of non-hydrodynamic modes at early times [101,109].
To identify the holographic attractor, it was posited in [101] that the relevant condition should be that the attractor experiences little oscillatory behavior and whose initial value tends to The attractor is then visible in the rightmost panel of Fig. 5 as a solid black line, while the various initial conditions produce a variety of oscillatory behavior. Such an ambiguity is not found in relativistic hydrodynamics and kinetic theory, where the approach to the attractor tends to be monotonic from a given set of initial conditions. In contrast to the Romatschke approach, another study was undertaken in [40], where the initial conditions of the anisotropic part of the metric was modelled instead in two ways (note: ρ = 1/r ) This particular choice of initial conditions was chosen such that the initial ratio of the longitudinal pressure to the energy density p L (τ 0 )/ε(τ 0 ) and its first derivative are equal. The meaningful difference between the two families of initial ior should approach the hydrodynamics, which are denoted by the dashed lines. Note that the agreement between the attractor and hydrodynamics at late times occur at different timescales, ω = τ T , for the different theories. Figure from [101] conditions is in their extent in the holographic direction, with the UV solution localized almost entirely on the boundary. Both sets of initial conditions display the expected oscillatory behavior due to the expected decay of quasinormal modes. Intriguingly, the more localized UV initial condition seems to approach the hydrodynamic behavior sooner than the corresponding IR family of initial conditions. For either family of initial conditions, the approach to the attractor takes longer than the corresponding set-up in kinetic theory or IS hydrodynamics, irrespective of initialization time. Thus, results of [40] seem to indicate that the attractor only exists in a meaningful way at late times for holographic set-ups with the early time dynamics dominated by transients, which is further supported by related work of higher-order fluid dynamics [110].

Borel summation
Whether the condition (116) is the "best" or unambigous identifier for a holographic attractor is up for debate. A perhaps more rigorous approach is due to [109,111]. Using the same metric as in (97), the hydrodynamic gradient expansion of the energy density in terms of the proper time was computed up to 240 th order [18]. The longitudinal and transverse pressure in the Bjorken case is given by (110) and (111), so that one can easily compute the anisotropy parameter (33), which has an expansion in terms of ω = T τ : Thus, the expansion of the energy density in terms of proper time (118) can be re-expressed as an expansion in the anisotropy parameter in terms of the dimensionless time, ε(τ ) ↔ A(ω). One then computes the Borel transform of the original power series, which is given by The intent of the Borel transform is to provide a method of handling divergent series. The gradient expansion (119) is a divergent expansion with larger order terms diverging factorially [16]. One can then perform the Borel resummation via the inverse Borel transform where the contour, C, connects x = 0 to x = ∞ in the complex plane. Practically, since only 240 terms are known due to the high technical cost of computing the coefficients, the analytic continuation is done using Padé approximants. The Padé approximants are, in some respects, a generalization of the Taylor series, with the improvement that singularities are explicitly included in the approximant. This is done by considering a ratio of polynomials, i.e. we can approximate a function via where the denominator can better capture the singularity structure of the function that is being approximated, compared to Taylor series, which does not contain such terms. For more on Padé approximants, see [112,113] for application to QCD sum rules. Does this method provide a good approximation to the attractor? In the case the attractor is known, such as in the previously mentioned HSJW model (20), the known attractor can be compared to the result of the integrating the series truncated at 240th order via (121), where there was quantitative agreement down to about ω ≈ 0.3, indicating the method provides meaningful results. In the case of holography, where the exact attractor is not known, it was found in [109] that there is good agreement with [101] up to ω ∼ 0.4, whereas the qualitative behavior is different for smaller ω.

Violation of energy conditions
As the main workhorse for the study of holographic attractors and the go-to holographic toy model for heavy ion collisions, it is important to understand the holographic Bjorken model in Sect. 3.1. For this reason, it is interesting to note the holographic Bjorken flow can evolve to violate the dominant and weak energy conditions [114], even with initial conditions that satisfy these bounds. For conformal Bjorken flow, the weak energy condition was given in the first holographic study of Bjorken flow [115] where t μ is any time-like vector, which implies for Bjorken flow Note that since we are considering a conformal system, the strong energy condition is equivalent to the weak energy condition. The dominant energy condition implies that It was found that the violations occur at times around τ ≈ 0.5 fm/c, which is within the prehydrodynamic regime in the QGP [28]. This can have implications for the validity of kinetic theory, the typical effective description used in this regime. Kinetic theory satisfies the energy conditions: the weak energy condition implies that the energy density is non-negative and the dominant energy condition is the requirement that matter does not exceed the speed of light [116]. If these violations were shown to exist in the energy momentum tensor during nuclear collisions, this could indicate that this transient regime would lie outside of the domain of validity of kinetic theory descriptions [114], indicating that strongly-coupled hydrodynamization as predicted by holography would be the preferred description.

Gauss-Bonnet attractor
Attractors have been studied in Gauss-Bonnet (GB) holography [22,117,118]. The key advantage of studying GB holographic systems is that one has access to intermediate coupling (in addition to infinitely strong and weak coupling) regimes. The only parameter is λ G B , which controls the size of the higher derivative corrections. The dual field theory is unknown, but finite and negative values of λ G B are at least qualitatively similar to what one would expect from finite coupling corrections to N = 4 SYM. The negative values of λ G B do lead to violations of causality in the ultraviolet [119], but the interest of these works is in the direction of the infrared dynamics and the system's approach to hydrodynamics.
Here, we will briefly outline some details of the construction and the lessons learnt. The GB action is where κ 5 is related to Newton's constant in 5 dimensions and L is the AdS radius for the λ G B = 0 theory. Note that the λ G B = 0 case corresponds to the example in Sect. 3.1 with the terms proportional to λ G B representing higher derivatives of the metric. The shear viscosity of this theory is given by Using a similar metric ansatz of the Eddington-Finkelstein form (97), one then solves the Einstein equations as outlined in Sect. 3.1.
For the values of λ G B studied in [22], the attractor is estimated using Borel summation as outlined in Sect. 3.1.1. It was found that the relaxation time of the non-hydrodynamic modes occurs when the gradients of the system are still large, i.e. for anisotropy parameter > 1, with first order viscous hydrodynamics comparable to the estimated attractor. Interestingly, it is found that by varying λ G B there is a qualitative interpolation in the singularities found in the complex plane between the holographic infinite coupling limit and the weakly coupled limit in kinetic theory, which can be inferred from Fig. 6.

D-dimensional Bjorken flow
The attractor in arbitrary D = n + 1 dimensions was studied in [120]. For a review of large D gravity, see [121]. The essential technical difference is that the metric (97) has the same form and dependence, but there are more transverse directions, i.e.
which is chosen such that the determinant of the term in brackets is equal to one. Operationally, the D−dimensional Einstein equations are rewritten in the same way as in (98)- (102). The equations are then expanded in a large n = D − 1 limit up to n −3 order. The key ingredient to compare holographic theories in different dimensions is the requirement that the late time behavior is the same, which is enforced by introducing the late time scale, Λ, which arises from considering the late time behavior of the ideal Bjorken expansion Interestingly in the case studied in [120], at the time scale τ ∼ 1/Λ the decoupling of non-hydrodynamic modes from the hydrodynamic modes occurs in an average sense: these modes are rapidly oscillating around the hydrodynamized evolution, averaging to zero over the typical expansion time scale. Coupled with the fact that the non-hydrodynamic excitations have amplitudes which are not parametrically suppressed, it seems like the typical approach to the attractor would be different to the usual D = 5 case.

Attractor in kinetic theory
In heavy ion collisions, the conceptual picture is that kinetic theory tracks the evolution of (quasi-)particles. The applicability of the theory is early times post collision up to the time when the build-up of the soft gluon bath is significant [122]. This corresponds to roughly the hydrodynamization time, τ ∼ 1 fm/c. The bridging between regimes where the dominant effective description is kinetic theory to hydrodynamics has been developed in KøMPøST [28,29].
As a workable theory when the heavy ion system is quite far-from-equilibrium, it is interesting to study which role attractors play in a kinetic description. To that end, we present the reader with an example kinetic computation in Sect. 4.1, before describing further progress in the following subsections.

Kinetic theory primer
Here we provide a short recap of key concepts regarding kinetic theory, following the discussion in [123]. Kinetic theory concerns itself with the microscopic detail of particles directly. The fundamental object of interest in kinetic theory, keeping track of the particles' positions and momenta, is the one particle distribution function, f = f (x μ , p μ ), where x μ are the spacetime coordinates and p μ is the four-momentum. Note that the distribution function should not be confused with the dimensionless quantity (30), used to characterize the hydrodynamic attractor. The evolution of the distribution function in phase space follows Liouville's theorem [124], namely where λ is an affine parameter. This equation is commonly known as the collisionless Boltzmann equation. It is common to consider the Boltzmann equation on-shell, e.g. p μ p μ = −m 2 in our sign convention. This means that the one of the components of the four-momentum is determined in terms of the other momenta with a common choice being p 0 = p 0 ( p i ). For a massive gas in a Minkowski background, this Macroscopic quantities can be computed straightforwardly by considering the moments of the distribution functions, i.e. integrating the distribution function with increasing powers of the four-momentum over momentum. For instance, the particle number current and energy momentum tensor are given by [123,125] respectively.
Noting that the four-momentum is p μ = dx μ dλ and the fourforce is F μ = dp μ dλ , we find that the distribution function evolves according to the Boltzmann equation in a curved spacetime in the absence of external forces where C[ f ] is the collision kernel, which encompasses all order scattering processes. In practice, it is more expedient to work with a simpler, truncated collision kernel. For analytic study, a helpful choice of the collision kernel is the relaxation time approximation (RTA) (also known as the Bhatnagar-Gross-Krook (BGK) approximation [126]), where the collision term is given by where τ R is the relaxation time, which in principle can depend on the proper time, τ R = τ R (τ ). Here, we will take the equilibrium distribution function to be given by where T is the temperature and will be a function of proper time, τ . The interpretation of the RTA collision term is that the system relaxes back to its equilibrium value, which can be easily seen in the following example. Assume that we are in flat Minkowski spacetime with only time dependence in the distribution function. Taking f eq to be a constant and shifting f − f eq → f to write the Boltzmann equation (134) in the RTA (135) as which has solution after shifting back as Thus, we see that the distribution function exponentially rapidly approaches its equilibrium value in a timescale, τ R . An important ingredient in the RTA is what is known as the matching condition. Since the energy momentum tensor (133) is conserved, ∇ μ T μν = 0, the first moment of the Boltzmann equation must vanish. This sets a condition on the collision kernel, namely that the moment of the collision kernel should vanish, leading to the dynamical matching conditions. In other words, the energy density, ε, should match the equilibrium energy density, ε eq , i.e.

Attractor via the method of moments
A straightforward way of seeing the attractor in kinetic theory is using the method of moments. In the context of attractors, the tower of moments was studied in the RTA [40,127] as well as with more realistic collision kernels [128], following prior work in studying the divergence of the Chapman-Enskog expansion [129]. Moreover, for an application of studying the fixed points, e.g. free streaming and hydrodynamic, found in kinetic theory via the method of moments, see [130][131][132].
Here we outline the basics of the construction. Instead of working with the distribution function directly, it is helpful to consider the generalized moments, ρ, of the distribution function [129] The moments contain all of the microscopic information of the distribution function. For instance, it is clear by inspection that the lower moments have a straightforward physical interpretation, namely the energy density is given by ε = ρ 1,0 and the longitudinal momentum is given by p L = ρ 1,1 .
The evolution equation for the moments is found by integrating the Boltzmann equation. One sees that where we integrated by parts and used that in the Bjorken massless case k 0 = k 2 x + k 2 y + k 2 η /τ 2 . Integrating over the equilibrium distribution function leads to ρ eq n,l = Putting together the pieces, we find that the general equation for moments is given by In the case considered by [40], namely the set of lowest lying equations (n, l) = (1, 0) and (n, l) = (1, 1), one finds the following set of equations As a reminder, the relaxation time is in principle a function of proper time. As above, we will consider the conformal case, where τ −1 R ∝ ε 1/4 . We will use the rescaled time t ≡ τ/τ R (not to be confused with the Minkowski time, t). A helpful identity to be used below is that where we introduced x ≡ p 1 /ε. Finally, calling y(t) ≡ p 2 /ε, we see that (146) and (147) can be combined in one equation as [40] In analyzing (149), a useful starting point is to consider is at t = 0. Requiring the derivative of x is finite at t = 0, we see that we have two cases At late times, it is clear that (149) becomes from which it follows that the generic late-time behavior is lim t→∞ x(t) → 1 3 . Depending on the initialization time, the approach to the attractor, which we call x A , has two qualitatively different regimes, as can be seen in Fig. 7. When the initial time is chosen to be t 0 1, the attractor is approached exponentially quickly, while for t 0 1 the approach to the attractor follows approximately a power law decay, x − x A ∼ t −8/3 . Similar behavior, namely a transition from power law to exponential behavior, is also seen in the case of when an approximate analytic attractor is known [133].
It is worth bearing in mind that the number of moments is infinite. It is beneficial to understand how the attractor might appear at higher moments in the distribution function. Work in the Bjorken flow in the RTA [127,134] indicates that the higher moments all approach the attractor at late times, although higher moments with l = 0 in (141) approach the exact attractor solution much slower than the other modes.

Exact solution in the RTA
Another method to study the attractor in kinetic theory is the use of the exact solution of the Boltzmann RTA, following e.g. [134]. Note that an analytic solution in terms of confluent hypergeometric functions was found in the RTA with massless particles with a relaxation time related to the proper time via τ R ∝ τ 1−Δ , where Δ is a constant [132,135]. As in the hydrodynamic case, we will limit our initial consideration to 0 + 1 dimensional Bjorken expansion. The RTA attractor has been further examined in [127].
We begin by taking the flow velocity to be Fig. 7 The RTA attractor. The ratio of the longitudinal pressure to energy density is plotted for various initializations of (149). Figure from [40] which is clearly spacelike z μ η μν z ν = (t 2 − z 2 )/τ 2 = 1. Due to our assumption of boost invariance, the distribution function will depend only on the proper time, τ , the transverse momenta, p T , and the boost-invariant variable [136,137] where E = p 0 . One can then define the quantity v(τ, ω, We will consider massless particles, setting m = 0 in what follows. The above relations lead to a straightforward identification of the energy and longitudinal momentum of a particle as Note that in this case, the integration measure reads With the simplifying assumptions and parameterizations outlined above, the RTA Boltzmann equation takes the simple form where we note that the relaxation time τ R is a function of proper time. The RTA Boltzmann equation has a formal solution [138,139] given by where the damping function is given by and f 0 is an arbitrary initial distribution function. Note that the first term in (160) carries much of the information of the initial conditions at early times, but at times greater than τ R , this first term decays exponentially. Thus the second term describes the late-time behavior of the system. Finally, note that for conformal RTA, we have [140,141] τ eq = 5η sT (τ ) . (162) To introduce anisotropy, one can use the Romatschke-Strickland distribution function [142] with anisotropy parameter, ξ , where without loss of generality, we aligned the anisotropy in the z direction. This was studied in the RTA and with a λφ 4 scalar collision kernel [52], which was extended further by considering dynamical fugacity in the RTA [127]. Interestingly, models with number-conserving kernels were found to approach the attractor more rapidly, especially the lowest lying mode, although that was not the case for the aHydro attractor [143]. One proceeds to solve the equations with the integral iterative method, first outlined in this context in [139,144]. It essentially boils down to integrating both sides of (160) over the momentum squared to arrive at an expression for the energy density. Using the matching conditions, one knows that ε = ε eq . Then one can read off the temperature from ε ∼ T 4 . With the temperature, the equilibrium distribution function is known at all times and so the integration can be done in (160) to find the distribution function, f .

Adiabatic hydrodynamization
Another interesting approach was to study a system with a gapped Hamiltonian, where the RTA Boltzmann equation of the momentum-weighted distribution function (integrating out the modulus p = p 2 T + p 2 z , but not performing the angular integration) is recast as a time dependent Schrödinger equation [145]. The lowest eigenmode is identified as the pre-hydrodynamic mode, which eventually evolves into the hydrodynamic mode, remaining gapped from the faster modes throughout the evolution. In this set-up, the approach to the attractor is mainly characterized by such prehydrodynamic modes, which the authors of [145] refer to as adiabatic hydrodynamization. This is due to adiabaticity, i.e. the excitations of the pre-hydrodynamic modes are suppressed. Future work involves relaxing the assumption of the RTA to a more realistic case, e.g. a Fokker-Planck collision kernel. To make better contact with the QGP, it is useful to go beyond the Bjorken expansion by considering transverse modes. A possible experimental signature could be found in the flow coefficients v n of small systems, where the pre-hydrodynamic modes might not have the chance to fully hydrodynamize.

Hydrodynamic generators
An intriguing notion that has been raised is the existence of hydrodynamic generators [146,147], where the exact solution of the RTA Boltzmann equation at late times generates the Borel-resummed Chapman-Enskog expansion. The Chapman-Enskog expansion is an expansion in the Knudsen number, K N , which is a dimensionless parameter given by the ratio of the mean free path and a characteristic macroscopic length scale, particular to the system at hand. For the RTA Boltzmann equation (159), the Chapman-Enskog expansion is [146] This series was found to be factorially divergent [129,148] even for small Knudsen number. One can try to Borel resum the Chapman-Enskog series as outlined in Sect. 3.1.1 to find Using the exact solution (160) in the case the relaxation time is kept constant, after some algebra one finds where z 0 = (τ − τ 0 )/τ R . Clearly, expanding (166) for large proper time reduces to (165). From here, the authors conjecture that even in the case the relaxation time is proper time dependent, the hydrodynamic generator is given by [146] where D(τ, τ ) is the damping factor (161). The hydrodynamic generator is a nice representation of the Chapman-Enskog expansion as it is finite, is amenable to generalizations to (3 + 1) dimensions and to systems without Bjorken flow, and satisfies the RTA Boltzmann equation Interestingly, the RTA hydrodynamic generator contains the Chapman-Enskog expansion with coupled non-hydrodynamic modes with different decay times. It is found that the higher Fig. 8 Evolution of the background energy density divided by the ideal hydrodynamic energy density (177) as a function of proper time (normalized by the kinetic relaxation time τ R ≡ (η/s)/T id ). Note that attractor curve goes smoothly from the early time free streaming limit to viscous hydrodynamics at late times. Figure from [28] order gradient corrections are suppressed during hydrodynamization, which indicates that these higher corrections do not necessarily play a large role, even with large gradients. This story seems to resemble the modes found via resurgence in [87,149]. Finally, it should be pointed out that the RTA has its limitations. The derivation above rested on the knowledge of an exact solution in the RTA, so incorporating more realistic collision kernels will require the damping function (161) to be replaced more generally with a Green's function, see e.g. [150].

An eye to phenomenology
We have seen that hydrodynamic attractors are found in many systems. The study of these attractors offer plenty of theoretical insights, e.g. to the applicability of hydrodynamics, resurgence and transseries, but what about about real world applications? In this section, we highlight some of the phenomenology in heavy ion physics that attractors played an important part.

KøMPøST
KøMPøST is a framework that matches the early time, nonequilibrium QCD kinetic theory to a later time when a hydrodynamic description is more appropriate [28,29]. The initial energy momentum tensor is described by kinetic theory beginning at the time scale, τ EKT . This is then propagated to later times when hydrodynamics is applicable, τ hydro , via nonequilibrium Green's functions, G μν αβ , which are computed in QCD kinetic theory. The linearized perturbations of the stress tensor evolve according to The barred quantities denote the background energy momentum tensor, which is computed by taking a spatial average over the causal circle. As observed in [28,87,151], the dynamics in these timescales becomes essentially independent of the microscopic kinetic coupling when measured in (macroscopic) dimensionless timē Even far from equilibrium forω < 1, it is observed that the energy momentum tensor is a function ofω only with hydrodynamization occuring afterω > 1. As such, the attractor can be viewed as a bridge between the kinetic and hydrodynamic regimes. The universal attractor behavior can be seen in Fig. 8, where two values of the kinetic theory coupling is plotted for a range of η/s. The nonequilibrium evolution of the system goes from early time free streaming to late time viscous hydrodynamics. The universal scaling function is due to the attractor behavior.
KøMPøST has been incorporated in a number of contexts, e.g. to pre-hydrodynamic evolution and its signatures in final-state heavy-ion observables [152] and more recently, to describe photon observables in small collision systems [153].

Entropy production
The fact that hydrodynamic attractors can serve as a bridge between late-time and earlier, far-from-equilibrium behavior was recently exploited to extract phenomenologically relevant information in the QGP context [24]. The central idea is that the multiplicity of final state particles measured in detectors is intimately linked to the entropy production during the pre-equilibrium phase. The total entropy per unit spacetime rapidity, η s , during the near-equilibrium transverse expansion of the QGP is which is related to the charged particle multiplicity at late times via where J denotes the Jacobian factor, y is the rapidity, dy/dη ≈ 1.1, A ⊥ = π R 2 and 2R is the transverse extent of the system. Then the initial energy density can be linked to the experimentally observed particle multiplicities via the neat pocket formula [24] (sτ ) hydro where C ∞ is a constant of order unity and ν eff is the number of effective degrees of freedom. For an ideal gas of gluons, ν eff = 16, while for a mixture of quarks and gluons, ν eff = 47.5.
To derive this result, the knowledge that there is an attractor can be put to work to find where the quantity on the right hand side is the attractor curve.
To see this, we follow the derivation of the link between the initial state energy density, ε 0 , and the near thermal system energy density, ε hydro , beginning from the boost-invariant conservation law The energy density is related to the effective temperature via It is then convenient to introduce the dimensionless scaling variableω given by (171). Finally, calling the anisotropy f (ω) ≡ P L /ε, we see that we can rewrite (176) as Using the definition ofω and the relation between the energy density and the effective temperature (177), we see that Thus, After integration, we arrive at From here, we arrive straightforwardly at (175). Note that the energy attractor (175) interpolates smoothly between free streaming and late time viscous hydrodynamics respectively. The pocket formula (174) and the charged particle multiplicity (173) can be used to estimate the initial energy density for central lead collisions at √ s N N = 2.76 TeV to be ε ≈ 270 GeV/fm 3 at the time τ 0 = 0.1 fm/c [24]. Furthermore, the study of such energy attractors has been used to characterize potential probes of the QGP at early times, such as intermediate mass dileptons [154]. Interestingly, the evolution of a system following an attractor solution has been shown to result in higher yields (in fact, maximum yields) in thermal dileption and photon particle production [155], compared to solutions away from the attractor.
The previous argument was turned on its head in [25], where the assumption of free streaming at early times was explored by positing the energy density scales with some to-be-determined parameter, β, namely ε ∼ τ −β (the free streaming case has β = 1). Using different initial state models and comparing the final state to experimental data, they find that the preferred expansion parameter β is slightly, but distinctly different from the free streaming case.

Non-zero chemical potential
Real world QGP has non-zero baryon chemical potential. As such, it is natural to include the effect of non-zero chemical potential on the attractor story described so far. This has been explored both in hydrodynamic [26,27] and kinetic setups [156,157]. In the hydrodynamic case, the 0 + 1d DNMR Bjorken case described in Sect. 2.2 is supplemented with an evolution of the baryon number density and relaxation of the bulk viscosity. The late-time attractor regime was not reached due to the short lifetime of the hydrodynamic runs, so it is more apt to talk of potential attractors, which the system looked like it was approaching. Although attractors were not systematically studied, there was the indication that near the critical point, the bulk attractor should have a pronounced effect. Also, the potential attractor line has a different behavior when comparing MIS to DNMR near the critical point, indicating care should be taken as there are phenomenological implications and differing results depending on which set of equations of motion is used [27].
In the kinetic approach of [156,157], quark/anti-quark and gluonic distribution functions are taken into account, which respectively take the form where ξ 0 is the initial anisotropy and Q 0 is the initial energy scale. The evolution of the anisotropic distribution functions with form similar to (163) is given by the QCD Boltzmann equation [158] with a collision kernel with leading order elastic (2 ↔ 2) and inelastic (1 ↔ 2) interactions taken into account, i.e. for a distribution function of species a, f a = f a (τ, p T , p ) The evolution of the system is given by the conservation of the energy momentum tensor and the conserved current, ΔJ μ f = (Δn f , 0, 0, 0), via The nonequilibrium system has a build up of longitudinal pressure due to interactions while the longitudinal expansion slows down, before finally at late times approaching local equilibrium with p L = p T = ε/3. Like in the case of zero chemical potential discussed previously, the validity of viscous hydrodynamics as a description of the non-equilibrium evolution in the case of non-zero chemical potential also begins at timescales roughly given by 4πη/(sT ) [156,157]. However, at early times non-zero chemical potential can have a large impact on the pressure for times earlier than 4πη/(sT ), although the late time approach is found to be affected by around 10%. This is to be contrasted with results from kinetic theory computations involving a single component distribution function, as discussed in Sect. 4, where the rapid memory loss of the initial conditions characteristic to attractor behavior allows for a universal approach to viscous hydrodynamics. As such, it is clear that the chemical composition of the QGP has an effect on the approach to the universal approach to viscous hydrodynamics at late time.

Improved freeze-out procedures
It is worthwhile to keep in mind that for all the discussion of hydrodynamic variables as a means to describe heavy ion collisions, we have left out a crucial ingredient. Experiments themselves measure the distribution of particles and not the fluid-dynamic variables directly as discussed previously. Following the collision and initial evolution of the quark gluon plasma, the distribution of particles "freeze-out" as color neutral hadrons and free stream towards the detectors at time scales around τ ∼ 10 fm/c. The translation of distribution of particles to hydrodynamic variables is known as the freeze-out procedure.
Recently, it was pointed out in [128] that possible improvements to freeze-out procedures can be made in light of the developments surrounding the attractor. Using Arnold-Moore-Yaffe (AMY) effective kinetic theory in one-dimensional Bjorken flow, where the collision kernel includes gluon scattering and inelastic gluon splitting using the pure glue code from [159], a comparison was made between aHydro and linearized freeze-out prescriptions. It was found that the freeze-out prescriptions studied tend to reproduce the low-order moments of the distribution at late times τ 5 fm/c. However, for higher moments, it was found that the linearized prescriptions do poorly compared to aHydro. The last fact indicates that the aHydro freeze-out prescription can be of use for the study of smaller systems, like p A or peripheral A A collisions.

Beyond conformal symmetry
The majority of the work discussed so far centered on conformal systems, which provides the advantage of theoretical clarity and sometimes the feasibility of analytical computations. However, in more realistic systems, one should not expect full conformal symmetry to hold, especially with a full three dimensional expansion and a realistic equation of state. There are a number of works that go beyond the status of the field, including the exploratory work of attractors with a non-conformal equation of state relevant for the QCD critical point considered in [26,27], see Sect. 5.3 for more details.
Here we discuss an initial survey that was made in [20] for a 0 + 1-dimensional Bjorken kinetic theory of a gas of particles of mass m and for inhomogeneous 2 + 1 dimensional BRSSS hydrodynamics. Two quantities were proposed to investigate whether attractor solutions are present, namely The first quantity represents a simple rewriting of the evolution of the energy density from the Navier-Stokes equations, while the second denotes the deviation of the equation of state from the conformal limit. When expressed as a function of the inverse gradient strength Γ , both A 1 and A 2 exhibit attractor behavior. However, the approach of different solutions is slower for A 2 than A 1 , which may indicate that at early times, the attractor is harder to discern in non-conformal systems. In the case of the 2 + 1 dimensional BRSSS hydrodynamics, the apparent attractor was seen by studying A 1 (Γ ), which seemed to greatly depart from the expected Navier-Stokes limit, although such a simulation proved too computationally costly to extend for much later times. In many respects, although the exploration of the nonconformal limit is just at beginning, there are a number of lessons that have been learned already. As mentioned above, moving away from the conformal limit also can provide insight into the nature of the early time attractor. Recent work involved studying a 0+1 dimensional gas of massive particles evolving under the Boltzmann equation in the RTA, which was compared to non-conformal hydrodynamics [160,161] . It is found that only the longitudinal pressure has an early time attractor, while that there is no early time hydrodynamic attractor for the bulk and shear viscous stresses, adding evidence to the related indication in [40]. Moreover, the comparison of the quantity (191) in [161] indicates that while for late times A 1 does exhibit universal attractor behavior, at early times there is a noticeable difference between the solutions depending on the choice of different masses. Another interesting observation in non-conformal Bjorken flow is that it is not a guarentee one even reaches local thermal equilibrium at late times. For instance, in the case of variable relaxation time scaling, namely τ R ∝ T −Δ (Δ = 1 corresponds to the conformal limit), the system can evolve to a late time regime out of local thermal equilibrium for Δ = 1 [149], in contrast to conformal Bjorken flow.

Conclusion and outlook
The story of hydrodynamic attractors in heavy ion collisions is still being written. We can take this point to reflect on where the field stands and where it might go. At this stage, it seems like attractors for Bjorken flow are well-established for a variety of effective theories relevant for heavy ion collisions, such as relativistic hydrodynamics, holography and kinetic theory. Attractors provide a conceptual mechanism to understand the applicability of hydrodynamics even in far-from-equilibrium systems, as well as providing some phenomenological insight. Notwithstanding the current successes, there is still much work to be done. Some open questions remain: -What are the theoretical properties of the QCD attractor and under which conditions does QCD exhibit attractor behaviour? -How do attractors behave in systems with relaxed symmetries? Studies of attractors have begun to go beyond the Bjorken and Gubser cases described here. Ultimately, the goal should be to consider the full symmetries in heavy ion collisions, moving away from the conformal limit as in [20,[160][161][162] and including transverse expansion in 3 + 1 dimensions, as well as the quark degrees of freedom as in [27,156,157,163,164]. There is some promising work to take transverse dynamics into account [165]. Recently, the attractor was studied in a Hubble expanding background [166], where the scale factor was assumed to follow a power law form a(t) ∼ t α instead of solving the Einstein equations to find a(t), an assumption informed by cosmology. -Beyond numerical techniques, how can we extend the determination of attractors to earlier times? Especially in the more complicated setting of relaxed symmetries, it is not clear if there is an unambiguous early time attractor. Also, in the holographic context, it seems like there is only the late-time attractor [40] (although the next question also can be relevant for the early time attractor). -How would finite coupling corrections affect the story of holographic attractors? The Gauss-Bonnet setup described in Sect. 3.3 can act as a laboratory to study higher derivative corrections. The observation that varying a single parameter allows the interpolation from the strongly coupled to the expectation from weakly coupled kinetic theory should be viewed as an encouraging step in this direction. It can be speculated that finite coupling corrections thus might provide an unambiguous method of determining the possible early time behavior of the attractor. -How do hydro attractors fit in with non-thermal fixed points? It would be interesting to see a study of the evolution of a kinetic model from a non-thermal fixed point to a late time hydrodynamic attractor. -Are there any signatures of attractors in smaller systems, e.g. pp collisions? In pp collisions, the number of particles produced is regularly quite low compared to typical heavy ion collisions. It is not clear if hydrodynamization would occur in such short-lived settings [167]. If so, there could still be the approach to the hydrodynamic attractor for generic initial conditions without reaching the attractor at all, potentially indicating some level of pre-hydrodynamization in the smaller systems.
-Which observables are sensitive to hydrodynamic attractors? Attractors seem to be a useful tool in using final state information to probe earlier times in heavy ion collisions, see Sect. 5. Furthermore, one can envisage using them to determine which final state observables depend directly on the initial state. -How universal are hydrodynamic attractors? The hydrodynamic limit exists in a number of neighboring fields, including cosmology and condensed matter systems. The lessons learned about attractors in the high energy community can surely be applied to expanding systems in a variety of fields.