Far-from-equilibrium attractors and nonlinear dynamical systems approach to the Gubser flow

The non-equilibrium attractors of systems undergoing Gubser flow within relativistic kinetic theory are studied. In doing so we employ well-established methods of nonlinear dynamical systems which rely on finding the fixed points, investigating the structure of the flow diagrams of the evolution equations, and characterizing the basin of attraction using a Lyapunov function near the stable fixed points. We obtain the attractors of anisotropic hydrodynamics, Israel-Stewart (IS) and transient fluid (DNMR) theories and show that they are indeed non-planar and the basin of attraction is essentially three dimensional. The attractors of each hydrodynamical model are compared with the one obtained from the exact Gubser solution of the Boltzmann equation within the relaxation time approximation. We observe that the anisotropic hydrodynamics is able to match up to high numerical accuracy the attractor of the exact solution while the second order hydrodynamical theories fail to describe it. We show that the IS and DNMR asymptotic series expansion diverge and use resurgence techniques to perform the resummation of these divergences. We also comment on a possible link between the manifold of steepest descent paths in path integrals and basin of attraction for the attractors via Lyapunov functions that opens a new horizon toward effective field theory description of hydrodynamics. Our findings indicate that anisotropic hydrodynamics is an effective theory for far-from-equilibrium fluid dynamics which resums the Knudsen and inverse Reynolds numbers to all orders.


Introduction and summary
Hydrodynamics is an effective theory which describes the long wavelength phenomena and/or small frequency of physical systems. The fluid dynamical equations of motion are derived by assuming that the mean free path is smaller than the typical size of the system [1]. The existence of a large separation between the microscopic to macroscopic scales can be reformulated as a small gradient expansion around a background (usually the thermal equilibrium state) which varies slowly. Thus, hydrodynamics would be invalid in far-from-equilibrium situations where the gradients of the hydrodynamical fields are large. This might be the case in systems of small size. However, recent experimental results of p−p collisions at high energies [2,3] have shown evidence of collective flow behavior similar to the one observed in ultra-relativistic heavy ion collisions. The experimental data measured in p − p collisions can be described quantitatively by using hydrodynamical models [4][5][6]. On the other hand, different theoretical toy models for both weak and strong coupling [7][8][9][10][11][12][13][14][15][16][17][18] have presented an overwhelming evidence that hydrodynamics becomes valid even in non-equilibrium situations where large gradients are present during the space-time evolution of the fluid. The success of hydrodynamical models to describe small systems as well as applying it to far-from-equilibrium situations calls for a better theoretical understanding of the foundations of hydrodynamics.
For different strongly coupled theories it has been shown that the hydrodynamical gradient series expansion has zero radius of convergence and thus, it diverges [19]. On the other hand, the divergent behavior of the hydrodynamical series expansion is also well known in weakly coupled systems based on the Boltzmann equation for relativistic and non-relativistic systems [20][21][22][23][24]. A more detailed mathematical analysis of the origin of this divergence has unveiled the existence of a unique universal solution, the so-called attractor [25]. The novel attractor solution is intrinsically related with the mathematical theory of resurgence [20,26,27] and the details of this solution depends on the particular theory under consideration [28][29][30][31][32][33][34][35][36][37]. In simple terms, the attractor is a set of points in the phase space of the dynamical variables to which a family of solutions of an evolution equation merge after transients have died out. In relativistic hydrodynamics it has been found in recent years that for far-from-equilibrium initial conditions the trajectories in the phase space merge quickly towards a non-thermal attractor before the system reaches the full thermal equilibrium. This type of non-equilibrium attractor can be fully determined by very few terms of the gradient series of relatively large size which involves transient non-hydrodynamical degrees of freedom [19,20]. This property of the attractor solution indicates that the system reaches its hydrodynamical behavior and thus, hydrodynamizes at scales of time shorter than the typical thermalization and isotropization scales. The fact that hydrodynamization happens in different size systems at short scales of time while exhibiting this degree of universality in far-from-equilibrium initial conditions might explain the unreasonable success of hydrodynamics in small systems such as p − p and heavy-light ion collisions [38].
In this work we continue exploring the properties of attractors for rapidly expanding systems within relativistic kinetic theory. Previous works have focused on fluids undergoing Bjorken flow [20,26,27,36,37] and non-homogeneous expanding plas-mas [39]. We expand these studies by investigating the properties of the attractors in the plasmas undergoing Gubser flow [40,41].
The Gubser flow describes a conformal system which expands azimuthally symmetrically in the transverse plane together with boost-invariant longitudinal expansion. The symmetry of the Gubser flow becomes manifest in the de Sitter space times a line dS 3 ⊗ R [40,41] and thus, the dynamics of this system is studied in this curved spacetime. The search for attractors in the Gubser flow poses new challenges due to the geometry and the symmetries associated to this velocity profile. We determine the location of the attractors with the help of well-known methods of nonlinear dynamical systems [42,43]. Our results bring new features and tools to the study of attractors which were not addressed previously in the context of relativistic hydrodynamics. For instance, in the 2d system of ordinary differential equations (ODE) derived from different hydrodynamical truncation schemes, we study the flow diagrams -streamline plots of the velocity vector fields in the space of state variables-and carefully examine the early-and late-time behavior of the flow lines near each fixed point and show that the system is exponentially asymptotically stable. We observe that the attractor is a 1d non-planar manifold only asymptotically because the 2d system cannot be dimensionally reduced to a 1d non-autonomous one by any reparametrization of the variables, which is a signature of Gubser flow geometry governed by a true 3d autonomous dynamical system as opposed to the Bjorken model. This is made more precise in the context of this 3d dynamical system, where the linearization problem is reviewed and a mathematically rigorous definition of attractor is given. We estimate the shape of the basin of attraction by giving an approximate Lyapunov function near the attracting fixed point. There we touch upon an important link between the basin of attraction and attractors with manifold of steepest descent paths in the path integral formalism of quantum field theories and briefly discuss the role of a Lyapunov function as an analogue of some effective action for the stable hydrodynamical and non-hydrodynamical modes.
The paper is structured as follows: In Sec. 2 we review the Gubser flow, its exact solution for the Boltzmann equation within the relaxation time approximation (RTA) as well as different hydrodynamical truncation schemes. In Sec. 3 we study extensively the flow lines of the IS theory in 2 and 3 dimensions. Our numerical studies of the attractors for different theories and their comparisons are discussed in Sec. 4. In this section we also address the issue of the divergence of the IS and DNMR hydrodynamical theories and how to fix it using resurgent asymptotics. Our findings are summarized in Sec. 5. Some technical details of our calculations are presented in the Appendices.

