Stable and causal relativistic Navier-Stokes equations

Relativistic Navier-Stokes equations express the conservation of the energy-momentum tensor and the particle number current in terms of the local hydrodynamic variables: temperature, fluid velocity, and the chemical potential. We show that the viscous-fluid equations are stable and causal if one adopts suitable non-equilibrium definitions of the hydrodynamic variables.


Introduction
Hydrodynamics is a classical effective description of macroscopic states of matter with small deviations from local thermal equilibrium. Hydrodynamics is conventionally formulated by starting from thermodynamics, and then promoting the constant parameters of global thermal equilibrium (temperature T , fluid velocity v, etc.) to slowly varying functions in space and time: T (t, x), v(t, x), etc. The evolution of these hydrodynamic variables is then determined by the local conservation laws of energy, momentum, and possibly other conserved quantities such as mass or particle number [1]. Intrinsic to hydrodynamics is thus a fundamental ambiguity: one must (somewhat arbitrarily) make a choice as to how to define "local temperature", "local fluid velocity", etc., out of equilibrium. For example, the non-equilibrium "fluid velocity" may be chosen to correspond to the flow of particles, or to the flow of energy, or to the flow of entropy, etc., with each choice resulting in different hydrodynamic equations.
In non-relativistic Navier-Stokes equations, the standard convention is to define the "fluid velocity" through the flow of mass [1]. Relativistic hydrodynamics was presented originally in two different formulations, one pioneered by Eckart [2] and one by Landau and Lifshitz [1]; both are still widely discussed. The non-equilibrium conventions differ between the two formulations. Consequently, the hydrodynamic equations of Eckart and of Landau and Lifshitz are different, mathematically inequivalent, equations. More generally, the arbitrariness in adopting different non-equilibrium definitions implies that there is simply no such thing as "the" equations of hydrodynamics. Still, one expects that conventions should not be physically relevant, and that different hydrodynamic equations must give rise to the same physical predictions within the domain of applicability of hydrodynamics.
The standard formulation of relativistic hydrodynamics (as in Ref. [1]) uses the following variables: temperature T , chemical potential µ, and fluid velocity u α . The first two are scalars, and the latter is a vector. The chemical potential is conjugate to a global conserved U (1) charge, such as the baryon number. The hydrodynamic equations are the conservation laws for the energy-momentum tensor T αβ and the corresponding U (1) current J α , where ∇ denotes the covariant derivative for the spacetime metric g αβ . The conservation laws have to be supplemented with the constitutive relations which express T αβ and J α in terms of T , µ, and u λ . We will take the constitutive relations T αβ = T αβ [T, µ, u λ , g ρσ ], J α = J α [T, µ, u λ , g ρσ ] to be local functions of the hydrodynamic variables, the metric, and their derivatives. The constitutive relations are then written as an expansion in derivatives, where the subscript denotes the number of derivatives. For example, T αβ (2) will have contributions proportional to (∂T )(∂µ), ∂ 2 u, as well as the purely geometric contributions proportional to the Ricci tensor R αβ , among others. The constitutive relations which only take into account T αβ (0) , J α (0) are said to correspond to "perfect fluid hydrodynamics", while the constitutive relations with terms up to T αβ (n) , J α (n) are said to describe "n-th order hydrodynamics". The standard physics of viscosity and heat conductivity is contained within first-order hydrodynamics. The equations of first-order hydrodynamics are often called the Navier-Stokes equations.
The first-order relativistic theories of Eckart, and of Landau and Lifshitz, suffer from two important pathologies: they both predict that the uniform thermal equilibrium state of a non-gravitating fluid in flat space is unstable [3], and they both predict that signals propagate faster than light [4]. The most popular remedy to these problems is provided by the Müller-Israel-Stewart (MIS) theories which introduce extra tensor variables besides T , u α , and µ into the hydrodynamic equations. See the recent book [5] for a modern perspective and references. The pathologies of the original hydrodynamic theories of [1,2] plus the successful practical applications of the MIS theories in describing the quark-gluon plasma produced in relativistic heavy-ion collisions [6,7] have led to a widespread belief which can be stated as "the relativistic Navier-Stokes equations are unstable and acausal".
Such a view, however, is misguided: as previously emphasized, there is no such thing as "the" relativistic Navier-Stokes equations. The freedom of convention 1 in defining the outof-equilibrium T , u α , and µ means that for the same physical fluid there are infinitely many different "Navier-Stokes equations"; the theories of [1] and [2] are just two examples. With some conventions, the Navier-Stokes equations are indeed unstable and acausal; with others, the Navier-Stokes equations may well be both stable and causal. One may as well choose a convention that makes physical sense.
In Refs. [8,9,10] (which we shall call BDNK) it was argued that there exist particularly convenient out-of-equilibrium definitions of the hydrodynamic variables, such that the equations of first-order hydrodynamics written in terms of these variables are causal, and the equilibrium state is stable. The BDNK conventions thus define stable and causal frames for relativistic fluids. The discussion in Refs. [9,10] was only concerned with "uncharged" fluids, i.e. fluids for which the only hydrodynamic variables are T and u α , and the only conservation laws are those of energy and momentum. In Ref. [8] general hydrodynamic field redefinitions were discussed for fluids with a global U (1) charge such as the baryon number (so called "charged" fluids); however, the stable and causal frames were only explicitly discussed for uncharged fluids.
The aim of the present paper is to extend the discussion of BDNK to fluids with a conserved U (1) charge, whose hydrodynamic equations are given by the conservation laws (1.1). For such fluids, we shall describe a class of stable and causal frames. In such frames, the relativistic Navier-Stokes equations are stable and causal, correcting the deficiencies of [1,2], but without introducing any extra variables besides the standard T , µ, and u α . The relation between the familiar Landau-Lifshitz frame and the general frames is described in the Appendix.

