Exact solutions and attractors of higher-order viscous fluid dynamics for Bjorken flow

TL;DR

This paper analyzes higher-order relativistic viscous hydrodynamics theories for Bjorken flow, demonstrating that third-order theories better approximate the exact kinetic theory attractor and characterizing the approach to equilibrium.

Contribution

It provides analytical solutions and a detailed comparison of second- and third-order hydrodynamics theories with exact kinetic theory for Bjorken flow.

Findings

Third-order theory outperforms MIS and DNMR in approximating the exact attractor.

Exponential approach to the attractor at small Knudsen numbers.

Power-law decay of deviations persists even at large Knudsen numbers.

Abstract

We consider causal higher order theories of relativistic viscous hydrodynamics in the limit of one-dimensional boost-invariant expansion and study the associated dynamical attractor. We obtain evolution equations for the inverse Reynolds number as a function of Knudsen number. The solutions of these equations exhibit attractor behavior which we analyze in terms of Lyapunov exponents using several different techniques. We compare the attractors of the second-order M\"uller-Israel-Stewart (MIS), transient DNMR, and third-order theories with the exact solution of the Boltzmann equation in the relaxation-time approximation. It is shown that for Bjorken flow the third-order theory provides a better approximation to the exact kinetic theory attractor than both MIS and DNMR theories. For three different choices of the time dependence of the shear relaxation rate we find analytical solutions for the energy density and shear stress and use these to study the attractors analytically. By studying these analytical solutions at both small and large Knudsen numbers we characterize and uniquely determine the attractors and Lyapunov exponents. While for small Knudsen numbers the approach to the attractor is exponential, strong power-law decay of deviations from the attractor and rapid loss of initial state memory is found even for large Knudsen numbers. Implications for the applicability of hydrodynamics in far-off-equilibrium situations are discussed.

Tables

Table 1. Table 1: Coefficients for the causal viscous hydrodynamic evolution of the shear stress in Bjorken flow for the three theories studied in this work.

	${\beta }_{\pi }$	$a$	$\lambda$	$\chi$	$\gamma$
MIS	$4P/5$	4/15	0	0	4/3
DNMR	$4P/5$	4/15	10/21	0	4/3
Third-order	$\mathrm{\hspace{0.17em}4}P/5$	4/15	10/21	72/245	412/147

Table 2. Table 2: Arguments and parameters of Eq. ( 49 ) for the analytic approximations studied in Secs.  V  A, B, and C, respectively.

$T\left(\tau \right)$	$w$	$\Lambda$	$k$	$m$
$|$ const. $|$	$\overline{\tau }$	$-1$	$-\frac{1}{2}\left(\lambda +1\right)$	$\frac{1}{2}\sqrt{4a\gamma +{\lambda }^2}$
$|$ ideal $|$	$\frac{3}{2}\overline{\tau }$	$-\frac{3}{2}$	$-\frac{3\lambda +2}{4}$	$\frac{3}{4}\sqrt{4a\gamma +{\lambda }^2}$
$|$ NS $|$	$\frac{3}{2}\left(\overline{\tau }+\frac{a}{2}\right)$	$-\frac{3}{2}$	$-\frac{6\lambda +4-3a}{8}$	$\frac{3}{8}\sqrt{16a\gamma +a^2-4a\lambda +4{\lambda }^2}$

Equations

$$T^{\mu \nu }=ϵu^{\mu }{u}^{\nu }-P{\Delta }^{\mu \nu }+{\pi }^{\mu \nu },$$

$$\dot{ϵ}+(ϵ+P)\theta -{\pi }^{\mu \nu }{\sigma }_{\mu \nu }$$

$$(ϵ+P){\dot{u}}^{\alpha }-{\nabla }^{\alpha }P+{\Delta }_{\nu }^{\alpha }{\partial }_{\mu }{\pi }^{\mu \nu }$$

$${\tau }_{\pi }{\dot{\pi }}^{⟨\mu \nu ⟩}+{\pi }^{\mu \nu }=2\eta {\sigma }^{\mu \nu }-\frac{4}{3}{\tau }_{\pi }{\pi }^{\mu \nu }\theta .$$

$${\dot{\pi }}^{⟨\mu \nu ⟩}+\frac{{\pi }^{\mu \nu }}{{\tau }_{\pi }}=2{\beta }_{\pi }{\sigma }^{\mu \nu }+2{\pi }_{\gamma }^{⟨\mu }{\omega }^{\nu ⟩\gamma }-\frac{10}{7}{\pi }_{\gamma }^{⟨\mu }{\sigma }^{\nu ⟩\gamma }-\frac{4}{3}{\pi }^{\mu \nu }\theta .$$

$${\dot{\pi }}^{⟨\mu \nu ⟩}=$$

$$-\frac{4}{3}{\pi }^{\mu \nu }\theta +\frac{25}{7{\beta }_{\pi }}{\pi }^{\rho ⟨\mu }{\omega }^{\nu ⟩\gamma }{\pi }_{\rho \gamma }-\frac{1}{3{\beta }_{\pi }}{\pi }_{\gamma }^{⟨\mu }{\pi }^{\nu ⟩\gamma }\theta$$

$$-\frac{38}{245{\beta }_{\pi }}{\pi }^{\mu \nu }{\pi }^{\rho \gamma }{\sigma }_{\rho \gamma }-\frac{22}{49{\beta }_{\pi }}{\pi }^{\rho ⟨\mu }{\pi }^{\nu ⟩\gamma }{\sigma }_{\rho \gamma }$$

$$-\frac{24}{35}{\nabla }^{⟨\mu }\left({\pi }^{\nu ⟩\gamma }{\dot{u}}_{\gamma }{\tau }_{\pi }\right)+\frac{4}{35}{\nabla }^{⟨\mu }\left({\tau }_{\pi }{\nabla }_{\gamma }{\pi }^{\nu ⟩\gamma }\right)$$

$$-\frac{2}{7}{\nabla }_{\gamma }\left({\tau }_{\pi }{\nabla }^{⟨\mu }{\pi }^{\nu ⟩\gamma }\right)+\frac{12}{7}{\nabla }_{\gamma }\left({\tau }_{\pi }{\dot{u}}^{⟨\mu }{\pi }^{\nu ⟩\gamma }\right)$$

$$-\frac{1}{7}{\nabla }_{\gamma }\left({\tau }_{\pi }{\nabla }^{\gamma }{\pi }^{⟨\mu \nu ⟩}\right)+\frac{6}{7}{\nabla }_{\gamma }\left({\tau }_{\pi }{\dot{u}}^{\gamma }{\pi }^{⟨\mu \nu ⟩}\right)$$

$$-\frac{2}{7}{\tau }_{\pi }{\omega }^{\rho ⟨\mu }{\omega }^{\nu ⟩\gamma }{\pi }_{\rho \gamma }-\frac{2}{7}{\tau }_{\pi }{\pi }^{\rho ⟨\mu }{\omega }^{\nu ⟩\gamma }{\omega }_{\rho \gamma }$$

$$-\frac{10}{63}{\tau }_{\pi }{\pi }^{\mu \nu }{\theta }^2+\frac{26}{21}{\tau }_{\pi }{\pi }_{\gamma }^{⟨\mu }{\omega }^{\nu ⟩\gamma }\theta .$$

$${\dot{\pi }}^{⟨\eta \eta ⟩}=-\frac{1}{{\tau }^2}\frac{d\pi }{d\tau },{\pi }_{\gamma }^{⟨\eta }{\sigma }^{\eta ⟩\gamma }=-\frac{\pi }{3{\tau }^3},{\pi }_{\gamma }^{⟨\eta }{\pi }^{\eta ⟩\gamma }=-\frac{{\pi }^2}{2{\tau }^2},$$

$${\pi }^{\rho ⟨\eta }{\pi }^{\eta ⟩\gamma }{\sigma }_{\rho \gamma }=-\frac{{\pi }^2}{2{\tau }^3},{\nabla }^{⟨\eta }{\nabla }_{\gamma }{\pi }^{\eta ⟩\gamma }=\frac{2\pi }{3{\tau }^4},$$

$${\nabla }_{\gamma }{\nabla }^{⟨\eta }{\pi }^{\eta ⟩\gamma }=\frac{4\pi }{3{\tau }^4},{\nabla }^2{\pi }^{⟨\eta \eta ⟩}=\frac{4\pi }{3{\tau }^4},{\pi }^{\rho \gamma }{\sigma }_{\rho \gamma }=\frac{\pi }{\tau },$$

$$\frac{dϵ}{d\tau }$$

$$\frac{d\pi }{d\tau }$$

$$\bar{\pi }=3\left(\frac{\tau }{\bar{\tau }}\right)\frac{d\bar{\tau }}{d\tau }-2.$$

$$\left(\frac{\bar{\pi }+2}{3}\right)\frac{d\bar{\pi }}{d\bar{\tau }}=-\bar{\pi }+\frac{1}{\bar{\tau }}\left(a-\lambda \bar{\pi }-\gamma {\bar{\pi }}^2\right).$$

$$|s(t)|\approx e^{\Lambda t}|s_0|,$$

$$\bar{\pi }(\bar{\tau })=\sum _{n=0}^{\infty }\frac{c_n}{{\bar{\tau }}^n},$$

$$c_n=a{\delta }_{n,1}$$

