(2) PHYSICAL REVIEW E 84, 026316 (2011). Relativistic viscoelastic fluid mechanics Masafumi Fukuma* and Yuho Sakatani† Department of Physics, Kyoto University, Kyoto 606-8502, Japan (Received 15 April 2011; published 15 August 2011) A detailed study is carried out for the relativistic theory of viscoelasticity which was recently constructed on the basis of Onsager’s linear nonequilibrium thermodynamics. After rederiving the theory using a local argument with the entropy current, we show that this theory universally reduces to the standard relativistic Navier-Stokes fluid mechanics in the long time limit. Since effects of elasticity are taken into account, the dynamics at short time scales is modified from that given by the Navier-Stokes equations, so that acausal problems intrinsic to relativistic Navier-Stokes fluids are significantly remedied. We in particular show that the wave equations for the propagation of disturbance around a hydrostatic equilibrium in Minkowski space-time become symmetric hyperbolic for some range of parameters, so that the model is free of acausality problems. This observation suggests that the relativistic viscoelastic model with such parameters can be regarded as a causal completion of relativistic Navier-Stokes fluid mechanics. By adjusting parameters to various values, this theory can treat a wide variety of materials including elastic materials, Maxwell materials, Kelvin-Voigt materials, and (a nonlinearly generalized version of) simplified Israel-Stewart fluids, and thus we expect the theory to be the most universal description of single-component relativistic continuum materials. We also show that the presence of strains and the corresponding change in temperature are naturally unified through the Tolman law in a generally covariant description of continuum mechanics. DOI: 10.1103/PhysRevE.84.026316. PACS number(s): 47.75.+f, 11.25.Tq, 24.10.Nz, 46.35.+z. I. INTRODUCTION. The dynamics of fluids at large scales is universally described by the Navier-Stokes equations, which represent the regression to a global equilibrium with transfers of conserved quantities (such as energy-momentum and particle number) among fluid particles [1]. This can be formulated in a generally covariant way, but it is known that there arises a problem of acausality. In fact, the obtained equations for the propagation of disturbance are basically parabolic and thus predict infinitely large speed of propagation for infinitely high frequency modes, leaving light cones. One should note here that this does not imply the breakdown of the internal consistency of the description because the Navier-Stokes equations are simply an effective description at large space-time scales and need not describe high-frequency modes correctly. However, this is still troublesome when adopting the equations in numerical simulations; the initial value problems are ill posed, and unacceptable numerical solutions can be obtained easily. To remedy the problem, M¨uller, Israel, and Stewart [2–4] extended the theory by treating the dissipative part of stress μν tensor, τ(d) , and the heat flux q μ (for the Eckart frame) or the particle diffusion current ν μ (for the Landau-Lifshitz frame) as additional thermodynamic variables on which the entropy density can depend. This prescription is based on the so-called extended thermodynamics and corresponds to taking into account higher derivative corrections to the effective theory. It has been shown that such modified theories have a good causal behavior and that linear perturbations around a hydrostatic equilibrium obey hyperbolic differential equations. This is now regarded as a fundamental. * †. fukuma@gauge.scphys.kyoto-u.ac.jp yuho@gauge.scphys.kyoto-u.ac.jp. 1539-3755/2011/84(2)/026316(25). framework for the numerical study of relativistic viscous fluids. Meanwhile, modifications of the Navier-Stokes equations have also been studied in the area of rheology, and the materials treated there are generically called viscoelastic materials or viscoelastic fluids. Historically, viscoelasticity was defined by Maxwell in the 19th century as the characteristic property of such continuum materials that behave as elastic solids at short time scales and as viscous fluids at long time scales [1,5]. In 1948, Eckart proposed a theory of elasticity and anelasticity [6], which describes the nonrelativistic dynamics of singlecomponent viscoelastic materials and was reinvented recently [7] in the light of the covariance under foliation preserving diffeomorphisms. In this description, elastic strains (or equivalently, the “intrinsic metric” defined below) are introduced as additional thermodynamic variables, as in the theory of elasticity. As explicitly shown in [7], this theory of viscoelasticity contains the theory of elasticity and the theory of fluids as special limiting cases, and correctly reproduces the NavierStokes equations in the fluid limit. Furthermore, as was pointed out in [8], since the dynamics at short time scales is dominated by elasticity, shear modes of linear perturbations around a hydrostatic equilibrium obey differential equations with second-order time derivatives (in contrast to the equations obtained from the Navier-Stokes equations that contain only a first-order time derivative), so that causal behaviors for large frequencies are significantly improved. Recently, on the basis of Onsager’s linear regression theory on nonequilibrium thermodynamics [9–12], the present authors proposed a relativistic theory of viscoelasticity [13] which generalizes the theory of elasticity and anelasticity [6,7] in a generally covariant form. In the present paper, after rederiving the theory relying on a local argument with the entropy current, we study the detailed properties of relativistic viscoelasticity. We show that fluidity is universally realized in the long time limit and also that acausal problems disappear. 026316-1. ©2011 American Physical Society.

(3) MASAFUMI FUKUMA AND YUHO SAKATANI. PHYSICAL REVIEW E 84, 026316 (2011). for a wide region of parameters. Thus, the relativistic theory of viscoelasticity with such parameters can be regarded as a causal completion of relativistic Navier-Stokes fluid mechanics, and we expect that it could be used as another basis in the numerical study of relativistic viscous fluids. This paper is organized as follows. In Sec. II we rederive the viscoelastic model of [13] using a local argument with the entropy current. We also show that the presence of strains and the corresponding change in temperature are naturally unified in a generally covariant description of continuum mechanics. In Sec. III we consider the long and short time limits of our viscoelastic model. We prove that the model universally gives relativistic Navier-Stokes fluids in the long time limit. In Sec. IV we show that when some parameters take specific values, our viscoelastic model reduces to (a higher-dimensional extension of) the nonlinear generalization of the simplified Israel-Stewart model [14]. In Sec. V we consider linear perturbations around a hydrostatic equilibrium in Minkowski space-time. The dispersion relations show that the evolutions are certainly stable. Although the wave equations for the linear perturbations are not always hyperbolic, if some parameters are chosen appropriately (including the parametrizations for the simplified Israel-Stewart model) they become symmetric hyperbolic and thus free of acausality problems. Section VI is devoted to conclusion and discussions. II. RELATIVISTIC VISCOELASTIC MECHANICS. In this section, we rederive the fundamental equations for relativistic viscoelastic mechanics using a local argument with the entropy current. In Appendix A we show that the present formulation is equivalent to the “entropic formulation” proposed in our previous paper [13] which is based on Onsager’s linear regression theory. A. Definitions. We start by giving a brief review on the generally covariant definitions of viscoelastic materials [13]. 1. Geometrical setup. We consider a single-component continuum material living in a (D + 1)-dimensional Lorentzian manifold M. The local coordinates are denoted by x μ (μ = 0,1, . . . ,D), and the background Lorentzian metric with signature (−, + , . . . ,+) by gμν (x). Following the convention of Landau and Lifshitz [1], we define the velocity field u = uμ (x)∂μ from the momentum (D + 1) vector pμ as uμ (x) ≡ g μν (x)pν (x)/e(x) = pμ (x)/e(x),. (1). where e(x) ≡ −g μν (x)pμ (x)pν (x) is the proper energy density. Note that uμ (x) is normalized as gμν (x)uμ (x)uν (x) = −1. Here and hereafter indices are subscripted (or superscripted) always with gμν (or with its inverse g μν ). Assuming that the velocity field is hypersurface orthogonal, we introduce a foliation of M consisting of spatial hypersurfaces (time slices) orthogonal to uμ . We parametrize the time slices with a real parameter t and denote them by t . We exclusively (except for Sec. V) use a coordinate system x = (x μ ) = (x 0 ,x) such that x 0 = t, x = (x i ) (i = 1, . . . ,D), for which the shape of the material at time t is given by the induced metric on t : hμν (x) ≡ gμν (x) + uμ (x)uν (x) pμ (x)pν (x) . = gμν (x) + e2 (x). (2). We also define the extrinsic curvature Kμν of the hypersurface as half the Lie derivative of hμν with respect to the velocity field u = uμ ∂μ : Kμν ≡ 12 £u hμν = 12 hρμ hσν (∇ρ uσ + ∇σ uρ ).. (3). This measures the rate of change in the induced metric hμν as material particles flow along uμ . Note that this tensor is symmetric and orthogonal to uμ , Kμν uν = 0. In the ADM parametrization, the metric and the velocity are represented with the lapse N (x) and the shifts N i (x) (i = 1, . . . ,D) as ds 2 = gμν (x)dx μ dx ν = −N 2 (x)dt 2 + hij (x)[dx i − N i (x)dt][dx j − N j (x)dt], (4) 1 N i (x) ∂0 + ∂i N (x) N (x) [⇔ uμ (x)dx μ = −N (x)dt]. u = uμ (x)∂μ =. (5). The volume element √ on the hypersurface √ is given by the D form hd D x ≡ det(hij )d D x = N −1 −gd D x. With a given foliation, we still have the symmetry of foliation preserving diffeomorphisms that give rise to transformations only among the points on each time slice. Using this residual gauge symmetry we can impose the synchronized gauge, N i (x) ≡ 0, so that the background metric and the velocity field are expressed as ds 2 = gμν (x)dx μ dx ν ≡ −N 2 (x)dt 2 + hij (x)dx i dx j , (6) ∂ 1 ∂ = , (7) u = uμ (x)∂μ = N (x) ∂t ∂τ where τ is the local proper time defined by dτ = N dt. In this gauge, due to the relation ∂/∂t = N (x)∂/∂τ , the proper. FIG. 1. Processes of deformation and stress relaxation [7,8]. 026316-2.