Constitutive relations
Following Ref. [2], we decompose the energy-momentum tensor and the current as where ∆ µν ≡ g µν + u µ u ν projects onto the space orthogonal to u, the vectors Q and J are transverse to u, and the tensor T is transverse to u, symmetric and traceless. For a given timelike vector u, the decomposition (2.1) defines the components E, Q µ , T µν , N , and J µ in terms of T µν and J µ . Using the notation of [8], the most general one-derivative constitutive relations in a charged fluid may be written as The dot signifies the derivative along the fluid velocity, i.e.Ṫ ≡ u λ ∂ λ T ,u α ≡ u λ ∇ λ u α . The shear tensor is where d is the number of spatial dimensions. We will refer to the coefficients ε i , π i , θ i , ν i , γ i as "transport parameters". There are fifteen transport parameters at one-derivative order, plus the shear viscosity η. We will work in the thermodynamic frame where the equilibrium constitutive relations follow from a partition function which is extensive in equilibrium. 2 Then the zero-derivative coefficients ǫ(T, µ), p(T, µ), n(T, µ) have the standard interpretations of the equilibrium energy density, pressure, and charge density, respectively, and moreover are related by n = ∂p/∂µ, ǫ + p − µn = T ∂p ∂T . Further, the same extensivity in the thermodynamic frame implies θ 1 = θ 2 , and γ 1 = γ 2 .
The bulk viscosity ζ and the charge conductivity σ are given by the following combinations of the transport parameters [8]: where the derivatives of the pressure are p ,ǫ ≡ (∂p/∂ǫ) n , p ,n ≡ (∂p/∂n) ǫ . Clearly, the transport parameters must be such that both ζ and σ are non-negative. A choice of "frame" corresponds to a choice of transport parameters. For example, the Landau-Lifshitz frame of [1] imposes that all one-derivative transport parameters vanish except π 2 and γ 3 . A stable and causal frame is a choice of the transport parameters such that the hydrodynamic equations (2.1), (2.2) are causal and predict that the thermal equilibrium state is stable.
In the grand canonical ensemble, the functions ǫ(T, µ) and n(T, µ) are not independent, as both are determined by the pressure p(T, µ). In particular, T (∂n/∂T ) + µ(∂n/∂µ) = ∂ǫ/∂µ. Further, we have the following thermodynamic inequalities: Denoting the Hamiltonian by H and the conserved particle number operator by N , the inequalities above follow by demanding that the connected equilibrium functions are non-negative: The equilibrium partition function is a functional of the external time-independent sources: the metric gµν and the gauge field Aµ that couples to the U (1) current. Equilibrium fluid velocity is defined to be aligned with the timelike Killing vector field V which specifies the direction of time, u µ = V µ / √ −V 2 , to all orders in the derivative expansion. Equilibrium temperature is defined so that the Tolman's law T = T0/ √ −V 2 holds, to all orders in the derivative expansion. See Ref. [11] for more details. Comparing T µν and J µ obtained by varying the equilibrium partition function with respect to gµν and Aµ with T µν and J µ given by the constitutive relations (2.1), (2.2) and evaluated in equilibrium, one obtains the constraints stated in the text.
where v s is the speed of sound, see e.g. [12].
The general constitutive relations (2.2) simplify if the underlying microscopic theory happens to be conformal. In d + 1 spacetime dimensions, conformal symmetry demands [8]: as well as θ 1 = θ 2 , γ 1 = γ 2 . The equation of state in a conformal theory is p(T, µ) = T d+1 f (µ/T ), with a dimensionless function f (µ/T ) which is determined by the microscopic dynamics. The inequalities (2.5) imply that the function f (x) must be such that The speed of sound in a conformal theory is v s = 1/ √ d, and the bulk viscosity ζ vanishes.