$$\displaystyle+\sum_{m=1}^{n}\Bigl{(}\frac{m{-}1}{3}-\gamma\Bigr{)}\,c_{n-m}\,c_{m-1},$$

$$\frac{2}{3}\frac{d\bar{\pi }}{d\bar{\tau }}=-\bar{\pi }+\frac{a}{\bar{\tau }}$$

$$\bar{\pi }=\alpha e^{-\frac{3}{2}\bar{\tau }}+\frac{3a}{2}{e}^{-\frac{3}{2}\bar{\tau }}\mathrm{Ei}\left[\frac{3\bar{\tau }}{2}\right]$$

$$\frac{\partial \bar{\pi }}{\partial \alpha }\sim \exp \left(-\frac{3}{2}\bar{\tau }\right).$$

$$\bar{\pi }(\bar{\tau })={\bar{\pi }}_0{e}^{-\frac{3}{2}\left(\bar{\tau }-{\bar{\tau }}_0\right)}+\frac{3a}{2}{e}^{-\frac{3}{2}\bar{\tau }}{\int }_{{\bar{\tau }}_0}^{\bar{\tau }}\frac{e^{\frac{3}{2}{\bar{\tau }}^{\prime }}}{{\bar{\tau }}^{\prime }}d{\bar{\tau }}^{\prime }$$

$${\int }_{{\bar{\tau }}_0}^{\bar{\tau }}\frac{e^{\frac{3}{2}{\bar{\tau }}^{\prime }}}{{\bar{\tau }}^{\prime }}d{\bar{\tau }}^{\prime }$$

$$+\frac{3}{2}m!{\int }_{{\bar{\tau }}_0}^{\bar{\tau }}d{\bar{\tau }}^{\prime }\frac{e^{\frac{3}{2}{\bar{\tau }}^{\prime }}}{(3{\bar{\tau }}^{\prime }/2{)}^{m+1}},$$

$$\bar{\pi }=$$

$$+a{\left(\frac{3}{2}\right)}^{2}m!e^{-3\bar{\tau }/2}{\int }_{{\bar{\tau }}_0}^{\bar{\tau }}d{\bar{\tau }}^{\prime }\frac{e^{\frac{3}{2}{\bar{\tau }}^{\prime }}}{(3{\bar{\tau }}^{\prime }/2{)}^{m+1}},$$

$$\bar{\pi }=$$

$$+\frac{3a}{2}{e}^{-\frac{3}{2}\bar{\tau }}\sum _{n=1}^{\infty }\frac{(3/2{)}^n}{n!n}\left({\bar{\tau }}^n-{\bar{\tau }}_0^n\right).$$

$$\bar{\pi }={\bar{\pi }}_0{e}^{-\frac{3}{2}\left(\bar{\tau }-{\bar{\tau }}_0\right)}+\frac{3a}{2}{e}^{-\frac{3}{2}\left[\bar{\tau }-\beta (\bar{\tau };{\bar{\tau }}_0)\right]}\frac{\bar{\tau }-{\bar{\tau }}_0}{\beta (\bar{\tau };{\bar{\tau }}_0)}.$$

$$\frac{d(\delta \bar{\pi })}{d\bar{\tau }}=-\frac{3}{2}\delta \bar{\pi }\left[1+\frac{2\lambda -a}{2\bar{\tau }}\right],$$

$$\delta \bar{\pi }\sim {\bar{\tau }}^{\frac{3}{4}(a-2\lambda )}\exp \left(-\frac{3}{2}\bar{\tau }\right)+\mathcal{O}\left(\frac{1}{{\bar{\tau }}^2}\right).$$

$$\mathcal{B}\left(\sum _n{a}_n\right)\equiv {\int }_0^{\infty }du{e}^{-u}\sum _n\frac{a_n}{n!}{u}^n.$$

Full text

I Introduction

Hydrodynamics is an effective macroscopic theory describing the large scale and slowly varying dynamical modes of multiparticle systems. The conventional formulation of hydrodynamic equations proceeds by assuming a separation of macroscopic and microscopic length and time scales such that gradients of local-equilibrium quantities, normalized to appropriate powers of the system’s energy density, are small and one can perform an order-by-order gradient expansion around local thermodynamic equilibrium Florkowski:2017olj . Therefore the remarkable success of relativistic hydrodynamics in describing the quark-gluon plasma (QGP) formed in ultra-relativistic heavy-ion collisions initially led to the belief that these collisions create a nearly thermalized medium close to local thermal equilibrium Heinz:2001xi (see also the reviews Heinz:2013th ; Jaiswal:2016hex ; Florkowski:2017olj ; Romatschke:2017ejr ). On the other hand, with the advent of numerical dissipative relativistic fluid dynamics Heinz:2005bw ; Romatschke:2007mq ; Song:2007fn ; Dusling:2007gi ; Song:2007ux ; Luzum:2008cw ; Song:2008si it became clear that the dynamical evolution of heavy-ion collisions is affected by persistent large dissipative corrections. In spite of this numerical evidence for large deviations from local thermal equilibrium, second-order dissipative relativistic fluid dynamics (especially when coupled with a hadronic cascade to describe the final dilute decoupling stage Song:2010aq ; Heinz:2011kt ) met with impressive phenomenological success and predictive power in the description of heavy-ion collision experiments Romatschke:2007mq ; Luzum:2008cw ; Song:2008si ; Song:2008hj ; Schenke:2010nt ; Song:2010mg ; Song:2011hk ; Song:2011qa ; Shen:2011eg ; Schenke:2011bn ; Gale:2012rq , and was even successfully applied to explain flow data in small collision systems formed in proton-proton (p-p) and proton-lead (p-Pb) collisions Bozek:2015swa ; Bozek:2016jhf . This “unreasonable effectiveness” of hydrodynamics as a dynamical description of high-energy hadronic collisions in situations that are even very far away from local thermal equilibrium has generated much recent interest in the very foundations of fluid dynamics Heller:2013fn ; Heller:2015dha ; Aniceto:2015mto ; Basar:2015ava ; Florkowski:2016zsi ; Heller:2016rtz ; Romatschke:2016hle ; Behtash:2017wqg ; Blaizot:2017lht ; Blaizot:2017ucy ; Denicol:2017lxn ; Romatschke:2017acs ; Romatschke:2017vte ; Strickland:2017kux ; Spalinski:2017mel ; Aniceto:2018uik ; Behtash:2018moe ; Chattopadhyay:2018apf ; Denicol:2018pak ; Heller:2018qvh ; Strickland:2018ayk ; Tinti:2018qfb ; Behtash:2019txb ; Blaizot:2019scw ; Grozdanov:2019uhi ; Grozdanov:2019kge ; Heinz:2019dbd ; Strickland:2019hff , culminating in the formulation of a new “far-from-local-equilibrium fluid dynamics” paradigm Romatschke:2017ejr ; Romatschke:2017vte . The present work is a contribution to this ongoing discussion, adding new analytic results for the heavily studied simple case of (0+1)-dimensional Bjorken expansion of a transversally homogeneous system with longitudinal boost-invariance Bjorken:1982qr .

The simplest relativistic dissipative theory, relativistic Navier-Stokes theory Eckart:1940zz ; Landau:1987 , imposes instantaneous constitutive relations between the dissipative flows and their generating forces, expressed through first-order gradients of equilibrium quantities. This approach was found to be plagued by acausality and intrinsic instability Hiscock:1983zz ; Hiscock:1985zz . The phenomenological second-order theory developed by Müller, Israel and Stewart (MIS) Muller:1967zza ; Israel:1976tn ; Israel:1979wp cures these problems by introducing a relaxation type equation for the dissipative flows and thus turning them into independent dynamical degrees of freedom of the system whose evolution is controlled by the competition between macroscopic expansion (driving the system away from local equilibrium) and microscopic scattering (driving it back towards local equilibrium). As discussed in Romatschke:2017ejr , even the minimal causal conformal theory given by MIS introduces new modes called non-hydrodynamic modes that were absent in Navier-Stokes theory. These non-hydrodynamic modes are now known to play an important role in the approach to the regime of applicability of hydrodynamics, also known as the “hydrodynamization” process Heller:2014wfa ; Heller:2015dha ; Kurkela:2015qoa ; Heller:2016rtz ; Heller:2016gbp . For conformal systems undergoing longitudinal Bjorken expansion, it was shown explicitly that hydrodynamization occurs at microscopic thermalization time scales which, in strongly and anisotropically expanding systems, are much shorter than both the time scales for local isotropization and thermalization Romatschke:2017vte . Similar conclusions were obtained in recent studies Behtash:2017wqg ; Chattopadhyay:2018apf of systems undergoing Gubser flow Gubser:2010ze where the fireball expands not only boost-invariantly in longitudinal, but also azimuthally symmetrically in transverse direction. These findings lead to the conclusion that the criterium of proximity to local thermodynamic equilibrium for hydrodynamics to be valid is (at least in these carefully analytically studied simplified situations) too strict and should be replaced by a different condition stating that contributions from the non-hydrodynamic modes can be neglected Romatschke:2016hle . In the present study, we will focus on yet another interesting feature that appears in a causal theory of relativistic dissipative hydrodynamics, “the hydrodynamic attractor” Heller:2015dha ; Romatschke:2017acs ; Denicol:2018pak ; Behtash:2017wqg ; Strickland:2017kux ; Chattopadhyay:2018apf , and its role in the hydrodynamization process.

