Hydrodynamic attractors for Gubser flow

The Boltzmann equation is solved in the relaxation time approximation using a hierarchy of angular moments of the distribution function. Our solution is obtained for an azimuthally symmetric radially expanding boost-invariant conformal system that is undergoing Gubser flow. The solution of moments that we get after truncating the infinite set of equations at various orders is compared to the exact kinetic solution. The dynamics of transition from an early time collisionless free streaming to the hydrodynamic regime at late times is described by the presence of fixed points. The attractor solution is found for various orders of moments as an interpolation between these fixed points. For arbitrary initial conditions, it is found that the solution approaches the attractor when the non-hydrodynamic modes decay. The relation of moments to various approximations of relativistic viscous hydrodynamics is investigated.


Introduction
Hydrodynamics is an effective theory for the description of long-wavelength phenomena of fluids, that can be expressed as a small gradient expansion relative to a thermal background [1]. Thus, hydrodynamics is expected to fail for systems which are far-from-equilibrium where the gradients are expected to be large. The medium produced in pp collisions at LHC and RHIC energies is an example of such a system. However, recent experimental results of high energy pp collision have shown evidence of collectivity similar to those observed in heavy-ion collisions [2][3][4][5][6]. The unprecedented success of hydrodynamics to describe collectivity in heavy-ion collisions, as well as small systems, can be attributed to the fact that there exists a stable universal attractor which makes the dynamical equations to quickly converge and enter a hydrodynamic regime, at a time scale much smaller than the typical isotropization time scales [7][8][9][10][11][12].
Previous works [8,[13][14][15][16][17] have mostly focused on studying the properties of attractors for rapidly expanding 1+1d boost invariant systems undergoing Bjorken flow using relativistic kinetic theory. However, the fireball produced in high energy heavy-ion collisions also expands in the transverse direction at late times. Therefore, it is natural to ask whether the systems undergoing simultaneous longitudinal and transverse expansion also shows the attractor nature as observed for 1+1d expansion. In the present work, we consider a system undergoing Gubser flow which has a simultaneous transverse and longitudinal expansion. In the context of Gubser flow, some recent works [18,19] study such dynamical properties using anisotropic hydrodynamics or DNMR hydrodynamical theories [20]. The Email addresses: ashutosh.dash@niser.ac.in (Ashutosh Dash), victor@niser.ac.in (Victor Roy) purpose of the present work is to investigate the dynamics of transition from free streaming regime to hydrodynamic regime using a moment method that translates the kinetic equation for distribution function in relaxation time approximation (RTA) to an infinite series of coupled ordinary differential equations advocated in [21][22][23].

Formulation
The hydrodynamic equations can be derived from the transport equation by taking appropriate moments of the distribution function. The collision term in the Boltzmann equation plays an important role to isotropize any arbitrary initially anisotropic out of equilibrium distribution function or the anisotropy generated by a strong expansion. The competition between the two effects is commonly investigated in terms of the ratio between longitudinal pressure P L and transverse P T , local equilibrium corresponds to P T /P L = 1 . In [21] it has been shown that a particular moment of the distribution function is very useful for studying an out of equilibrium system undergoing longitudinally boost-invariant expansion. The details of the moment of the distribution function will be discussed later. The Boltzmann equation in an arbitrary coordinate system for on-shell particles is where Γ λ µ j are the Christoffel symbols and f = f (x µ , p j ) is the one particle distribution function.
Following [24,25] we look for solutions of the hydrodynamic equations with S O(3) q ⊗ S O(1, 1) ⊗ Z 2 symmetry in flat spacetime. This can by implemented by a Weyl transformation of Minkowski spacetime in Milne coordinates x µ = (τ, x, y, η) to dS 3 ⊗ R spacetime. Here 3 stands for the 3-dimensional de Sitter spacetime assuming that the fluid is homogeneous. This spacetime is described by the line element dŝ 2 = −dρ 2 + cosh 2 ρ dθ 2 + sin 2 θdφ 2 + dη 2 (2) Here we have introduced the de Sitter coordinatesx µ = (ρ, θ, φ, η), with q −1 is an arbitrary length scale and sets the size of the system. In these coordinates, the Gubser flow profile takes the form: u τ = − cosh κ(τ, r), u r = sinh κ(τ, r), with transverse flow rapidity κ = tanh −1 2q 2 τr/ 1 + q 2 τ 2 + q 2 r 2 . In this paper, all quantities in de Sitter coordinates are denoted with a hat. Gubser symmetry requires massless degrees of freedom of the fluid, i.e.,p 2 = 0 and thus the spatial momentum can be written aŝ The S O(3) q symmetry implies the distribution function to be independent of (θ, φ) and that the distribution function depends only on the following combination of momentum components: Similarly the S O(1, 1) symmetry subgroup implements longitudinal boost invariance as a result of which the distribution function is independent of η. Using the above constraints, the RTA Boltzmann equation Eq. (1) in de Sitter space takes the following form [26,27] Herep ρ = (p Ω / cosh ρ) 2 +p 2 η andT = τT . c is a dimensionless parameter which, in RTA, can be expressed in terms of the shear viscosity (η s ) to entropy density (s) ratioη s ≡ η s /s as c = 5η s .

Evolution of moments in Gubser flow
Although Eq. (7) can be solved numerically, much insight can be gained by taking various moments of the distribution function as elucidated in [21,22]. These moments capture the deviation of distribution function from isotropy and for an expanding system they approach local equilibrium at late times due to collisions. The n th order moment L n of the distribution function is defined as : where dP = dp η dp θ dp φ / (2π) 3pρ cosh 2 ρ sin θ is the phase space measure and P 2n is a Legendre polynomial of order 2n.  Except for the moment n = 0 which corresponds to energy densityˆ , all higher order n quantify the details of longitudinal momentum anisotropy. For example, the first order moment L 1 =P L −P T , describes the anisotropy of pressure components in longitudinal and transverse directions.
Taking the derivative of Eq. (8) with respect to de Sitter time ρ and using Eq. (7) one arrives at the following infinite set of coupled differential equations, where the coefficients a n , b n , and c n are given as As one can clearly see from Eq. (9), there is a competition between collisions (the term containing relaxation time) which washes out the effect of anisotropy and expansion (determined from the dimensionless coefficients a n , b n , and c n ) which drive the system out of equilibrium. One should note that a key difference between the present work and that of [22] is that the relaxation time τ R in a conformal setting is not a constant but is related to the only available scale in the system i.e., temperature through the relationτ R = c/T (ρ). A comparison of the coefficients a n , b n , and c n for Gubser flow against the coefficients a n ,b n . andc n for Bjorken flow of [22] gives us the following relations: a n +ã n = 4, b n = −b n and c n = −c n . The system of equations Eq. (9) can be easily solved by truncating the series at a given n and ignoring all higher-order terms. In Fig. 1 we compare the solution of these equations with the exact numerical solution of Eq. (7) by using a method given in [27]. First, we see that even for n < 2, the solution captures the qualitative details of the exact solution. Second, one notices that the approach to the exact solution is alternating for even and odd n. We would like to also point that, the exact kinetic solution shows numerical instability at negative ρ values when the system is initialized with non-zero shear stresses at ρ 0 = 0. The same problem also persists in the moment method when we try to extrapolate the solution to negative de Sitter coordinates.

The free streaming fixed point
The RTA Boltzmann equation Eq. (7) has the following exact solution [26,27]: where tion, f 0 ρ 0 ;p 2 Ω , p η is the initial distribution function at de Sitter time ρ 0 and f eq ρ ;p 2 Ω ,p η is the equilibrium distribution function. The moments can be calculated from the distribution function Eq. (13) using the definition given in Eq. (8).
Assuming that the initial distribution f 0 ρ 0 ;p 2 Ω ,p η at ρ 0 is an isotropic Boltzmann distribution, we find the following free streaming solution (taking c → ∞ in Eq. (13)) for the moments where x = cosh ρ 0 cosh ρ and F n (x) as HereT 0 is the initial temperature at ρ 0 . The function F n (x) has the following limits: F n (x) → 1/2 as x → 0 and F n (x) → sin 2nπ 2nπ(1+2n) as x → 1 (i.e, F 0 = 1 for n = 0 and F 0 = 0 for n 0). Consequently, for asymptotically large de Sitter time, i.e., |ρ| |ρ 0 |, F n (x) is a constant and L FS n (ρ) decays as 1/cosh 2 ρ. However, the scaled moments, saturate in the free streaming limit, which would have decayed to zero in the presence of collisions.
In Fig. 2 we show the analytical solution of the evolution of scaled moments by solving Eq. (15). As one can see, all the scaled moments approach unity at large de Sitter times. The free streaming regime can also be described by Eq. (9) by taking the limit c → ∞. We have also checked (not shown here) that even for arbitrary initial condition, the moments approach unity both from above and below at large de Sitter times.
Since the moments continuously evolve to a constant value at large times, it is natural to define the quantity g n (ρ) [21,22], For Eq. (9), g n (ρ) turns out to be, Taking the limit c → ∞, and using the expression Eq. (17) for the ratio of the moments (for large ρ) one can obtain the free streaming fixed point. We find that for all n the solution yields g n (ρ) = −2. For the above initial condition, one finds that solution of g n (ρ) from Eq. (9) does not evolve with ρ and hence is indeed a fixed point. One may also verify that for arbitrary initial conditions, the series of equations Eq. (9) also gives the approximate result g n (ρ) ≈ −2 with error due to finite truncation of the series. Alternatively, one may define an equation of motion for the quantity g n (ρ), The fixed points correspond to the zeros of β(g n ). For n = 0 and assuming vanishing L 2 and higher order moments, β(g 0 ) turns out to be The plot of the function β(g 0 ) is shown in Fig. 2(b). We find that there is yet another fixed point g 0 ≈ −3 apart from the previously found point g 0 ≈ −2. Applying linear perturbation around the fixed point g 0 , i.e., g 0 (ρ) = g 0 + f (ρ), and substituting this in Eq. (20) we find that the solution around g 0 = −2 yields decaying modes and hence a stable fixed point where as the solution at g 0 = −3 yields growing modes and therefore an unstable fixed point. If one needs to keep the second moment L 2 in the expression for β(g 0 ) in Eq. (21), we get For the stable fixed point g 0 = −2, we get the value L 2 /L 0 = 1 as expected. For the unstable fixed point g 0 = −3, L 2 /L 0 = 3/8. The correction to β(g 0 ) because of the term L 2 is also shown in Fig. (2b).
It is interesting to note here that the late time behavior of the ratio of moments for the unstable fixed point g n = −3 (for any arbitrary n) turns out to be . One notices that the late time behavior of the ratio of moments for stable and unstable fixed points in the Bjorken scenario is exactly the opposite of Gubser flow. This is because if we substitute x → 1/x in the expression for F n (x) in Eq. (16) and then take the ratio of moments, they happen to be exactly like Eq. (23).

The hydrodynamic fixed point
Previously, we have seen that for an initially isotropic equilibrium distribution, the collisionless free streaming solution drives the system to anisotropic distribution with a wealth of moments being generated as time goes on. Owing, to the presence of fixed point that we have already seen, the time dependence can be qualitatively captured by the two lowest moments.
As we have discussed earlier, the effect of collision is to wash out the effect of initial anisotropy and the anisotropy generated due to expansion. For a 1+1 dimensional expansion, the late time behavior of the system can be described as a series expansion in the power of expansion scalar 1/τ. The late-time behavior in such a system has been attributed to the presence of a different kind of fixed point, called hydrodynamic fixed point. To proceed, we shall assume that all the moments L n admit a gradient expansion in powers of tanh ρ [19]: where γ (m) n 's are coefficients of the expansion. Due to Landau matching condition, for all orders the zeroth moment corresponds to the energy density i.e., γ obtained from Eq. (9) by ignoring the contribution from L 1 . We find that the ideal hydrodynamics limit of the time evolution of (ρ) is governed by the term a 0 = 8/3, i.e.,ˆ (ρ) ∼ (cosh ρ) −8/3 . The leading terms in the expansion can be determined by demanding the cancellation of relaxation terms in Eq. (9), since other terms can be ignored at large times. The above condition can be met if we have Starting from n = 1, the above equation leads to a recursion relation for L n from which the coefficients γ (n) n can be deduced as Since,ˆ (ρ) ∼T 4 (ρ), the time dependence of leading order is given as: γ (n) n ∼ (cosh ρ) −2(4−n)/3 . Substituting, this in Eq. (24) gives us the time dependence of L n ∼ (tanh ρ) n (cosh ρ) −2(4−n)/3 . The fixed points can be calculated using the definition of g n (ρ) in Eq. (18), Eq. (27) shows that unlike the 1+1d Bjorken case where the g n (ρ) is a constant, for 2+1d Gubser flow the hydrodynamic fixed point is a function of time. However, in the asymptotic limit, we have g n (ρ → ±∞) = 2(n − 4)/3. It is also interesting to note that for n = 1 at large time, the hydrodynamic and free streaming stable fixed point are exactly same. The next to leading order correction to γ (n) n can be determined by ignoring the contribution from the term L n+1 /L n but still keeping the constant contribution in Eq. (18). This gives The hydrodynamic attractor Following [22], we define the 'attractor' as the solution of Eq. (9) that joins smoothly the free streaming and hydrodynamic fixed point. In Fig. 3, the numerical solution of attractor is shown for n < 3 in terms of the quantity g n . While solving numerically we used ρ 0 = −10, because this value avoids unphysical behavior e.g. negative temperatures [27] if one initialise at ρ 0 = 0. The solution is obtained, for a particular c = 25/(2π) by truncating Eq. (9) after 30 terms. Thin solid black lines correspond to arbitrary initial conditions, while the thick solid green line and red dotted line corresponds to initial condition pertaining to stable free streaming fixed point Eq. (17) and unstable free streaming fixed point Eq. (23) respectively. It is really interesting to observe that solutions for arbitrary initial conditions merge with the hydrodynamic attractor after the non-hydrodynamic modes have decayed around de Sitter time |ρ − ρ 0 | ∼ c.
At this point, it is worth emphasizing that for modes n 0, the hydrodynamic fixed points diverge when ρ = 0 owing to the presence of a pole in Eq. (27).

Matching to hydrodynamics
Here we pause and compare the various orders of moments generated from truncation of Eq. (9) to successive orders of viscous hydrodynamic corrections. In the last section we have already shown that truncation at lowest order with vanishing L 1 results in ideal hydrodynamic equation of motion. Truncation at n = 1 yields two coupled set of equations: whereπ = −c 0 L 1 . The above set of equations constitute the second order viscous hydrodynamic equations for Gubser flow in the Chapman-Enskog expansion [28]. The term b 1 , a 1 are related to the the termβ π ,λ π through the relation:β π = −b 1ˆ /2, λ π = a 1 respectively and as we have already mentionedτ R = c/T (ρ). In the first order Navier-Stokes approximation the Eq.(29) reduces for small c:π = (4/3)τ Rβπ tanh ρ.
For third and higher-order viscous correction from the moment equation by keeping the term L 2 in Eq. (29) and also considering the time evolution of moment L 2 : In the above equation, L 2 appears as a new dynamical variable with its own evolution equation which in contrast to thirdorder hydrodynamics is related toπ by a constitutive relation. To make a connection with hydrodynamics it suffices to take a small c limit in which, where we have used the Navier-Stokes limit to express tanh ρ in terms of momentπ. Eq. (32) can be matched to the coefficient χ appearing in third-order hydrodynamic [28] through the relation:χ = −(3/4)c 1 b 2 . For higher-order corrections, the series γ (n) n is divergent since it grows as n! for large n as has already been seen previously in Bjorken flow.

Conclusion
To summarize, the moment method which has been formulated in [22] for 1+1d boost invariant system acts both as a practical tool for solving the kinetic equation and to address dynamical questions like the presence of fixed points and attractors, etc. for the system under consideration. In the present work, we applied the moment method for a system undergoing Gubser flow which has a simultaneous transverse and longitudinal expansion. Our study suggests that the presence of an attractor, to which the solution of the dynamical equations quickly converges before eventually reaching the viscous hydrodynamic regime is a feature not limited to the 1+1d system with an overwhelming amount of symmetries like Bjorken flow. We believe that these results offer detailed insights into the dynamics of longitudinal/transverse momentum isotropization in relativistic systems undergoing simultaneous transverse and longitudinal expansion. In passing, we would also mention that although Gubser flow has a transverse expansion, in itself is a highly idealized model when confronted with dynamics of matter produced in heavy-ion collisions. Therefore we think it is important to investigate the appearance of attractors further, by relaxing certain symmetries e.g. conformality and homogeneity.