Small fluctuations in equilibrium
Let us now look at small fluctuations of the equilibrium state with constant T = T 0 , µ = µ 0 , and v = v 0 . Taking the fluctuations δT , δµ, δv i proportional to e −iωt+ik·x , the linearized hydrodynamic equations give rise to polynomial equations in ω and k which we schematically write as F (ω, k) = 0. Their solutions determine the dispersion relations ω = ω a (k). We take k to be real; the ω a (k) will be (in general, complex) functions of k ≡ |k| and (k·v 0 ). The modes with ω a (k→0) = 0 are "gapless", while the modes with ω a (k→0) = 0 are "gapped". All genuine hydrodynamic modes (sound waves, shear waves, heat diffusion) are gapless, reflecting the existence of conserved densities. On the other hand, the gapped modes, if they are present, should be viewed as parametrizations of non-hydrodynamic physics.

Hydrodynamic modes
There are d+2 hydrodynamic (gapless) modes: d−1 transverse shear modes, two sound modes, and one heat diffusion mode. Their dispersion relations for the fluid at rest (v 0 = 0) take the following form at small k, as described for example in Ref. [12]: The speed of sound is expressed in terms of the equilibrium thermodynamic quantities as which can also be written as v 2 s = (∂p/∂ǫ) S . The damping coefficient of the sound waves is determined by the viscosities and the charge conductivity. The heat diffusion coefficient is Thermodynamic inequality (2.6) implies that D is positive for positive σ. For the state with n = 0 in equilibrium, the diffusion constant is related to the conductivity by σ = (∂n/∂µ)D.
The hydrodynamic (gapless) dispersion relations in a moving fluid with v 0 = 0 can be found by applying a Lorentz boost to the spectral function F (ω, k, v 0 =0), as described for example in Ref. [8]. If a mode has a quadratic dispersion relation ω = −iDk 2 + . . . at small k in the fluid at rest, then in a moving fluid one finds at small k: The corresponding formulas for the sound mode at v 0 = 0 can be found in Ref. [8].