Attractor behavior was first identified by considering the hydrodynamic formulation as a gradient series expansion. Recently, this gradient series was shown, for several highly symmetric flow configurations that are amenable to analytic treatment, to have zero radius of convergence. This suggests that hydrodynamic theory cannot be systematically improved by taking into account higher order-terms in the gradient series Heller:2015dha . The gradient expansion generates an asymptotic series which exhibits initial signs of convergence for a few terms before eventually diverging Heller:2015dha ; Basar:2015ava . The initial appearance of convergence may explain the observed remarkable phenomenological success of hydrodynamic formulations based on truncations of the gradient expansion at second or third order but, due to the ultimate divergence of the series, the theory cannot be improved beyond a certain order by keeping additional terms. Fortunately, the (diverging) gradient expansion series can be Borel resummed, giving rise to a unique hydrodynamic attractor solution (hydrodynamic mode), which is well defined even for large gradients, and a series of rapidly decaying non-hydrodynamic modes that describe the approach towards this attractor from arbitrary initial conditions Heller:2013fn ; Heller:2015dha . This suggests that hydrodynamics displays a novel type of universality even far from local equilibrium which is independent of the initial state of the system, indicating the existence of a new, far-from-local-equilibrium hydrodynamic theory Romatschke:2017vte . Recently, it was discovered Strickland:2019jut that this phenomenon extends even beyond hydrodynamics: in kinetic theory also the evolution of non-hydrodynamic higher-order momentum moments of the distribution function is controlled by attractors.

The theory of dynamical attractors can be complex, and hence it is instructive to have examples were the attractor is quantitatively understood using analytical methods. Perhaps the most useful quantities associated with an attractor are Lyapunov exponents which characterize the rate of separation of infinitesimally close trajectories in the phase-space evolution of dynamical systems Goldstein . While negative Lyapunov exponents are associated with dissipative systems and indicate the existence of attractors, positive values are usually associated with chaotic systems. For a conservative system one obtains vanishing Lyapunov exponents. In this article, we employ the theory of Lyapunov exponents to study the attractors for MIS theory and two other, improved versions of causal relativistic dissipative hydrodynamics (see Refs. Behtash:2017wqg ; Behtash:2019txb for earlier related work).

Motivated by its application to QGP evolution, several improvements over MIS theory were proposed during the last decade Martinez:2010sc ; Florkowski:2010cf ; Denicol:2012cn ; Strickland:2014pga ; Bazow:2013ifa ; Jaiswal:2013vta ; Jaiswal:2014isa ; Tinti:2015xwa ; Molnar:2016vvu ; Molnar:2016gwq ; McNelis:2018jho ; Nopoush:2019vqc . For Bjorken and Gubser flows, some of these theories lead to very good agreement with the exact solution of the Boltzmann equation even in situations where the deviations from thermal equilibrium are large Florkowski:2013lya ; Martinez:2017ibh ; Jaiswal:2013vta ; Chattopadhyay:2018apf . In this article, we consider the Denicol-Niemi-Molnar-Rischke (DNMR) theory Denicol:2012cn which is an improved version of second-order MIS theory, as well as a third-order hydrodynamic theory derived form relativistic kinetic theory by going to third order in the Chapman-Enskog expansion Jaiswal:2013vta . To understand the emergent attractor behavior we go beyond previous numerical studies of these theories by finding analytical solutions of the corresponding hydrodynamic evolution equations Chattopadhyay:2018pwe . Analytical solutions of higher-order dissipative hydrodynamics were found previously for a few very special cases Hatta:2014gqa ; Denicol:2017lxn . We here expand this portfolio and use the new analytical solutions to study their dependence on initial conditions, their attractors, and their late-time behavior.

The rest of this article is structured as follows. In Sec. II we briefly review causal theories of second-order (MIS and DNMR) and third-order relativistic viscous hydrodynamics. In Sec. III we simplify these theories for Bjorken flow and obtain a generic evolution equation for the inverse Reynolds number as a function of the inverse Knudsen number, which can be adapted to all three hydrodynamic theories by adjusting a set of two parameters. This ordinary differential equation is decoupled from the hydrodynamic evolution of the energy density (but feeds back into it), and thus it can be solved independently. In Sec. IV we demonstrate numerically that the solutions of these differential equations exhibit attractor-like behavior which we analyze in terms of Lyapunov exponents, using several different approaches. We also compare the exact numerical attractors of the second-order MIS and DNMR theories, as well as our third-order theory, with the exact solution of the relativistic Boltzmann equation in the relaxation-time approximation (RTA). We show that the third-order theory provides a better approximation to the exact kinetic theory attractor than both the MIS and DNMR attractors. In Sec. V we provide further clarification of the numerical behavior found in the preceding sections by working out explicit analytical solutions for the energy density and shear stress as a function of the inverse Knudsen number, for three different forms of the shear relaxation time. These analytical solutions allow us to also study the attractors analytically. In Sec. VI, finally, we explore the universal behavior of these solutions at both small and large Knudsen numbers and use the results to characterize and uniquely determine the hydrodynamic attractor and its associated Lyapunov exponents in a new way. Our results are summarized and some conclusions offered in Sec. VII.

II Relativistic dissipative hydrodynamics

In this section, we will briefly review the equations for relativistic dissipative hydrodynamics. We consider a conformal system, which corresponds in kinetic theory to a system of massless particles. The energy-momentum tensor for such a system, in the Landau frame, has the form

[TABLE]

where $\epsilon$ and $P$ are the local energy density and pressure. Conformal symmetry implies an equation of state (EoS) $\epsilon=3P$ and zero bulk viscous pressure, $\Pi=0$. We also define $\Delta^{\mu\nu}\equiv g^{\mu\nu}{-}u^{\mu}u^{\nu}$ which serves as a projection operator to the space orthogonal to $u^{\mu}$ (i.e. onto the spatial directions in the local rest frame (LRF)). The shear stress tensor, $\pi^{\mu\nu}$, is traceless and orthogonal to $u^{\mu}$. The metric convention used here is $g^{\mu\nu}={\rm diag}(+\,-\,-\,-)$.

Evolution equations for $\epsilon$ and $u^{\mu}$ are obtained from energy-momentum conservation, $D_{\mu}T^{\mu\nu}=0$:

[TABLE]

Here we use the notations $D_{\mu}$ for the covariant derivative, $\dot{A}\equiv u^{\mu}D_{\mu}A$ for the co-moving time derivative, $\nabla^{\alpha}\equiv\Delta^{\mu\alpha}D_{\mu}$ for space-like derivative, $\theta\equiv D_{\mu}u^{\mu}$ for the expansion scalar, and $\sigma_{\mu\nu}\equiv\frac{1}{2}(\nabla_{\mu}u_{\nu}{+}\nabla_{\nu}u_{\mu})-\frac{1}{3}\theta\Delta_{\mu\nu}$ for the velocity shear tensor.

To close Eqs. (2),(3) we need additional equations for the shear stress $\pi^{\mu\nu}$. The simplest form of $\pi^{\mu\nu}$ is the Navier-Stokes form, which is first order in velocity gradients, $\pi^{\mu\nu}_{\mathrm{NS}}=2\eta\sigma^{\mu\nu}$, where $\eta$ is the shear viscosity coefficient. As already mentioned in the Introduction, relativistic Navier-Stokes theory violates causality and is unstable. The simplest way to restore causality is by introducing a relaxation-type equation for $\pi^{\mu\nu}$. This prescription, also known as the “Maxwell-Cattaneo law”, requires that the dissipative forces relax to their Navier-Stokes values in some finite relaxation time, i.e., $\tau_{\pi}\dot{\pi}^{\langle\mu\nu\rangle}+\pi^{\mu\nu}=2\eta\sigma^{\mu\nu}$,111Angular brackets around pairs of Lorentz indices indicate projection of the tensor onto its traceless and locally spatial part, e.g., $\dot{\pi}^{\langle\mu\nu\rangle}=\Delta^{\mu\nu}_{\alpha\beta}\dot{\pi}^{\alpha\beta}$, where $\Delta^{\mu\nu}_{\alpha\beta}=\frac{1}{2}(\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}{+}\Delta^{\mu}_{\beta}\Delta^{\nu}_{\alpha})-\frac{1}{3}\Delta^{\mu\nu}\Delta_{\alpha\beta}$. where $\tau_{\pi}$ is the shear relaxation time. For conformally symmetric systems one more term must be added in the evolution of shear stress Baier:2006um ; Baier:2007ix :

[TABLE]

This equation is a close variant Baier:2006um of the one first derived by Müller, Israel and Stewart Israel:1979wp ; Muller:1967zza ; Israel:1976tn , and we will therefore refer to it as the “MIS” theory. Its derivation was based on an analysis of the entropy current and the second law of thermodynamics, without recourse to a specific theory for the underlying microscopic dynamics.222It should be noted that Eq. (4) does not include all possible second-order terms Baier:2007ix , and that (unlike the theories discussed below) its transport coefficients are not matched to an underlying kinetic theory.

A systematic derivation of second-order (“transient”) relativistic fluid dynamics from relativistic kinetic theory, using an expansion of the dissipative flows in momentum-moments of the distribution function, was performed in Denicol:2012cn . For conformal systems and an RTA collision term, the result obtained in the 14-moment approximation differs from Eq. (4) by two additional terms that are of second order in gradients:

[TABLE]

Here $\beta_{\pi}\equiv\eta/\tau_{\pi}=4P/5$, while $\omega^{\mu\nu}\equiv\frac{1}{2}(\nabla^{\mu}u^{\nu}{-}\nabla^{\nu}u^{\mu})$ is the vorticity tensor. This “DNMR” theory Denicol:2012cn can also be derived from a Chapman-Enskog like iterative solution of the RTA Boltzmann equation Jaiswal:2013npa .

Carrying the Chapman-Enskog expansion to one additional order, a third-order evolution equation for the shear stress was derived for the same system in Jaiswal:2013vta :

[TABLE]

For conformal systems the complete set of possible third order terms was given in Grozdanov:2015kqa . The specific form (II) implements transport coefficients obtained from the RTA Boltzmann equation Jaiswal:2013vta . We will here refer to Eq. (II) as the “third-order” theory. Let us now simplify all of these equations for Bjorken flow.

III Bjorken Flow

Milne coordinates $x^{\mu}=(\tau,x,y,\eta_{s})$ (with $\tau=\sqrt{t^{2}{-}z^{2}}$ and $\eta_{s}=\tanh^{-1}(z/t)$) are the natural choice for describing ultra-relativistic heavy-ion collisions where the colliding nuclei approach each other approximately following light-cone trajectories. For Bjorken flow Bjorken:1982qr of transversally homogeneous and longitudinally boost-invariant systems, macroscopic fields such as the energy density, pressure and shear stress, can depend neither on the transverse coordinates $(x,y)$ nor on the space-time rapidity $\eta_{s}$, but only on the longitudinal proper time $\tau$. The hydrodynamic evolution equations thus reduce to a set of coupled ordinary differential equations (ODEs) in $\tau$. The flow is irrotational ($\omega^{\mu\nu}=0$) and unaccelerated ($u^{\mu}=(1,0,0,0)$, $\dot{u}^{\mu}=0$), but (owing to the curvilinear nature of Milne coordinates) it has a non-zero local expansion rate, $\theta=1/\tau$, and velocity shear, e.g. $\sigma^{\eta\eta}=-2/(3\tau^{3})$.333To avoid clutter we drop the subscript on $\eta_{s}$ whenever we use it as a sub- or superscript. Symmetries further constrain the shear tensor to be diagonal and space-like in Milne coordinates, leaving only one independent component which we take to be the $\eta\eta$ component: $\pi^{xx}=\pi^{yy}=-\tau^{2}\pi^{\eta\eta}/2\equiv\pi/2$.444From here on we will use $\pi$ (not to be confused with the mathematical constant denoted by the same symbol), or its normalized version $\bar{\pi}\equiv\pi/(\epsilon{+}P)=3\pi/(4\epsilon)$, as the independent shear stress component. Note that for Bjorken flow the Navier-Stokes value for $\pi$ is positive, $\pi_{{}_{\mathrm{NS}}}\geq 0$. Using the following relations that hold for Bjorken flow,

[TABLE]

the shear evolution equations (4)-(II) can be brought into the following generic form:

[TABLE]

The coefficients $\beta_{\pi}$, $a$, $\lambda$, $\chi$, and $\gamma$ appearing in Eq. (9) above and in Eq. (11) below are tabulated in Table 1 for the three theories studied in this work.

Since $\beta_{\pi}=4P/5=4\epsilon/15$, Eqs. (8) and (9) are mutually coupled. Eq. (9) for the shear stress can be completely decoupled from the evolution of the energy density by rewriting it in terms of the normalized shear stress (inverse Reynolds number) $\bar{\pi}=\pi/(\epsilon{+}P)=\pi/(4P)$. Introducing at the same time the rescaled time variable Heller:2016rtz $\bar{\tau}\equiv\tau/\tau_{\pi}$ (which is the inverse Knudsen number for Bjorken flow), Eq. (8) can be used to obtain the relation

[TABLE]

Here we also used that for a conformal system $\epsilon{\,\propto\,}T^{4}$ and $T\tau_{\pi}=5\bar{\eta}=\mathrm{const.}$ where $\bar{\eta}\equiv\eta/s$ is the specific shear viscosity. Equations (9) and (10) can now be combined to obtain a first-order nonlinear ordinary differential equation (ODE) for the inverse Reynolds number that is completely decoupled555This decoupling works as long as the relaxation time has a power law temperature dependence, $\tau_{\pi}\propto T^{-\Delta}$, in which case in Eq. (11) the coefficient multiplying $d\bar{\pi}/d\bar{\tau}$ must be generalized to $\displaystyle{\frac{\Delta(\bar{\pi}{-}1){+}3}{3}}$. from the evolution of the energy density:666Even if the relaxation time $\tau_{\pi}$ depends on temperature (as it does, for example, in systems with conformal symmetry), the temperature evolution has been completely absorbed into the rescaled time variable (inverse Knudsen number) $\bar{\tau}$.

[TABLE]

The three hydrodynamic theories studied here can be selected by choosing for $\lambda$ and $\gamma$ the appropriate combinations of constants given in Table 1. All three theories share the same constant $a=4/15$, but we can solve Eq. (11) numerically777Equation (11) has the form of an Abel differential equation of the second kind for which, to the best of our knowledge, an analytical solution does not exist. (and the very closely related Eq. (28) to be discussed in Sec. V even analytically) for general $a$ and will therefore keep it as a free parameter until the end.

An important feature that we observe in Eq. (11) is that the derivative $d\bar{\pi}/d\bar{\tau}$ diverges for $\bar{\pi}\to-2$, indicating a discontinuity in $\bar{\pi}$ at $-2$. This feature is not present in Eq. (9) and is merely an artifact of changing the evolution parameter from $\tau$ to $\bar{\tau}$. Moreover, we notice that the transverse pressure, $P_{T}\equiv P+\pi/2=P(1+2\bar{\pi})$, becomes negative for $\bar{\pi}<-1/2$ indicating cavitation in the transverse direction which may lead to mechanical instability. Therefore one may conclude that $\bar{\pi}=-2$ already lies in the physically excluded region. In the next section, we study several other interesting properties inherent in the solutions of Eq. (11).

IV Gradient expansion, Lyapunov exponents and attractors

Consistent formulations of relativistic dissipative hydrodynamics involve short-lived non-hydrodynamic modes. These cannot be captured by a standard gradient expansion in terms of the Bjorken expansion rate $\theta=1/\tau$, causing such an expansion to be asymptotic, with zero radius of convergence Heller:2013fn . Borel resummation of this divergent series leads to a hydrodynamic attractor which is well defined even for large gradients. The existence of an attractor in a theory is indicated by negative Lyapunov exponents888Lyapunov exponents characterize the rate of separation of initially infinitesimally close trajectories in a dynamical system. A Lyapunov exponent is defined as

$|s(t)|\approx e^{\Lambda t}\,|s_{0}|,$

where $s_{0}$ and $s(t)$ are the separation between two trajectories at initial time $t_{0}$ and at a later time $t$, respectively. Negative $\Lambda$ indicates the rate at which the system approaches a regular attractor. which govern the rate at which the system loses information about its initial conditions and evolves towards the attractor. In this section, we study the gradient expansion, Lyapunov exponents and their implications for attractor behavior in solutions of hydrodynamic equations for Bjorken flow.

IV.1 Gradient expansion

To start with we consider the late proper-time expansion of $\bar{\pi}$. Substituting a power series ansatz for $\bar{\pi}$ in terms of powers of $1/\bar{\tau}$,

[TABLE]

into the non-linear ODE (11), one finds a recursion relation for the coefficients $c_{n}$ ($n\geq 1$)

[TABLE]

with initial value $c_{0}{\,=\,}0$. The first non-vanishing coefficient $c_{1}{\,=\,}a$ corresponds to the Navier-Stokes term in the gradient expansion.

For large $n$ the behavior of the coefficients $c_{n}$ is dominated by the term proportional to $c_{n-1}$ in Eq. (IV.1) and shows factorial growth. This is shown in Fig. 1 where the ratio of consecutive coefficients, $|c_{n+1}/c_{n}|$, is plotted against $n$; as shown in the inset this ratio is proportional to $n$ for $n\gtrsim 5$. Divergence of the hydrodynamic gradient expansion was found before in a variety of hydrodynamic theories Heller:2015dha ; Florkowski:2016zsi ; Denicol:2017lxn . A new observation from Fig. 1 is that, while the coefficients of the DNMR and third-order theories are very similar, those for MIS theory feature a constant upward shift which can be attributed to $\lambda=0$ in Eq. (IV.1). The inset of Fig. 1 shows the ratio $|c_{n+1}/c_{n}|$ for MIS theory ($\lambda{\,=\,}0$, red dots) compared with the same ratio obtained from Eq. (IV.1) without the nonlinear last term (black line). One concludes that this nonlinear term plays no role in the asymptotic factorial growth of the expansion coefficients. The series solution (12) is thus dominated by the term proportional to $c_{n-1}$ in Eq. (IV.1). We verified that the same statement holds for the DNMR and third-order theories.

IV.2 “Effective MIS” and Lyapunov exponents

