On higher order and anisotropic hydrodynamics for Bjorken and Gubser flows

We study the evolution of hydrodynamic and non-hydrodynamic moments of the distribution function using anisotropic and third-order Chapman-Enskog hydrodynamics for systems undergoing Bjorken and Gubser flows. The hydrodynamic results are compared with the exact solution of the Boltzmann equation with a collision term in relaxation time approximation. While the evolution of the hydrodynamic moments of the distribution function (i.e. of the energy momentum tensor) can be described with high accuracy by both hydrodynamic approximation schemes, their description of the evolution of the entropy of the system is much less precise. We attribute this to large contributions from non-hydrodynamic modes coupling into the entropy evolution which are not well captured by the hydrodynamic approximations. The differences between the exact solution and the hydrodynamic approximations are larger for the third-order Chapman-Enskog hydrodynamics than for anisotropic hydrodynamics, which effectively resums some of the dissipative effects from anisotropic expansion to all orders in the anisotropy, and are larger for Gubser flow than for Bjorken flow. Overall, anisotropic hydrodynamics provides the most precise macroscopic description for these highly anisotropically expanding systems.


I. INTRODUCTION
A remarkable property of the hot and dense matter formed in ultra-relativistic heavy ion collisions at RHIC and LHC is a strong collective motion which has been successfully modeled using relativistic hydrodynamics (see [1] for a recent review). Dissipative hydrodynamics is formulated as an expansion in gradients of the fluid four-velocity, the simplest of them being the first-order Navier-Stokes theory due to Eckart [2] and Landau and Lifshitz [3]. Second-order dissipative theories developed later by Grad [4], Müller [5] and Israel and Stewart [6] cure an undesirable feature of relativistic Navier-Stokes theory, its acausality and instability [7,8]. These theories, based on the principle of non-negative entropy production, are formulated by assuming an algebraic form for the entropy-four current in terms of dissipative quantities. Unfortunately, this method does not provide a unique set of higher-order viscous evolution equations. This has motivated a broad spectrum of attempts to derive dissipative relativistic hydrodynamics from a more fundamental framework.
Hydrodynamics may be regarded as a macroscopic effective theory of a many-body system in which the complex interactions occurring over short distance and time scales are averaged out, and the effective degrees of freedom are a small number of conserved charge currents coupled to dissipative fluxes. For sufficiently weak coupling among its microscopic constituents, such a system can be described statistically by a more involved kinetic theory, based on a single particle phase-space distribution function f (x, p) whose evolution is typically governed by some generalized form of Boltzmann equation. The macro-scopic conserved currents and dissipative fluxes can be formulated in terms of momentum moments of this distribution function for which equations of motion are then derived from the Boltzmann equation. Closing the set of moment equations requires approximations to truncate the resulting moment hierarchy. Different such approximation schemes result in different sets of hydrodynamic equations. The validity and accuracy of the applied approximations can be judged by comparing, for specific highly symmetric situations in which the underlying kinetic theory can be solved exactly, the solutions of the different hydrodynamic approximations to the corresponding momentum moments of the exact microscopic solution [9][10][11][12]. This idea has generated increased interest for the search of new exact solutions of the relativistic Boltzmann equation [9,[13][14][15][16][17][18][19][20].
The equilibrium distribution function in the local rest frame (LRF) of a fluid is, by definition, isotropic in the momentum space, irrespective of the macroscopic motion of the fluid. Using it as an approximation for the true LRF distribution function in rapidly expanding systems is justified only in the limit of vanishing mean free path, i.e. instantaneous local thermalization. For realistic systems with small but non-zero mean free paths this approximation fails to properly account for the competition between microscopic scattering processes driving the system towards local momentum isotropy (and eventually into local thermal equilibrium) and the macroscopic expansion rate which drives the local phase-space distribution away from local thermal equilibrium (and, in the case of anisotropic expansion, also away from local momentum isotropy). This leads to a deviation δf (x, p) of the distribution function from its local equilibrium form, f (x, p) = f eq (x, p) + δf (x, p), with the relative size of δf increasing with the Knudsen number, i.e. with the product of the microscopic mean free time between collisions and the macroscopic expansion rate θ = ∂ ·u where u µ (x) is the fluid's flow four-velocity.
The first attempt to include such non-equilibrium δf effects in the distribution function f (x, p) was based on Grad's 14-moment approximation [4,6,21]. However, this moment expansion does not follow systematically from the underlying kinetic theory, such as the Boltzmann equation. A systematic approach of obtaining viscous hydrodynamics, to any given order in gradients of the macroscopic flow velocity, is based on a Chapman-Enskog-like iterative solution of Boltzmann equation [22][23][24][25][26][27][28][29][30][31]. Recently, this method was employed to derive higher order dissipative hydrodynamic equations [29]. Another novel way of formulating hydrodynamics from kinetic theory is based on an expansion controlled by the Knudsen number and the inverse Reynold's number [32]. For a conformal system and using the relaxation time approximation (RTA) for the collision term of Boltzmann equation, this approach leads to identical viscous evolution equations as obtained in [29], up to second order in gradients.
All these formulations, however, assume that the local deviations of f (x, p) from equilibrium are small, and that an expansion of f (x, p) about its equilibrium value to a few low orders in derivatives should suffice. Anisotropic hydrodynamics [11,20,[33][34][35][36][37][38][39][40] aims to extend the domain of applicability of traditional hydrodynamics, i.e., it attempts to better describe physical situations where the deviation of f (x, p) from local momentum isotropy is non-perturbatively large. This is achieved by explicitly including in the leading-order LRF distribution function an anisotropy parameter ξ describing the momentum-space deformation along the direction of largest anisotropy of the local expansion rate, and then expanding perturbatively the dynamical equations for the residual dissipative effects caused by the residual deviation δf , defined by writing f (x, p) ≡ f a (x, p; ξ) + δf (x, p). The non-trivial additional task in this approach is to determine the time evolution of the anisotropy parameter ξ non-perturbatively such that the residual dissipative effects encoded in δf are minimized and can again be described perturbatively. The recent works [12,39,41] have made significant progress in this direction.
It is necessary that the different macroscopic hydrodynamic formalisms described above are tested in scenarios where the microscopic dynamics can be solved exactly. We here study expanding systems with longitudinal boost-invariance and reflection symmetry, and either transverse homogeneity ((0+1)-dimensional Bjorken flow [42]) or azimuthally symmetric transverse density and flow gradients dictated by Gubser symmetry [43] ((1+1)dimensional Gubser flow [43,44]). For these highly symmetric flow patterns in each case a convenient system of coordinates can be found in which the macroscopic hydrodynamic flow appears static and the microscopic rel-ativistic Boltzmann equation, using the relaxation time approximation (RTA) for the collision term [45], reduces to an ordinary differential equation in longitudinal proper time τ [42] or de Sitter time ρ [43], respectively, and can be easily solved analytically [9,10,13].
In this work, we compare with these exact solutions the evolution of various macroscopic variables obtained using hydrodynamic equations obtained from the (perturbative) third-order Chapman Enskog (CE) approach [29,31] and the (non-perturbative) anisotropic hydrodynamic approach in the P L matching scheme [12,41]. The present work goes beyond similar earlier comparisons [10,12,31,41] by presenting for the first time the solution of third-order CE evolution equations for Gubser flow and a detailed analysis of the evolution of the systems' entropy in the various approximations (see [20] for an earlier study of entropy production in the isotropic FLRW universe). We find that entropy production is a sensitive discriminator between different hydrodynamic approximations and exhibits generically much larger deviations from the exact solution of the Boltzmann equation than all of the hydrodynamic observables. This reflects a significant contribution to entropy production by nonhydrodynamic modes whose dynamics is not constrained by macroscopic conservation laws.
The paper is organized as follows. In Section II we briefly describe the Bjorken and Gubser flow profiles and the coordinates we use to describe them. Section III reviews the exact solution of the Boltzmann equation in relaxation time approximation for the two flow profiles. In Sec. IV we elaborate on the Chapman-Enskog formalism and derive third-order dissipative hydrodynamics for Gubser flow. This is followed in Sec. V by a brief review of anisotropic hydrodynamics in the P L matching scheme for the Bjorken and Gubser flows. Numerical results for the comparison of the different approaches are presented and discussed in Sec. VI. We close with conclusions and an outlook in Sec VII.
In this paper, any quantity expressed in Gubser coordinatesx µ is denoted by a hat and made unitless by scaling with appropriate powers of the Weyl rescaling parameter (longitudinal proper time in Milne coordinates) τ [43,44]. For example,

III. EXACT SOLUTION OF THE BOLTZMANN EQUATION FOR BJORKEN AND GUBSER FLOWS
In this section, we review the central idea common to deriving from microscopic dynamics the dissipative hydrodynamic equations considered in this article. We consider a conformally symmetric system of weakly interacting massless Boltzmann particles without conserved charges whose phase-space distribution function f (x, p) evolves according to the Boltzmann equation. In the absence of external forces, and with a relaxation-time approximation for the collisional kernel, the Boltzmann equation has the form [45] p µ ∂ µ f = (u · p) δf τ r , where τ r (x) is the momentum-independent relaxation time and δf ≡ f − f eq is the deviation of the distribution function from its local equilibrium form Here β(x) ≡ 1/T (x) is the inverse local temperature and u µ (x) is the velocity of the local rest frame, defined as the velocity associated with the local energy flow (LRF = Landau frame). Conformal symmetry requires τ r = 5η/T ≡ c/T , where specific shear viscosityη ≡ η/s is defined the ratio of shear viscosity η to entropy density s.
On obtaining a solution of Eq. (7), either exact or in some approximation, the macroscopic hydrodynamic variables are constructed from the momentum moments of f (x, p). Specifically, the conserved energy momentum tensor T µν is the second moment of f (x, p) [46]: where we use the shorthand notation is the invariant momentum-space integration measure, with g being the determinant of the metric tensor.
In the remainder of this section we discuss the exact solutions of Eq. (7) for Bjorken and Gubser flows. Two other, more general methods of obtaining from Eq. (7) approximate solutions for f (x, p), namely the Chapman-Enskog iterative scheme and anisotropic hydrodynamics with P L matching, are presented in the next two sections, including again their specific forms for Bjorken and Gubser flows.
For massless systems with Bjorken symmetry the single particle phase-space distribution f (x, p) = f (τ ; p T , w) can only depend on the longitudinal proper time τ , the magnitude of the transverse momentum p T , and the longitudinally boost-invariant variable w = tp z − zp 0 = p T τ sinh(y−η) (where y is the kinematic rapidity of a particle) [9,13]. With this simplification, Eq. (7) reduces, at every point (p T , w) in momentum space, to an ordinary differential equation in τ , with the integral solution [9] f (τ ; p T , w) =D(τ, τ 0 )f (τ 0 ; p T , w) Here D(τ 2 , τ 1 ) = exp − τ2 τ1 dτ /τ r (τ ) is the so-called damping function, and the energy p τ is obtained from p T and w through the mass-shell constraint. The temperature defining the local equilibrium distribution under the integral in the last term is obtained from the Landau matching condition = (u · p) 2 = eq = (u · p) 2 eq = (3/π 2 )T 4 . This condition involves an integral over all momenta (p T , w) and renders the solution (10) highly nonlinear, in spite of its apparent simplicity.
For Gubser flow the symmetries constrain the dependence of the phase-space distribution f as follows: [10,15], andp η = w is the same boost-invariant longitudinal momentum variable as in the Bjorken case. Again, the symmetry constraints reduce the RTA Boltzmann equation (7) to an ordinary differential equation at each point (p 2 Ω ,p η ) in momentum space, with the solution [10,15] f (ρ; Again, the temperature in the equilibrium distribution on the right hand side is obtained by Landau matching, andp ρ is obtained from (p 2 Ω ,p η ) through the mass-shell constraint. The exact solutions (10,12) can be evaluated numerically [9,10], and the exact evolution of any macroscopic quantity (in particular of all the components of the energy momentum tensor (8)) can then be obtained by taking appropriate momentum moments of the exact f (x, p).

IV. DISSIPATIVE HYDRODYNAMICS FROM THE CHAPMAN-ENSKOG METHOD
This method is based on the assumption that the deviation of f (x, p) from its local equilibrium value is small, such that the RTA Boltzmann equation, Eq. (7), can be solved iteratively to obtain a Chapman-Enskog-like expansion for the non-equilibrium part of the distribution function in powers of space-time gradients [22,47]: where δf (1) is first-order in derivatives, δf (2) is secondorder, and so on. To first and second order in derivatives one obtains The above expansion may also be seen as a perturbation series in powers of the expansion parameter τ r . The energy-momentum tensor has the general form where and P are local energy density and pressure, respectively, related to the temperature by = 3P = 3/(π 2 β 4 ) through the Landau matching condition. ∆ µν ≡ g µν + u µ u ν projects a tensor to the space orthogonal to u µ , and the shear stress tensor π µν is traceless and orthogonal to u µ .
To close the equations (17,18) we need additional equations for the shear stress π µν . To obtain them we express π µν in terms of δf , If one substitutes on the r.h.s. for δf the first-order term (14) of the expansion (13) and uses the energy-momentum conservation laws (17)(18) together with ∝ β −4 to eliminate all temperature derivatives on the r.h.s. of Eq. (14) in terms of velocity gradients, one obtains the well-know Navier-Stokes result π µν = −2τ r β π σ µν . Here β π is a thermodynamic integral over the local equilibrium distribution, related to the relaxation time τ r and shear viscosity η by τ r = η/β π .
To obtain higher order approximations for π µν we take the co-moving time derivative of Eq. (19), where A µν ≡ ∆ µν αβ A αβ denotes the traceless symmetric projection orthogonal to u µ of the tensor A µν , and express δḟ through the Boltzmann equation (7), by rewriting it as Inserting this back into Eq. (20) one obtainṡ From this equation it is clear that the shear relaxation time τ π is equal to the Boltzmann relaxation time τ r . Now we can substitute f = f eq + δf (1) , with δf (1) from Eq. (14), on the r.h.s. of Eq. (22) to obtain the secondorder evolution equation [25] (see also [32]) To go to third-order, δf is required up to second-order in velocity gradients, where φ 1 and φ 2 are first-and second-order corrections, respectively. They are found to be [29] Please note the change of sign of several terms in the above equation compared to [29] where a "mostly minus" signature for the metric was used. We note that φ 1 and φ 2 in Eqs. (25,26) satisfy the Landau matching conditions u ν T µν = u µ and = eq [28].
Substituting f = f eq (1 + φ 1 + φ 2 ) into Eq. (22), some algebra yields the third-order evolution equation [29] The right-hand side of this equation contains three second-order and fourteen third-order terms.
An expression for the entropy four-current is derived using the kinetic theory definition for particles with Boltzmann statistics [46] S µ = − dp p µ f ln f −1 . (28) Assuming small deviations from local thermodynamic equilibrium, f = f eq (1 + φ), where φ 1, we obtain an expression for the non-equilibrium entropy four-current up to third-order in φ as where s eq = β( + P ) is the equilibrium definition of the entropy density.
where we ignore terms higher than third-order in the derivative expansion. Substituting φ 1 and φ 2 from Eqs. (25) and (26) and performing the integrations, we obtain [31] recalling that β π = 4P/5. The LRF entropy density, s ≡ −u µ S µ , is given by whereas the entropy flux in the LRF, S µ ≡ ∆ µ ν S ν , reduces to We observe that, beginning at third order in the derivative expansion, the Chapman-Enskog method leads to a non-vanishing entropy flux in the LRF.
Using the above results the evolution equations forˆ andπ take the form As in the Bjorken case some third-order contributions were brought into the formπ 2θ by using the first order (Navier-Stokes) relationπ = (4/3)τ πβπ tanh ρ. The transport coefficients in Eq. (40) are given bŷ For Gubser flow the expression for the LRF entropy densityŝ(ρ) is given in terms ofπ aŝ Similar to the Bjorken result in Milne coordinates, we find that for Gubser flow the entropy flux vanishes in de Sitter (i.e. LRF) coordinates.

V. ANISOTROPIC HYDRODYNAMICS
Anisotropic hydrodynamics makes the single particle phase space distribution function f (x, p) explicitly dependent on a spacelike four-vector l µ , which denotes the local anisotropy direction, and a momentum anisotropy parameter β l which controls the amount of deformation from the usual isotropic form. The leading-order part of the distribution function f (x, p) ≡ f a (x, p) + δf (x, p) is in general written as f a β u (−u·p), β l (l ·p) [40] such that lim β l →0 f a β u (−u · p), β l (l · p) = f eq β u (−u · p) . (43) Different from conventional hydrodynamics, the parameter β l can be arbitrarily large, enabling anisotropic hydrodynamics to handle large deviations of the system from local momentum isotropy and equilibrium.
In this work, the vector l µ is taken to point in the longitudinal η direction in the LRF, i.e. l µ = (0, 0, 0, 1) in LRF coordinates, and we consider the widely used Romatschke-Strickland (RS) [33] ansatz for the anisotropic distribution function: where Note that for this choice β u ≡ β RS and β l = β RS √ ξ. The parameter β RS is related to the inverse temperature β = 1/T through the Landau matching condition, as we shall see later.
Owing to the presence of an intrinsic directionality l µ in the system, the energy-momentum tensor T µν RS corresponding to the leading-order distribution f a only has the general decomposition in the Landau frame [40] T µν RS = RS u µ u ν + P L,RS l µ l ν − P T,RS Ξ µν .
Here Ξ µν ≡ g µν + u µ u ν − l µ l ν projects onto the space orthogonal to both u µ and l µ . The local energy density RS , longitudinal pressure P L,RS , and transverse pressure P T,RS can be expressed as moments of f RS [12]: For massless systems they are related by conformal invariance, RS = (2P T,RS + P L,RS ), and one has (β RS ) = 3P (β RS ) = 3/(π 2 β 4 RS ). The Landau matching condition RS (β RS ) = (β) yields β = β RS /R 1/4 200 . For a massless Boltzmann gas the anisotropic integrals R nrq (ξ) in Eq. (47) can be calculated analytically [12,41]: R 220 (ξ) = − 1 ξ The residual deviation δf of the distribution function generates, in principle, additional contributions to the longitudinal and transverse pressures, δP L = (−u·p) 2 δf and δP T = (l · p) 2 δf . We here use the P L matching scheme [12,41] in which the anisotropy parameter ξ(x) is chosen such that these contributions vanish exactly. This is a dynamical matching scheme similar to Landau matching which defines the local temperature T (x) in such a way that the deviation δf from local equilibrium makes no contribution to the energy density (x). With this matching scheme we can drop the subscripts RS on P L and P T .
For massless systems with Bjorken or Gubser symmetry it can be shown that there are no other dissipative contributions from δf to the energy momentum tensor [12,41]. The bulk viscous pressure vanishes by conformal symmetry, and the shear stress tensor is fully specified by the difference between the longitudinal and transverse pressures, It can thus be reduced to a single independent component for which we choose π ≡ −π η η = −τ 2 π ηη = 2 3 (P L −P T ) in the Bjorken case andπ ≡π ηη = 2 3 (P L −P T ) in the Gubser case. Evolution equations for π andπ are obtained from the Boltzmann equation following Refs. [12,41].
For Bjorken flow one finds [41] that Eqs. (34,35) are in anisotropic hydrodynamics replaced by Here we followed [41] and expressed the shear stress π through the longitudinal pressure P L via π = P −P L = 1 3 −P L . The thermodynamic integral I RS 240 over the RS distribution function is given in terms of the momentum deformation parameter ξ as I RS 240 (β, ξ) = (β) R 240 (ξ)/R 200 (ξ), with Eq. (54) agrees with Eq. (34) in Sec. IV B while the evolution equations for the shear stress π = P −P L , Eqs. (35) and (55), differ. We solve Eqs. (54,55) by using the relations (47,48) to write P L = (β)R 220 (ξ)/R 200 (ξ) and convert Eq. (55) into an evolution equation for ξ. For Gubser flow one obtains [12] the energy conservation law (39) and, instead of Eq. (40), the shear stress evolution withÎ 240 (β, ξ) =ˆ (β) R 240 (ξ)/R 200 (ξ) and the modified transport coefficientλ a = 4 3 . For the definition of the out-of-equilibrium entropy current we substitute f (x, p) = f RS (x, p)+δf in Eq. (28): Using f RS = exp(−β RS (1+ξ)w 2 /τ 2 + p 2 T ) for Bjorken flow, together with the integration measure dp = dw d 2 p T /[(2π) 3 τ (−u·p)], and applying the transformation w → w = w √ 1 + ξ for which f RS → f eq (τ, p T , w ; β RS ), the leading contribution s a can be evaluated exactly: To linear order the δf correction to the entropy density is given by To evaluate it an approximate solution of the Boltzmann equation for δf is needed. We here use the moments method [32] in the 14-moment approximation, δf ≈ δf 14 . This shows that in the P L -matching scheme only nonhydrodynamic moments of the distribution function contribute to the residual non-equilibrium entropy density δs. We leave a detailed study of such non-hydrodynamic mode contributions to entropy production to future work. For Gubser flow a similar calculation yields for the leading contribution The correction δŝ, at linear order in δf , again vanishes in the 14-moment approximation.
tems with Gubser symmetry is given in Eq. (48) of Ref. [12]. For systems with Bjorken symmetry the same expression holds without the hats. It is straightforward to see that in the P L -matching scheme this expression is zero in both cases.

VI. NUMERICAL RESULTS AND DISCUSSION
We compare the numerical results obtained from three different formalisms: Chapman-Enskog third-order viscous hydrodynamics, anisotropic hydrodynamics with P L matching, and the exact solution of the RTA Boltzmann equation. Although in principle each set of evolution equations can be solved for any initial condition, we here only show results evolving from local thermal equilibrium (with vanishing initial momentum-space deformation ξ 0 = 0 and shear stress π 0 = 0) at some initial time.

A. Bjorken flow
For our Bjorken flow results we initialize the system at longitudinal proper time τ 0 = 0.25 fm/c with initial temperature T 0 = 300 MeV. Figure 1 shows the resulting proper time evolution of the normalised shear stress π/( + P ) (panel a), the pressure anisotropy P L /P T ≡ (P − π)/(P +π/2) (panel b), the entropy density per unit rapidity and transverse area sτ (panel c), and the normalised entropy density s/s eq (panel d), for the three theories listed above and three choices of the specific shear viscosity as indicated in the figure.
For small specific viscosityη ≡ η/s = 1/4π, all three formalisms yield very similar results. Except for the normalized entropy density s/s eq in Fig. 1d, the three curves agree within line thickness. As the specific shear viscosity increases, increasing differences between the three formalisms become visible. Generally the differences remain small for the hydrodynamic moments of the distribution function, i.e. for the evolution of the normalized shear stressπ ≡ π/( +P ) and pressure anisotropy (which, because of P L /P T = (1−4π)/(1 + 2π), are basically the same quantity). Third-order Chapman-Enskog hydrodynamics performs somewhat better at late times whereas anisotropic hydrodynamics reproduces the exact solution more accurately at earlier times; the exact solution lies between these two hydrodynamic approximations.
As seen in Figs. 1c and 1d, the differences between macroscopic hydrodynamic and exact microscopic kinetic evolution are larger for the entropy. In ideal fluid dynamics with Bjorken flow, sτ is a constant of motion. The increase of sτ with time shown in Fig. 1c thus illustrates the rate of entropy production by dissipative effects in the different approaches. One sees that nonequilibrium effects on the rate of entropy production are not as well described by the hydrodynamic models as is the non-equilibrium evolution of the energymomentum tensor shown in panels a and b. Fig. 1d shows the non-equilibrium deviation of the entropy from the value expected from the first law of thermodynamics, s eq = ( +P )/T = 4P/T (where both P and T evolve according to viscous fluid dynamics). For η/s = 10 times the "minimal" KSS value of 1/(4π) [48], the entropy differs from the "equilibrium" value s eq by up to 10% for third-order Chapman-Enskog hydrodynamics, and even for anisotropic hydrodynamics (where, as discussed at the end of the previous section, only non-hydrodynamic moments of the distribution function contribute to the residual entropy) the deviation is still 5-7% over most of the evolution history. This indicates that, while the coupling of non-hydrodynamic modes into the evolution of the hydrodynamic moments of the distribution func-tion is rather weak, the same is not true for the entropy density.
A remarkable feature of the entropy evolution predicted by the exact solution of the RTA Boltzmann equation is the crossing of the three green curves in Fig. 1c corresponding to different values ofη: As the value of η increases, the initial rate of entropy production decreases, but entropy is produced over a longer time period such that its eventual saturation value increases withη. This feature, which is shared by aHydro, but not by the third-order Chapman-Enskog approach, appears counterintuitive at first sight: In first-order Navier-Stokes theory, the slope of sτ as a function of τ is proportional tō η: d(sτ )/dτ = 4ηs/3τ T . However, this argument implicitly assumes the equilibrium definition of the entropy density, s ≡ s eq = ( + P )/T , and the substantial deviation of the exact result from this lowest-order expectation illustrated in Fig. 1d (which shows that the deviation increases with increasingη) demonstrates the importance of higher order terms in the definition of the entropy density. (Note that both aHydro and the third-order Chapman-Enskog approach have trouble accounting for this non-equilibrium deviation of the entropy from the first law of thermodynamics.) Microscopically, the increasing deviation from naive Navier-Stokes expecations is related to the growth of the relaxation time with increasingη, resulting in a slower response to the expansion driving the system away from equilbrium. We have checked that the curve crossing disappears when plotting s eq τ instead of sτ ; in this case the initial slope of the curves is directly proportional toη.
Larger values of the specific shear viscosityη lead to stronger viscous heating, thereby delaying the cooling by expansion of the fireball. At a given (sufficiently late) proper time τ the more viscous fluid thus has a higher temperature than the less viscous one if both started out with the same initial temperature T 0 . The authors of Ref. [49] showed that this effect can be scaled out of the evolution plots for dimensionless ratios such as π/P or P L /P T if one plots them as a function of the dimensionless scaling variablew = τ T /(4πη/s) instead of τ . Fig. 2 shows this for the four quantities plotted in Fig. 1. The dimensionless ratios π/( +P ), P L /P T , and even the non-equilibrium entropy ratio s/s eq exhibit clear scaling behavior, converging at aroundw 1 to a universal late-time attractor given by relativistic Navier-Stokes theory. That the aHydro attractor, whose equation involves a resummation of terms in powers of inverse Reynolds number, closely matches with the exact attractor has already been demonstrated in Ref. [50], albeit with a slightly different version of aHydro that did not implement P L -matching. We note that the dimensionful quantity sτ does not scale, but the crossing of the curves seen in Fig. 1c is removed by rescaling the time evolution variable. The scaling plots shown in Fig. 2 reinforce the observation made in Fig. 1 that the hydrodynamic approximations reproduce the exact evolution of the energy momentum tensor, in particular the nor-malized shear stress and pressure anisotropy, much more accurately than that of the entropy ratio s/s eq . Eventually, however, even this latter ratio approaches a universal Navier-Stokes attractor, albeit only atw > ∼ 2, i.e. twice later than the hydrodynamic moments.

B. Gubser flow
Gubser flow is interesting because of its very strong transverse expansion which asymptotically (i.e. for very large de Sitter times) drives the system arbitrarily far away from local thermal equilibrium, into a state of freestreaming [10,15]. This is in contrast to Bjorken flow where there is no transverse flow and the longitudinal expansion rate decreases for late longitudinal proper times,  For each theory three sets of curves are shown, corresponding to three different values for the specific shear viscosity, 4πη/s = 1 (solid), 3 (dashed), and 10 (dash-dotted). Thermal equilibrium initial conditions (π = 0) with initial temperatureT0 = 0.002 were implemented at ρ0 = −10.
allowing the system to settle into a state of approximate local thermal equilibrium. The dramatic transverse expansion encoded in Gubser flow thus provides a testbed for the performance of hydrodynamic approximations in situations very far from equilibrium.
For Gubser flow we initialize the system in equilibrium (i.e. withπ 0 = ξ 0 = 0) at de Sitter time ρ 0 = −10 with initial normalized temperatureT = 0.002. through a transient state of approximate local momentum isotropy (π = ξ = 0) before again being driven away from it by increasingly strong transverse expansion (resulting in positive pressure anisotropyπ ∼P L −P T > 0), eventually leading to free-streaming withπ/(ˆ +P ) → 0.5 at late de Sitter times. Fig. 3 shows that at late de Sitter times anisotropic hydrodynamics slightly overpredicts the temperature corresponding (by Landau matching) to the energy density of the exact solution of the RTA Boltzmann equation, by a constant factor. For the third-order Chapman-Enskog approach, the asymptotic temperature is seen to keep falling further and further below that of the exact solution, indicating (as the following figures will show more clearly) that this hydrodynamic model does not correctly approach the asymptotic free-streaming state and underpredicts the shear stress and viscous heating at late de Sitter times. For both hydrodynamic approximations the asymptotic deviation from the exact solution increases with the specific shear viscosityη.
In Fig. 4 we show the de Sitter time evolution of the Gubser analogues of the quantities plotted in Fig. 1 above for Bjorken flow. As already reported in [12], anisotropic hydrodynamics with P L matching provides a very accurate approximation to the exact solution of the RTA Boltzmann equation for the evolution of the shear stress and pressure anisotropy (panels a and b). In particular, it approaches the correct free-streaming limit at large de Sitter times. This approach is faster (in ρ) for larger specific shear viscosity η/s. However, as was the case for Bjorken flow, the ability of aHydro to describe the evolution of the entropy content of the system (panel c) and of the non-equilibrium correction to the first law of thermodynamics (shown in panel d) is much more limited.
Especially at late de Sitter times, the aHydro curves appear to move farther and farther away from the exact solution.
For the third-order Chapman-Enskog approach, large deviations from the exact solution at late de Sitter times are even observed for the hydrodynamic moments shown in Figs. 4a and b: Instead of saturating at the freestreaming limitπ ≡π/(ˆ +P ) = 0.5, the normalized shear stress in Fig. 4a saturates at 0.4. As a result, the pressure anisotropyP L /P T = (1+4π)/(1−2π) shown in Fig. 4b saturates at large de Sitter times in the thirdorder Chapman-Enskog approach instead of continuing to grow as dictated by the exact solution of the RTA Boltzmann equation and is correctly reproduced by aHydro with P L -matching. This failure is similar to the one observed in DNMR theory [32] (which is a second-order viscous hydrodynamic approach based on an expansion around a locally isotropic momentum distribution function) except that in DNMR theoryπ saturates at a value > 0.5, corresponding to negative transverse pressure and instability against cavitation [12].
As far as the de Sitter time evolution of the entropy content of the system (Fig. 4c) and of the nonequilibrium correction to the first law of thermodynamics (Fig. 4d) are concerned, the discrepancies between third-order Chapman-Enskog hydrodynamics and the exact solution of the Boltzmann equation are even larger than those observed for anisotropic hydrodynamics. Although all approaches correctly predict that the entropy densityŝ cosh 2 ρ increases as dictated by the second law of thermodynamics (Fig. 4c), the rate of increase is overpredicted by the hydrodynamic models at late de Sitter times. The rate of viscous entropy production is controlled by the normalized shear stressπ shown in Fig. 4a; near ρ = 0 it is small in all three approaches becauseπ passes through zero. Finally, Fig. 4d shows that thirdorder Chapman-Enskog hydrodynamics predicts a saturation of the ratio s/s eq at large de Sitter times whereas the exact solution shows that this ratio should continue to decrease as ρ keeps increasing. aHydro reproduces this continued decrease, but at an incorrect rate.
It is worth noting that for both hydrodynamic approximations studied here, the ratio s/s eq shown in Fig. 4d passes through 1 near ρ = 0 whereπ passes through zero. This is not the case for the exact solution which shows non-vanishing deviations of this ratio from unity (whose magnitude increases with η/s) even whenπ = 0. As similar observation was made before in Ref. [20], it shows that the exact solution of the Boltzmann equation includes contributions to the non-equilibrium entropy from nonhydrodynamic moments [20] that are not captured by the hydrodynamic approximations studied here.
We close this section by replotting Fig. 4 as a function of the scaling variablew = (4πη/s)(2 tanh ρ)/T [51] in Fig. 5. 4 Our findings are consistent with the detailed study of the Gubser flow fixed point presented in Ref. [51]. As for the case of Bjorken flow, one observes convergence of the curves describing the evolution of the normalized shear stress (Fig. 5a) and pressure anisotropy (Fig. 5b) for different specific shear viscosities to a common attractor at large values ofw. 5 In this case, however, the attractor for the normalized shear stressπ differs for third-order Chapman-Enskog hydrodynamics from the shared "free-streaming attractor" for aHydro and the exact solution of the RTA Boltzmann equation. This reflects the above observation that third-order Chapman-Enskog hydrodynamics does not approach the correct free-streaming limit at large de Sitter times. Fig. 5b additionally shows that the rate at which the trajectories for aHydro and the RTA Boltzmann equation approach the asymptotic valueπ = 0.5 is slightly different for the two theories, but insensitive to the value of η/s in each case.
In contrast to the dimensionless ratios shown in panels a, b, and d, the evolution of the dimensionful entropy density shown in Fig. 5c exhibits no clear scaling behavior. For the non-equilibrium entropy ratio s/s eq in Fig. 5d one observes different attractors for all three dynamical approaches: whereas in each case the curves corresponding to different specific shear viscosity converge at largew, the attractors they converge to are very different for the exact solution, aHydro and third-order Chapman-Enskog. The difference between the aHydro and exact attractors is smaller than between third-order Chapman-Enskog and the exact result, but still large. Clearly, the hydrodynamic approximations are having difficulties reproducing the non-equilibrium contributions to the entropy density at largew, i.e. deep in the free-streaming region of the exact solution.

VII. CONCLUSIONS
In this work, we have considered two different formalisms for deriving macroscopic descriptions of the non-equilibrium dynamics of a system, namely, dissi-pative hydrodynamics using the Chapman-Enskog iterative scheme to third order and anisotropic hydrodynamics (aHydro) with P L -matching. The performance of these different hydrodynamic schemes was tested by comparing their predictions with the exact solution of the RTA Boltzmann equation in two situations where such an exact solution is available, namely for the Bjorken and Gubser flows. Both situations are effectively onedimensional such that the energy-momentum tensor can be characterized by just two hydrodynamic moments of the microscopic distribution function, the energy density (or, equivalently, the temperature) and a single shear stress component. The shear stress also defines the phenomenologically important longitudinal-transverse pressure anisotropy P L /P T . Bjorken and Gubser flows describe two extreme situations that bracket realistic situations: while both share boost-invariant longitudinal expansion, Bjorken flow lacks any transverse expansion (and correspondingly allows the system to approach a state of local thermal equilibrium at late times) whereas Gubser flow features very strong radial expansion in the transverse directions which at late times drives the system completely away from local equilibrium into an asymptotic state of free-streaming. Both flows start out with strong longitudinal expansion in which dissipative effects deform the local rest frame momentum distribution by making it narrower in the longitudinal momentum p η than in transverse momentum p T (such that P L −P T < 0), but for Bjorken flow the local momentum distribution becomes asymptotically isotropic whereas for Gubser flow it eventually becomes narrower in p T than p η (leading to P L −P T > 0). The two flows thus present a testbed for macroscopic hydrodynamic approximations of the microscopic dynamics under very different conditions of anisotropic expansion, with opposite signs of the pressure anisotropy P L − P T at late times.
In addition to the evolution of the abovementioned hydrodynamic moments (whose dynamics has been studied before) we also explored here the evolution of the entropy density of the system (which, in practical situations such as relativistic heavy-ion collisions, controls the multiplicity of finally emitted hadrons). Our interest in the entropy arises from previous observations [20] that suggested that the entropy evolution is more strongly influenced by dynamical couplings to non-hydrodynamic moments of the distribution function, and we wanted to know how well these couplings can be captured in macroscopic hydrodynamic treatments.
In all cases (i.e. for both anisotropic flow patterns and for all the observables studied) we found that anisotropic hydrodynamics with P L -matching provides a more accurate approximation to the exact evolution obtained from the exact solution of the Boltzmann equation than does the dissipative hydrodynamics derived from a thirdorder Chapman-Enskog expansion of the distribution function. The latter is found to consistently underpredict the deviation from local equilibrium even when terms up to third-order in velocity gradients are kept in the expansion of the shear stress tensor. As a consequence, the results of third-order Chapman-Enskog hydrodynamics deviate substantially from the exact solution whenever momentum-space anisotropies become non-perturbatively large. This feature is most apparent for Gubser flow, both during early times when longitudinal expansion dominates the pressure anisotropy and at late times when the strong transverse expansion drives the pressure anisotropy and pushes the system towards free-streaming.
Our findings can be understood most intuitively when plotting them against a dimensionless time variablew (defined in the text) that is scaled by the microscopic relaxation time (which increases with increasing specific shear viscosity η/s). For Bjorken flow one finds that the exact solution and the two hydrodynamic approximation schemes studied in this paper share a common attractor to which all solutions converge at late times, irrespective of initial conditions. For Gubser flow, aHydro shares a common attractor with the exact solution for the normalized shear stress and pressure anisotropy, whereas these quantities approach a different attractor for third-order Chapman-Enskog hydrodynamics. For the non-equilibrium entropy, the asymptotic evolution in Gubser flow is controlled by three different attractors for the exact solution and the two hydrodynamic approximation schemes, with the differences between the exact and aHydro attractors being smaller than between the exact solution and third-order Chapman-Enskog hydrodynamics.