(4) RELATIVISTIC VISCOELASTIC FLUID MECHANICS. PHYSICAL REVIEW E 84, 026316 (2011). energy density e(x) measured with the proper time τ is related to the energy density e(x) measured with time t as e(x) = N (x)e(x).. (8). Note that e(x) includes the gravitational potential through the factor N(x). Accordingly, the local temperature T measured with τ is related to the temperature T measured with t through the following Tolman law: T(x) = N (x)T (x).. (9). 2. Definition of (relativistic) viscoelastic materials. According to the definition of Maxwell, viscoelastic materials behave as elastic solids at short time scales and as viscous fluids at long time scales (see, e.g., Sec. 36 in [5]). In order to understand how such materials evolve in time, we consider a material consisting of many molecules bonding each other and assume that the molecules first stay at their equilibrium positions in the absence of strains (as in the leftmost illustration of Fig. 1) [7,8]. We now suppose that an external force is applied to deform the material. An internal strain is then produced in the body, and according to the definition, the accompanied internal stress can be treated as an elastic force at least during short intervals of time. However, if we keep the deformation much longer than the relaxation times (characteristic to each material), then the bonding structure changes to maximize the entropy, and the internal strain vanishes eventually as in the rightmost of Fig. 1. The point is that two figures (the central and the rightmost) have the same shape (same induced metric) hμν , but different bonding structures. The internal bonding structure can be specified by the intrinsic metric h¯ μν , which measures the shape that the material would take when all the internal strains are removed virtually [6,7]. For the example given in Fig. 1, the intrinsic metric for the center illustration is given by the induced metric for the leftmost illustration, while the intrinsic metric for the rightmost illustration agrees with the induced metric for itself. Thus, the plastic (i.e., nonelastic) deformation from the center illustration to the rightmost illustration is described as the evolution of the intrinsic metric.1 Its generally covariant generalization can be defined in the following way. Suppose that we have two adjacent, spatially separated space-time points P and Q, each of which represents a point on the trajectory of a material particle (see Fig. 2). By denoting their coordinates by x = (x μ ) and x + dx = (x μ + dx μ ), respectively, the distance between P and Q in the real configuration is of course given with the metric gμν as (the square root of) ds 2 = gμν (x)dx μ dx ν .. (10). We now virtually remove all the strains in a sufficiently small space-time region including the two points. Then P and Q. FIG. 2. Real [x μ (τ )] and virtual [x¯ μ (τ )] trajectories of material particles. The distance between P¯ and Q¯ gives the definition of the intrinsic metric g¯ μν .. ¯ whose coordinates we would move to other positions P¯ and Q, denote by x¯ = (x¯ μ ) and x¯ + d x¯ = (x¯ μ + d x¯ μ ), respectively. ¯ This correspondence defines a local map f : x → x¯ = x(x), with which we define the intrinsic metric g¯ μν (x) as the metric ¯ (or, as the measuring the virtual distance between P¯ and Q ∗ pullback of the metric gμν for the map; g¯ μν ≡ f gμν ):2 ¯ x¯ ρ d x¯ σ = gρσ [x(x)] ¯ d s¯ 2 ≡ gρσ (x)d. ∂ x¯ ρ ∂ x¯ σ μ ν dx dx ∂x μ ∂x ν. ≡ g¯ μν (x)dx μ dx ν .. (11). With the velocity vector u = u (x)∂μ , we parametrize g¯ μν as μ. g¯ μν = −(1 + 2θ )uμ uν − εμ uν − εν uμ + h¯ μν (εμ uμ = 0, hμν uν = 0, h¯ μν uν = 0).. (12). The strain tensor is then introduced as Eμν (x) ≡ 12 [gμν (x) − g¯ μν (x)] = θ uμ uν + 12 (εμ uν + εν uμ ) + εμν ,. (13). where εμν (x) ≡ 12 [hμν (x) − h¯ μν (x)]. (14). is the spatial strain tensor. Note that if we define the extrinsic curvature associated with the spatial intrinsic metric h¯ μν as K¯ μν ≡ 12 £u h¯ μν = 12 (uλ ∂λ h¯ μν + ∂μ uλ h¯ λν + ∂ν uλ h¯ μλ ),. (15). the following identity holds: £u εμν = Kμν − K¯ μν .. (16). A viscoelastic material is a thermodynamic system consisting of material particles as its subsystems. While the system regresses to a thermodynamic equilibrium, one can imagine that the virtual trajectory of each material particle approaches its real trajectory, so that the strain tensor Eμν approaches zero. Such an irreversible process is called plastic (i.e., nonelastic), and thus we see that the dynamics of Eμν includes plastic evolutions (in addition to reversible, elastic evolutions). In the following discussions, we assume that Eμν = (εμν ,εμ ,θ ) are all small quantities, such that their nonlinear effects can be neglected. We shall denote. 1¯. hμν is also called the “strain metric” and was first introduced by Eckart to embody “the principle of relaxability-in-the-small” in anelasticity [6]. Some examples of the explicit form of hμν and h¯ μν under various deformations can be found in [7,8].. 2 As in the standard theory of elasticity [5], there may be an arbitrariness in defining x¯ μ , but the intrinsic metric g¯ μν can still be defined uniquely.. 026316-3.