The nonlinear last term in Eq. (IV.1) arises from the quadratic terms $\propto\bar{\pi}^{2},d\bar{\pi}^{2}/d\bar{\tau}$ in Eq. (11). At late times $\bar{\pi}{\,\ll\,}1$, hence $\bar{\pi}^{2}{\,\ll\,}\bar{\pi}$, and the nonlinear terms can be ignored in (IV.1). The late-time behavior in MIS theory is thus controlled by

[TABLE]

which we call “effective MIS theory”. Its analytic solution is

[TABLE]

where $\mathrm{Ei}[z]=-\int_{-z}^{\infty}\left(e^{-t}/t\right)dt$ is the exponential integral and $\alpha$ is the integration constant. Equation (15) implies that the separation between two solutions for $\bar{\pi}$ that are initialized with different initial conditions is damped exponentially:

[TABLE]

The negative Lyapunov exponent $\Lambda=-3/2$ in this equation confirms the existence of an attractor. This attractor, given by Eq. (15) with $\alpha{\,=\,}0$, is shown in Fig. 2 as the solid black line, together with a swarm of particular solutions of Eq. (14) with different integration constants $\alpha{\,\neq\,}0$ (dashed lines) and the asymptotic Navier-Stokes solution (red solid line) for comparison. One sees that the attractor, as well as the other solutions with different initial conditions, join the asymptotic Navier-Stokes behavior after $\bar{\tau}{\,\gtrsim\,}3$, i.e. for Knudsen number Kn${\,\lesssim\,}0.3$.

The solution (15) can also be written as

[TABLE]

where $\bar{\pi}_{0}\equiv\bar{\pi}(\bar{\tau}_{0})$ is the initial condition at time $\bar{\tau}_{0}$. From Eq. (17) one can extract the gradient expansion for $\bar{\pi}$ by integrating the last term by parts,

[TABLE]

resulting in

[TABLE]

where $\tilde{c}_{n}=(n{-}1)!/(3/2)^{n{-}1}$. The first term within parentheses contains all dependence on the initial conditions memory of which is here seen to decay exponentially with a decay time of $\frac{2}{3}\tau_{\pi}$. The second term is the divergent gradient series (12) up to order $m$. The third term can be thought of as the error in approximating the late time solution for $\bar{\pi}$ using a truncated gradient series of order $m$. Minimizing the error term with respect to $m$ would give an estimate of the optimal truncation for the gradient series at a given value of $\bar{\tau}$.

We point out that Eq. (17) also admits a convergent series in positive powers of $\bar{\tau}$ of the form

[TABLE]

It is obtained by expanding the exponential term $e^{\frac{3}{2}\bar{\tau}^{\prime}}$ in the integral solution for $\bar{\pi}$ into a power series. This series is, however, of limited use at late times as it would require including a large number of terms for convergence.

We conclude this subsection by offering another simple, yet interesting, form of the “effective MIS” solution for the inverse Reynolds number:

[TABLE]

Here the function $\beta(\bar{\tau};\bar{\tau}_{0})$ has units of time, and for any given pair $(\bar{\tau}_{0},\bar{\tau})$ its value can be shown to lie in the interval $[\bar{\tau}_{0},\bar{\tau}]$. To obtain this form we used the mean value theorem for the integral in Eq. (17). It would be illuminating to derive an approximate analytical functional form for $\beta(\bar{\tau};\bar{\tau}_{0})$; this is left for future work.

IV.3 Lyapunov exponents from linear perturbation

The exercise in the preceding subsection demonstrated exponential loss of memory of initial conditions in the “modified MIS” theory, on a time scale $\sim\frac{2}{3}\tau_{\pi}$. One arrives at similar conclusions without the assumption of ignoring in Eq. (11) the terms quadratic in $\bar{\pi}$, by studying the fate of a linear perturbation around late-time solutions of the full equation (11). This is important because, as we will see further below, the effects on the evolution of $\bar{\pi}$ of the non-linear terms that were ignored in “modified” MIS theory last much longer than the memory-decay time over which initial-state information is erased.

We employ linear perturbation theory around the solution $\bar{\pi}$, i.e. we set $\bar{\pi}\to\bar{\pi}+\delta\bar{\pi}$ and insert this into the generic evolution equation (11). To simplify the nonlinear terms $\sim\bar{\pi}\delta\bar{\pi}$ we expand the solution $\bar{\pi}$ for late times according to Eq. (12), keeping only the first non-vanishing term which yields the Navier-Stokes approximation: $\bar{\pi}\approx a/\bar{\tau}$. Keeping terms up to first order in $1/\bar{\tau}$ we find

[TABLE]

with the solution

[TABLE]

The prefactor $\bar{\tau}^{\frac{3}{4}a}$ in front of the exponential decay factor is easily traced back to the quadratic term $\bar{\pi}\left(d\bar{\pi}/d\bar{\tau}\right)$ on the left hand side of Eq. (11). The $\lambda$-dependence of this prefactor was not seen in the preceding subsection because $\lambda{\,=\,}0$ for MIS theory. The solution (23) exhibits the same negative Lyapunov exponent $\Lambda=-3/2$ as Eq. (16). This shows that keeping the non-linear terms in Eq. (11) does not change the initial state memory loss rate with which the system approaches the hydrodynamic attractor at late times. This rate is entirely controlled by microscopic physics and independent of the precise macroscopic dynamical state of the expanding system.

IV.4 Lyapunov exponents from Borel resummation

The results in the preceding subsection may also be looked at from the perspective of Borel resummation, which replaces a divergent series $\sum_{n}a_{n}$ by Kleinert:2001ax

[TABLE]

Borel resummation interchanges the order in which the sum and integral in $\sum_{n}a_{n}=\sum_{n}(a_{n}/n!)\int_{0}^{\infty}du\,e^{-u}\,u^{n}$ are performed. While for divergent series this yields inequivalent results, Borel resummation of asymptotic series can help with their interpretation, as shown below.

Starting from the recursion relation (IV.1) and ignoring the non-linearities therein one finds the coefficients of the gradient series to be

[TABLE]

where the normalisation ${\cal C}\equiv a(3/2)^{1{-}3\lambda/2}/\Gamma(1{-}3\lambda/2)$ ensures $c_{1}{\,=\,}a$. The Borel resummed version of Eq. (12) reads

[TABLE]

The Borel integrand has a pole at $u{\,=\,}1$ which results in an ambiguous imaginary part of the sum, $\mathrm{Im}\left[{\cal B}(\bar{\pi})\right]=\pm\pi\,{\cal C}\,e^{-3\bar{\tau}/2}\bar{\tau}^{-3\lambda/2}$, arising from non-unique choices of deforming the contour around the pole. In order to remove the ambiguity in the Borel resummed gradient series one must include in Eq. (12) a “non-hydrodynamic” term $\delta\pi\sim e^{-3\bar{\tau}/2}\bar{\tau}^{-3\lambda/2}$.999For a detailed study of this issue see the discussion of BRSSS theory Baier:2007ix in Ref. Basar:2015ava . Notice the similarity between this non-hydrodynamic term and that obtained via the method of linear perturbation in Eq. (23).101010The extra factor of $\bar{\tau}^{3a/4}$ in the latter stems from non-linear terms in Eq. (IV.1) which were here ignored. It should be emphasized that when considering the full non-linear version of Eq. (11), not just one, but an entire series of exponentially damped terms must be added to the gradient series, turning it into a trans-series Heller:2015dha ; Basar:2015ava . In this case the above-mentioned non-hydrodynamic term plays the role of the leading order correction.

IV.5 Numerical attractors for various theories

We close this section with a numerical study of the attractor solution of Eq. (11) towards which specific solutions with generic initial conditions decay exponentially. To identify the attractor we follow the prescription outlined in Ref. Heller:2015dha : The initial condition for the attractor solution is obtained by imposing the boundary conditions that both $\bar{\pi}$ and $d\bar{\pi}/d\bar{\tau}$ remain finite as $\bar{\tau}\to 0$. As shown in Heller:2015dha this results in the quadratic equation $\gamma\,\bar{\pi}^{2}+\lambda\,\bar{\pi}-a=0$ for the initial value of $\bar{\pi}$. One finds two solutions, one of which (the positive root) is stable and corresponds to the attractor solution whereas the negative root corresponds to a repulsor. With this initial condition we then solve Eq. (11) numerically, for the different parameter combinations listed in Table 1. The resulting attractors for the MIS, DNMR, and third-order theories are shown in Fig. 3 where we also compare them with the exact numerical attractor of the RTA Boltzmann equation Romatschke:2017vte . Please note that, unlike any particular solution of Eq. (11) which depends on both the initial condition $\bar{\pi}(\bar{\tau}_{0})$ and the initial time $\bar{\tau}_{0}$ at which it is imposed, the attractors are universal, i.e. they attract any particular solution with initial conditions within their basin of attraction, irrespective of its starting time.