Stability and causality
We will call the a-th mode "stable" if Im ω a (k) 0 , and we will call the a-th mode "causal" if See Ref. [13] for a discussion of causality and large-k dispersion relations. The modes are as follows. The shear-channel fluctuations decouple from the sound-channel fluctuations as a consequence of rotation invariance, so that F (ω, k) = F shear (ω, k) d−1 F sound (ω, k). The function F shear (ω, k) is a second-order polynomial in ω, which gives rise to one gapless mode (3.1) and one gapped mode. The function F sound (ω, k) is a sixth-order polynomial in ω which gives rise to three gapless modes (3.2), (3.3) and three gapped modes. For shear-channel fluctuations, the modes of a charged fluid are identical to those of an uncharged fluid. Demanding stability and causality then gives rise to the constraint [8] (3.10) Let us now look at sound-channel fluctuations. The spectral function F sound (ω, k) is lengthy, and while one could in principle derive the constraints by applying Eqs. (3.8), (3.9) to the roots of F sound (ω, k) = 0, the constraints are unwieldy, and depend on the equation of state. The spectral function can be analyzed in various limiting cases in order to derive various necessary or sufficient conditions for stability and causality. One set of necessary conditions follows by requiring that the gaps (obtained by solving F sound (ω, k=0) = 0) have negative imaginary parts. Demanding that the sound-channel gaps are stable at v 0 = 0 gives the following necessary conditions for the stability of equilibrium: Note that T (∂ǫ/∂T ) µ/T = T (∂ǫ/∂T ) + µ(∂ǫ/∂µ) 0, thanks to the thermodynamic inequalities (2.5). For example, in a frame with ε 3 = ν 1 = 0, stability conditions (3.11), (3.12) will be satisfied for ε 1 > 0, ν 3 > 0. The constraints on the transport parameters arising from the stability of the gaps in a moving fluid are more involved. With ∆ ≡ −iω, the gaps satisfy a cubic equation a 0 ∆ 3 + a 1 ∆ 2 + a 2 ∆ + a 3 = 0, with coefficients a n that depend on the transport parameters and on the equation of state. The coefficient a 3 is proportional to ( ∂ǫ ∂T ∂n ∂µ − ∂n ∂T ∂ǫ ∂µ )(1−v 2 0 v 2 s ), hence one can always take a 3 > 0. Demanding that the sound-channel gaps are stable, i.e. Re ∆ < 0, by the Routh-Hurwitz criterion [14] then amounts to a 0 > 0, a 1 > 0, a 1 a 2 > a 0 a 3 .
Another simple set of constraints comes from the high-momentum modes. As k → ∞, the six sound-channel modes have linear dispersion relations ω = ±c s k, where c 2 s satisfies a cubic equation. The transport coefficients must be constrained by demanding c 2 s > 0 (stability) and c 2 s < 1 (causality). If we choose a frame with ε 3 = π 3 = θ 3 = 0, the cubic equation for c 2 s factorizes into a product of a linear equation and a quadratic equation: 14) The first one implies The constraints on the transport coefficients from Eq. (3.14) can be obtained as follows. For a quadratic equation ax 2 + bx + c = 0 with a > 0, the conditions that the roots are real and fall between 0 and 1 amount to the following: Applying these to Eq. (3.14) with gives a set of non-linear constraints among the coefficients ε 1,2 , π 1,2 , θ 1 , and η. Note that Eq. (3.14) is exactly the same equation that determines the propagation speed of the large-k eigenmodes in an uncharged fluid [8]. Thus in the frame with ε 3 = π 3 = θ 3 = 0, the causality constraints (3.16), (3.17) on the coefficients ε 1,2 , π 1,2 , θ 1 , and η will be exactly the same as in uncharged fluids. In order to express the causality constraints in terms of physical transport coefficients, we need the bulk viscosity (2.3) and the conductivity (2.4) in the frame with ε 3 = π 3 = θ 3 = 0. From Eq. (3.15) one immediately finds For example, if we choose a frame in which γ 1 = nθ 1 /(ǫ+p), then γ 3 = −σT , and the causality constraint (3.15) becomes simply ν 3 > σT . Similarly, in order to simplify the bulk viscosity one could further choose a frame in which ν 2 = (p ,ǫ )ν 1 + (p ,n )ν 3 /T . Then the transport parameters ν i drop out from the bulk viscosity, and the large-k causality constraints (3.16), (3.17) take exactly the same form as the corresponding constraints in uncharged fluids [8].