(5) MASAFUMI FUKUMA AND YUHO SAKATANI. PHYSICAL REVIEW E 84, 026316 (2011). the contraction of a spatial tensor3 Aμν with g μν by trA, so that trε ≡ g μν εμν = hμν εμν , trK ≡ g μν Kμν = hμν Kμν .. (17) (18). We close this section by explaining the physical meaning of the strain tensor Eμν = (εμν ,εμ ,θ ). The spatial strain tensor εμν stands for the standard strains, measuring the difference between the induced metric hμν and the spatial induced metric h¯ μν . One can easily see that the quantity εμ represents the relative velocity of a material particle in its real trajectory with respect to that in its virtual trajectory, εμ = uμ − u¯ μ ≡ dx μ /dτ − d x¯ μ /dτ , where τ is a common proper time (see Fig. 2). In order to understand the meaning of θ , we first recall that the covariant vector uμ is expressed as uμ dx μ = −N dx 0 . We can then rewrite ds 2 and d s¯ 2 as ds 2 = −N 2 (x)(dx 0 )2 + hμν (x) dx μ dx ν ,. (19). d s¯ = −[1 + 2θ (x)]N (x)(dx ) + 2N(x)εμ (x) dx dx 0 + h¯ μν (x) dx μ dx ν , (20) 2. 2. 0 2. μ. with hμν dx μ dx ν = hij (dx i − N i dx 0 )(dx j − N j dx 0 ) and a similar (but a bit more complicated) expression for μ 0 μ ν ¯ ¯ 2Nε √ μ dx dx + hμν dx dx . These equations mean that N ≡ 1 + 2θ N (1 + θ )N represents the lapse function for the intrinsic metric. Then, through the Tolman law, we can relate the virtual temperature T¯ observed in the absence of strains to the actual temperature T as N T = N¯ T¯ (= T). We thus obtain the relation θ = (N¯ 2 /N 2 − 1)/2 = (T 2 /T¯ 2 − 1)/2 (T − T¯ )/T¯ , and conclude that the scalar θ expresses the increase of the temperature due to strains. This conclusion shows that the presence of strains and the corresponding change in temperature are naturally unified in a generally covariant description of continuum mechanics. B. Entropy production rate. As was adopted in [13], in order to develop thermodynamics in a generally covariant manner, it is convenient to distinguish density quantities from other intensive quantities, √ and, by multiplying them with the spatial volume element h, we construct new quantities which are spatial densities on each time slice. For example, the entropy density s, the energymomentum density pμ , and the number density n are density quantities, and for them we introduce the following spatial densities: √ √ √ s˜ ≡ hs, p˜ μ ≡ hpμ , n˜ ≡ hn. (21) We assume that each material particle is in its local thermodynamic equilibrium, and that the local entropy s˜ is a function of ˜ and gμν as well as of the strain tensor Eμν = (εμν ,εμ ,θ ): p˜ μ , n, ˜ s˜ (x) = s˜ (Eμν (x),p˜ μ (x),n(x),g μν (x)).. (22). We further assume that s˜ depends on p˜ μ only through the local ˜ p˜ μ ,gμν ) ≡ −g μν p˜ μ p˜ ν , so that s˜ can also be proper energy e( expressed as ˜ ˜ s˜ (x) = σ˜ (Eμν (x),e(x), n(x),g μν (x)) ˜ p˜ μ (x),gμν (x)),n(x),g ˜ = σ˜ (εμν (x),εμ (x),θ (x),e( μν (x)). (23) Since we are only interested in linear nonequilibrium thermodynamics, we only need to expand s˜ in Eμν to second order:4 1 [2λ1 ε μν ε μν s˜ = (terms independent of Eμν ) − 2T + λ2 εμ εμ + γ1 (trε)2 + 2γ2 (trε)θ + γ3 θ 2 ]. (24) We require the stability of the system under the change in strains Eμν , so that the constants λ1 and λ2 are non-negative, and the matrix γ = (γγ12 γγ23 ) is positive semidefinite. Then the fundamental thermodynamic relation can be written as √ uν h μν μ δ s˜ = − δ p˜ ν − δ n˜ + T δgμν T T 2T (q) √ √ h h. μν 2λ1 ε δε μν − (γ1 trε + γ2 θ )δ(trε) − T T √ √ h h λ2 εμ δεμ − (γ3 θ + γ2 trε)δθ. − (25) T T Here the temperature T , the chemical potential μ and the μν quasiconservative part of the stress tensor, τ(q) , are defined as5 √ 1 ∂ σ˜ μ ∂ σ˜ ∂ σ˜ h μν = , =− , τ , = (26) ∂ e˜ T ∂ n˜ T ∂gμν 2T (q) μν. μν μν T(q) ≡ euμ uν + τ(q) .. (27). In deriving Eq. (25), we have used the relations ˜ p˜ μ ,gμν ) ∂ e( = −uν , ∂ p˜ ν. ˜ p˜ μ ,gμν ) ∂ e( e˜ = uμ uν . ∂gμν 2. (28). We now set the variation in Eq. (25) to be δ = £u . We then obtain √ h∇μ (suμ ) ν √ u μ 1 μν = h − ∇μ (pν uμ ) − ∇μ (nuμ ) + τ(q) Kμν T T T 2λ1 μν 1 ε £u ε μν − (γ1 trε + γ2 θ )£u (trε) − T T 1 λ2 (29) − εμ £u εμ − (γ3 θ + γ2 trε)£u θ . T T. For a tensor Aμν , we define A μν ≡ (1/2)hρμ hσν [Aρσ + Aσρ − (2/D)hαβ Aαβ hρσ ]. 5 We here use a convention that the quasiconservative stress tensor μν does not include stresses originated from strains. τ(q) 4. 3 By spatial we mean that Aμν is orthogonal to uμ , Aμν uν = 0 = Aμν uμ . Recall that g μν = −uμ uν + hμν .. μν. where we require that τ(q) be orthogonal to uμ , τ(q) uν = 0. The quasiconservative part of the energy-momentum tensor is then defined as. 026316-4.

(6) RELATIVISTIC VISCOELASTIC FLUID MECHANICS. PHYSICAL REVIEW E 84, 026316 (2011). Here we have used the identities for Lie derivatives: √ √ £u s˜ = h∇μ (suμ ), £u p˜ ν = h∇μ (pν uμ ), (30) √ £u n˜ = h∇μ (nuμ ), (31) √ √ which can be shown by using the identities £u h = h∇μ uμ and pμ ∇ν uμ = 0. Note that tr(£u εμν ) = hμν £u εμν can be replaced by £u (trε) in our approximation because the difference £u (trε) − tr(£u εμν ) = (£u h¯ μν )εμν = −2K μν εμν is of higher orders. The full energy-momentum tensor T μν and the full number current nμ are given by T μν ≡ euμ uν + τ μν ,. nμ ≡ nuμ + ν μ ,. (τ μν uν = 0 = ν μ uμ ),. (32). where τ μν and ν μ are the stress tensor and the diffusion current, respectively. Then, by introducing the entropy current μ (33) s μ ≡ suμ − ν μ , T and by using Eq. (29) together with the current conservation laws ∇ν T μν = 0, ∇μ nμ = 0,. (34). the local entropy production rate can be evaluated as μ 1 μν μν μ μ ∇μ s = − τ − τ(q) Kμν + ν ∂μ − T T 2λ1 μν 1 ε £u εμν − (γ1 trε + γ2 θ )£u (trε) T T 1 λ2 − εμ £u εμ − (γ3 θ + γ2 trε)£u θ T T. 2λ1 £ ε − 1 u. μν T = ε μν − K μν (q) T τ μν − τ μν. λ2 − T £u εμ μ μ μ + ε ∇ − T νμ ⎞ ⎛ £u (trε) 1 1 − γ + tr ε θ − trK ⎝ T £u θ ⎠ . T 1 (trτ − trτ(q) ) D. Here G, H, and K are antisymmetric matrices, 0 G 0 H G= , H= , −G 0 −H 0 ⎞ ⎛ 0 K

(7) K ⎟ ⎜ K = ⎝−K

(8) 0 −Ka ⎠ , −K Ka 0. (39). (40). and η, σ , and ζ are positive semidefinite symmetric matrices, η1 η2 σ1 σ2 η= , σ = , (41) η2 η3 σ2 σ3 ⎛ ⎞ ζ 1 ζ2 ζ4 ⎜ ⎟ (42) ζ = ⎝ζ2 ζ3 ζ5 ⎠ . ζ4 ζ5 ζ6 Note that only the symmetric matrices contribute when substituted to the entropy production rate (35). This means that the matrices η, σ , and ζ are associated with irreversible processes, while the matrices G, H, and K are with reversible ones. The relationship between the equations given above and the corresponding ones given in [13] is summarized in Appendix A. C. Fundamental equations. Using Eqs. (36)–(42) at each point x = (x 0 = t,x) on time slice t , we can express (A) the currents τ μν and ν μ and (B) the evolution of strains, £u ε μν , £u εμ , £u trε, and £u θ , only in terms of local thermodynamic quantities on t . We thus conclude that the dynamics of relativistic viscoelastic materials is described by the following two sets of equations [7,13]: (A) Current conservation laws:. −. ∇μ T μν = ∇μ (euμ uν + τ μν ) = 0, ∇μ nμ = ∇μ (nuμ + ν μ ) = 0,. (43) (44). with the constitutive equations (35). 2η3 μν K T ζ6 − (K − ζ4 )trε − (Ka + ζ5 )θ + trK hμν , T μ . ν μ = −(H − σ2 )εμ + σ3 hμν ∂ν − T. τ μν = τ(q) − 2(G − η2 )ε μν − μν. Thus, if we require that each term be separately positive definite, we obtain the following equations: 2λ. − T 1 £u ε μν ε μν = 2(G + η) , (36) (q) − T1 K μν τ μν − τ μν. εμ − λT2 £u εμ , (37) = (H + σ ) hνμ ∂ν − Tμ νμ ⎞ ⎛ ⎛ ⎞ trε £u (trε) 1 γ − ⎝ T £u θ ⎠ = (K + ζ ) ⎝ θ ⎠ . (38) 1 1 (trτ − trτ ) (q) trK − D T. (45) (46). (B) Rheology equations:. 026316-5. £u ε μν = −. η1 T G + η2 ε μν + K μν , λ1 λ1. σ1 T (H + σ2 )T ν μ £u εμ = − εμ − h μ ∂ν − , λ2 λ2 T. (47) (48).

(9) MASAFUMI FUKUMA AND YUHO SAKATANI. £u (trε) £u θ. PHYSICAL REVIEW E 84, 026316 (2011). −1 −ζ1 T trε − (K

(10) + ζ2 )T θ + (K + ζ4 )trK γ1 γ2 (K

(11) − ζ2 )T trε − ζ3 T θ − (Ka − ζ5 )trK γ2 γ3 ⎞ ⎛

(12)

(13) )+γ2 (Ka−ζ5 ) 2 (K −ζ2 )]T 3 (K +ζ2 )]T − [γ3 ζ1 +γdet trε + [γ2 ζ3 −γdet θ + γ3 (K+ζ4det trK γ γ γ ⎠. = ⎝ [γ ζ +γ (K

(14) −ζ )]T [γ1 ζ3 −γ2 (K

(15) +ζ2 )]T γ2 (K+ζ4 )+γ1 (Ka−ζ5 ) 2 1 1 2 trε − θ − trK det γ det γ det γ. =. The former set of equations describes the dynamics of D + 2 conserved quantities (pμ = euμ ,n), while the latter that of D(D + 1)/2 dynamical variables Eμν = (εμν ,εμ ,θ ). It is convenient to introduce the following parameters: τs ≡. λ1 , η1 T. τ± ≡. , T Pζ γ ∓ Pζ2γ − 4 det γ (det ζ s + K

(16) 2 ). a± ≡. τσ ≡. λ2 , σ1 T 2 det γ. (50) (51). −2[ζ3 γ2 − (K

(17) + ζ2 )γ3 ] , ζ3 γ1 − ζ1 γ3 − 2K

(18) γ2 ± Pζ2γ − 4 det γ (det ζ s + K

(19) 2 ) (52). where Pζ γ ≡ ζ3 γ1 + ζ1 γ3 − 2ζ2 γ2 0, and ζ s is the principal submatrix of ζ defined by ζ s ≡ ( ζζ12 ζζ23 ). Since ζ s is positive semidefinite, det ζ s is non-negative. Note that τs , τσ , and Reτ± are all non-negative. We further introduce the scalar variables ε± ≡ 12 (trε − a± θ ).. (53). Then the rheology equations (47)–(49) can be rewritten in a more compact form. (B

(20) ) Rheology equations: 1 G + η2 £u ε μν = − ε μν + K μν , (54) τs λ1 μ 1 (H + σ2 )T ν , (55) £u εμ = − εμ − h μ ∂ν − τσ λ2 T (Ka − ζ5 )(a± γ1 + γ2 ) + (K + ζ4 )(a± γ2 + γ3 ) £u ε± = trK 2 det γ 1 − ε± . (56) τ± From these, we see that τs , τσ , and Reτ± give the typical time scales for the relaxation of strains. The relation between the viscoelastic models and a few well-known rheological models (such as the Maxwell model and the Kelvin-Voigt model) is discussed in Appendix B.. (49). III. FLUID AND ELASTIC LIMITS. In this section, we discuss the limits of elasticity and fluidity in the relativistic theory of viscoelasticity. We first identify the properties that characterize a given material as a fluid or as an elastic material. We then consider the long-time and short-time limits of our dynamical equations and show that fluidity is universally realized in the long time limit. We also make a comment on the subtlety existing in Maxwell’s definition of viscoelasticity. A. Fluidity and elasticity. Fluidity is characterized by the property that the relaxation of the strains Eμν = (εμν ,εμ ,θ ) proceeds instantaneously. Thus, their rheology equations are expressed as £u εμν = 0,. £u εμ = 0,. £u θ = 0, (fluids). (57). £u ε± = 0.. (58). or equivalently, £u ε μν = 0, £u εμ = 0,. (fluids). This situation can also be realized in the long time limit, and we show in the next section that the constitutive equations for our viscoelastic model universally reduces to those for the Navier-Stokes fluids in the long time limit. On the other hand, elastic materials by definition do not undergo any plastic deformations, and thus their intrinsic metric h¯ μν does not evolve for any processes. Thus, a given viscoelastic material is regarded as being elastic when its rheology equations are expressed as [6,7,15] K¯ μν = 12 £u h¯ μν = 0.. (elastics). (59). B. Long time limit as a fluid limit. Let the time scale of observation be Tobs . If the observation is made much longer than the relaxation times (i.e., Tobs τs ,τσ ,Reτ± ), then we can neglect the terms £u ε μν , £u εμ , and £u ε± in Eqs. (54)–(56) because, for example, £u ε μν ∼ −1 ε μν τs−1 ε μν . We thus obtain Tobs. G + η2 G + η2 K μν , K μν = λ1 η1 T μ μ H + σ2 ν (H + σ2 )T ν = , h μ ∂ν − h μ ∂ν − εμ τσ λ2 T σ1 T. ε μν τs. 026316-6. (60) (61).

(21) RELATIVISTIC VISCOELASTIC FLUID MECHANICS. PHYSICAL REVIEW E 84, 026316 (2011). ⎡. or (Ka − ζ5 )(a± γ1 + γ2 ) + (K + ζ4 )(a± γ2 + γ3 ) ⎢ ε± τ± trK ⎣ trε. 2 det γ θ. ⎤. ζ3 (ζ4 + K) − (ζ2 + K

(22) )(ζ5 − Ka) ⎥ ⎦. −(ζ2 − K

(23) )(ζ4 + K) + ζ1 (ζ5 − Ka). trK (det ζ s +K

(24) 2 )T. (62). By substituting these equations to Eqs. (45) and (46), the constitutive equations take the following form: (long) (q) τμν = τμν − 2ηNS K μν − ζNS (trK)hμν , μ (long) ν , νμ = σNS hμ ∂ν − T. By substituting Eqs. (71)–(73) into Eq. (45), the stress tensor can be rewritten in the following form: 2Gη3 (short) (q) = τμν − 2(G − η2 )ε μν − £u ε μν τμν (G + η2 )T − (K − ζ4 )trε − (Ka + ζ5 )θ ζ6 det γ + £u (trε) hμν . [(Ka − ζ5 )γ2 + (K + ζ4 )γ3 ]T (74). (63) (64). where we have defined viscosity and diffusion coefficients by det η + G 2 , (65) η1 T det ζ +K2 (a 2 ζ1 +2aζ2 +ζ3 )−2KK

(25) (aζ4 +ζ5 )+K

(26) 2 ζ6 , ζNS ≡ (det ζ s + K

(27) 2 )T (66) 2 det σ + H σNS ≡ . (67) σ1 ηNS ≡. Note that they are always non-negative, as can be seen from the inequality K2 (a 2 ζ1 + 2aζ2 + ζ3 ) − 2KK

(28) (aζ4 + ζ5 ) + K

(29) 2 ζ6 ⎞ ⎛ Ka ⎟ ⎜ (68) = (Ka K −K

(30) )ζ ⎝ K ⎠ 0.

(31) −K In particular, when the material is locally isotropic, we can take μν τ(q) = P hμν , with P the pressure, and thus the stress tensor certainly gives the constitutive equations for a relativistic Navier-Stokes fluid: τ(long) = −2ηNS K μν + (P − ζNS trK)hμν . μν. (69). These constitutive equations have the same form as those of a Kelvin-Voigt material (see Appendix B). However, one cannot yet identify the material at short time scales with a KelvinVoigt material, because they generically obey a different type of rheology equations. As discussed in the first subsection, elasticity is characterized by the condition that the intrinsic metric h¯ μν does not evolve, and the rheology equations for elastic materials are given by K¯ μν = 0, or equivalently by £u εμν = Kμν [6,7,15]. However, this is realized only when the conditions G + η2 = λ1 and (Ka − ζ5 )γ2 + (K + ζ4 )γ3 = det γ are satisfied. That is, for generic values of parameters, even if the observation time is sufficiently shorter than the relaxation times, the intrinsic metric h¯ μν evolves when the induced metric hμν does (i.e., K¯ μν = 0 if Kμν = 0). Thus, Maxwell’s original definition of viscoelasticity [considered only for the situations where the induced metric is static, Kμν = (1/2)£u hμν = 0] needs to be modified for generic values of parameters, such that h¯ μν is allowed to evolve when hμν does.. We thus confirm that our viscoelastic model always exhibits fluidity in the long time limit. C. Short time limit as an elastic limit. In contrast, at short time scales (Tobs τs ,Reτ± ), we have 1 1 £u ε μν − ε μν , £u ε± − ε± , τs τ±. (70). IV. SIMPLIFIED ISRAEL-STEWART FLUIDS. In this section, as an interesting example, we consider μν the case where K

(32) = η3 = σ3 = ζ6 = 0 and τ(q) = P hμν . In this case, from the positivity of matrices η, σ , and ζ , the conditions η2 = σ2 = ζ4 = ζ5 = 0 also must be imposed. Then the conserved currents take the following form:6. so that Eqs. (54)–(56) can be approximated as G + η2 £u ε μν K μν , (71) λ1 (Ka − ζ5 )(a± γ1 + γ2 ) + (K + ζ4 )(a± γ2 + γ3 ) £u ε± trK, 2 det γ (72) (Ka − ζ5 )γ2 + (K + ζ4 )γ3 ⇒£u (trε) trK . (73) det γ. T μν = euμ uν + P hμν − 2Gε μν − K(trε − aθ )hμν , nμ = nuμ − Hεμ ,. (75) (76). From this form of the bulk stress and the relation θ (T − T¯ )/T¯ , we see that a/T¯ can be identified with the thermal expansion coefficient. 6. 026316-7.

(33) MASAFUMI FUKUMA AND YUHO SAKATANI. PHYSICAL REVIEW E 84, 026316 (2011). and the rheology equations become 1 G £u ε μν = − ε μν + K μν , τs λ1 1 HT ν μ £u εμ = − εμ − h μ ∂ν − , τσ λ2 T £u (trε) = −. (77) (78). (γ2 ζ3 − γ3 ζ2 )T (γ3 ζ1 − γ2 ζ2 )T trε + θ det γ det γ. K(aγ2 + γ3 ) trK, det γ (γ2 ζ1 − γ1 ζ2 )T (γ1 ζ3 − γ2 ζ2 )T £u θ = trε − θ det γ det γ K(aγ1 + γ2 ) trK. − det γ +. (79). Here we have introduced τb ≡ γ1 /(ζ1 T ), and the viscosity and diffusion coefficients are given in this case by ηNS = τs G 2 /λ1 = G 2 /(η1 T ), ζNS = τb K2 /γ1 = K2 /(ζ1 T ), and σNS = H2 /σ1 . These equations look like the nonlinear causal dissipative hydrodynamics proposed in [14]. Although the nonlinear terms in [14] (e.g., hρμ νν ∇ρ uν ) are important for numerical simulations of ultrarelativistic dynamics, these terms, in principle, cannot be treated properly in our first-order formalism. However, if we do not make the approximation £u (trε) tr(£u εμν ), then Eq. (88) becomes −τb Ktr(£u εμν ) = τb [£u + (1/D)trK − KK μν ε μν ] = − − ζNS trK and coincides with Eq. (14) in [14] where the spatial dimension is set to be D = 1. If we neglect the nonlinear terms, we then get relations of Maxwell-Cattaneo type:. (80). π μν = −2ηNS K μν − τs hμγ hνδ uρ ∇ρ πγ δ , = −ζNS trK − τb uγ ∇γ , μ μ ν − τσ hμν uγ ∇γ ν ν , ν = σNS hμ ∂ν − T. By using the relations τ μν = −2Gε μν , ν μ = −Hεμ , 1 ≡ (trτ − trτ(q) ) = −K(trε − aθ ), D. (81). μν. the rheology equations can be rewritten as 1 2G 2 K μν , £u τ μν = − τ μν − τs λ1 μ 1 H2 T ν , h ∂ν − £u νμ = − νμ + τσ λ2 μ T. (82) (83). [(aγ2 + γ3 )ζ1 − (aγ1 + γ2 )ζ2 ]T det γ K2 (a 2 γ1 + 2aγ2 + γ3 ) trK − det γ KT [aζ1 (aγ2 + γ3 ) − ζ2 (a 2 γ1 − γ3 ) − ζ3 (aγ1 + γ2 )] + θ, det γ (84). where π μν ≡ τ μν − τ(q) . They are the constitutive equations for the simplified version of the Israel-Stewart model.7 Thus, in this case the rheology equations are equivalent to the constitutive equations for the simplified Israel-Stewart model (90), and the [D + 1 + 1 + D(D + 1)/2 + D] dynamical variables (excluding θ ) can be determined from the D + 2 conservation laws (∇μ nμ = ∇μ T μν = 0) and the D(D + 1)/2 + D equations (90).. £u = −. £u θ =. K(aγ1 + γ2 ) (γ1 ζ2 − γ2 ζ1 )T − trK K det γ det γ [γ1 (aζ2 + ζ3 ) − γ2 (aζ1 + ζ2 )]T − θ. det γ. (85). This model gives hyperbolic differential equations for small perturbations around a hydrostatic equilibrium, as is shown in Sec. V. For brevity, we here consider the case when θ is decoupled from other variables. This can be realized by setting a = γ2 = ζ2 = 0 in the above equations, and the rheology equations become τs £u τ μν = −τ μν − ηNS K μν , μ ν τσ £u νμ = −νμ + σNS hμ ∂ν − , T τb £u = − − ζNS trK, ζ3 T £u θ = − θ. γ3. (90). V. HYPERBOLICITY AND DISPERSION RELATIONS. In this section, we study linear perturbations around a hydrostatic equilibrium in Minkowski space-time. We exclusively take a coordinate system (x μ ) = (x 0 ,x i ) in which the background metric is written as gμν = ημν ≡ diag(−1,1, . . . ,1). A hydrostatic equilibrium is then specified by the velocity u(0) = uμ(0) ∂μ ≡ ∂0 (i.e., uμ(0) = δ0μ ), the proper energy density e(0) , the number density n(0) , and the (0) ≡ 0. The induced metric is then vanishing strain tensor Eμν (0) (0) (0) given by hμν = ημν + uμ uν = diag(0,1, . . . ,1). Note that from the fundamental relation for the hydrostatic equilib rium, s˜(0) = σ˜ (0) (e˜(0) ,n˜ (0) , h(0) ) ≡ h(0) s(0) (e(0) ,n(0) ), other thermodynamic quantities such as the temperature T(0) , the chemical potential μ(0) and the pressure P(0) are determined as δ s˜(0) =. 1 μ(0) P(0) δ e˜(0) − δ n˜ (0) + δ h(0) , T(0) T(0) T(0). (91). 1 μ(0) δe(0) − δn(0) , T(0) T(0). (92). or δs(0) =. (86) (87) 7. (88) (89). The constitutive equations for a simplified Israel-Stewart fluid is obtained by setting the viscous-heat coupling coefficients to be zero in those for an Israel-Stewart fluid (i.e., α0 = α1 = 0 in Eqs. (8a)–(8c) in [3]).. 026316-8.

(34) RELATIVISTIC VISCOELASTIC FLUID MECHANICS. PHYSICAL REVIEW E 84, 026316 (2011). ∂μ (T μν uν ) − T μν ∂μ uν = −∂μ (euμ ) − τ μν ∂μ uν . From this we obtain. with the Euler-Gibbs-Duhem relation s(0) =. e(0) μ(0) P(0) − + . T(0) T(0) T(0). (93). ∂μ (euμ ) = ∂0 δe + e(0) ∂i δui μν. = −τ μν ∂μ uν = −τ(0) ∂μ δuν = −P(0) ∂i δui , (104) or. A. Linear perturbations around a hydrostatic equilibrium. We now consider linear perturbations around the hydrostatic equilibrium, gμν = ημν + 0,. μ. uμ = δ0 + δuμ ,. hμν = h(0) μν + η0μ δuν + δuμ η0ν , e = e(0) + δe, n = n(0) + δn,. (94) Eμν = 0 + Eμν ,. (95) (96). ∂0 δe = −w(0) ∂i δui .. Here w(0) ≡ e(0) + P(0) is the enthalpy density at the hydrostatic equilibrium. As for the components orthogonal to uμ , from the equations 0 = hλν ∂μ T μν = hλν ∂μ (euμ uν ) + hλν ∂μ τ μν = euμ ∂μ uλ + hλν ∂μ τ μν = eaλ + hνλ ∂ μ τμν , we obtain eai = e(0) ∂0 δui = −hνi ∂ μ τμν = −∂ μ δτμi. and denote their conjugate thermodynamic variables by T = T(0) + δT , μ = μ(0) + δμ, P = P(0) + δP .. (97) (98). = −∂ 0 δτ0i − ∂ k δτik = P(0) ∂ 0 δui − ∂ k δτik ,. w(0) ∂0 δui = −∂ k δτik = −∂i δP + 2(G − η2 )∂ k ε ik (D − 2)η3 ζ6 η3 + + δui ∂i ∂k δuk + DT(0) T(0) T(0) + (K − ζ4 )∂i (trε) − (Ka + ζ5 )∂i θ, (107) where is the spatial Laplacian, ≡ δ ij ∂i ∂j . The conservation law of particle number current becomes. (99). 0 = ∂μ (nuμ + ν μ ). . μ = ∂0 δn + n(0) ∂i δu + σ3 δ − T. or. i. K ij . trK = ∂i δui ,. 1 2 k (0) ∂i δuj + ∂j δui − (∂k δu )hij . = 2 D . (100) (101). As for the stress tensor (45), by decomposing it as τμν = (0) (0) τμν + δτμν , the zeroth part is given by τμi = P(0) h(0) μi , and (0) from 0 = τμν uν = τμi δui + δτμ0 we can show that δτ00 = 0, δτi0 = −τij(0) δuj = −P(0) δui , and the spatial components are written as. ν = −(H − σ2 )ε + μ. μ. μν σ3 h(0) ∂ν δ. μ . − T. − (H − σ2 )∂i εi .. (B) The rheology equations are linearized as 1 G + η2 2 ∂0 ε ij = ∂i δuj + ∂j δui − (∂k δuk )h(0) − ε ij , ij 2λ1 D τs (109) 1 (H + σ2 )T(0) μ , (110) ∂0 εi = − εi − ∂i δ − τσ λ2 T [γ3 ζ1 + γ2 (K

(35) − ζ2 )]T(0) trε det γ [γ2 ζ3 − γ3 (K

(36) + ζ2 )]T(0) θ + det γ γ2 (Ka − ζ5 ) + γ3 (K + ζ4 ) + ∂i δui , det γ. ∂0 (trε) = −. ∂0 θ =. (103). We now substitute the above expressions to the set of fundamental equations, consisting of (A) the conservation laws (43)–(46) and (B) the rheology equations (47)–(49) [or (54)–(56)]. (A) As for the conservation laws of energy-momentum tensor, the component along uμ is given by 0 = uν ∂μ T μν =. (108). δτij = δP h(0) ij − 2(G − η2 )ε ij η3 2 k (0) ∂i δuj + ∂j δui − (∂k δu )hij − [(K − ζ4 )trε − T(0) D ζ 6 − (Ka + ζ5 )θ ]h(0) (∂k δuk )h(0) (102) ij − ij . T(0) The diffusion current is written as. (106). or. μν. We only consider the locally isotropic case: τ(q) = P hμν . Using the identity −1 = uμ uμ = −1 + 2δu0 = −1 − 2δu0 , we can show that δu0 = δu0 = 0, and the acceleration vector a μ = uν ∂ν uμ = ∂0 δuμ has only spatial components: a 0 = ∂0 δu0 = 0 and a i = ∂0 δui . Moreover, from 0 = εμν uν = εμ0 , εμν also has only spatial components, εij , in this linear approximation. Similarly, since 0 = Kμν uν = Kμ0 , the extrinsic curvature also has only spatial components, which are expressed as Kij = 12 hμi hνj (∂μ uν + ∂ν uμ ) = 12 (∂i δuj + ∂i δuj ),. (105). [γ1 (K

(37) − ζ2 ) + γ2 ζ1 ]T(0) trε det γ [γ1 ζ3 − γ2 (K

(38) + ζ2 )]T(0) θ − det γ γ2 (K + ζ4 ) + γ1 (Ka − ζ5 ) ∂i δui , − det γ. (111). (112). where we have used the approximation £u εij ∂0 εij , £u εi ∂0 εi , £u (trε) ∂0 (trε), and £u θ ∂0 θ .. 026316-9.

(39) MASAFUMI FUKUMA AND YUHO SAKATANI. PHYSICAL REVIEW E 84, 026316 (2011). Since we are considering locally isotropic materials, the fundamental thermodynamic relation (25) can be rewritten with the use of the Euler relation (C6) as μ 1 1 1 δe − δn − 2λ1 ε μν δε μν − (γ1 trε + γ2 θ )δ(trε) T T T T 1 1 (113) − λ2 εμ δεμ − (γ3 θ + γ2 trε)δθ. T T If we denote the thermodynamic variables collectively by (a r ) = (e,n,ε μν ,εμ ,trε,θ ), the matrix A ≡ −(∂ 2 s/∂a r ∂a s )|(0) is positive definite from the convexity of entropy. Here |(0) means that the matrix is evaluated at the hydrostatic state. In the following discussions, we assume for brevity that the matrix takes the following form: ⎛ ⎞ A1 A2 0 0 0 0 ⎜ ⎟ 0 0 0 0⎟ ⎜A2 A3 ⎜ ⎟. μν , ρσ ⎜0 0 A4 0 0 0⎟ ⎜ ⎟ (114) A=⎜ ⎟, μν ⎜0 ⎟ 0 0 A 0 0 5 ⎜ ⎟ ⎜ ⎟ 0 0 0 A6 A7 ⎠ ⎝0 0 0 0 0 A7 A8 δs =. where the principal submatrix . ∂ 2 s ∂2s − ∂e∂n − ∂e2 (0) A1 A2 (0) As ≡ = ∂2s ∂2s A2 A3 − ∂e∂n − ∂n2 (0) (0) . ∂(1/T ) ∂(μ/T ) − ∂e (0) ∂e (0) = ∂(μ/T ) − ∂(1/T ) ∂n. (0). ∂n. ∂0 δn = −n(0) ∂i δui + σ3 (A2 δe + A3 δn) + (H − σ2 )∂i εi , (120) 1 G + η2 2 ∂i δuj + ∂j δui − (∂k δuk )h(0) − ε ij , ∂0 ε ij = ij 2λ1 D τs (121) (H + σ2 )T(0) 1 ∂0 εi = (A2 ∂i δe + A3 ∂i δn) − εi , (122) λ2 τσ [γ3 ζ1 + γ2 (K

(40) − ζ2 )]T(0) ∂0 (trε) = − trε det γ [γ2 ζ3 − γ3 (K

(41) + ζ2 )]T(0) θ + det γ γ2 (Ka − ζ5 ) + γ3 (K + ζ4 ) + (123) ∂i δui , det γ [γ1 (K

(42) − ζ2 ) + γ2 ζ1 ]T(0) trε ∂0 θ = det γ [γ1 ζ3 − γ2 (K

(43) + ζ2 )]T(0) θ − det γ γ2 (K + ζ4 ) + γ1 (Ka − ζ5 ) ∂i δui . − (124) det γ We now consider wave propagations in the x D direction, demanding that perturbations depend only on x 0 and x D :. (115). (0). = s(0) ∂i δ[(1/T )−1 ] + n(0) ∂i δ[(1/T )−1 (μ/T )]. ∂i δ −. μ T. δe = δe(x ,x ), δn = δn(x ,x ).. (126). (116). = −(A2 ∂i δe + A3 ∂i δn),. (117). w(0) ∂0 δui = 2(G − η2 )∂ k ε ik +. (D − 2)η3 ζ6 + DT(0) T(0). (118) ∂i ∂k δuk. η3 δui + (K − ζ4 )∂i (trε) − (Ka + ζ5 )∂i θ T(0) − (w(0) A1 + n(0) A2 )T(0) ∂i δe − (w(0) A2 + n(0) A3 )T(0) ∂i δn, (119) +. 8 Note that the right-hand side of Eq. (C7) can be set to zero for the linear perturbations around a hydrostatic equilibrium.. 026316-10. G + η2 1 ∂D δuD − ε I I , Dλ1 τs 1 − ε I J (I = J ), τs G + η2 1 ∂D δuI − ε I D , 2λ1 τs η3 2 2(G − η2 )∂D ε I D + ∂ δuI , T(0) D 1 − εI , τσ −w(0) ∂D δuD ,. (127). ∂0 ε I J =. (128). w(0) ∂0 δuI = ∂0 εI =. and we finally obtain the following set of linearized equations of motion: ∂0 δe = −w(0) ∂i δui ,. D. ∂0 ε I I = −. ∂0 ε I D =. = (w(0) A1 + n(0) A2 )T(0) ∂i δe + (w(0) A2 + n(0) A3 )T(0) ∂i δn,. 0. D. Then the above equations can be rewritten as follows:. ∂i δP = s(0) ∂i δT + n(0) ∂i δμ. . (125). 0. is positive definite. Then the Gibbs-Duhem equation (C7) can be written as8. . δui = δui (x 0 ,x D ), εij = εij (x 0 ,x D ),. ∂0 δe = w(0) ∂0 δuD. (129) (130) (131) (132). 2(D − 1)η3 ζ6 ∂ 2 δuD = 2(G − η2 )∂D ε DD + + DT(0) T(0) D + (K − ζ4 )∂D (trε) − (Ka + ζ5 )∂D θ − (w(0) A1 + n(0) A2 )T(0) ∂D δe − (w(0) A2 + n(0) A3 )T(0) ∂D δn, (133) ∂0 δn = −n(0) ∂D δuD + σ3 A2 ∂D2 δe + A3 ∂D2 δn . + (H − σ2 )∂D εD , (134) (D − 1)(G + η2 ) 1 ∂0 ε DD = ∂D δuD − ε DD , (135) Dλ1 τs (H + σ2 )T(0) 1 (A2 ∂D δe + A3 ∂D δn) − εD , ∂0 εD = λ2 τσ (136).

(44) RELATIVISTIC VISCOELASTIC FLUID MECHANICS. γ2 (Ka − ζ5 ) + γ3 (K + ζ4 ) ∂D δuD det γ [γ3 ζ1 + γ2 (K

(45) − ζ2 )]T(0) trε − det γ [γ2 ζ3 − γ3 (K

(46) + ζ2 )]T(0) θ, + det γ γ2 (K + ζ4 ) + γ1 (Ka − ζ5 ) ∂D δuD ∂0 θ = − det γ [γ1 (K

(47) − ζ2 ) + γ2 ζ1 ]T(0) trε + det γ [γ1 ζ3 − γ2 (K

(48) + ζ2 )]T(0) θ, − det γ. PHYSICAL REVIEW E 84, 026316 (2011). scales (Tobs ∼ τs ) and will disappear at hydrodynamic time scales (Tobs τs ).. ∂0 (trε) =. C. Shear modes. (137). (138). where I,J = 1, . . . ,D − 1. This set of equations can be further decomposed according to the transformation properties under the little group SO(D − 1): (1) tensor modes: (ε I I ,ε I J ); (2) shear modes: (ε I D ,δuI ,εI ); (3) sound modes: (trε,ε DD ,δuD ,δe,δn,εD ,θ ). In the remainder of this section, we study hyperbolicity and dispersion relations for each type of perturbation modes. B. Tensor modes. For tensor modes, the set of equations can be written as G + η2 1 ∂D δuD − ε I I , (139) Dλ1 τs 1 ∂0 ε I J = − ε I J . (140) τs From the identity I ε I I + ε DD = 0, the number of independent variables of ε I I is D − 2. If we define the variables EI J by ∂0 ε I I = −. EI J ≡. ε I I − ε (D−1)(D−1). (for I = J ),. ε I J . (for I = J ),. (141). then the number of independent EI J is D(D − 1)/2 − 1 = (D − 2)(D + 1)/2 because E(D−1)(D−1) = 0, and the equations for EI J become 0 = ∂0 EI J +. 1 EI J . τs. (142). For shear modes, we have the equations G + η2 ∂D δuI , 0 = ∂0 + τs−1 ε I D − 2λ1 2(G − η2 ) η3 ∂D ε I D − ∂ 2 δuI , 0 = ∂0 δuI − w(0) w(0) T(0) D 0 = ∂0 + τσ−1 εI .. εij = εij (ω,k)eikx. D. −iωx 0. ,. (143) (144). we obtain the dispersion relation ω = −i/τs which represents nonpropagating, purely dissipating modes. Since τs is positive, the imaginary part of ω is always negative, and thus we find that the tensor modes are always stable. Such relaxation modes correspond to stress relaxations observed at rheological time. (146) (147). Note that εI is decoupled from the other variables, and Eq. (147) represents its pure relaxation with relaxation time τσ (0). If we set η3 = 0, by redefining the variables by sI ± ≡ ±. λ1 1 ε I D + δuI , w(0) 2. (148). the set of linearized equations for sI ± can be written as sI + − sI − = 0, (149) (∂0 ∓ cshear ∂D )sI ± ± 2τs for I = 1,2, . . . ,D − 2. These are hyperbolic equations and the characteristic velocity is given by ! ηNS cshear ≡ . (150) w(0) τs For generic cases, from Eqs. (145) and (146), we obtain telegrapher’s equations with Kelvin-Voigt damping, ε I D 1 η3 2 ∂D2 ∂0 − cshear ∂D2 ∂02 + ∂0 − = 0. (151) δuI τs w(0) T(0) Although they are generically nonhyperbolic and have infinite wave-front velocity as in the standard relativistic fluid mechanics, they can be made into hyperbolic telegrapher’s equations by setting η3 = 0.9 If we consider the short time limit (τs → ∞), the differential equations become η3 ε I D 2 = 0. (152) ∂02 − ∂D2 ∂0 − cshear ∂D2 δuI w(0) T(0) The wave equations in this form also appear for viscous solids such as Kelvin-Voigt materials, and reduce to the wave equations when η3 = 0. Finally, for plane waves i (ω,k)eikx D −iωx 0 , δui = δu. Thus, if we consider plane waves propagating in the x D direction, i (ω,k)eikx D −iωx 0 , δui = δu. (145). εij = εij (ω,k)e. ikx D −iωx 0. ,. from (151), we obtain the dispersion relation 1 η3 2 2 + + k 2 + cshear k 2 = 0, τs w(0) T(0). (153) (154). (155). 9 In this case, from the non-negativeness of the matrix η, η2 (and thus det η) must vanish. However, this still gives a positive shear viscosity if G = 0, as can be seen from Eq. (65).. 026316-11.

(49) MASAFUMI FUKUMA AND YUHO SAKATANI. PHYSICAL REVIEW E 84, 026316 (2011). where ≡ −iω. Since all the coefficients are positive, the real part of (or the imaginary part of ω) always takes negative values, and thus we see that there are no unstable growing modes in the shear modes. Equation (155) has two solutions, which are expanded around k = 0 as i 2 4 τs k 2 + i(1−rs )cshear τs3 k 4 + O(k 6 ) − τs + i(1−rs )cshear ω= 2 4 τs k 2 − i(1 − rs )cshear τs3 k 4 + O(k 6 ), −icshear. relation of Maxwell-Cattaneo type, the effective relaxation time associated with the hydrodynamic modes is read off from the coefficients of k 4 as (1 − rs )τs . Indeed, if we set rs = 1, the effective relaxation time becomes zero and the dispersion relation becomes purely diffusive; ω = −i(ηNS /w(0) )k 2 . If we are interested only in the hydrodynamic modes, the dispersion relation coincides with that of the Israel-Stewart model up to O(k 4 ) by identifying (1 − rs )τs with the relaxation time τπ in the Israel-Stewart model. However, if the relaxation modes are also taken into account, our viscoelastic model has a richer structure than the Israel-Stewart model, which is the special case (rs = 0) of the viscoelastic model.. (156) with rs ≡ η3 /(ηNS T(0) ). The former represents the relaxation modes which are not observed at hydrodynamic time scales (Tobs τs ). The latter represents the hydrodynamic modes where ω → 0 in the limit k 2 → 0, and from the coefficients 2 τs = of k 2 , the diffusion coefficient is found to be cshear ηNS /w(0) . Moreover, by the comparison with the dispersion. D. Sound modes. Finally, for sound modes, we have the following set of differential equations:. 0 = ∂0 δe + w(0) ∂D δuD , 0 = ∂0 δuD −. (157). 2(G − η2 ) ∂D ε DD − w(0). . 2(D − 1)η3 ζ6 + Dw(0) T(0) w(0) T(0). ∂D2 δuD −. K − ζ4 Ka + ζ5 ∂D (trε) + ∂D θ w(0) w(0). (w(0) A1 + n(0) A2 )T(0) (w(0) A2 + n(0) A3 )T(0) ∂D δe + ∂D δn, w(0) w(0) 0 = ∂0 δn + n(0) ∂D δuD − σ3 A2 ∂D2 δe + A3 ∂D2 δn − (H − σ2 )∂D εD , +. 0 = ∂0 ε DD − 0 = ∂0 εD +. (159). (D − 1)(G + η2 ) 1 ∂D δuD + ε DD , Dλ1 τs. (160). 1 (H + σ2 )T(0) εD − (A2 ∂D δe + A3 ∂D δn), τσ λ2. 0 = ∂0 (trε) + 0 = ∂0 θ +. (158). (161). γ2 (Ka − ζ5 ) + γ3 (K + ζ4 ) [γ3 ζ1 + γ2 (K

(50) − ζ2 )]T(0) [γ2 ζ3 − γ3 (K

(51) + ζ2 )]T(0) trε − ∂D δuD − θ, det γ det γ det γ. [γ1 ζ3 − γ2 (K

(52) + ζ2 )]T(0) γ2 (K + ζ4 ) + γ1 (Ka − ζ5 ) [γ1 (K

(53) − ζ2 ) + γ2 ζ1 ]T(0) ∂D δuD − trε + θ. det γ det γ det γ. (162). (163). In particular, if we consider the case where η3 = ζ6 = σ3 = 0, the set of equations reduces to the following linear differential equations: (∂0 +B0 ∂D + B1 )Y = 0, ⎛. 0. ⎜ ⎜ ⎜ 0 ⎜ ⎜ K ⎜− ⎜ w(0) ⎜ B0 ≡ ⎜ ⎜ 0 ⎜ ⎜ 0 ⎜ ⎜ ⎜ 0 ⎜ ⎝ 0. (164). 0. 2 +γ3 ) − K(aγ det γ. 0. 0. 0. 0. − (D−1)G Dλ1. 0. 0. 0. − w2G(0). 0. 0. w(0). 0. 0. 0. 0. n(0). 0. 0. −H. 0. 0. 0. K(aγ1 +γ2 ) det γ. T(0) (w(0) A1 w(0). −. + n(0) A2 ). A2 T(0) H λ2. 0. 026316-12. T(0) (w(0) A2 w(0). −. + n(0) A3 ). A3 T(0) H λ2. 0. 0. 0 0. 0. ⎞. ⎟ ⎟ 0 ⎟ ⎟ ⎟ Ka ⎟ w(0) ⎟ ⎟ 0 ⎟ ⎟, ⎟ 0 ⎟ ⎟ ⎟ 0 ⎟ ⎟ ⎠ 0. (165).

(54) RELATIVISTIC VISCOELASTIC FLUID MECHANICS. ⎛. [γ3 ζ1 +γ2 (K

(55) −ζ2 )]T(0) det γ. ⎜ ⎜ 0 ⎜ ⎜ ⎜ 0 ⎜ ⎜ 0 B1 ≡ ⎜ ⎜ ⎜ 0 ⎜ ⎜ ⎜ 0 ⎝

(56) [γ (K −ζ2 )+γ2 ζ1 ]T(0) − 1 det γ Here we have defined cs2 ≡. PHYSICAL REVIEW E 84, 026316 (2011). −. [γ2 ζ3 −γ3 (K

(57) +ζ2 )]T(0) ⎞ det γ. 0. 0. 0 0. 0. 1 τs. 0. 0. 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0 0 [γ1 ζ3 −γ2 (K

(58) +ζ2 )]T(0) det γ. 0. 0. 0 0. 1 τσ. 0. 0. 0 0. 0. ⎛ ⎞ trε ⎟ ⎟ ⎜ε DD ⎟ ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ ⎟ ⎜ δuD ⎟ ⎟ ⎜ ⎟ ⎟ δe ⎟ ⎟ , Y ≡ ⎜ ⎜ ⎟. ⎟ ⎜ δn ⎟ ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ ⎟ ⎝ ε D ⎠ ⎟ ⎠ θ. ∂p T(0) 2 , w(0) A1 + 2w(0) n(0) A2 + n2(0) A3 = w(0) ∂e s. (166). (167). n. and. ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ M≡⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝. w(0) (aγ2 +γ3 ) det γ. ⎞ . 0. 0 w(0) (D−1) 2Dλ1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. w(0) cs. w A2 +n(0) A3 − (0)√det As cs w(0) A1 +n(0) A2 √ det A(1) cs. 0. 0. 0. n(0) cs. 0. 0. 0. 0. 0. w(0) λ2. 0. 0. 0. 0. 0. 0. . . 0 w(0) (aγ1 +γ2 ) a det γ. Y

(59) ≡ M −1 Y .. ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠. (168). (169). We then have (∂0 + M −1 B0 M∂D + M −1 B1 M)Y

(60) = 0, with. ⎛. 0 ⎜ 0 ⎜ ⎜ ⎜−M1 ⎜ ⎜ M −1 B0 M = ⎜ 0 ⎜ ⎜ 0 ⎜ ⎜ 0 ⎝ 0 ⎛ M6 ⎜ 0 ⎜ ⎜ ⎜ 0 ⎜ −1 M B1 M = ⎜ ⎜ 0 ⎜ 0 ⎜ ⎜ ⎝ 0 M8. −M1 −M2 0 cs 0 0 M3. 0 0 −M2 0 0 0 0. 0 0 cs 0 0 −M4 0. 0. 0. 0 0. 0. 1 τs. 0 0 0. 0 0 0 0. 0 0 0 0. 0 0 0 0. 0 0 0 0. 0 0. 0 0. 0 0 0 0. 1 τσ. 0. M7. 0 0 0 0 0 −M5 0 ⎞. ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎟ 0 ⎠ M9 0 0 0 0. (170) 0 0 0 −M4 −M5 0 0. ⎞ 0 0 ⎟ ⎟ ⎟ M3 ⎟ ⎟ 0 ⎟ ⎟, ⎟ 0 ⎟ ⎟ 0 ⎟ ⎠ 0. (171). (172). M1 ≡. K2 (aγ2 + γ3 ) , w(0) det γ. (173). M2 ≡. 2(D − 1)ηNS , Dw(0) τs. (174). 026316-13.

(61) MASAFUMI FUKUMA AND YUHO SAKATANI. PHYSICAL REVIEW E 84, 026316 (2011). K2 a(aγ1 + γ2 ) , w(0) det γ H w(0) A2 + n(0) A3 M4 ≡ cs M3 ≡. (175) T(0) , τσ w(0) σ1. H det As T(0) w(0) , cs τσ σ1 [γ3 ζ1 + γ2 (K

(62) − ζ2 )]T(0) M6 ≡ , det γ M5 ≡. M7 ≡. [(K

(63) + ζ2 )γ3 − ζ3 γ2 ]T(0) det γ. [(K

(64) − ζ2 )γ1 + ζ1 γ2 ]T(0) det γ [γ1 ζ3 − γ2 (K

(65) + ζ2 )]T(0) M9 ≡ . det γ. M8 ≡ −. (176) (177) (178). aγ1 + γ2 , a(aγ2 + γ3 ) a(aγ2 + γ3 ) , aγ1 + γ2. (179) (180) (181). The real matrix M −1 B0 M is symmetric and can be diagonalized. The eigenvalues are calculated to be {0,0,0, ± v± }, where 1 2 ζNS 2(D − 1)ηNS A3 σNS 2 cs + ≡ + + v± 2 w(0) τb Dw(0) τs τσ ζNS 2(D − 1)ηNS A3 σNS 2 4A3 σNS ζNS 2(D − 1)ηNS det As w(0) T(0) 1 cs2 + (182) + + − + + ± 2 w(0) τb Dw(0) τs τσ τσ w(0) τb Dw(0) τs A3 give the characteristic velocities. Since all the eigenvalues are real, we see that the system of differential equations (164) is hyperbolic. If we particularly set H = 0 (and thus σNS = 0), the characteristic velocity reduces to v± =. cs2 +. ζNS 2(D − 1)ηNS + w(0) τb Dw(0) τs. (183). and agrees with the large wave-number limit of the group velocity (which in our case coincides with the front velocity and the characteristic velocity) in the M¨uller-Israel-Stewart theory (see, e.g., Eq. (49) in [16]). If we take the long time limit, τb ,τs → 0, the characteristic velocity becomes infinitely large, and thus causality gets violated. For generic cases (i.e., when we do not impose the conditions η3 = ζ6 = σ3 = 0), from Eqs. (157)–(163), the dispersion relation for the plane wave i (ω,k)eikx D −iωx 0 , δui = δu D 0 εij = εij (ω,k)eikx −iωx , D 0 δe = δe(ω,k)eikx −iωx , ikx D −iωx 0 δn = δn(ω,k)e ,. (184) (185) (186) (187). is obtained as 7 + (c60 + c62 k 2 ) 6 + (c50 + c52 k 2 + c54 k 4 ) 5 + (c40 + c42 k 2 + c44 k 4 ) 4 + (c30 + c32 k 2 + c34 k 4 ) 3 + (c22 k 2 + c24 k 4 ) 2 + (c12 k 2 + c14 k 4 ) + c04 k 4 = 0,. (188). where = −iω and c60 = τs−1 + τσ−1 + τ+−1 + τ−−1 , ζ6 2(D − 1)η3 + + A3 σ3 , c62 = T(0) w(0) DT(0) w(0) τs τσ + (τs + τσ )(τ+ + τ− ) + τ+ τ− c50 = , τs τσ τ+ τ− ζ6 τs−1 + τσ−1 2(D − 1)η3 τˆs−1 + τσ−1 + τ+−1 + τ−−1 ζNS τb−1 c52 = cs2 + + + + A3 σ3 τs−1 + τˆσ−1 + τ+−1 + τ−−1 , T(0) w(0) w(0) Dw(0) T(0) 026316-14. (189) (190) (191) (192).

(66) RELATIVISTIC VISCOELASTIC FLUID MECHANICS. . c54 c40 c42. c44. ζ6 2(D − 1)η3 = A3 σ3 + w(0) T(0) Dw(0) T(0) τs + τσ + τ+ + τ− = , τs τσ τ+ τ−. PHYSICAL REVIEW E 84, 026316 (2011). (193). ,. (194). ζNS τs τσ + (τs + τσ )τb−1 τ+ τ− ζ6 2 −1 −1 −1 −1 = cs τs + τσ + τ+ + τ− + + w(0) T(0) τs τσ w(0) τs τσ τ+ τ− 2(D − 1)ηNS [τˆs τσ + (τˆs + τσ )(τ+ + τ− ) + τ+ τ− ] A3 σNS [τs τˆσ + (τs + τˆσ )(τ+ + τ− ) + τ+ τ− ] + + , Dw(0) τs τσ τ+ τ− τs τσ τ+ τ− . ζ6 τˆσ−1 + τs−1 2(D − 1)η3 τˆs−1 + τˆσ−1 + τ+−1 + τ−−1 ζNS τb−1 det AT(0) w(0) = A3 σ3 + + + , T(0) w(0) w(0) DT(0) w(0) A3. c30 =. 1 , τs τσ τ+ τ−. (195) (196) (197). cs2 [τs τσ + (τs + τσ )(τ+ + τ− ) + τ+ τ− ] ζNS τs + τσ + τb−1 τ+ τ− + τs τσ τ+ τ− w(0) τs τσ τ+ τ− 2(D − 1)ηNS (τˆs + τσ + τ+ + τ− ) A3 σNS (τs + τˆσ + τ+ + τ− ) + + , Dw(0) τs τσ τ+ τ− τs τσ τ+ τ− ζNS [τs τˆσ + (τs + τˆσ )τb−1 τ+ τ− ] ζ6 = A3 σNS + T(0) w(0) τs τσ w(0) τs τσ τ+ τ−. c32 =. c34. c22 c24. 2(D − 1)ηNS [τˆs τˆσ + (τˆs + τˆσ )(τ+ + τ− ) + τ+ τ− ] det Aσ3 T(0) w(0) τˆσ−1 + τs−1 + τ+−1 + τ−−1 , + + Dw(0) τs τσ τ+ τ− A3 σNS 1 ζNS 2(D − 1)ηNS cs2 (τs + τσ + τ+ + τ− ) + = + + A3 σNS , τs τσ τ+ τ− w(0) Dw(0) ζNS (τs + τˆσ + τb−1 τ+ τ− ) 2(D − 1)ηNS (τˆs + τˆσ + τ+ + τ− ) A3 σNS = + τs τσ τ+ τ− w(0) Dw(0) det Aw(0) T(0) [τs τˆσ + (τs + τˆσ )(τ+ + τ− ) + τ+ τ− ] , + A3 cs2 , τs τσ τ+ τ− ζNS A3 σNS 2(D − 1)ηNS det Aw(0) T(0) (τs + τˆσ + τ+ + τ− ) , = + + τs τσ τ+ τ− w(0) Dw(0) A3 det AσNS T(0) w(0) = . τs τσ τ+ τ−. (198). (199) (200). (201). c12 =. (202). c14. (203). c04. Here we have defined non-negative constants η3 σ3 τs , τˆσ ≡ τσ , τˆs ≡ rs τs = ηNS T(0) σNS. (205). and redefined τb as τb ≡. ζNS det γ , K2 γ+ + Pζ ζ γ. γ+ ≡ a 2 γ1 + 2aγ2 + γ3 0, Pζ ζ γ ≡ ζ3 ζ6 − ζ52 γ1 + 2(ζ4 ζ5 − ζ2 ζ6 )γ2 + ζ1 ζ6 − ζ42 γ3 0,. (204). that all the coefficients are positive, and thus at least the necessary condition for the stability is satisfied. For a full analysis to be performed, one should further check the RouthHurwitz stability criterion, which we have not carried out yet. The dispersion relation around k = 0 gives seven solutions, and four of the seven take the following form:. (206) i i + O(k 2 ), ω = − + O(k 2 ), τ± τs i ω = − + O(k 2 ). τσ. ω=−. (207) (208). which becomes γ1 /(ζ1 T ) when the parameters are taken as in Sec. IV. Note that complex parameters τ± appear always through the combinations τ+ + τ− = 2Reτ+ (0), τ+ τ− = |τ+ |2 (0), or τ+−1 + τ−−1 = 2Reτ+ /|τ+ |2 (0). One can check. (209) (210). They correspond to the relaxation modes, and as the observation time becomes much longer than the relaxation times Re τ± , τs , and τσ , these modes fade away in time and will not be observed eventually.. 026316-15.