In Fig. 3, we show the attractor solutions for $\bar{\pi}$ for the MIS, DNMR, and third-order theories, as well as for the exact solution of RTA Boltzmann equation Strickland:2018ayk and the Navier-Stokes solution. We see that of these the MIS attractor approaches the exact attractor most slowly,111111In light of footnote 2 this should perhaps not be too surprising. while the attractor of the third-order theory exhibits the best agreement with the exact RTA BE attractor, almost as good as the anisotropic hydrodynamic (aHydro) attractor studied in Strickland:2017kux ; Chattopadhyay:2018apf . This adds to the evidence of the superior performance of the third-order theory over different variants of second-order theories that are based on expansions around a locally isotropic momentum distribution.121212aHydro, a second-order approach that is based on an expansion around a self-consistently adjusted ellipsoidally deformed local momentum distribution Florkowski:2013lya ; Tinti:2015xwa ; Molnar:2016gwq ; Strickland:2017kux ; Martinez:2017ibh ; McNelis:2018jho , performs even better than the third-order theory Chattopadhyay:2018apf .

V Approximate analytical solutions

Up to this point we focused our discussion of the evolution of the inverse Reynolds number $\bar{\pi}$ on Eq. (11) which holds for conformal systems where $T\tau_{\pi}{\,=\,}$const. It has the advantage of completely absorbing any dependence of the shear relaxation time $\tau_{\pi}$ on the energy density or temperature into the rescaled time variable (inverse Knudsen number) $\bar{\tau}$, but at the expense of not being able to solve this ODE analytically. In this section we derive analytical solutions for the evolution of $\bar{\pi}$ for Bjorken flow, at the expense of not being able to ensure the conformal relation $T\tau_{\pi}{\,=\,}$const. consistently with the evolution of the energy density. Instead, we find three separate classes of analytical solutions, corresponding to three different approximations of $\tau_{\pi}$ as a function of time.

Starting from Eqs. (8),(9), we decouple them as before by rewriting them in terms of the inverse Reynolds number $\bar{\pi}$ but without rescaling the time Chattopadhyay:2018pwe :

[TABLE]

In the following, we find analytical solutions of Eq. (28), using different approximations for the form of shear relaxation time $\tau_{\pi}$.

V.1 Constant relaxation time

In this subsection we ignore the scaling of $\tau_{\pi}$ with temperature, by simply setting it constant. This constitutes a rather drastic violation of conformal symmetry by introducing, in addition to the inverse temperature $1/T$, a second, independent length scale $\tau_{\pi}$. In the following two subsections we will successively improve on this approximation.

The first analytical solution for the evolution in second-order hydrodynamics with Bjorken flow of the energy density and inverse Reynolds number for a constant shear relaxation time $\tau_{\pi}$ was found by Denicol and Noronha Denicol:2017lxn . For completeness we briefly review this solution, generalizing it to the generic form (27),(28) of the evolution equations which also includes third-order hydrodynamics. Introducing again the rescaled time $\bar{\tau}=\tau/\tau_{\pi}$, for constant $\tau_{\pi}$ Eq. (28) turns directly into

[TABLE]

which is similar to Eq (11) but without the nonlinearity on the left hand side.131313Equation (29) can be readily derived by setting $\Delta=0$ in footnote 5. As will be discussed below, this difference has important consequences for the attractor solutions and Lyapunov exponents.

Equation (29) is a first-order nonlinear ODE of Riccati type which can be written as a second-order linear ODE with the help of the following transformation of variables,

[TABLE]

which turns Eq. (29) into

[TABLE]

The general solution of this linear ODE can be expressed in terms of Whittaker functions $M_{k,m}(\bar{\tau})$ and $W_{k,m}(\bar{\tau})$ Denicol:2017lxn 141414Note that the authors of Denicol:2017lxn used $\gamma=4/3$ and a different definition of $a$.:

[TABLE]

Here $k=-\frac{1}{2}(\lambda{+}1)$ and $m=\frac{1}{2}\sqrt{4a\gamma+\lambda^{2}}$ while $A$ and $\alpha$ are arbitrary constants. Substituting this solution in Eqs. (30) and (27) one finds

[TABLE]

Here $\epsilon_{0}$ is the initial energy density at time ${\bar{\tau}}_{0}$, and the constant $\alpha$ encodes the initial normalized shear stress $\bar{\pi}_{0}$. Note that $\alpha$ can only take values for which the energy density is positive-definite for $\bar{\tau}>0$.

It is easy to see from Eq. (33) that the solution for $\bar{\pi}(\bar{\tau})$ loses all memory about initial conditions at late times when $M_{k,m}(\bar{\tau})$ dominates over $W_{k,m}(\bar{\tau})$: for large arguments $\bar{\tau}\to\infty$ the ratio

[TABLE]

decays exponentially. Thus at late times the terms proportional to $\alpha$ in (33), (34) decay like $e^{-\bar{\tau}}$, corresponding to a Lyapunov exponent of $\Lambda{\,=\,}-1$. This obviously differs from $\Lambda{\,=\,}-\frac{3}{2}$ for the conformally invariant theories described by Eq. (11); the difference is a direct consequence of the breaking of conformal symmetry by setting $\tau_{\pi}$ constant instead of $\propto 1/T$.

V.2 Relaxation time from ideal hydrodynamics

A better approximation to Eq. (11) can be obtained by setting $T\tau_{\pi}{\,=\,}$const. but, instead of using the exact time dependence of the temperature $T$, approximating the latter such that the evolution equations can still be integrated analytically. Seeing that for Bjorken flow the system asymptotically approaches local thermal equilibrium, we may approximate the time-dependence of $T$ at late times by the ideal fluid law Bjorken:1982qr

[TABLE]

where $T_{0}$ is the temperature at initial time $\tau_{0}$. For $T\tau_{\pi}=5\bar{\eta}$ this yields

[TABLE]

Using this to define the scaled time variable $\bar{\tau}\equiv\tau/\tau_{\pi}$ one finds

[TABLE]

and Eq. (28) turns into

[TABLE]

independent of $b$. This equation again misses the nonlinear term on the left hand side (l.h.s.) of Eq. (11) and, except for the factor $2/3$ on the l.h.s., has the same structure as Eq. (29). Its analytical solution is therefore very similar to Eq. (33), except for a change of the argument of the Whittaker functions by a factor 2/3:

[TABLE]

Here $k=-\frac{3\lambda+2}{4}$ and $m=\frac{3}{4}\sqrt{4a\gamma+\lambda^{2}}$, different from Eqs. (33)-(34). The asymptotic behavior (35) of the Whittaker functions now tells us that, for the choice (37), memory of the initial conditions is lost exponentially according to $e^{-\frac{3}{2}\bar{\tau}}$, corresponding to the same Lyapunov exponent $\Lambda{\,=\,}-\frac{3}{2}$ as for the conformally invariant theories described by Eq. (11).

V.3 Relaxation time from Navier-Stokes evolution

We can further improve our approximation by accounting for first-order gradient effects in the evolution of the temperature, by replacing the ideal fluid law (36) by the Navier-Stokes result Kouno:1989ps ; Muronga:2001zk ; Baier:2006um

[TABLE]

For $\bar{\eta}{\,=\,}0$ this reduces to (36). Substituting this into $T\tau_{\pi}=5\bar{\eta}$ we find

[TABLE]

For the scaled time variable $\bar{\tau}\equiv\tau/\tau_{\pi}$ we now have

[TABLE]

Using this in Eq. (28) one obtains

[TABLE]

independent of the constant $d$. This shares with Eq. (38) the factor 2/3 on the l.h.s. which, as we saw in the preceding subsection, leads to the correct Lyapunov exponent for conformally symmetric systems. Comparing with Eq. (11) one sees that they are identical up to the substitution $\bar{\pi}\mapsto a/\bar{\tau}$ (which is the first non-zero term in the series expansion (12)). This should not be surprising as the term within parenthesis on the l.h.s. of Eq. (11) stems solely from the energy (or, equivalently, temperature) evolution equation (8) which, when making the replacement $\bar{\pi}\mapsto a/\bar{\tau}$, leads to the Navier-Stokes solution $T_{{}_{\rm NS}}$.

Comparing with the preceding subsection, this suggests that it might be possible to account for the time evolution of the temperature in the relation $\tau_{\pi}(\tau)\sim 1/T(\tau)$ with ever increasing precision by substituting the gradient series (12) for $\bar{\pi}$ in the prefactor on the l.h.s. of Eq. (11) and truncating it at increasingly higher order: zeroth order for ideal hydrodynamics, first order for Navier-Stokes dynamics, and so on. Unfortunately, this is not justified as Eq. (12) is an asymptotic (i.e. divergent) series. The consequences of this on the hydrodynamic attractor will be discussed in the next section.

Equation (44) can be recast in a form similar to Eq. (29) and solved again analytically:

[TABLE]

This looks formally identical to Eqs. (39)-(40), except for the substitution $\frac{3}{2}\bar{\tau}\mapsto w\equiv\frac{3}{2}(\bar{\tau}+\frac{a}{2})$ in the arguments of the Whittaker functions on the right hand side (r.h.s.), together with modified definitions for the indices $k$ and $m$:

[TABLE]

Using Eq. (35) for the asymptotic behavior of the Whittaker functions we see that once again memory of the initial conditions is lost exponentially according to $e^{-\frac{3}{2}\bar{\tau}}$, with no effect from the constant shift of the time variable in the arguments of the Whittaker functions. This corresponds to the same Lyapunov exponent $\Lambda{\,=\,}-\frac{3}{2}$ as in the preceding subsection and, more generally, for all the conformally invariant theories described by Eq. (11).