Stable and causal frames for conformal fluids
The discussion above simplifies for conformal fluids, where the one-derivative transport parameters satisfy the constraints of Eq. (2.7) due to conformal symmetry. Let us further choose a frame in which π 3 , θ 3 , ν 1 , and γ 1 all vanish, leaving one with five independent transport parameters π 2 , θ 1 , η, ν 3 , and γ 3 . The five parameters can be thought of as two genuine transport coefficients η and σ, and three "relaxation times" corresponding to the relaxation of the energy density (π 1 ), momentum density (θ 1 ), and charge density (ν 3 ). As k → 0, the stability conditions (3.11), (3.12) for the gapped modes become ν 3 π 2 > 0, and T ν 3 +dλπ 2 > 0 where λ ≡ T 2 (∂n/∂µ)/(ǫ+p) > 0 is the dimensionless charge susceptibility. These are satisfied for ν 3 > 0, π 2 > 0. The large-k constraint (3.18) becomes where κ ≡ nT /(ǫ+p) is a dimensionless measure of the equilibrium charge density. Note that the thermodynamic inequalities (2.5) imply λ dκ 2 . In fact, the only information about the equation of state that is relevant for the linearized analysis of stability and causality is contained in κ and λ. The large-k constraints (3.20) -(3.23) reduce simply to Restricting to d = 3 space dimensions, the inequalities (3.25) become exactly the same constraints found earlier for uncharged conformal fluids in Refs. [9,8]. In order to satisfy them, it is sufficient to demand π 2 > 4 3 η, θ 1 > 4η. It is a straightforward exercise to check that the large-k constraints (3.24), (3.25) in d = 3 also ensure that the sound-channel gaps are stable at all v 2 0 < 1, as long as the equation of state obeys the standard thermodynamic inequalities (2.5). In other words, conditions (3.24), (3.25) ensure that in our chosen frame all small-k modes are stable, and all large-k modes are stable and causal. Finally, we point out that the large-k causality alone does not guarantee stability at all k: for example, the causality conditions (3.9) allow for negaive π 2 , which is ruled out by the small-k stability conditions.

Real-space causality
So far, we have looked at the linearized stability of the equilibrium state, and the linearized causality of the near-equilibrium perturbations in momentum space. In fact, similar to what was done in Refs. [9,10], one can study the non-linear causality of charged first-order hydrodynamics in real space using the quasi-linear character of the hydrodynamic equations. The hydrodynamic conservation laws are partial differential equations for the variables T , µ, and u α which satisfies u α u α = −1. Rather than working with the vector u α which is constrained by u α u α = −1, we find it more convenient to work with the unconstrained vector β α = βu α , where β ≡ (−β α β α ) 1/2 > 0 is the inverse temperature, β = 1/T . The hydrodynamic equations (1.1) are second-order quasilinear partial differential equations for U A = (β α , µ/T ) that can schematically be written as where the indices A, B, C range from 1 to d+2 (again, d is the number of spatial dimensions), and the coefficients M µν , N µν , P µ depend on U , but not on the derivatives of U . We thus have d+2 differential equations for d+2 unconstrained variables. 3 The hyperbolicity of the equations and the causality of the solutions are determined by the principal part (M µν ) AB , see e.g. Ref. [15], Ch. VI. The characteristic surfaces φ(x) = const. are found from where the vectors ξ µ ≡ ∂ µ φ(x) are normal to the characteristic surfaces at the point x. For the hydrodynamic equations (4.1) to be hyperbolic, Eq. (4.2) must only have non-zero real solutions ξ 0 = ξ 0 (ξ i ). For the equations to be causal, the surfaces swept out by these normals must lie either outside or on the lightcone ξ µ ξ µ = 0.
3 Ref. [10] chooses to work with u α which is not constrained by uαu α = −1, and instead considers the projected conservation laws u β ∇αT αβ = 0, ∆ µ β ∇αT αβ = 0 as independent equations. Even though the norm of u α is not preserved under time evolution in this approach, the causality can be argued by noting that the causal structure of the projected equations is the same as the causal structure of the original equations ∇αT αβ = 0. We find it conceptually cleaner to work with β µ = u µ /T whose norm is not fixed. The causality constraints one finds by projecting the energy-momentum conservation laws are the same as those in Sec. 3. Let us work in a frame with ε 3 = π 3 = θ 3 = 0. The conservation laws can be written as as well as The first factor in Eq. (4.5) gives At a given point in spacetime, passing to a local coordinate system in which u α = (1, 0) at that point, the solutions are ξ 0 = ± −γ 3 /ν 3 |ξ i |. Demanding that these are real and lie outside the lightcone gives the same constraint (3.15) we found earlier from the linearized analysis in momentum space. In d+1 dimensions, the determinant in Eq. (4.5) can be computed by using the following easily derived identity: (4.7) Applying this to det (M ρσ ) α β ξ ρ ξ σ , we find B = θ 1 (ξ·u) 2 − η(ξ·∆·ξ), hence there are characteristic surfaces determined by The conditions of hyperbolicity and causality will be satisfied for (ξ·u) = c s (ξ·∆·ξ), with c s real and 0 < c 2 s < 1. This gives exactly the same equation (3.14) for c s , and consequently the same causality conditions (3.16), (3.17) we found earlier from the linearized analysis in momentum space (recall that the formulas in Sec. 3 are written in a thermodynamic frame with θ 2 = θ 1 ). In other words, the real-space analysis of hyperbolicity and causality gives the same constraints on the transport parameters as the linearized analysis of Sec. 3.
As in Refs. [9,10], it is straightforward to couple Eq. (4.1) to Einstein's equations. The latter have no terms with second derivatives of β α or µ/T . Thus, extending the variables to U A = (β α , µ/T, g µν ), the part of the determinant (4.2) describing the metric degrees of freedom decouples, and the conditions of Sec. 3 give rise to causal Navier-Stokes equations coupled to dynamical gravity.