Setup
Before starting our discussion we first introduce the notation used in this paper. We work in natural units where = c = k B = 1. The metric signature is taken to be "mostly plus" (−, +, +, +). In Minkowski space with Milne coordinates x µ = (τ, r, φ, ς) the line element is given by where the longitudinal proper time τ , spacetime rapidity ς and polar coordinates r and φ are related with the usual Cartesian coordinates (t, x, y, z) through the following expressions It is better to study the Gubser flow in de Sitter space times a line (dS 3 ⊗ R) that is a curved spacetime in which the flow is static and the symmetries are manifest [40,41]. This is obtained by applying a conformal map between dS 3 ⊗ R and Minkowski space, which consists of rescaling the metric ds 2 → dŝ 2 = e Ω ds 2 with Ω = log τ −2 . Afterwards one performs the following coordinate transformation Here,τ = qτ andr = qr with q being an arbitrary energy scale that sets the transverse size of the system [40,41]. The time-like coordinate is the variable ρ ∈ (−∞, ∞) and the polar coordinate is θ ∈ [0, 2π). Therefore, the line element in dS 3 ⊗ R reads as so the metric isĝ µν = diag −1, cosh 2 ρ, cosh 2 ρ sin 2 θ, 1 . In dS 3 ⊗ R the flow velocity is the normalized time-like vectorû µ = (1, 0, 0, 0), which is invariant under the q-deformed symmetry group SO(3) q ⊗ SO(1, 1) ⊗ Z 2 [40,41]. We choose the fluid velocity to be defined in the Landau frame, i.e.,T µνû ν ≡ˆ û µ .

Relativistic kinetic theory for the Gubser flow
In kinetic theory the macroscopic hydrodynamic variables are calculated as momentum moments of the distribution function. The symmetry of the Gubser flow restricts the phase-space distribution f (x,p) = f (ρ,p 2 Ω ,p ς ), i.e., the distribution function depends only on invariants of the SO(3) q ⊗ SO(1, 1) ⊗ Z 2 symmetry group: the de Sitter time ρ, the combination of momentum componentsp 2 Ω =p 2 θ +p 2 φ / sin 2 θ wherep θ and the longitudinal momentum componentp φ are conjugate to the coordinates θ, φ, respectively, as well asp ς , that is conjugate to the coordinate ς [9, 13] 2 . Thus, in dS 3 ⊗ R the RTA Boltzmann equation reduces to a one-dimensional relaxation type equation [9,13] whereT is the temperature. In this work we take f eq (z) = e −z as the local thermal equilibrium distribution. The conformal symmetry demandsτ r (ρ) = c/T (ρ) with c = 5 η/s where η and s are the shear viscosity and entropy density, respectively. Given a solution for the Boltzmann equation (2.5) one calculates the energy-momentum tensor as a second-rank tensor moment defined as 3 The relaxation equations ofT µν are obtained by either solving Eq. (2.5) exactly or by finding an approximate perturbative solution. In the next section we review different approximate methods to obtain the fluid dynamical equations within kinetic theory as well as the exact solution of the kinetic Eq. (2.5).

Fluid dynamical theories
For weakly coupled systems the equations of the macroscopic fluid dynamical variables are derived from a microscopic underlying kinetic theory based on the Boltzmann equation. The derivation of these equations assumes that the Boltzmann equation admits a generic solution of the form where f b is a distribution function that describes the evolution of an existing background, M α are orthogonal polynomials of degree l, whose orthogonality properties depend on f 0 , and a α are moments of the full distribution function f . The leading order background distribution function is chosen based on the problem at hand [62]. In the rest of the section we shall briefly describe two choices of the leading order background distribution function which leads to different relaxation equations for the components of the energy-momentum tensor.

Expansion around an equilibrium background
The canonical derivation of fluid dynamics assumes an expansion around a local equilibrium background For relativistic systems f eq is taken to be the equilibrium Jüttner distribution function f (x µ , p i ) = 1/ e β[(u·p)−µ] + a where β = 1/T with T being the temperature, µ the chemical potential and a = +1, −1, 0 for particles following Fermi-Dirac, Bose-Einstein or Maxwell-Boltzmann statistics, respectively. The energy of the particle −û ·p is taken to be isotropic in the local rest frame (LRF) while δf encodes the deviations from the thermal equilibrium state. It is implicitly assumed that δf f eq . For systems close to equilibrium the four-momentum is decomposed asp µ =Êpûû µ +p µ whereÊpû = −(û ·p) (in LRFÊpû ≡p ρ ) whilep µ =∆ µνp ν projects over the spatial momentum component orthogonal to the flow velocity. For the Gubser flow this vectorial decomposition of the four-momentum allows to write the most general conformal isotropic energy-momentum tensor [63] T µν =ˆ û µûµ +P∆ µν +π µν , (2.9) whereˆ is the energy density,P is the isotropic pressure andπ µν is the shear stress tensor, which are the momentum moments of the distribution function given bŷ  [40,41]. For conformal systems the equation of state reads asP =ˆ /3.
The local temperature of the system introduced in Eq. (2.8) is found from the Landau matching conditionˆ =ˆ eq. (T ) [1]. This phenomenological constraint ensures that the parameterT in f eq is adjusted such that the dissipative corrections encoded in δf do not shift the value of the energy density. As a consequence, deviations from the local thermal equilibrium are captured entirely by the shear stress tensorπ µν ≡ p µpν δ where · · · δ indicates a momentum moment weighted by δf .
The energy-momentum conservation lawD µT µν = 0 for the energy-momentum tensor (2.10) gives us the following evolution equation for the energy densityˆ [40,41] ∂ ρˆ ˆ + 8 3 tanh ρ =π tanh ρ , (2.11) which can be rewritten in terms of the temperature since for conformal systemsˆ ∼T 4 and thus, ∂ ρT /T = 4∂ ρˆ /ˆ . As a result one finds [40,41] ∂ ρT whereπ :=π ςς /(ˆ +P). We point out that for the Gubser flow there is only one independent component of the shear viscous tensor [40,41]. It is needed an additional evolution equation forπ obtained by expanding δf within some approximation. The form of δf is found by expanding it in terms of a set of orthogonal polynomials in energyÊpû, and a set of irreducible tensors in momentum invariant under the little group SO(3) of the Lorentz transformations, i.e., 1, p µ , p µ p ν , etc. Afterwards, the relaxation equations of the dissipative macroscopic quantities are obtained by applying a systematic truncation method based on a power counting in the Knudsen Kn and inverse Reynolds Re −1 numbers [45]. The lowest truncation order provides the IS equations while the inclusion of all other second order terms, i.e. in Kn 2 , Re −2 and Kn · Re −1 , gives the DNMR equations. A detailed discussion of this power counting schemes can be found in Ref. [45]. For the normalized shear stressπ =π/(ˆ +P) one finds the following evolution equation for the IS and DNMR theory respectively [13,64] For IS theory:τπ dπ dρ (2.13b) The conformal symmetry impliesτπ = c/T with c = 5η/s.