We note in passing that including higher order terms $c_{n}/{\bar{\tau}}^{n}$ for $n{\,>\,}1$ in the gradient series while approximating $\bar{\pi}$ on the l.h.s. of Eq. (44) spoils its reducibility to the analytically solvable form explored in this work. Whether for such approximations Eq. (11) can still be solved in terms of known functions remains to be seen.

We summarize this section by observing that, while the assumption of a fixed shear relaxation time $\tau_{\pi}$ (i.e. of a fixed temperature $T$ when writing $\tau_{\pi}\propto 1/T$) leads to an incorrect Lyapunov exponent describing the rate of approach towards the hydrodynamic attractor, the correct decay rate is recovered as soon as one allows the temperature $T$ to vary with time hydrodynamically even if the exact time dependence is replaced by an approximation based on a truncated hydrodynamic gradient expansion.

VI Analytical attractors

In this section we investigate the hydrodynamic attractors associated with the analytic approximate solutions of the evolution equation for the inverse Reynolds number $\bar{\pi}$ found in the preceding section. In Section IV.5 we saw that, as $\bar{\tau}\to 0$, $\bar{\pi}$ can take one of only two finite values of opposite sign, and we identified the attractor as the unique solution which connects to the positive value. All other solutions were found to connect in the limit $\bar{\tau}\to 0$ to the negative value. An illustration of this generic behavior is shown in Fig. 4 for the case of the attractor corresponding to the analytical solution (39) for the third-order theory.

In this section we introduce the following procedure for identifying the hydrodynamic attractor Chattopadhyay:2018pwe : In terms of the parameter $\alpha$ encoding the initial condition for $\bar{\pi}$, we search for the value $\alpha_{0}$ at which the quantity

[TABLE]

diverges at the time $\bar{\tau}_{0}$ where the two fixed points of the evolution trajectories (see Fig. 4) are located Behtash:2017wqg . For Bjorken flow, this time is usually $\tau_{0}{\,=\,}0$. For the numerical solution of Eq. (11) in Sec. IV.5 this can be seen in Fig. 3, and for the analytical solutions Eqs. (33) and (39) in Secs. V.1 and V.2, respectively, this can be seen by studying their behavior near $\bar{\tau}=0$. The structural similarity of Eqs. (45) and (39) shows that for the solution (45) the fixed points are instead located at $w_{0}{\,=\,}0$, which corresponds to a negative (i.e. unphysical) longitudinal proper time $\bar{\tau}_{0}=-a/2$. This arises from the Navier-Stokes substitution $\bar{\pi}\mapsto a/\bar{\tau}$ on the l.h.s. of Eq. (44) which breaks down (and thus renders the analytic solution (45) unreliable) in the region $\bar{\tau}\ll 1$.

We can use the same approach of studying the sensitivity to initial conditions to obtain the Lyapunov exponent $\Lambda$ from the formula

[TABLE]

For the analytical solutions (33), (39), (45) this prescription reproduces the same results as obtained from the late-time behavior (35) of the Whittaker functions but its advantage is that it can also be used numerically where exact solutions for $\bar{\pi}$ are not available (such as for the numerical solutions of the generic equation (11)).

For the case shown in Fig. 4 we have verified explicitly that indeed there are two fixed points at $\bar{\tau}_{0}=0$ and that only for the attractor solution, characterized by $\alpha_{0}{\,=\,}0$, $\psi(0)$ (defined in Eq. (47)) diverges. For all other solutions $\alpha\neq 0$ we found that $\psi(\alpha)=0$ in the limit $\bar{\tau}\to 0$, indicating that they all converge to the negative branch, as shown in Fig. 4.

The approximate analytical solutions (33), (39), (45) discussed in Sec. V can be written in the generic form

[TABLE]

with arguments and parameters for the three cases (i.e. for $\tau_{\pi}\propto 1/T$ with $T$ either constant or with time dependence taken from ideal or Navier-Stokes hydrodynamics) compiled for convenience in Table 2.

The attractor solutions are obtained from Eq. (49) by setting the initial condition parameter $\alpha{\,=\,}0$:

[TABLE]

They are shown in Fig. 5 for the three different hydrodynamic theories discussed in this paper (MIS (a), DNMR (b), and third-order (c)) and compared with the corresponding exact numerical attractors as well as with the attractor for the exact analytical solution of the RTA Boltzmann equation Romatschke:2017vte (the latter is, of course, the same in all three subpanels). Comparison of these attractors provides insights not only about the performance of the three different hydrodynamic theories as approximations to the underlying kinetic theory, but also about the relative accuracy of the additional approximations made in Sec. V in order to obtain analytical results.

For all three hydrodynamic theories, we note (especially at early times when the system is farthest away from local equilibrium) that the differences between the exact numerical attractors and their analytical approximations from Sec. V are significantly smaller than their discrepancy from the attractor of the underlying kinetic theory. At late times the breaking of conformal symmetry by choosing a constant relaxation time leads generally to the largest difference between the exact numerical and analytically approximated hydrodynamic attractors; this may be attributed to the fact that the $\tau_{\pi}{\,=\,}$const. approximation underestimates the rate of approach towards the attractor by a factor 2/3. More surprisingly, the analytic approximation that performs best at late times (which uses $\tau_{\pi}(\tau)\sim 1/T_{{}_{\mathrm{NS}}}(\tau)$) performs worst at early times when compared with the exact numerical result. This reflects a different value for the fixed point of the inverse Reynolds number compared to the other analytic approximations and may be related to the fact that in this case the fixed point is shifted outside the physical region to $\bar{\tau}=-a/2$. As already seen in Fig. 3, in comparison with the exact solution of the RTA Boltzmann equation, the third-order theory performs much better than both the MIS and DNMR theories, both when evaluated exactly numerically or approximately analytically. The only known theory that performs even better than the third-order theory studied here is second-order anisotropic hydrodynamics Bazow:2013ifa ; Strickland:2017kux ; Chattopadhyay:2018apf which effectively resums terms of all orders in the inverse Reynolds number.

Numerical studies, such as those shown in Figs. 2 and 4, show that at late times $\bar{\tau}{\,>\,}1$ any initial deviation from the attractor approaches the attractor exponentially, with the Lyapunov exponents discussed before. The approximate analytical result (49) allows to understand this approach analytically over the entire range of $\bar{\tau}$, i.e. also for large Knudsen numbers $\bar{\tau}\ll 1$. Following the analysis Kurkela:2019set we write151515Note that $\delta$ as defined in Kurkela:2019set differs from ours by a factor 4/3.

[TABLE]

where in the last expression we expanded the deviation $\delta$ to first order in the deviation of the initial value parameter $\alpha$ from the value $\alpha{\,=\,}0$ characterizing the attractor solution. The decay of the inverse Reynolds number $\bar{\pi}$ towards its attractor value is for small deviations161616It is worth pointing out that, for fixed initial deviation $\delta_{0}$ at initial time $\bar{\tau}_{0}$, the corresponding initial state parameter $\alpha$ approaches 0 as $\bar{\tau}_{0}\to 0$. Small $\delta_{0}$ thus implies small $\alpha$ (especially for small values of $\bar{\tau}_{0}$), but not vice versa. $\delta$ given by

[TABLE]

This is a function of the scaling variable $w\propto|\Lambda|\bar{\tau}$ which is proportional to the inverse Knudsen number $\bar{\tau}$. At first sight this suggests that the competition between the global expansion rate $1/\tau$ and the microscopic relaxation rate $1/\tau_{\pi}$ not only rules the evolution of the hydrodynamic attractor $\bar{\pi}_{\mathrm{attr}}$ itself, but also the way initial excursions of the inverse Reynolds number from this attractor decay as the full dynamical solution approaches the attractor. However, as first pointed out in Kurkela:2019set , a deeper analysis exhibits that the dependence of the lifetime of such excursions on the interaction rate $1/\tau_{\pi}$ changes dramatically between early times $\tau{\,\ll\,}\tau_{\pi}/|\Lambda|$ and late times $\tau{\,\gg\,}\tau_{\pi}/|\Lambda|$.

Using the following Whittaker function properties for small arguments,

[TABLE]

where the approximation holds at leading order, we see from Eq. (52) that, in contrast to the exponential decay (35) at late times $\bar{\tau}\gg 1$, the decay of excursions away from the attractor decay with a power law at early times $\bar{\tau}\ll 1$ Kurkela:2019set . The transition from power law to exponential decay around $w\propto|\Lambda|\bar{\tau}=1$ is illustrated in Fig. 6.171717Note that the “effective MIS” approximation (15) discussed in Sec. IV.2 misses the early-time power-law decay of excursions away from the hydrodynamic attractor. We found that early-time power-law decay of initial deviations from the attractor requires that in the generic evolution equation (11) for the inverse Reynolds number at least one of the two coefficients $\lambda$ and $\gamma$ must be nonzero. This is not the case for the “effective MIS” theory.,181818A few comments regarding the domain of validity of the linearised approach are in order. Using Eq. (49) to determine the deviation $\delta(w)$ from the attractor solution, one readily checks that, for $w\ll 1$, terms beyond the linear order in Eq. (51) are proportional to $\alpha^{n}(w^{-2m})^{n}$ for $n\geq 2$. Even for fixed small values of $\alpha$, this series diverges for sufficiently early “times” $w^{2m}\ll\alpha$, where linearization in $\alpha$ must break down. This is closely related to footnote 16: in order to keep $\delta_{0}\equiv\delta(w_{0})$ fixed for small values of $w_{0}$, $\alpha$ must be tuned down to ensure $\alpha w_{0}^{-2m}\ll 1$ such that higher-order contributions may be safely neglected. The breakdown of the linearised method at very early times owes itself to the presence of the repulsive fixed point at $w=0$, i.e., as long as $\alpha$ is not strictly equal to 0, all solutions, $\bar{\pi}(\alpha,w)$ will ultimately hit the repulsor at $w=0$, where the deviation from the attractor no longer stays infinitesimal. This mathematical feature of $\bar{\pi}(\alpha,w=0)$ being discontinuous at $\alpha=0$ was, in fact, used to uniquely determine the attractor solution via Eq. (47) Chattopadhyay:2018pwe .