Conclusions
We have proposed a class of stable and causal frames for relativistic hydrodynamics with conservation laws given by Eqs. (1.1), and whose constitutive relations contain up to one derivative of the hydrodynamic variables T , u α , and µ. The stable and causal frames generalize the frames proposed in BDNK [9,8,10] to fluids with a global U (1) charge such as the baryon number. A choice of frame ultimately amounts to a convention specifying how one chooses to fix the arbitrariness of defining T , u α , and µ beyond the perfect-fluid approximation. In a causal frame, the Navier-Stokes equations are causal both in flat and in curved space.
The entropy production for hydrodynamics in the stable and causal frames is exactly the same as the entropy production in the classic Eckart or Landau-Lifshitz frames [1]. This is because in first-order hydrodynamics the divergence of the entropy current has to be evaluated on-shell (when the hydrodynamic equations are satisfied), and in the derivative expansion. Just like in Ref. [8], the divergence of the entropy current on-shell and to first order in derivatives is the same as in [1], and is only determined by η, ζ, and σ.
The procedure that gives rise to the stable and causal Navier-Stokes equations is quite straightforward, and embodies the spirit of effective field theory. In the standard quantum field theory, the effective description is constructed by writing down the action in terms of all possible operators consistent with the symmetry, up to a given dimension, and then constraining the coefficients of these operators based on the stability of the vacuum and unitarity. Similarly, in hydrodynamics we write down all possible terms in the constitutive relations up to a given derivative order, and then constrain the coefficients of these terms based on the stability of equilibrium and causality.
The most general one-derivative constitutive relations are given by Eqs. (2.1) and (2.2). Using a frame in which ε 3 = π 3 = θ 3 = 0, the hyperbolicity and causality of the equations are easily demonstrated, provided the remaining coefficients obey the inequalities discussed in Sec. 3. One can further choose a frame in which the causality constraints look exactly like the constraints in uncharged fluids, plus a lower bound on the coefficient ν 3 , see Eqs. (3.19)-(3.23).
We have not performed an exhaustive analysis of stability. While the stability of the gaps requires the inequalities (3.11) to hold, a full study of stability would require working with a specific equation of state. Performing such a study would be straightforward. In conformal theories, it appears that the conditions of stability and causality at large k combined with the stability at small k also ensure stability at all k, though at the moment we do not have proof of that. We hope that our observations will stimulate further work on stable and causal relativistic Navier-Stokes equations, including their applications to heavy-ion physics and astrophysics.
Consider now the energy-momentum tensor and the current (2.1). Following Ref. [12], we write E = ǫ(T, µ) + f E , N = n(T, µ) + f N , where f E and f N are the derivative corrections, as well as The hydrodynamic variables T , u α , µ are related to T L , u α L , µ L by where δT , δu α , δµ are O(∂). Demanding that T αβ and J α are frame-independent, one finds [12] δu α = Q α ǫ + p , δT = f E One can use the above expressions in order to find T , u α , µ in the general frame, if T µν and J µ happen to be known. Given only T µν (x, t=t 0 ) and J µ (x, t=t 0 ), one can not tell whether these single-time values correspond to a physical state that is describable by hydrodynamics. However, if one has access to T µν (x, t) and J µ (x, t) for a range of times around t 0 , then the time derivatives of the Landau-Lifshitz variables can be evaluated, the importance of the derivative corrections at t = t 0 can be estimated, and relations such as (A.8), (A.9) can be used in order to find the general-frame T , u α , µ in that state, within the derivative expansion. See also Ref. [9], sec. VII C for related comments.