P L matching
When introducing the leading order anisotropic distribution function defined in Eq. (2.14) we did not comment on how to determine the parametersΛ and ξ. The different variants of aHydro have been proposed in recent years in order to perform this matching procedure [46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61]. In a previous study of the Gubser flow, it was shown that thê P L matching scheme [48] provides the most accurate macroscopic description when comparing its predictions with the ones obtained from an exact solution of the RTA Boltzmann equation. Thus, we shall consider in this work this matching procedure. We start by introducing first the anisotropic integrals The first term on the RHS of the previous equation comes from the leading order contribution associated to the anisotropic distribution function f a (2.14) while the second term corresponds to the subleading contribution from δf in Eq. (2.14). In Appendix A we show how to perform the integrals of the leading order anisotropic distribution function for the massless case. The functional form of f a is taken to be the RS ansatz which in the LRF looks like [69] f a (x,p;Λ, ξ) = f eq E RS (ξ)/Λ (2.19) where f eq (z) = e −z is a Maxwellian distribution function evaluated for the momentumanisotropic argument The leading order anisotropic variables contributing to the anisotropic energy-momentum tensor (2.15) areT The traceless condition of the conformal anisotropic energy-momentum tensor implies thatˆ RS = 2P RS ⊥ +P RS L . The conservation law gives the evolution equation for the energy density The matching condition for the energy densityˆ RS (Λ, ξ) =ˆ eq (T ) leads to the following relation between the temperatureT and the momentum scaleΛ The functionR 200 (ξ) is given in Eq. (A.3). This relation ensures that the energy density does not receive any contribution from the residual deviation δf in Eq. (2.14). Now, for the Gubser flow the conservation law (2.23) indicates that the pressure anisotropy is the force that drives the system far from equilibrium. The microscopic origin of this force is the momentum anisotropy created by the expansion of the system which in the LRF is measured by the anisotropy parameter ξ. Hence, one simply adjusts the value of ξ such that the leading order anisotropic distribution function f a fully captures the information of the full pressure anisotropy. For the Gubser flow this means that the longitudinal pressureP L (2.16c) does not receive contributions from the residual deviation δf of the distribution function (2.14), i.e.P L =P RS L . Now the effective shear viscous tensor is related with the total pressure anisotropy via the identity (2.17) which in this case gives [48] The functionsR 240 andR 200 are listed in App. A. Furthermore, ξ(π) in Eq. (2.29) is the inverse of the function An advantage of using aHydro is that by construction the transverse and longitudinal pressures remain positive during the entire evolution of the system as it is expected in kinetic theory 5 . These constraints imply −1/4 <π < 1/2 which is in agreement with the anisotropy parameter ξ ∈ (−1, ∞). On the other hand, the LHS of Eq. (2.30) is identified as a term proportional to the inverse Reynolds number Re −1 = √ π µν π µν /P 0 (with P 0 =ˆ /3 for the conformal case) [45]. Eq. (2.30) therefore indicates that the anisotropy parameter resums non-perturbatively not only large gradients (i.e., Knudsen number) but also large Re −1 numbers. An analogous relation between the effective shear and the anisotropy parameter was found in the Bjorken case [46]. The dissipative corrections O(Re −2 ) arise in general from the most nonlinear sector of the collisional kernel [45] and thus, the calculation of these terms is cumbersome even for the simplest kernels [70]. Moreover, some of the nonlinear terms O(Re −2 ) calculated in the DNMR theory lead to violations of causality [71].