An immediate consequence of this discussion is that for early starting times $\tau_{0}{\,\ll\,}\tau_{\pi}/|\Lambda|$ (where, for sufficiently small $\delta$, $\delta(\bar{\tau})\propto(|\Lambda|\bar{\tau})^{-2m}$) the decay time $\tau_{1/e}$ (defined as the time over which the magnitude of $\delta$ decreases by a factor $1/e$) is given by $\tau_{1/e}=\rho\,\tau_{0}$ (where, to leading order in $\delta$ and $\bar{\tau}$, $\rho$ is a constant given in terms of the attractive fixed point $\bar{\pi}_{\mathrm{attr},0}$ by $\rho=(e^{1/2m}{-}1)^{\frac{3}{\bar{\pi}_{\mathrm{attr},0}{+}2}}$) and thus independent of $\tau_{\pi}$, whereas for late starting times $\tau_{0}{\,\gg\,}\tau_{\pi}/|\Lambda|$ it is given by $\tau_{1/e}=\tau_{\pi}/|\Lambda|$, independent of $\tau_{0}$. $\tau_{1/e}$ thus smoothly increases from 0 to $\tau_{\pi}/|\Lambda|$ (where it saturates) as $\tau_{0}$ increases from 0 to values much larger than $\bar{\tau}/|\Lambda|$. Sufficiently small excursions of the inverse Reynolds number from its attractor thus decay faster at early times $\tau{\,\ll\,}\tau_{\pi}/|\Lambda|$ than at late times $\tau{\,\gg\,}\tau_{\pi}/|\Lambda|$.

Based on the observation that for $\tau\ll\tau_{\pi}/|\Lambda|$ the decay rate is controlled by the initial expansion rate $1/\tau_{0}$ and independent of the scattering rate $1/\tau_{\pi}$ whereas the opposite holds for $\tau\gg\tau_{\pi}/|\Lambda|$, the authors of Kurkela:2019set attribute the transition from power-law decay at early times to exponential decay at later times to a transition from the “pre-hydrodynamic” to the “hydrodynamic” stage, i.e. they see it as associated with the process of hydrodynamization. We have a different view of this matter: It is known from earlier numerical studies Strickland:2017kux ; Behtash:2017wqg ; Chattopadhyay:2018apf of both Bjorken and Gubser flows that some hydrodynamic approximations (in particular anisotropic and third-order hydrodynamics, but to a somewhat lesser degree also DNMR theory) evolve the inverse Reynolds number very accurately (when compared with the evolution predicted by the underlying RTA Boltzmann equation) already at very early times $\bar{\tau}\ll 1$, even when initialized far away from the attractor. That is, not only do the attractors for these hydrodynamic approximations agree well with the exact attractor even for large Knudsen numbers, as shown for DNMR and third-order theories here in Fig. 3 and for anisotropic hydrodynamics in Fig. 3 of Strickland:2017kux , but the theories also describe accurately the evolution of the system towards the attractor even when initialized far away from it. To us this implies that the system hydrodynamizes well before $\bar{\tau}=1/|\Lambda|$, and power-law rather than exponential decay of deviations from the attractor are not a tell-tale signature for “pre-hydrodynamic” behavior.

VII Summary and conclusions

We studied analytically and numerically the evolution of the inverse Reynolds number in causal theories of second- and third-order relativistic viscous fluid dynamics for Bjorken flow. In this situation there is only a single non-vanishing component of the shear stress, describing an anisotropy between longitudinal and transverse pressure, which is generically very large at early times. For Bjorken flow the evolution of the associated inverse Reynolds number (i.e. the ratio of the shear stress to the enthalpy density of the system) decouples from that of the energy density and temperature and thus can be solved independently. When expressed as a function of time in units of the microscopic shear relaxation time (which measures the inverse Knudsen number of Bjorken flow), the solution is universal, i.e. independent of the specific shear viscosity of the medium. For three different macroscopic hydrodynamic theories, we studied these solutions, numerically and with various analytical approximations, their hydrodynamic attractors, the rate of initial state memory loss and approach to the attractor (expressed through Lyapunov exponents), and compared all these with the corresponding solution of the relativistic Boltzmann equation in RTA approximation which describes the underlying microscopic dynamics and can, for Bjorken flow, be solved exactly.

When comparing the exact numerical solutions for the attractor of the inverse Reynolds number for the three different hydrodynamic theories studied here with the exact solution from the Boltzmann equation we find significant differences at early times, i.e. at large Knudsen numbers where the dynamics happens far away from equilibrium. While for the third-order theory the discrepancy remains always below 10%,191919It has been shown elsewhere Strickland:2017kux ; Chattopadhyay:2018apf that the attractor for anisotropic hydrodynamics (aHydro) is even closer to the exact Boltzmann attractor than the one for the third-order theory. it increases to ${\sim\,}30\%$ for DNMR theory and to almost a factor of 2 for MIS theory (both being second-order hydrodynamic theories). Compared to these, the additional discrepancies caused by the various approximations we made to arrive at analytical solutions for the dynamics of $\bar{\pi}$ are small. At late times (small Knudsen numbers) all attractors approach the Navier-Stokes solution. Again third-order hydrodynamics is closest to the exact solution of the underlying Boltzmann kinetics, but the excellent agreement is spoiled somewhat if conformal symmetry is broken by an approximation that sets the relaxation time $\tau_{\pi}$ to a constant rather than allowing it to vary inversely with the temperature.

As $\tau\to 0$, the Navier-Stokes value of the inverse Reynolds number diverges while its attractor value approaches the finite value $\bar{\pi}=0.25$, corresponding to a shear stress over thermal pressure ratio $\pi/P{\,=\,}1$. The Lyapunov exponent associated with the evolution of $\bar{\pi}$, $\Lambda{\,=\,}-\frac{3}{2}$, indicates that even far-from-equilibrium initial conditions, $\pi_{0}/P_{0}\gg 1$, relax exponentially to the attractor value with a decay time of order $\frac{2}{3}\tau_{\pi}=\frac{10}{3}\frac{\bar{\eta}}{T}$. For minimal specific shear viscosity $\bar{\eta}=\frac{1}{4\pi}$ and a medium temperature of, say, $T{\,=\,}0.5$ GeV, this corresponds to a decay time of $\approx 0.1$ fm/$c$. Subsequently, the hydrodynamic evolution follows essentially the hydrodynamic attractor of the theory which agrees, within the precision stated above, with the exact attractor associated with of the underlying microscopic Boltzmann kinetics.

While the ODE describing the evolution of the inverse Reynolds number for Bjorken flow is easily solved on a computer, with arbitrary precision, the analytic approximations studied here are surprisingly accurate, and they yield valuable insights into the details of initial state memory loss and the approach to attractor dynamics in Bjorken flow. Similar methods may be applicable to different situations (for example, Gubser flow, which is physically quite different from Bjorken flow but shares with it many mathematical similarities) where they can lead to similarly valuable qualitative insights. Generalization to Gubser flow should be particularly interesting because it does not thermalize at late times but rather approaches an asymptotic free-streaming state. In such a situation the questions of initial state memory loss and the approach to a hydrodynamic attractor Behtash:2017wqg ; Chattopadhyay:2018apf ; Denicol:2018pak have not yet been fully understood.

Acknowledgements.

S.J. and A.J. acknowledge kind hospitality of IISER Pune where parts of this work were completed. We gratefully acknowledge interesting discussions with A. Behtash, S. Kamata, M. Martinez and H. Shi as well as the stimulating effect of Ref. Kurkela:2019set which appeared simultaneously when the present work was first circulated and which led to several additional insights reported in Sec. VI of this final version of our work. U.H. thanks Urs Wiedemann for an important clarifying discussion of Ref. Kurkela:2019set . A.J. is supported in part by the DST-INSPIRE faculty award from the Indian Department of Science and Technology under Grant No. DST/INSPIRE/04/2017/000038. The research of U.H. and C.C. was supported in part by the U.S. Department of Energy (DOE), Office of Science, Office for Nuclear Physics under Award No. DE-SC0004286 and in part by the National Science Foundation (NSF) within the framework of the JETSCAPE Collaboration under Award No. ACI-1550223. S.J. thanks Pranav S. and C.C. thanks Yuri Kovchegov for illuminating discussions. U.H. acknowledges the kind hospitality of the Institut für Theoretische Physik at the J. W. Goethe-Universität during a research visit supported by a Research Prize from the Alexander von Humboldt Foundation. He thanks Michał Heller and Viktor Svensson for fruitful discussions that motivated some of the research reported here.