Exact solution to the Boltzmann equation in the RTA approximation
The RTA Boltzmann equation (2.5) admits the following exact solution [9,13] is the damping function. For the initial condition of the distribution function f 0 at ρ 0 we shall consider the RS ansatz [69] whereΛ 0 is the initial temperature and ξ 0 is the initial momentum anisotropy along the ς direction. From the exact solution for f one gets the energy density and the only independent component of the shear stress [9,13]: These integral equations are solved numerically by means of the method described in Refs. [11][12][13]72]. The temperature of the system is obtained from the energy density through the Landau matching condition, i.e.,ˆ RS (Λ,

Flow diagrams of the IS theory
In dS 3 ⊗ R the expansion rate of the Gubser flowθ = 2 tanh ρ becomes a constant when ρ → ±∞. As we shall see below this geometrical property of the velocity profile poses a challenge to find the attractors of different hydrodynamical theories. In this section we explain a method for finding the attractors based on the mathematical theory of non-autonomous systems [42,43]. We discuss extensively the case of IS theory. The same method can be used in different models.
The gist of what we are about to do in this section is to study the flow diagrams of the Gubser flow from the perspective of nonlinear dynamical systems. The flow lines at early times do dictate the far-from-equilibrium behavior of any system governed by a set of differential equations and at late times they show what happens to the matter distribution until it evolves to a steady state (thermally non-equilibrium) at ρ → ∞ 6 .
As a minor digression, we recall that for the Bjorken flow, the second order hydrodynamical equations can be entirely reduced to one single explicitly time-dependent ODE (i.e. 1d non-autonomous parametrized by a new time w) due to scaling symmetry [19,25]. This suggests that its attractor is a planar 1-manifold characterized by w and the basin of attraction is indeed 2d. But what is the proper definition of attractor? It is an invariant set of points and every flow line starting from a point within the field of attraction of the attracting fixed point, will always limit to this set of points at late times. Given a set of initial values for Bjorken time, the flow lines always lie on a plane and as long as they belong to a special set of numbers, they tend to the attractor at late times. This set defines the basin of attraction: Given an initial time τ 0 , there is an associatedT (τ 0 ) andπ(τ 0 ), and the basin of attraction is elaborated as the set of all pairs of (T (τ 0 ),π(τ 0 )) such that the Bjorken flow lines are doomed to limit to the equilibrium. We remind that (T ,π) is the actual phase space 7 of the Bjorken flow which fixes the dimensionality of the basin of attraction as well.
In the case of Gubser flow, we simply identify the fixed points and discuss whether or not they are stable by analyzing the eigenvalues of the Jacobian matrix of the linearized flow equations near each fixed point. We find that there are two unstable saddle points and one stable fixed point which determines the steady state of the system in any hydro theory. This task in carried out by analyzing the nulllines of the 2d nonautonomous system (3.1a)-(3.1b). The intersection of these lines gives the fixed points. By analyzing the asymptotics of flow lines close to the steady state, it is determined the fixed point is indeed exponentially asymptotically stable. One observation that we make is that the flow equations cannot be reduced to one single non-autonomous ODE by any kind of reparametrization unless at very late times such that ρ is large enough that we can take tanh 2 ρ ∼ 1. But this is fine from the viewpoint of attractors since they are asymptotic or in other words limit sets. This does not affect the dimensionality of the attractor and hence Gubser flow lines still evolve into a 1-manifold but in reality it is no longer a planar curve since the time direction is not fully decoupled from the ODE in terms of the new time w, which will be clarified in subsection 3.2.1. What about the basin of attraction? It is indeed what we go at great length into its details in subsection 3.2.2 where the impossibility of this dimensional reduction to a 1d system is a symptom of the peculiarity of the q-deformed conformal symmetry SO(3) q which is not a simple scaling of variables and time as in the Bjorken case. The variable τ = tanh ρ is promoted to be independent which lifts the phase space dimensionality by 1, now being labeled by (T ,π, τ ). Therefore, we lift the dynamical system describing the flow equations to a 3d autonomous system where ρ-dependency is now implicit. So the basin of attraction is three dimensional, being identical to the dimensions of the autonomous system (or phase space).
We finally emphasize on the importance of the basin of attraction and the way to define it locally and globally. As was said before, the basin of attraction defines an effective attractive potential field for the flow lines, thus forcing them to evolve into the stable steady state at late times. This is an important problem in the farfrom-equilibrium hydrodynamics to map out all the divergent flows and only probe the relevant stable ones. Similar to the path integral where a steepest descent path is one that counts as contributing to the path integral manifold, in hydrodynamics the space of all flow lines taking initial values in the basin of attraction would determine this manifold. Effectively, it is just enough to consider all the paths starting at ρ → −∞ somewhere on the boundary of basin where the attractive force field of the fixed point is the weakest. This field is determined by an effective potential function(al) known as Lyapunov function V that satisfies two key properties: V > 0, and dV/dρ ≤ 0. The existence of V hints at the stability of steady state and that the linearized flow equations are asymptotically stable. The local function V loc can always be computed in the nearequilibrium region by solving Lyapunov equation given in 3.2.2 and we hence solve it to estimate the shape of the local basin of attraction. The global function is hard to be precisely built but there are numerical and analytical optimization techniques that are mentioned in the same subsection. The global Lyapunov function mimics the behavior of all stable fluctuations both at early and late times which in this regard might be a very useful tool to construct an effective partition function for hydrodynamics. We will close this section by commenting on this issue.

From the perspective of 2d non-autonomous dynamical system
Using a secondary time parameter τ = tanh ρ ∈ [−1, 1], Eqs. (2.12) and (2.13a) are put into the form These two equations can be combined in such a way that we are left with the following second-order nonlinear non-autonomous differential equation for the temperatureT (τ ) in the IS theory To find the fixed points of the system, we solve (3.1) for theπ N andT N along the respective nulllines (An A-nullline is a trajectory in the phase space along which dA/ds = 0 where s is the flow time variable.). We first consider the case where theT -nullline is given by the solutionπ N (τ ) = 2 to Eq. (3.1a) with vanishing temperature derivative, andT (τ ) is a nonzero constant. The Gubserπ-nullline equation, then fixes the value of temperature atT N (τ ) = −38cτ /15, which at large time τ ∼ 1 is simplyT c = −38c/15. The intersection of these nulllines are not generally reached by all flow lines since near τ ∼ 1 there are flow lines, for instance, that either diverge from or converge to this point. Furthermore, knowing that the system converges to a fixed point sometime at late times τ 0, it is certainly not physical to consider the point to be reached has a negative temperature (−38c/15) and thus we are led to a situation where we look for the stable steady state in the rangeπ =π c < 2. Fig. 1 indicates explicitly why this point cannot define a stable steady state for the system. Taking this to be the case, one may solve for temperature from (3.1a), to only obtain for a constantT 0 > 0. Note that this represents the evolution of temperature near τ = 1 and it is always positive for all times. We can check that Inserting this temperature into (3.1b) gives the shear along theπ-nullline toward the stable steady state, We also notice that lim τ →1 − (1 − τ 2 ) dπ ± N dτ vanishes as expected. At the limit τ → 1 − , this yields the value ofπ ± c = ±1/ √ 5. Similarly the rayT = 0 is aT -nullline, along which Eq. (3.1a) can be solved subject to conditionsπ ± (∞) = ±1/ √ 5 to givē which demonstrates the evolution of the shear component of the trajectories near the fixed points (0, π ± c ). The fixed points in the IS theory of the Gubser flow are therefore determined to be As defines the initial condition in Eq. (3.11). Physically, this is the situation where the distribution of matter is most anisotropic, thus farthest from the stable steady state. In the panel (b) the flow diagram shows that the stable fixed point (3.8) is approached by the flow lines at late times (here τ = 0.9) Note that, in terms of the new variable τ the time flows from τ = −1 to τ = 1 and therefore an observer having a light of sight along the time direction would see the flow lines coming out of the point (0, 1/ √ 5) at early times and going in at late times, which is suggestive of the basin of attraction being three dimensional. At late-times the flow diagram for Gubser flow in the IS theory at ∼ τ = 0.9 shows that there are one attracting fixed point e.g. a sink, and two repelling fixed points e.g. saddle points. The area trapped between the green lines on the top right represents the physical portion of the basin of attraction at a given time slice.
which is consistent with the flow diagrams plotted in Fig. 2. We can also read the Lyapunov exponent 8 ofT from the asymptotic solution (3.4), for both fixed points , thus hinting at the exponentially fast running up of the flow lines at ρ → ∞ to the stable steady state and even faster convergence to the repelling point only along the attracting direction. Along the shear component of the flow trajectories close to (0, ±1/ √ 5) the Lyapunov exponent is given by λπ = −8/(3 √ 5), being a sign of seemingly drastic convergence (∼ e λ + π ρ ) to the steady state or divergence from the repeller. This provides a de facto evidence of the asymptotically exponential stability of the steady state in Gubser model which can be rigorously derived by building Lyapunov functions of the linearized system.
We will be mainly interested in the late-time behavior of the Gubser flow near the attracting fixed point. We refer to this as "attractor" that points to the existence of some bounded set, B A e , namely an absorbing set, in the phase space X of all independent state variables (T ,π) and possibly time τ . B A e is indeed invariant under the forward evolution of state equations in time such that every solution of the system of ordinary differential equations (ODEs) has to ultimately enter B A e provided that the chosen initial conditions for the flow trajectories put them in a bounded set B e known as "basin of attraction" for the attracting fixed point. Notice that by way of definition, B A e ⊂ B e . Let φ τ be the flow defined by the Gubser dynamical system, Eqs. (3.1a)-(3.1b), then we can define an attractor for this system by which will be compact and invariant, and subject to the condition every flow trajectory will approach this set as τ → 1 − . Since there are two other fixed points in the Gubser flow geometry, this condition guarantees that the flow trajectories will not start in the other basins ⊂ X \ B e otherwise φ τ >−1 would always lie in X \ B e . Hence any flow line starting in the basin of attraction B e independently of their whereabouts will always converge asymptotically to the attractor A at late times. 8 The Gubser flow as a dynamical system is deterministic in the sense that the stability of its fixed points is quantified by the eigenvalues of the linearization matrix (so-called Jacobian matrix, cf. Eq. (3.17)) at these points. Therefore, the Lyapunov exponents are the real parts of the eigenvalues of this matrix computed along a flow trajectory that satisfies the ordinary differential equations (3.15). The more interesting case of deterministic chaos occurs when the maximal Lyapunov exponent is positive. In the current system under study, or in any known hydrodynamical system, the maximal Lyapunov exponent is always negative near a stable fixed point which does not give rise to chaos.
In general an attracting fixed point defining a basin of attraction for any (nonchaotic) model is referred to as the dynamical equilibrium state of the system and for the Gubser system this is denoted by (T e ,π e ). From a differentiable geometry point of view, in the phase space of a 2d autonomous dynamical system, one would expect the attractor A to be a manifold of co-dimension 1. For a 2d non-autonomous dynamical system such as the one at hand, this may not be the case unless it could be reduced to a single non-autonomous ODE, by a reparametrization of state variables and τ . In the case such a reduction is possible, the attractor A can be formally defined using, instead of the flows φ τ , the map A(w) : U → R for some U = [w min , w max ] ⊆ R, where for instance, w := w(T , τ ), and A(w) is an algebraic function ofπ subject to the boundary conditions A(w → w max ) = A e and A(w → w min ) = A s (3.11) such that (T e , A −1 e , τ max ) ∈ B A e and (T s , A −1 s , τ min ) ∈ ∂B e 9 . We keep in mind that (T s , A −1 s , τ min ) is basically a saddle point in second order hydrodynamical theories in which the underlying dynamical system entails a term proportional toπ 2 . Finally, the one dimensional manifold of A can be represented by (3.12) Let us introduce the parametrization w = tanh ρ/T of Gubser time and the function A(w) defined as . whereπ = 3A(w)+2 from the conservation law (2.12). This reduces the dimensionality of the problem only asymptotically (e.g. ρ → ±∞) to one as opposed to the Bjorken model where a truly one dimensional ODE was achieved via a similar trick in [25]. Therefore it is expected that the attractors for the Gubser flow in all the hydrodynamical schemes are 1d non-planar manifolds and correspondingly the basins of attraction are three dimensional. We analyze this below briefly in the context of dynamical systems for the IS theory but bear in mind that a similar line of thought can be applied to any other theory as well.

From the perspective of 3d autonomous dynamical system
In this section we study the linearization problem by enhancing the 2d system to a 3d autonomous one, determine the stability conditions, and show that the attractor is a non-planar manifold of co-dimension 2. We will sketch the proof of the exponential asymptotic stability of the steady state in Gubser flow and create a local Lyapunov function to estimate the basin of attraction that is conjectured to be crucial in the study of an effective field theory of stable hydrodynamical and non-hydrodynamical modes.

Linearization and exponential asymptotic stability
The non-autonomous system (3.1a)-(3.1b) can easily be extended to an autonomous system by considering the time τ as an additional variable. Therefore the IS theory of Gubser flow is a truly 3d autonomous system of ODEs given by This is a polynomial vector field whose fixed points are given by Because of the symmetry (T ,π, τ ) → (−T ,π, −τ ) of the 3d problem (3.15), we are left with only three fixed points A = (0, 1/ √ 5, 1), B = (0, −1/ √ 5, −1), and C = (−38c/15, 2, 1). We now solve the linearized system around any fixed point, namely .
One can then immediately see that A is a sink (all eigenvalues are negative), and B, C have two positive eigenvalues, thus making them saddle points, as expected and since the eigenvalues all are nonzero, the fixed points are hyperbolic. Hence, we can apply the Hartman-Grobman (HG) theorem [73,74] 10 , that allows the local flow structure (phase space portrait) near a hyperbolic fixed point to be topologically equivalent to the flow diagram of its linearized system. For the Gubser flow in the IS theory, the flow diagram shown in Fig. 2 confirms the HG theorem in the vicinity of all the fixed points. This figure also portrays the state of flow lines on the time slice τ = 0.7 of the three dimensional phase space. Using the HG and the fact that A is an exponentially asymptotically stable fixed point of Gubser flow, it is then easy to write down the 10 This mathematical theorem states the behavior of a nonlinear system of differential equations in a domain near its hyperbolic equilibrium points is qualitatively the same as the behavior of the linearized version of these differential equations near this fixed point. Hyperbolicity means that no eigenvalue of the Jacobian matrix associated with the linearized dynamical system has real part equal to zero. Therefore, one can use the linearization of the original dynamical system to analyze its behavior around equilibria.
Lyapunov exponents of phase space variables along all the trajectories converging to the attractor by the eigenvalues given above, that are which matches the result obtained from the 2d analysis previously. The attractor may therefore be parametrized vectorially in the phase space spanned by (u 1 , u 2 , u 3 ) as for some constant ρ 0 . Evidently only for large ρ, τ (ρ) tends to one exponentially faster than do the other two converge to the fixed point i.e., λ τ < λπ < λT such that u 3 is just a point in the phase space and Eq. (3.19) at ρ → ∞ lies on the plane spanned by {u 1 , u 2 }, and is well-approximated bŷ given the initial conditions recovered from the asymptotics.
Since the non-autonomous system (3.1a)-(3.1b) is exponentially asymptotically autonomous for the attracting fixed point, with the limiting functions g i (T ,π) being the RHS of the two equations resulted from tanh 2 ρ ∼ 1, Theorem 2.2 of [75] applies where solving for the attractor of either 2d or 3d systems would yield the same result. The main assumption is to suppose a 1-to-2 transformation of the form t = ±e −2ρ where for each value of ρ there is a pair of values for t, such that the steady state is not on the boundary of the basin of attraction, B e . Then ρ = − 1 2 log |t| in both cases. Also, let F 3 (T ,π, t) := −2t where the index 3 denotes the RHS of the 3rd equation in the 3d problem. We define for i = 1, 2 (the RHS of ith state equation) where we have shifted theπ for later purposes. With this new time parametrization, Eqs. (3.15) may be cast into the form We note that the Jacobian matrix of this 3d system, Jac(F i ) has a block form at (T ,π, 0) and the origin is a fixed point, i.e. F i (0, 0, 0) = 0, and Jac(g i ) (0,0) has negative eigenvalues. Since the original system is exponentially asymptotically autonomous, and g i (0, 0) = 0 the exponential asymptotic stability immediately follows.
Lastly, a few remarks are due. There is a sense in which one must consider the aforementioned coordinate transformation inspired by the asymptotics of the series solutions to the flow equations (3.1a)-(3.1b) mainly the study of resurgence properties and transseries [26,27]. Focusing on just a simple series solution is problematic on its own because of divergence issues but above all else lies the fact that there is a challenge to pick a good expansion variable for the Gubser flow due to the peculiarities of de Sitter time. The natural choice for such a series expansion should at first glance be either tanh ρ or coth ρ, but they come with a caveat; they never grow bigger than 1 as ρ → ∞ thus not useful from the perspective of series asymptotics. It turns out that the exponential asymptotic stability of the fixed point offers a suitable and rather natural choice of the expansion variable which will be briefly touched upon in Sec. 4.2.2. It suffices to state that the rate at which the converging flow trajectories approach the sink is known from the analysis above to be exponentially fast and therefore one can expand around this fixed point by considering a variable of the type 1/t or any arbitrary positive real power of it.

Estimate of basin of attraction locally and globally
A path integral analogy: Before going into any details, we seek to motivate the reader to think about the question of why the basin of attraction and Lyapunov functions are extremely important concepts from a physical standpoint.
In the Feynman path integral formalism of quantum field theory and quantum mechanics, one often encounters integrals of the sort (3.23) where S[φ] is the action functional of the underlying theory that involves some field φ and the manifold M defines the space of fields or paths over which the integral is performed. For keeping the generality of the problem, we complexify φ and thus S[φ] is complex-valued and M is a middle-dimensional manifold in the space of complex fields. We will come to the dynamical system interpretation of this integral in a moment. But before, we have to note that the main burden is to understand what M is made of and how to find it in a general field theory is extremely difficult. Qualitatively, M is considered to be built out of the union of all the paths each connected to a critical point (stable fixed point of S[φ] that satisfies some equation of motion subject to stationary action principle) which are the physically relevant paths (steepest descent paths) in the phase space or configuration space over which the field φ takes values. So the most favored configuration is the one along which Re(S[φ]) increases as one approaches the critical point attached to it by some dynamical flow involving the action functional and an appropriate flow time t f defined over the real line, for instance. This is described by an ODE of the form where the RHS is nothing but the variation of action functional with respect to the dynamical field φ which now depends on t f and the bar represents complex conjugation.
There is in principle no universally well-established picture where we have an effective field theory approach to hydrodynamics, we rather choose to go to the phase space of state variables and flow time to make our analogies with the picture given above. But what about the action functional in the phase space? It is kind of obvious that the Lyapunov function(al) V that depends on the phase space parameters of the underlying hydrodynamical system acts like an effective action for all the stable thus relevant hydro modes. Once defined, V is always positive definite, and more importantly its derivative with respect to the flow time has to be always non-positive i.e. dV dρ ≤ 0, (3.26) similarly to the negative of real part of action functional decreasing on M due to stability conditions set by the properties of (3.24) 13 . Given a global Lyapunov function, we therefore can write an effective action functional for the Gubser flow in any theory where x(ρ) = (T ,π, t), and c is a 'coupling constant'. Note that this effectively describes the underlying hydrodynamical system as a theory of phase space variablesπ,T , t taking values in B e -the basin of attraction for the stable fixed point. Finally the partition function can be formulated as where M was defined in (3.25). This can be generalized to any theory of hydrodynamics which probes far-equilibrium aspects as well. In what follows, we create a local Lyapunov function and review a few techniques of obtaining a global one that practically speaking is the ultimate goal of taking this approach. In an upcoming work, we explore this effective partition function and to what extent it captures the right properties of the original (relativistic) hydrodynamics. Local Lyapunov function: We showed in the previous subsection that the dynamical ODEs in (3.22) limit to a 2d autonomous system close to t = 0,

29)
13 One can think of the effective action as a Morse-Bott functional, which satisfies dRe(−S)/dρ ≤ 0. 14 For a similar discussion on the relation between the effective potentials and Lyapunov potential functions based on gradient flow equations, check [80].
where (0, 0) is an exponentially asymptotically stable fixed point similarly to (0, 0, 0) being the same for the 3d autonomous system.
A local Lyapunov function V loc (T ,π, t) can estimate the basin of attraction near the fixed point (0, 0, 0) of (3.22) that contains the origin inside its basin. Like the global Lyapunov function, V loc (T ,π, t) is a continuous positive-definite scalar function V loc (T ,π, t) defined on the set D = {|t| > 0,T ∈ R,π ∈ R}. V loc has continuous firstorder partial derivatives at every point of D. Here we first focus on a local construction of the basin of attraction, where a local Lyapunov function V loc = x T P x should exist such that x T = (T ,π, t) and P satisfies the Lyapunov equation [81] Jac(F i ) T (0,0,0) P + P Jac(F i ) (0,0,0) = −diag. (1, 1, 1), (3.30) Solving this equation for P and inserting the resulting matrix in the formula for V loc , we obtain (3.31) As discussed before, along every flow line in the basin of attraction, V loc should effectively decrease with time, therefore dV loc (T ,π, t) dρ = 1 2204c 2 (|t| + 1) − 58c 2 38t 2 (|t| + 1) + (|t| − 1) 2((10 + This provides the main restriction on the shape of the basin. We note that the function V loc defines a symmetric basin of attraction under t → −t which is consistent with Proposition 2.9 of [75]. The two-time coordinate patch (T ,π, ±e −ρ ) does in fact tell us that from the eyes of an observer sitting at the origin, there is a mirror symmetric copy of the attractor, say A R , with respect to the t = 0 plane to which all the flow lines starting at t 0 will be seen to converge as well. This is nothing new other than a symmetry of the new coordinate system. The left copy A L is depicted in Fig. 3 along with the basin of attraction corresponding to the Lyapunov function in (3.31). Finding a Lyapunov function that captures the global basin of  (3.31). The boundary of the basin is topologically isomorphic to a deformed two-sphere which is locally a ρ = const hypersurface. For better visibility, the basin shape and size are set by the condition −100 < dV loc /dρ ≤ 0 and V loc < 5 respectively, which can be set, for instance, based on a near-equilibrium energy scale of the hydrodynamical modes. In the new coordinate system the origin (red dot) is the stable steady state that corresponds to the point (0, 1/ √ 5, 1) in the original coordinate system and, locally speaking, all the initial conditions chosen from the inside of this basin will yield converging trajectories towards the origin at late times. We remark that the heart-shaped gap mimics the asymptotics of boundary flow lines converging to the stable fixed point from below on the boundary of the basin of attraction in Fig. 1. The black line represents the IS asymptotic attractor solved numerically.
attraction for non-autonomous systems of higher dimensions (> 1) is in general a hard problem and in most cases an analytic result cannot be obtained unless the dynamical system entails some nice properties due to hidden symmetries that could give rise to further simplifications.
Sum of squares (SOS) polynomials is another recent method developed based on global optimization for many applications including Lyapunov stability analysis and control theory [82]. This method attempts at seeking a sum of the form solving an n dimensional optimization problem for some positive integer N where g i (x) are all polynomials. Here all x i depend on a time parameter ρ. If there is a λ ∈ R that allows a quadratic expansion of f (x) = x T P (λ)x in terms of monomials for some n × n matrix P (λ) such that P factorizes as P (λ) = QQ(λ) T for a rank-N matrix Q, then f is said to have an SOS decomposition of the form (3.32). For f to be a Lyapunov function V, one more condition comes from the fact that −dV/dρ has to be SOS as well.
Another method that is similar in spirit to SOS programming calls for a solution of the PDE [75,83] V where ||..|| shows the usual distance from the origin assuming that it defines a stable fixed point for the dynamical system. The differentiation is implemented with respect to a flow time ρ. It has been comprehensively studied in literature that in various (non)linear dynamical systems, global Lyapunov functions can be well approximated by solving Eq. (3.33) using radial basis functions, see for instance [81]. We leave this to a future work and rather show a numerical plot of a 2d surfaces in the basin of attraction for different initial values of τ in Fig. 4. Let us summarize the key take-home lessons learned in this section: 1. For any value of flow time, as long as the initial values initiate the flows in the basin of attraction, the flow lines will always go toward the attractor; 2. The attractor is an invariant set of numbers toward which the flow lines evolve inside the basin of attraction at late times. If there are no repelling fixed points around, with an underlying regular geometry, the rate of convergence to the attractor for all the flows should be uniform; 3. The flow line that begins at the saddle point located on the boundary of the basin, goes to the attractor at the fastest possible rate among all other flow lines because it gets propelled by the repulsive force of this point. For a flow starting from (T (ρ 0 ),π(ρ 0 ), τ (ρ 0 )), this rate can be quantified by which is literally the length of the curve of velocity vector field or flow line.

4.
Out of all the flows starting at the same time-slice, say τ = −1, the ones closer to the saddle points around, will have more repelling force and thus a faster convergence; 5. Less often it is mentioned that the role of basin of attraction is important in determining the (multi)stability and strength of the attractors as well as the usefulness of the dynamical systems in consideration even in the regimes far from (thermal) equilibrium. So knowing the basin of attraction, its topology and size, is necessary for the study of hydrodynamical flows and the search for attractors per se in a dynamical system with a stable fixed point and some unstable fixed points would not be illuminating. Things evolve toward stability, fast or slow as long as they are initiated with values in the basin of attraction; 6. The attractor line has the fastest rate of convergence due to the propulsion of saddle point it starts from and therefore naturally there is always a fast convergence of modes nearby even in the far-from-equilibrium regime of hydrodynamics; 7. Any effective theory of hydrodynamics can be written as a simple kinetic term with a Lyapunov potential functional that describes the relevant physical modes that only belong to the basin of attraction. In the effective field theory language, anything outside of the basin of attraction is completely integrated out.

Universal asymptotic attractors for different dynamical models
The relaxation equations of the different components of the energy-momentum tensor derived from a particular microscopic theory do not lead necessarily to the same attractor. In this section we present the numerical results of calculating the attractors for different hydrodynamical models and exact solution of the RTA Boltzmann equation.
In order to determine the attractors we follow closely the procedure outlined in the work of Heller and Spalinski [25] which is based on the slow roll-down approximation [84]. The reader should bear in mind that the attractors calculated in this section are asymptotical since the basin of attraction is three dimensional as we argue in the previous section.
The equations of IS, DNMR and aHydro can be combined into an unique ODE for the function A(w) (3.13) which reads as These expressions were determined by using the conservation law (2.12) together with Eqs (2.13a) and (2.13b) for IS and DNMR respectively while for anisotropic hydro we consider Eqs. (2.23) and (2.25). We point out that for aHydro it is necessary to rewrite the function F(π) (2.29) in terms of A(w). When applying the slow roll down approximation dA/dw = 0 in Eq. (4.1) one needs to find the roots of H(A(w), w) ≡ 0. We take the late time or asymptotic limit tanh 2 ρ ∼ 1 which was explained before and it poses no problem in the universality of the tor due to exponentially fast convergence of flow lines to the attractor at late times. In the case of IS and DNMR this constraint gives two different solutions so one chooses only the stable one A + (w) [25,26]. For these theories we get

IS:
A (4.3b) The initial condition for solving the differential equation (4.1) for IS and DNMR is obtained by evaluating lim w→−∞ A + (w) ≡ A i 15 . Afterwards, as mentioned before, one simply solves Eq. (4.1) while taking ρ → −∞ (and thus coth 2 ρ − 1 → 0). For the case of aHydro the initial condition is determined by finding numerically the roots of H(A(w), w) (4.2c) which gives A i ≈ −0.75. In the next subsection we discuss the numerical solutions of Eq. (4.1) for each approximation scheme.

Numerical results
We are now ready to discuss the numerical solutions of the attractors (late-time asymptotic behavior) for each truncation scheme together with the exact attractor of the Gubser solution (2.31). The exact location of the attractor was found by following the technique explained in Ref. [35]. For completeness we briefly explain it in App. B. When solving numerically we used ρ 0 = −10, which is good enough for our purposes because this value avoids unphysical behaviour e.g. negative temperatures [13,85].
In Fig In each plot we chose c = 5η/s with η/s = 3/(4π). The set of initial conditions were chosen by allowing the initial condition of the effective shear to be either prolate (π < 0 and thusP L <P ⊥ ), or oblate (π > 0 and thuŝ P L >P ⊥ ). The former configuration corresponds to A i < −2/3 while the latter indicates A i > −2/3. The initial values A i 's of the different asymptotic attractors (late-time tail of solid red line) in Fig. 5 are always below A i < −2/3 so that the steady state attractor of fluids undergoing Gubser flow corresponds to a prolate configuration. This result is valid independently of the hydrodynamical model and in agreement with 15 Alternatively, one can also predict the value of the initial condition for A + (w → −∞) from the method discussed in Sec. 3. In this case one equates the value ofπ(ρ → −∞) at the stable fixed point into the conservation law (2.12). For instance, for c = 15/(4π) the values ofπ(ρ → −∞) at the stable fixed points for IS and DNMR are 1/  the methods discussed in Sec. 3 since the saddle fix points are located precisely when π ρ→−∞ < 0. We also notice that independently of the theory, prolate configurations merge faster to their corresponding attractors than the oblate ones. Furthermore, the positive definite condition of the transverse and longitudinal pressures implies that −3/4 < A(w) < −1/2. The positivity of the pressures is satisfied only by aHydro and the exact Gubser solution and thus one cannot initialize A(w) below these values. This statement was tested numerically and explains why there are no gray lines below the early time line of attractor in the bottom panels of Fig. 5. However, for the IS and DNMR evolution equations this condition can be slightly broken by initializing A(w) below the attractor as it is shown in the top panels of Fig. 5. This indicates that the basin of attraction of IS and DNMR theories is larger than aHydro and the exact solution while being physically invalid. If one restricts IS and DNMR to satisfy the positivity condition of the pressures, it would physically smallen the phase space of initial conditions and thus, the basin of attraction [86]. A natural question which arises in our context is: What is the best model that matches the underlying exact microscopic kinetic theory? We answer this question by plotting in Fig. 6 the attractors of different hydrodynamical theories together with the one obtained for the exact Gubser solution (2.31). First, we observe that none of the truncated approximated schemes -DNMR and IS -are able to be in good agreement with the exact attractor over the entire w regime studied here. One might be at first surprised that none of the hydrodynamical truncation schemes do work even at large w when the system supposedly reaches its thermal state. However, the Gubser flow does not reach this state asymptotically since the expansion rateθ = 2 tanh ρ saturates at large ±ρ without vanishing exactly. Among these two hydrodynamical truncation schemes, we find the IS to be closer to the numerical values of the exact attractor than DNMR albeit still unable to match it to high numerical accuracy. Now, the best theory to describe the exact attractor shown Fig. 6 is aHydro. A closer look shows that all the hydrodynamical models are not able to match the exact result in the small w region (−0.8 w 2) as it is shown in the inset of Fig. 6. However, the numerical difference of IS and aHydro with respect to the exact result is no larger than 4% in this w interval while DNMR deviates entirely in this regime of w. In the large or intermediate regime, on the other hand, we verify numerically that the largest numerical deviation between the aHydro attractor and the exact one does not exceed 0.06%. The numerical results presented here provide a conclusive prove that aHydro resums effectively the Knudsen and inverse Reynolds numbers to all orders independent of the initial conditions.
We point out again that the notion of "attractor solution" is ill-defined and one should not care about what occurs in mid-range or early time regime of w because the attractor is actually an statement about late-time asymptotics of the flow lines. We can only say that the attractor solution is just a solution to some system of ODEs with a given initial value that is located exactly at the saddle point on the boundary of the basin of attraction. But because there are two stable and one unstable directions at this point (one of the eigenvalues of Jacobian matrix being positive means unstable direction), then being a metastable, the saddle point will turn out to initiate a flow inside the basin of attraction and finally the attractor will absorb it. We point out that this flow line is the fastest to converge to the stable fixed point and thus the attractor. In general this is an invariant set of numbers that flow lines converge to at very large w and as can be seen in Fig. 6, even for w > 0 there is no separation between flow lines due to the exponentially fast convergence. This is way better than what we could have asked for from the approximation tanh 2 ρ ∼ 1 that was applied to get A(w) in (4.1).

Asymptotic perturbative series expansion
The numerical comparisons between different models discussed previously demonstrated that the IS and DNMR theories cannot describe the exact steady state attractor. This failure is somehow expected since both theories are derived from an asymptotic series expansion of the distribution function. Here we show strong evidence of the divergence of this series and we briefly comment on how to fix it using resurgence techniques.

Divergence of IS and DNMR theories
Since w ∈ (−∞, ∞), one can propose an asymptotic series ansatz of the form that is claimed to solve Eq. (4.1) in the asymptotic limit ρ → ±∞ with H(A(w), w) given by Eqs. (4.2a) and (4.2b). After equating this ansatz to Eq. (4.1) one obtains the following recursive relation for the kth coefficient (k ≥ 2) IS: . The truncation of this power series at n = 30 is sufficient for the purpose of convergence tests in general. In Fig. 7 we present the numerical solution of the recursive equations (4.5) for both IS and DNMR theories when c = 15/(4π). From this figure we observe that for both of these hydrodynamical models so the root test [87] fails for the ansatz of A(w) (4.4) and thus, the asymptotic series expansions of the IS and DNMR theories are divergent. Fig. 7 also indicates that in this asymptotic regime the divergence in the DNMR case is more severe. This is due to the new term in Eq. (4.2b) which survives in the large w limit instead of converging to zero and thus any divergence in A(w) in the IS theory will be magnified in the DNMR theory. This result also confirms earlier numerical studies of the Gubser flow (see for instance Fig. 3 in Ref. [13]) where it was found that in the large ρ + limit the DNMR solution does a poorer job than IS when comparing the predictions of both theories with the RTA exact solution (2.31).

Resurgence to the rescue
The true nature of the attractor from the perspective of solving for a series solution would be unveiled by studying the resurgent aysmptotic expansions around this solution which mimics the instanton corrections on top of perturbation theory in quantum field theories. This can be packaged into a formal exponential series solution known as transseries accounting for all sorts of corrections such as exponential, logarithmic and so on [88]. For the system of equations (3.1a)-(3.1a) with the substitution τ = (1 − e 2ρ )/(1 + e 2ρ ) an analytic solution of this sort can be written as [89] T (e ρ ) =T (0,0) (e ρ ) + σ 1 1 σ 0 2 e m 1 ρ e −|λT |e ρT (1,0) (e ρ ) + σ 0 1 σ 1 2 e m 2 ρ e −|λπ|e ρT (0,1) (e ρ ) + . . .

Conclusions
The properties of the attractors of different hydrodynammical systems undergoing Gubser flow were studied within relativistic kinetic theory. Our work extended previous studies of attractors by incorporating techniques and tools of nonlinear dynamical systems. We extensively discuss the application of these methods to the IS evolution equations to investigate the stability properties of the fixed points, the flow lines around those and the Lyapunov exponents. These techniques were proven to be strikingly useful in the understanding of the reason behind the impossibility to reduce the 2d system of ODEs of each hydrodynamical scheme into a single one universally as it was done in the Bjorken case [25]. It was shown that this non-reduction is intrinsically related with the dimensionality of the basin of attraction and that the tanh ρ also acts as an independent variable of ρ that is responsible for the extra dimension of the phase space compared to the Bjorken case. Inspired by the exponential asymptotic stability and resurgent transseries type arguments, we defined a natural time variable, t = ±e −2ρ , which was expounded in subsection 3.2.2. In this coordinate system, the observer sees two copies of the same flow started at t → ±∞ that approach to the point (0, 0, 0) which is now the stable fixed point of the system, hence making the basin symmetric under t → −t. This method is a very well-known trick used in usually exponentially asymptotically stable systems to connect the basin of attraction in non-autonomous systems to that of autonomous ones as well as compactifying the basin in case the time-dependency was in the form of an unbounded function. It also brings the fixed point to the origin that is crucial for the construction of Lyapunov functions. This helped us build such a function for the IS theory and discuss the shape of the basin of attraction at least locally near the steady state equilibrium. This notion of the Lyapunov function from dynamical systems was shown to be related to an effective-action description of hydrodynamics motivated by Picard-Lefschetz theory. The attractors of the IS, DNMR and aHydro models were obtained via the slow roll-down approximation while the attractor associated to the Gubser RTA solution was determined through Romatschke's method [35]. From the numerical comparisons between the attractors of different theories with the exact solution we conclude: (a) Hydrodynamical models based on an asymptotic series expansion of the distribution function, aka IS and DNMR, are unable to provide a quantitative description of the attractor of the exact solution and, (b) the best agreement with the exact attractor is achieved by aHydro up to high numerical accuracy. For the IS and DNMR approximation schemes we showed that their corresponding asymptotic solutions to the respective differential equations are divergent. The origin of the divergences and how to resum them was briefly discussed by applying resurgence techniques.
The success of aHydro to describe the asymptotic attractor of the exact Gubser solution demonstrates that this theory is able to take into account both, large inhomogeneities in the fluid due to collisions -quantified by the Knudsen number-as well as big spacetime inhomogeneities of the macroscopic fluid variables -quantified by the inverse Reynolds number-. As a matter of fact, the match between aHydro attractor and the exact one demonstrates that aHydro resums the Kn and Re −1 to all orders. This resummation is carried out in a non-perturbative way by including the largest momentum-space anisotropies present in the plasma into the leading order anisotropic distribution function.
The very common techniques in the context of nonlinear dynamical systems presented in this work and applied to hydrodynamics, can be extended to study a long list of properties of other relevant more complex physical systems than the one studied here. This list may include an investigation of hydrodynamization processes and dynamics of attractors for more general nonlinear collisional kernels [7,22,90], holographic models [8,[14][15][16][17][18], spatially non-homogeneous expanding fluids [39]. On the other hand, it is also interesting to investigate the rich structure and topology of the basin of attractors in turbulent flows and other chaotic systems of interest. On a more theoretical subject, the possibility of formulating effective actions for hydrodynamics by exploring the analogy between the steepest descent directions in the path integrals and the flow lines starting at the boundary of the basin of attraction for the attractors opens a new perspective to re-formulate this old problem. Moreover, the issues addressed in this work shed more light on new questions that remain to be answered within the aHydro framework. For instance, in the resurgence program, the attractor is understood as the leading-order asymptotic transseries [43,91]. This mathematical statement together with our results suggest a highly non-trivial relation between the nonperturbative resummation of large Knudsen and inverse Reynolds numbers carried out by aHydro and a certain class of solutions of the nonlinear Boltzmann equation that can be expressed as a transseries. We leave these matters to future works.

A Anisotropic integrals
Here we calculate the anisotropic integralsÎ nlq that appear in this paper. The explicit forms of the functionsR nlq used in this paper are given bŷ R 200 (ξ) = 1 2 The momentsÎ eq nlq associated with the equilibrium distribution function are obtained by considering the ξ → 0 limit of Eq. The attractor for the exact solution (2.31) was found by considering Romatschke's technique [35]. For the Gubser case the slow roll condition, cf. Refs. [25,35,84], dA(w)/dw = 0 gives the following condition dA(w) dw = 0 ⇒ 1 4