Linear stability of Israel-Stewart theory in the presence of net-charge diffusion
LLinear stability of Israel-Stewart theory in the presence of net-charge diffusion
C. V. Brito ∗ and G. S. Denicol † Instituto de F´ısica, Universidade Federal Fluminense, UFF, Niter´oi, 24210-346, RJ, Brazil
In this paper, we perform a linear stability analysis of Israel-Stewart theory around a global equi-librium state, including the effects of shear-stress tensor, net-baryon diffusion current and diffusion-viscous coupling. We find all the relevant modes of this theory and derive necessary conditions thatthese modes must satisfy in order to be stable and subluminal. With these conditions, we thenderive constraints for the shear and diffusion relaxation times and the transport coefficients relatedto diffusion-viscous coupling.
I. INTRODUCTION
Relativistic fluid-dynamics is an effective theory constructed to describe the long-wavelength and long-time dynamicsof a many-body system, with applications that range from astrophysics [1] to high-energy nuclear physics [2]. In thepast 20 years, there has been a growing interest in understanding fundamental aspects of relativistic hydrodynamics[3], mostly due to the fluid-dynamical modeling of the quark-gluon plasma (QGP), produced in modern particleaccelerators via ultra-relativistic heavy-ion collisions [3–6].Relativistic generalizations of Navier-Stokes theory, derived by Landau and Lifshitz [7] and, independently, byEckart [8], are known to be ill-defined, containing intrinsic instabilities when perturbed around an arbitrary globalequilibrium state [9–11]. Such instabilities are intrinsically related to the acausal nature of Navier-Stokes theory,which allows perturbations to propagate with an infinite speed. These fundamental problems prohibit the applicationof Navier-Stokes theory to describe any practical fluid-dynamical problem, may it be in the description of neutronstar mergers or in the description of the quark-gluon plasma produced in heavy ion collisions.Stable and causal theories of relativistic fluid dynamics were later derived by Israel and Stewart, following theprocedure initially developed by H. Grad [12] for nonrelativistic systems. Israel and Stewart performed this task intwo distinct ways: (i) the first being a phenomenological derivation, based on the second law of thermodynamics [13],and (ii) the second being a microscopic derivation starting from the relativistic Boltzmann equation [14]. Similartheories have been widely developed in the past decades [15–24], but all carry the same fundamental aspects: incontrast to Navier-Stokes theory, such causal theories of fluid dynamics include in their description the transientdynamics of the non-conserved dissipative currents. For this reason, they were initially named by Israel and Stewartas transient fluid dynamics (nowadays, they are often referred to as second order theories ).However, it is important to remark that the theory formulated by Israel and Stewart is not guaranteed to be causaland stable. As was first shown by Hiscock, Lindblom and, later, by Olson, such transient theories of fluid dynamicsare only linearly causal and stable if their transport coefficients satisfy certain conditions [25, 26]. Such conclusionswere obtained by analysing the properties of the theory in the linear regime and by imposing that the perturbationsaround a global equilibrium state are stable and propagate subluminally. More recent analysis were developed inRef. [10], including only the effects of bulk viscosity, and, later, in Ref. [11], which included the effects of both shearand bulk viscosity. In both these papers, constraints for the shear and bulk relaxation times were explicitly derived.Nowadays, causality analysis have even been performed also in the nonlinear regime [27, 28] (in this case, includingthe effects of shear and bulk viscosity), where more general inequalities required to ensure the causal propagation ofthe theory were derived. In the latter case, the inequalities constrain not only the tranport coefficients, but also thevalues of the dissipative currents (in the linear regime, the inequalities derived in Ref. [28] reduce to those derivedin Refs.[10, 11]). Such constraints are relevant for, e.g. fluid-dynamical applications in heavy ion collisions, sincethe transport coefficients of QCD matter are not precisely known (often, they are completely unknown) and suchfundamental constraints on transport coefficients (and the values of the dissipative currents) can be extremely useful.Recently, several programs to experimentally study QCD matter at finite net-baryon density have been put in motionat the Relativistic Heavy Ion Collider (RHIC), in Brookhaven National Lab (Upton, USA), and at the Nuclotron-basedIon Collider fAcility (NICA), in the Joint Institute for Nuclear Research (Dubna, Russia), and, will be starting soon atthe Facility for Antiproton and Ion Research (FAIR), in GSI Helmholtzzentrum fur Schwerionenforschung (Darmstadt, ∗ Electronic address: caio [email protected]ff.br † Electronic address: [email protected] a r X i v : . [ nu c l - t h ] J u l Germany). Nevertheless, the more recent investigations on the stability and causality of fluid-dynamical descriptions[10, 11, 27, 28], have not yet considered the complete set of the Israel-Stewart equations, usually neglecting anydissipation by net-baryon diffusion and, also, possible diffusion-viscous coupling terms . In this paper, we actuallyperform a linear stability analysis around global equilibrium of Israel-Stewart theory, including the effects of shear-stress tensor and net-baryon diffusion 4-current (all effects of bulk viscous pressure are neglected). We find all therelevant modes of this theory and derive the conditions that these modes must satisfy in order to be stable andsubluminal. With this result, we obtain new conditions that the shear and diffusion relaxation times must satisfy sothat Israel-Stewart theory remains linearly causal and stable. We further find constraints for the transport coefficientsthat couple the shear stress tensor and the net-baryon diffusion current (diffusion-viscous coupling). In other words,we show that the inclusion of diffusion-viscous coupling in the equations of motion drives the theory unstable, ifthese transport coefficients do not satisfy certain bounds. Such novel constraints may be useful when such transportcoefficients are included in the current fluid dynamical simulations of the quark-gluon plasma.This paper is divided as follows. In Sec. II we review the fundamentals of relativistic hydrodynamics, as proposed byIsrael and Stewart [13]. Then, in Sec. III, we linearize the Israel-Stewart equations around a global equilibrium state,expressing the resulting equations in Fourier space. In particular, we demonstrate how to decompose the linearizedequations of motion in Fourier space into independent equations of motion for their longitudinal and transversecomponents – a procedure that considerably simplifies the calculations. Next, Secs. IV and V are dedicated tothe study of the theory’s linear stability in the absence and in the presence of diffusion-viscous coupling terms,respectively. All our conclusions are summarized in Sec. VI. In this paper, we use natural units c = (cid:126) = k B = 1, andthe mostly-minus convention for the metric tensor g µν = diag(+ , − , − , − ). II. RELATIVISTIC FLUID DYNAMICS
In this section, we briefly review the equations of relativistic dissipative fluid dynamics. The main equations arethe continuity equations associated to the conservation of net-charge, energy, and momentum ∂ µ T µν = 0 , (1) ∂ µ N µ = 0 , (2)where T µν is the energy-momentum tensor and N µ is the net-charge 4-current. In our case, the only conserved chargeconsidered is baryon number, with electric charge and strangeness being neglected.It is convenient to write the conserved currents in a fluid-dynamical form, T µν = εu µ u ν − P ∆ µν + π µν , (3) N µ = n B u µ + n µ , (4)where ε is the energy density, n B is the net-baryon number density, P ( n B , ε ) is the thermodynamic pressure, u µ isthe normalized 4-velocity ( u µ u µ = 1), n µ is the net-baryon diffusion current, and π µν is the shear-stress tensor. Inthis work, we employ Landau’s picture for the velocity field [7], where it is defined as an Eigenvector of the energy-momentum tensor, i.e., T µν u ν ≡ εu µ . We also introduced the projection operator onto the 3-space orthogonal to u µ ,∆ µν ≡ g µν − u µ u ν . All effects of dissipation due to bulk viscous pressure are neglected in our calculations.When considering ideal fluids, the conservation laws, supplemented by an equation of state, are sufficient to describethe time evolution of the system. On the other hand, when considering viscous fluids, the conservation laws mustalso be complemented by dynamical equations for the dissipative currents. Here, we employ relaxation-type equationsderived from kinetic theory as our baseline [21, 29, 30], τ n ˙ n (cid:104) µ (cid:105) + n µ = κ n ∇ µ α B − n ν ω νµ − δ nn n µ θ + (cid:96) nπ ∆ µν ∇ λ π λν − τ nπ π µν ∇ ν P − λ nn n ν σ µν − λ nπ π µν ∇ ν α B , (5) τ π ˙ π (cid:104) µν (cid:105) + π µν = 2 ησ µν + 2 τ π π (cid:104) µλ ω ν (cid:105) λ − δ ππ π µν θ − τ ππ π λ (cid:104) µ σ ν (cid:105) λ − τ πn n (cid:104) µ ∇ ν (cid:105) P + (cid:96) πn ∇ (cid:104) µ n ν (cid:105) + λ πn n (cid:104) µ ∇ ν (cid:105) α B , (6) The work by Olson [26] considered the complete Israel-Stewart equations, including all sources of fluctuations. Nevertheless, in thiswork, the dispersion relations for the hydrodynamic modes were not explicitly evaluated. where α B ≡ µ B /T , with µ B being the baryon chemical potential and T the temperature. Above, we defined thecomoving time derivative, ˙ A ≡ u µ ∂ µ A , the expansion rate, θ ≡ ∂ µ u µ , the shear tensor σ µν ≡ ∂ (cid:104) µ u ν (cid:105) , the vorticitytensor, ω µν ≡ ( ∇ µ u ν − ∇ ν u µ ) /
2, and the projected derivative ∇ µ ≡ ∂ (cid:104) µ (cid:105) . We further employ the notation, A (cid:104) µ (cid:105) ≡ ∆ µν A ν , and A (cid:104) µν (cid:105) ≡ ∆ µναβ A αβ , where we make use of the double, symmetric, and traceless projection operator ∆ µναβ = (cid:16) ∆ µα ∆ νβ + ∆ µβ ∆ να (cid:17) / − ∆ µν ∆ αβ /
3. Note that all the dissipative currents are constructed to be orthogonal to the4-velocity, u µ n µ = 0 , u µ π µν = 0. (7)In the equations of motion for the dissipative currents, we introduced a wide set of transport coefficients, allbeing functions of temperature and chemical potential. The most relevant for our work are the net-charge diffusioncoefficient, κ n , the shear viscosity coefficient, η , and the diffusion and shear relaxation times, τ n and τ π , respectively.The effects of the relaxation times have already been extensively investigated in the linear regime in [11, 25, 26] andwere shown to be essential to render the theory causal and stable. So far, the remaining coefficients have not beenwidely investigated, even in the linear regime. Such coefficients can be relevant as they determine the strength of thesecond-order terms that appear in the fluid-dynamical equations. In particular, we are interested in the effects of thecoupling terms ∆ µν ∇ λ π λν and ∇ (cid:104) µ n ν (cid:105) , which are associated with the transport coefficients (cid:96) nπ and (cid:96) πn , respectively.The effect of these transport coefficients can be investigated in the linear regime. We note that such coupling termsalso appear in the phenomenological derivation of Israel-Stewart theory [13, 14].For the sake of convenience, we further re-express the conservation laws, Eqs. (1) and (2), in the following way, u ν ∂ µ T µν = ˙ ε + ( ε + P ) θ − π αβ σ αβ = 0 , (8)∆ λν ∂ µ T µν = ( ε + P ) ˙ u λ − ∇ λ P − π λβ ˙ u β + ∆ λν ∇ µ π µν = 0 , (9) ∂ µ N µ = ˙ n B + n B θ − n µ ˙ u µ + ∇ µ n µ = 0 . (10)The fluid-dynamical equations will be linearized in the above form. III. LINEARIZED FLUID DYNAMICS
In this section, we linearize the fluid-dynamical equations described in the previous section around a global equilib-rium state. For this purpose, we consider small fluid-dynamical perturbations around a hydrostatic equilibrium state,with an energy density, ε , a vanishing net-baryon number density, n B, = 0, a constant 4-velocity, u µ , and vanishingdissipative currents, n µ = π µν = 0.Perturbations around such an equilibrium state can be expressed in the following form, ε = ε + δε, n B = δn B , u µ = u µ + δu µ , (11) n µ = δn µ , π µν = δπ µν . (12)As already stated, in this analysis we neglect any effects from bulk viscous pressure, i.e., δ Π = 0.Since the fluid 4-velocity is normalized, i.e., u µ u µ = 1, it is straightforward to demonstrate that the perturbationsof the fluid velocity satisfy u µ δu µ = O (2) , (13)where O (2) denote all possible terms that are second order or higher in perturbations of the fluid-dynamical fields.That is, up to first order in perturbations, the fluctuations of the fluid 4-velocity are orthogonal to the background4-velocity. Similarly, due to the orthogonality relations satisfied by the dissipative currents, Eq. (7), one can showthat all perturbations of the dissipative currents are also orthogonal to the background 4-velocity, u µ δπ µν = O (2) , (14) u µ δn µ = O (2) . (15)Such linear orthogonality relations satisfied with the background 4-velocity, motivate us to also introduce projectionoperators that are constructed from the background 4-velocity,∆ µν ≡ g µν − u µ u ν , (16)∆ µναβ ≡ (cid:16) ∆ µα ∆ νβ + ∆ µβ ∆ να (cid:17) −
13 ∆ µν ∆ αβ , (17)and define the following projected derivatives D ≡ u µ ∂ µ , ∇ µ ≡ ∆ µν ∂ µ . (18)Retaining only contributions which are linear in the perturbations, the fluid-dynamical equations of motion reduceto D (cid:18) δεw (cid:19) + ∇ µ δu µ = O (2) , (19) D δu µ − ∇ µ (cid:18) δPw (cid:19) + ∇ ν δχ µν = O (2) , (20) D (cid:18) δn B n (cid:19) + ∇ µ δξ µ = O (2) , (21)where w = ε + P is the enthalpy and n ≡ ( ε + P ) / T is the particle number density at vanishing chemical potential.We further defined the dimensionless fields associated with the hydrodynamic fluctuations, δχ µν ≡ δπ µν / ( ε + P )and δξ µ ≡ δn µ /n . In this notation, the linearized equations of motion for the dissipative currents become τ n D δξ µ + δξ µ = ¯ n B n τ κ ∇ µ δα B + L nπ ∇ ν δχ µν , (22) τ π D δχ µν + δχ µν = 2 τ η ∆ µναβ ∂ α δu β + L πn ∆ µναβ ∂ α δξ β . (23)All transport coefficients in the above equations are functions only of the temperature, since we assumed the chemicalpotential of the unperturbed system is zero. We also introduced several time scales associated to the transportcoefficients that appear in the fluid-dynamical equations, τ η ≡ ηε + P , τ κ ≡ κ n ¯ n B , L πn ≡ (cid:96) πn T , L nπ ≡ T (cid:96) nπ , (24)with ¯ n B being the baryon number density, not to be confused with the net-baryon number density. The first twoscales appear already in Navier-Stokes theory, while the remaining two are new scales related to the coupling termsinvestigated in this work. Naturally, there are also intrinsic time scales in Israel-Stewart theory – the relaxation times– which are required to render the theory causal and stable. It would be interesting to see whether the appearance ofsuch new time scales can affect the stability of the linearized theory (investigating the causality of the full nonlineartheory is a complex task, with recent progress on this topic being developed in Refs. [27, 28], but without the effectsof net-charge diffusion). Linearized equations of motion in Fourier space
It is practical to express these equations in Fourier space. We adopt the following convention for the Fouriertransformation ˜ M ( k µ ) = (cid:90) d x exp ( − ix µ k µ ) M ( x µ ) , (25) M ( x µ ) = (cid:90) d k (2 π ) exp ( ix µ k µ ) M ( k µ ) , (26)where k µ = ( ω, k ), with ω being the frequency and k the wavevector. It is convenient to introduce the covariantvariables, Ω ≡ u µ k µ , (27) κ µ ≡ ∆ µν k ν , (28)which correspond to the frequency and wavevector in the local rest frame of the unperturbed system. We furtherintroduce a covariant wavenumber, κ ≡ (cid:112) − κ µ κ µ . (29)Using this notation, the conservation laws are written asΩ δ ˜ εw + κ µ δ ˜ u µ = 0 , (30)Ω δ ˜ u µ − κ µ δ ˜ Pw + κ ν δ ˜ χ µν = 0 , (31)Ω δ ˜ n B n + κ µ δ ˜ ξ µ = 0 , (32)while the equations for the dissipative currents become( iτ n Ω + 1) δ ˜ ξ µ = iτ κ κ µ δ ˜ α B + i L nπ κ ν δ ˜ χ µν , (33)( iτ π Ω + 1) δ ˜ χ µν = 2 iτ η (cid:20) κ ( µ δ ˜ u ν ) −
13 ∆ µν κ λ δ ˜ u λ (cid:21) + i L πn (cid:20) κ ( µ δ ˜ ξ ν ) −
13 ∆ µν κ λ δ ˜ ξ λ (cid:21) , (34)where the parentheses in the indices denote the symmetrised tensor A ( µν ) ≡ ( A µν + A νµ ) / κ µ ( longitudinal degrees of freedom ) and those that are orthogonal to it ( transverse degrees of freedom ).In the linear regime, the transverse and longitudinal degrees of freedom are no longer coupled and can be solvedindependently. For instance, the decomposition of an arbitrary 4-vector orthogonal to u µ , A µ , can be implemented as A µ = A (cid:107) κ µ κ + A µ ⊥ , (35)where we define A (cid:107) ≡ − κ µ A µ /κ and A µ ⊥ ≡ ∆ µνκ A ν . Here,∆ µνκ ≡ g µν + κ µ κ ν κ − u µ u ν , (36)is the orthogonal projector onto the 2-space orthogonal to κ µ and u µ . A symmetric traceless rank two tensor orthogonalto u µ can be decomposed in a similar manner A µν = A (cid:107) κ µ κ ν κ + 12 A (cid:107) ∆ µνκ + A µ ⊥ κ ν κ + A ν ⊥ κ µ κ + A µν ⊥ , (37)where we define the projections A (cid:107) ≡ κ µ κ ν A µν /κ , A µ ⊥ ≡ − κ λ ∆ µνκ A λν /κ , and A µν ⊥ ≡ ∆ µναβκ A αβ , with∆ µναβκ ≡ (cid:0) ∆ µακ ∆ νβκ + ∆ µβκ ∆ νακ − ∆ µνκ ∆ αβκ (cid:1) . (38)The equations of motion satisfied by the longitudinal modes are obtained by contracting Eqs. (31) and (33) withthe tensor κ µ /κ and by contracting Eq. (34) with the tensor κ µ κ ν /κ . Note that Eqs. (30) and (32) are already interms of the longitudinal components of the fluctuations. This leads to the equationsΩ δ ˜ εw − κδ ˜ u (cid:107) = 0 , (39)Ω δ ˜ u (cid:107) − κ δ ˜ ε w − κδ ˜ χ (cid:107) = 0 , (40)Ω δ ˜ n B n − κδ ˜ ξ (cid:107) = 0 , (41) (cid:16) i ˆ τ n ˆΩ + 1 (cid:17) δ ˜ ξ (cid:107) + i ˆ L nπ ˆ κδ ˜ χ (cid:107) = i ˆ τ κ ˆ κ δ ˜ n B n , (42) (cid:16) i ˆ τ π ˆΩ + 1 (cid:17) δ ˜ χ (cid:107) − i ˆ L πn ˆ κδ ˜ ξ (cid:107) = 43 i ˆ κδ ˜ u (cid:107) , (43)Above, we expressed all dimensionful quantities in terms of the viscous time scale, τ η . That is, we defined,ˆΩ ≡ τ η Ω , ˆ κ ≡ τ η κ, ˆ τ n,π,κ ≡ τ n,π,κ τ η , (44)ˆ L nπ ≡ L nπ τ η , ˆ L πn ≡ L πn τ η . (45)In deriving the above equations, we have already made assumptions regarding the equation of state. We assume anequation of state of a gas composed solely of a massless particle and its corresponding anti-particle. In this case, theperturbation of pressure and chemical potential can be expressed as δ ˜ P = 13 δ ˜ ε, (46) δ ˜ α B = δ ˜ n B ¯ n B . (47)The transverse equations are obtained projecting Eqs. (31) and (33) with ∆ µλκ and Eq. (34) with ∆ µλκ κ ν . Then, weobtain the following set of equations ˆΩ δ ˜ u λ ⊥ − ˆ κδ ˜ χ λ ⊥ = 0 , (48)( i ˆ τ n ˆΩ + 1) δ ˜ ξ λ ⊥ + i ˆ L nπ ˆ κδ ˜ χ λ ⊥ = 0 , (49)( i ˆ τ π ˆΩ + 1) δ ˜ χ λ ⊥ − i ˆ κδ ˜ u λ ⊥ − i ˆ L πn κδ ˜ ξ λ ⊥ = 0 . (50)The equation for the fully transverse component of the shear-stress tensor, δ ˜ χ µν ⊥ , decouples from energy density andvelocity fluctuations and will not be considered in this analysis.When required, we use the transport coefficients calculated in Ref. [29] for an ultra-relativistic gas of hard spheres.In this case, one finds for the re-scaled transport coefficients the following values: ˆ τ π = 5, ˆ τ n = 27 /
4, ˆ τ κ = 9 /
16. Thecoefficients related to the coupling terms, L nπ and L πn , will not be fixed to a particular value. Nevertheless, we shallassume in this work, unless stated otherwise, that L nπ L nπ < . (51)This assumption is supported by kinetic theory calculations [14, 29, 30]. Furthermore, this constraint is obtained inthe phenomenological derivation of Israel-Stewart theory from the second law of thermodynamics [14, 31]. IV. STABILITY ANALYSIS WITHOUT COUPLING TERMS
We begin our analysis considering the case where coupling terms between shear-stress tensor and diffusion 4-currentare not present. This analysis was performed before, without the presence of diffusion fluctuations in Ref. [11].Preliminary calculations on the effect of diffusion fluctuations have been presented in Ref. [32].
A. Transverse modes
We start our discussion with the transverse degrees of freedom. In this case, the equations of motion, takingˆ L nπ = ˆ L πn = 0, can be expressed in the following form ˆΩ − ˆ κ
00 0 i ˆΩˆ τ n + 1 − i ˆ κ i ˆ τ π ˆΩ + 1 0 δ ˜ u µ ⊥ δ ˜ χ µ ⊥ δ ˜ ξ µ ⊥ = 0 . (52)The non-trivial solutions are obtained when the determinant vanishes, leading to the dispersion relation (cid:104)(cid:16) i ˆ τ π ˆΩ (cid:17) ˆΩ − i ˆ κ (cid:105) (cid:16) i ˆ τ n ˆΩ (cid:17) = 0 . Note that, in the absence of coupling terms, the transverse dispersion relations related to the shear-stress tensor andnet-baryon diffusion current decouple, and can be expressed as (cid:16) i ˆ τ π ˆΩ (cid:17) ˆΩ − i ˆ κ = 0 , (53)1 + i ˆ τ n ˆΩ = 0 . (54)Such decoupling of the modes would not necessarily occur for other choices of equation of state. Also, this will notoccur if we consider fluctuations around an equilibrium state with a finite net-baryon number density.Initially, we shall consider the case where the unperturbed fluid is at rest, i.e., u µ = (1 , , , ω and ˆ κ = ˆ k . In this case, the solutions can be expressed analytically asˆ ω diff T = i ˆ τ n , ˆ ω shear T, ± = i ± (cid:112) − τ π ˆ k τ π . (55)Clearly, the modes ˆ ω diff T and ˆ ω shear T, + are non-hydrodynamic, i.e., they do not vanish when the wavenumber is takento zero, while the mode ˆ ω shear T, − is hydrodynamic. Also, it is straightforward to see that, when taking the Navier-Stokes limit, i.e., when the relaxation times are sent to zero, ˆ τ π →
0, ˆ τ n →
0, the non-hydrodynamic modes goesto infinity, lim ˆ τ n → ˆ ω diff T ∼ i/ ˆ τ n and lim ˆ τ π → ˆ ω shear T, + ∼ i/ ˆ τ π , while the hydrodynamic mode goes to the usual Navier-Stokes solution, lim ˆ τ π → ˆ ω shear T, − = i ˆ k . We also see that the shear modes become propagating, i.e., they have a realcomponent, for wavenumbers that are sufficiently large, ˆ k > / (cid:0) √ ˆ τ π (cid:1) . All these results for the shear modes werepreviously obtained in Ref. [11].Since they are the only ones that carry any dependence on the wavenumber ˆ k , it is useful to look at the modesˆ ω shear T, ± in the small, ˆ k (cid:28)
1, and large, ˆ k (cid:29)
1, wavenumber limits. In the first case, we haveˆ ω shear T, + = i ˆ τ π − i ˆ k + O (cid:16) ˆ k (cid:17) , (56)ˆ ω shear T, − = i ˆ k + i ˆ τ π ˆ k + O (cid:16) ˆ k (cid:17) , (57)while, in the second case, we obtainˆ ω shear T, ± = i τ π ± (cid:32) ˆ k √ ˆ τ π − τ / π ˆ k (cid:33) + O (cid:18) k (cid:19) , (58)We can see that the modes obtained above are stable as long as the relaxation times are positiveˆ τ n ≥ , ˆ τ π ≥ . (59)Furthermore, since the shear modes become propagating when ˆ k (cid:29)
1, causality imposes the following constraint tothe asymptotic group velocity [11, 33] lim ˆ k →∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ Re(ˆ ω ) ∂ ˆ k (cid:12)(cid:12)(cid:12)(cid:12) ≤ ⇒ ˆ τ π ≥ . (60)The modes ˆ ω shear T, ± and ˆ ω diff T are plotted in Fig. 1, considering the relaxation times calculated from the Boltzmannequation, using the 14-moment approximation in the ultra-relativistic limit, i.e., ˆ τ π = 5 and ˆ τ n = 27 / τ π and ˆ τ n are employed. The transverse modes related to theshear-stress tensor are usually referred to as shear modes. These modes have been calculated before in Ref. [11] andare identical to the results presented here, since the inclusion of net-baryon diffusion does not affect them.We now analyse the case in which the unperturbed fluid is moving. Such problem involves intrinsically relativisticphenomena, since velocities that are close to the velocity of light are possible even in such perturbative scenario. Wedefine the Lorentz factor, γ = 1 / √ − V , and, for the sake of simplicity, assume that the 3-velocity and the wavevectorare taken in the same direction, e.g., the x –axis. This case corresponds to u µ = γ (1 , V, ,
0) and k µ = ( ω, k, , γ (ˆ ω − V ˆ k ) , (61)ˆ κ = γ (ˆ ωV − ˆ k ) . (62)In this case, the dispersion relations given by Eqs. (53) and (54), become i ˆ τ π ( γ ˆ ω − γV ˆ k ) + ( γ ˆ ω − γV ˆ k ) − i ( γ ˆ ωV − γ ˆ k ) = 0 , (63)1 + iγ ˆ τ n (ˆ ω − V ˆ k ) = 0 . (64)The solution for the diffusion transverse mode is, once more, the simplest one,ˆ ω diff T = V ˆ k + iγ ˆ τ n , (65) k I m ( ω ) - - k R e ( ω ) FIG. 1: The imaginary and real parts of the transverse modes in the absence of coupling terms, i.e., ˆ L nπ ˆ L nπ = 0. given by a propagating part, with velocity identical to that of the unperturbed fluid, and a non-propagating part thatrelaxes to equilibrium within a Lorentz dilated relaxation time, γ ˆ τ n (this corresponds to a non-hydrodynamic mode).This mode is always stable as long as the diffusion relaxation time is positive. Furthermore, as expected, causalityis satisfied as long as the background fluid has a velocity smaller than the speed of light. The dispersion relationsatisfied by the remaining transverse modes is, i (cid:0) ˆ τ π − V (cid:1) ( γ ˆ ω ) + (cid:104) − i (ˆ τ π − γV ˆ k (cid:105) γ ˆ ω − γV ˆ k + i (ˆ τ π V − γ ˆ k ) = 0 , (66)In order to understand the stability of such modes, we first look at their behavior when ˆ k = 0, i.e., in the homogeneouslimit. In this limiting case, the hydrodynamic mode will naturally disappear, while the non-hydrodynamic modeachieves the following value ˆ ω shear T, + (ˆ k = 0) = iγ (ˆ τ π − V ) . (67)We note that this mode does not vanish when the shear relaxation time is taken to zero. This implies that suchnon-hydrodynamic mode also appears in Navier-Stokes theory.The non-hydrodynamic mode should be stable for all possible values of the background velocity V . In order toguarantee this at least for k = 0, the shear relaxation time must satisfy the following conditionˆ τ π ≥ ⇒ τ π ≥ τ η . (68)This constraint is identical to the causality condition obtained for the same modes when the unperturbed fluid wasat rest, see Eq. (60).In Fig. 2, we plot the solutions of Eq. (66) for V = 0 . V = 0 .
4, and V = 0 .
9, using the same values of relaxationtimes as before: ˆ τ π = 5 and ˆ τ n = 27 /
4. In this case, we see that the modes are indeed stable for all values of ˆ k . Anexample of an unstable fluid configuration is shown in Fig. 3, by taking ˆ τ π = 0 . τ n = 27 /
4, values that do notsatisfy the causality and stability conditions derived in this section for the transverse modes.
B. Longitudinal modes
We now discuss the longitudinal projections of the fluid-dynamical equations. The equations of motion for thelongitudinal degrees of freedom, taking ˆ L nπ = ˆ L nπ = 0, can be expressed in the following matrix form Ω 0 0 − κ
00 Ω − κ − κ Ω 0 − κ − i ˆ τ κ ˆ κ i ˆ τ n ˆΩ + 1 00 0 − i ˆ κ i ˆ τ π ˆΩ + 1 δ ˜ n B /n δ ˜ ε/w δ ˜ u (cid:107) δ ˜ ξ (cid:107) δ ˜ χ (cid:107) = 0 . (69) k I m ( ω ) V = k I m ( ω ) V = k I m ( ω ) V = - - k R e ( ω ) V = k R e ( ω ) V = k R e ( ω ) V = FIG. 2: The imaginary and real parts of the transverse modes for perturbations around a moving fluid, considering V = 0 . V = 0 .
4, and V = 0 .
9, in the absence of coupling terms, i.e., ˆ L nπ ˆ L nπ = 0. - - - k I m ( ω ) FIG. 3: The imaginary part of the unstable shear mode for ˆ τ π = 0 . V = 0 .
9, in the absence ofcoupling terms, i.e., ˆ L nπ ˆ L nπ = 0. Non-trivial solutions for these equations are obtained when the determinant is zero, leading to the dispersion relation (cid:20)(cid:18) ˆΩ −
13 ˆ κ (cid:19) ( i ˆ τ π ˆΩ + 1) − i ˆ κ ˆΩ (cid:21) (cid:104) ˆΩ( i ˆ τ n ˆΩ + 1) − i ˆ τ κ ˆ κ (cid:105) = 0 . (70)Once more, the dispersion relation related to net-baryon current perturbations decouples from those related to energy-momentum tensor fluctuations and can be solved independently,ˆΩ( i ˆ τ n ˆΩ + 1) − i ˆ τ κ ˆ κ = 0 , (71) (cid:18) ˆΩ −
13 ˆ κ (cid:19) ( i ˆ τ π ˆΩ + 1) − i ˆ κ ˆΩ = 0 . (72)As before, we first consider the case where the unperturbed fluid is at rest, i.e., u µ = (1 , , , ω and ˆ κ = ˆ k . In this case, the solution for the modes related to the net-baryon diffusion current fluctuations can be0cast in a simple form ω B L, ± = i ± (cid:112) − τ n ˆ τ κ ˆ k τ n . (73)The mode ω B L, − is clearly hydrodynamic, as can be seen from the small wavenumber limit ω B L, − = i ˆ τ κ ˆ k + i ˆ τ n ˆ τ κ ˆ k + O (cid:16) ˆ k (cid:17) , (74)with the leading term being the dispersion related usually obtained in Navier-Stoke theory. On the other hand, ω B L, + is a non-hydrodynamic mode that simply does not exist in Navier-Stokes theory, ω B L, + = i ˆ τ n − i ˆ τ κ ˆ k + O (cid:16) ˆ k (cid:17) . (75)For asymptotically large values of wavenumber, these modes become ω B L, ± = i τ n ± (cid:114) ˆ τ κ ˆ τ n (cid:18) ˆ k − τ n ˆ τ κ ˆ k (cid:19) + O (cid:18) k (cid:19) . (76)The modes ω diff L, ± become propagating when ˆ k > / (cid:0) √ ˆ τ n ˆ τ κ (cid:1) and causality dictates that [11, 33]lim k →∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ Re( ω ) ∂k (cid:12)(cid:12)(cid:12)(cid:12) ≤ ⇒ ˆ τ n ≥ ˆ τ κ . (77)The remaining longitudinal modes (which include the sound modes) are a solution of a cubic equation and, thus,their analytical solution cannot be expressed in a simple form. For the sake of completeness, we plot the imaginaryand real parts of the longitudinal modes in Fig. 4, taking ˆ τ π = 5, ˆ τ n = 27 / τ κ = 9 /
16. In the following, werestrict our discussion to the asymptotic form of these modes. We take a look at the behavior of such modes for smallvalues of wavenumber, ω sound ± = ±
13 ˆ k + 23 i ˆ k + O (cid:16) ˆ k (cid:17) , (78) ω shear L = i ˆ τ π − i ˆ k + O (cid:16) ˆ k (cid:17) , (79)while at large values of wavenumber, the same modes behave as ω sound ± = ± (cid:114) τ π τ π ˆ k + 2 i ˆ τ π (4 + ˆ τ π ) + O (cid:18) k (cid:19) , (80) ω shear L = i τ π + O (cid:18) k (cid:19) . (81)Clearly, the modes ˆ ω sound ± are hydrodynamic (describing the propagation of sound waves) and, as expected, theyreduce to the modes obtained from Navier-Stokes theory when the wavenumber is sufficiently small. The other modeis non-hydrodynamic and does not exist in Navier-Stokes theory. The stability of these modes is guaranteed for smallor large values of wavenumber as long as ˆ τ π >
0. Furthermore, causality imposes that the asymptotic group velocitymust satisfy [11, 33] lim k →∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ Re( ω ) ∂k (cid:12)(cid:12)(cid:12)(cid:12) ≤ ⇒ ˆ τ π ≥ . (82)We now perform the same analysis considering that the unperturbed fluid is moving. As already done in theprevious section, we assume that the background four-velocity and the perturbations are in the same direction, e.g.,the x axis. Therefore, we have u µ = γ (1 , V, ,
0) and k µ = ( ω, k, , (cid:16) γ ˆ ω − γV ˆ k (cid:17) (cid:104) i ˆ τ n (cid:16) γ ˆ ω − γV ˆ k (cid:17) + 1 (cid:105) − i ˆ τ κ (cid:16) ˆ ωγV − γ ˆ k (cid:17) = 0 , (83) (cid:20)(cid:16) γ ˆ ω − γV ˆ k (cid:17) − (cid:16) ˆ ωγV − γ ˆ k (cid:17) (cid:21) (cid:104) i ˆ τ π (cid:16) γ ˆ ω − γV ˆ k (cid:17) + 1 (cid:105) − i (cid:16) ˆ ωγV − γ ˆ k (cid:17) (cid:16) γ ˆ ω − γV ˆ k (cid:17) = 0 . (84)1 k I m ( ω ) - k R e ( ω ) FIG. 4: Real and imaginary part of the longitudinal modes in the absence of coupling terms, ˆ L nπ ˆ L nπ = 0, for a staticbackground, V = 0. Since analysing the stability of the modes for any value of wavenumber would be extremely complicated, we shallinitially restrict our focus to the vanishing wavenumber limit, as before. In this limit, we have that the longitudinalmodes related to net-baryon diffusion fluctuations can be cast in the form ω B L, − = 0 , ω B L, + = iγ (ˆ τ n − ˆ τ κ V ) , (85)while the remaining longitudinal modes become ω sound ± = 0 , ω shear0 = i (cid:0) − V (cid:1) γ [3ˆ τ π − (ˆ τ π + 4) V ] . (86)It is essential that the modes are stable at k = 0. Furthermore, the perturbations must be stable for any value ofthe background fluid velocity. In order to ensure that the imaginary part of the mode has the correct sign (that leadsto modes that are damped at late times) for any background velocity, we must have that,ˆ τ n ≥ ˆ τ κ , ˆ τ π ≥ , (87)which are the same conditions that are obtained when imposing causality for perturbations on top of a fluid thatis at rest, see Eqs. (77) and (82). These stability conditions are stronger than the ones found using the transversemodes and, therefore, are sufficient to ensure that the system is linearly causal and stable. The condition for the shearrelaxation time was previously obtained in [11]. The stability condition that was obtained for the diffusion relaxationtime is, to the best of our knowledge, new.The modes that are obtained as the solutions of Eq. (70) for a moving background fluid are displayed in Fig. 5,considering ˆ τ π = 5, ˆ τ n = 27 /
4, and ˆ τ κ = 9 /
16. One can easily see that all modes are stable for all values of ˆ k . Further,in Fig. 6, we consider two examples of unstable fluid configurations, i.e., fluids that do not satisfy the conditionsderived in Eq. (87). On the left panel we analyse an unstable case in which ˆ τ n < ˆ τ κ . Here, we considered ˆ τ n = 3 / τ κ = 9 /
16 for an unperturbed system with velocity V = 0 .
9. In this scenario, the longitudinal non-hydrodynamicmode related to net-baryon current fluctuations is unstable. On the right panel, we analyse the case where ˆ τ π = 0 . V = 0 .
9. Again, there is the occurrence of an unstable non-hydrodynamicmode, related to fluctuations of the shear-stress tensor.So far, we have derived the stability conditions of the Israel-Stewart equations in the absence of coupling terms.These conditions, which are constraints for the relaxation times, were shown to be identical to the causality conditionsobtained in the case where the unperturbed system is at rest. We concluded that the stability and causality conditionsimposed by analysing longitudinal modes are stronger and supersede the ones obtained investigating transverse modes.In the next section, we will be investigating the stability conditions in the presence of the coupling terms.2 k I m ( ω ) V = k I m ( ω ) V = k I m ( ω ) V = - k R e ( ω ) V = - - - k R e ( ω ) V = k R e ( ω ) V = FIG. 5: The imaginary and real parts of the longitudinal modes for a moving background fluid, for V = 0 . V = 0 .
4, and V = 0 . L nπ ˆ L nπ = 0. - - - k I m ( ω ) - - - - - k I m ( ω ) FIG. 6: Imaginary part of the unstable longitudinal mode related to baryon-number fluctuations for ˆ τ κ = 9 /
16 and ˆ τ n = 3 / V = 0 . τ π = 0 . V = 0 . V. STABILITY ANALYSIS WITH COUPLING TERMS
We now consider the case where coupling terms between shear-stress tensor and diffusion diffusion 4-current arepresent. The following stability analysis will be very similar to the one performed in the previous section.3
A. Transverse modes
As before, we start our discussion with the linearized equations of motion satisfied by the transverse degrees offreedom. The equations of motion can be cast in the following form ˆΩ − ˆ κ i ˆ L nπ ˆ κ i ˆΩˆ τ n + 1 − i ˆ κ i ˆ τ π ˆΩ + 1 − i ˆ L πn ˆ κ δ ˜ u µ ⊥ δ ˜ χ µ ⊥ δ ˜ ξ µ ⊥ = 0 . (88)These equations lead to the following dispersion relation, obtained from the determinant of the matrix above − A ˆΩ + i B ˆΩ + (cid:0) C ˆ κ (cid:1) ˆΩ − i ˆ κ = 0 , (89)where we defined the quantities A ≡ ˆ τ π ˆ τ n , B ≡ ˆ τ n + ˆ τ π , (90) C ≡ ˆ τ n −
12 ˆ L πn ˆ L nπ . (91)Naturally, we can see that the addition of coupling terms in the Israel-Stewart equations does not add a new modeto the dispersion relation for the transverse modes, since it continues to be a third-order polynomial. However, nowthat the coupling terms are included, the modes related to energy-momentum tensor fluctuations and net-baryoncurrent fluctuations no longer factorize. This renders the solution of the problem considerably more complicated andanalytical solutions can no longer be cast in a simple form. Furthermore, we also see that the modes only depend onthe product of the coupling terms, ˆ L πn ˆ L nπ , that is contained within our definitions of variables in the parameter C .Note that, if the relaxation times are set to zero, one no longer recovers a simple Navier-Stokes dispersion relation.Instead, one obtains the solution ˆΩ = i ˆ κ − ˆ L πn ˆ L nπ ˆ κ . (92)This solution is purely imaginary and is only stable if ˆ L πn ˆ L nπ is negative. If the relaxation times are not zero, we willsee that stable solutions with positive values of ˆ L πn ˆ L nπ are still possible, even though we will not investigate suchcases thoroughly. ℒ n π ℒ π n = ℒ n π ℒ π n = - ℒ n π ℒ π n = - - - k I m ( ω ) FIG. 7: The imaginary part of the dispersion relation of the transverse modes for an unperturbed system at rest, V = 0, forthree different values of the product of the coupling terms ˆ L nπ ˆ L nπ = − , ,
1. In this plot the relaxation times are set to zero,ˆ τ π = ˆ τ n = 0. We plot this mode assuming a static background (Ω = ω and κ = k ) in Fig. 7 for the following choices of couplingˆ L πn ˆ L nπ = − , , and 1. Even though the modes show a familiar behavior in the small-wavenumber region, as thewavenumber increases, they become rather different. Positive values of the effective coupling term render the theory4unstable even in the case where the unperturbed system is at rest. On the other hand, for negative values, theimaginary part of the mode is positive-definite, and therefore, it is stable.In the following stability analysis, all the relaxation times are once again taken into account. As before, we initiallyconsider the case where the unperturbed system is at rest, u µ = (1 , , , − A ˆ ω + i B ˆ ω + (cid:16) C ˆ k (cid:17) ˆ ω − i ˆ k = 0 . (93)We first study these solutions in two different limits: for small and large wavenumber.In the small wavenumber limit, ˆ k (cid:28)
1, one obtainsˆ ω diff T = i ˆ τ n + i ˆ L πn ˆ L nπ τ π − ˆ τ n ) ˆ k + O (cid:16) ˆ k (cid:17) , (94)ˆ ω shear T, + = i ˆ τ π + i τ n − ˆ τ π ) − ˆ L πn ˆ L nπ ˆ τ π − ˆ τ n ˆ k + O (cid:16) ˆ k (cid:17) , (95)ˆ ω shear T, − = i ˆ k + O (cid:16) ˆ k (cid:17) . (96)We see that, when either ˆ L πn or ˆ L nπ is set to zero, we recover the result from the previous section. It is also interestingto notice that the coupling terms lead to corrections that are of higher order in ˆ k .In the limit of large wavenumbers, one can demonstrate that an expansion of the following form existsˆ ω = c − ˆ k + ∞ (cid:88) n =0 c n ˆ k − n . (97)The expansion coefficients can be identified order by order by replacing the above ansatz into the dispersion relation.We then obtain the following solutionsˆ ω shear T, ± = ± (cid:114) CA ˆ k + i BC − A AC + O (cid:18) k (cid:19) , (98)ˆ ω diff T = i C + O (cid:18) k (cid:19) . (99)The solutions derived above, obtained from the asymptotic expansion of ˆ ω , already carry information on the linearstability of the theory. In order for the system to be stable, we must have that (i) (cid:112) C / A is real, (ii) C ≥
0, and (iii)
BC − A ≥
0. Note that condition (i) automatically guarantees condition (ii), since A is guaranteed to be positive aslong as the relaxation times are also positive. The stability conditions (ii) and (iii) lead to, C > ⇒ ˆ L πn ˆ L nπ < τ n , (100) BC − A > ⇒ ˆ L πn ˆ L nπ < τ n ˆ τ n + ˆ τ π . (101)These conditions can be shown to be equivalent to those obtained using the Routh-Hurwitz stability criterion [34–36].In this sense, they are more general, and guarantee the stability conditions for any value of wavenumber k (and notjust for asymptotically large values of k ). The last condition is clearly stronger and imposes restrictions on the valuesthat the coupling terms coefficients can have. It is interesting that, with the inclusion of the coupling terms, we alreadyhave to impose stability conditions even for perturbations around a background that is at rest – something that wasnot required in Navier-Stokes theory or in simplified versions of Israel-Stewart theory. Note that both conditions areautomatically satisfied if the product of the coupling terms has a negative sign. Furthermore, a causality conditioncan be extracted from the expansion of the modes in the large wavenumber limit,lim k →∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ Re( ω ) ∂k (cid:12)(cid:12)(cid:12)(cid:12) = CA ≤ ⇒ ˆ L πn ˆ L nπ ≥ − τ n (ˆ τ π − . (102)For the sake of illustration, we display, in Fig. 8, the real and imaginary parts of the modes ˆ ω shear T, ± and ˆ ω diff T for thefollowing negative values of the product of the coupling terms, ˆ L πn ˆ L nπ = − . , − , −
4. In Fig. 9, we display suchmodes for positive values of ˆ L πn ˆ L nπ , chosen to be ˆ L πn ˆ L nπ = 0 . , ,
6. As before, these plots were made using the5transport coefficients calculated from the Boltzmann equation, ˆ τ π = 5 and ˆ τ n = 27 /
4. For these values of transportcoefficients, the stability condition given by Eq. (101) becomes, ˆ L πn ˆ L nπ (cid:46) .
75. We note that it is possible to obtainstable modes even for positive values of ˆ L πn ˆ L nπ . As the value of ˆ L πn ˆ L nπ becomes negative, the coupling terms renderthe imaginary parts of the non-hydrodynamic modes degenerate at larger values of wavenumbers. On the other hand,if ˆ L πn ˆ L nπ is positive, the non-hydrodynamic mode related to the diffusion 4-current, ˆ ω diff T , becomes degenerate withthe hydrodynamic mode ˆ ω shear T, − instead, when the wavenumber increases. In Fig. 10, we show a case in which therelaxation times are chosen in such a way to ensure the causality and stability conditions derived in the previoussections, but the modes are driven unstable by the coupling term, with ˆ L πn ˆ L nπ = 10. k I m ( ω ) ℒ n π ℒ π n = - k I m ( ω ) ℒ n π ℒ π n = - k I m ( ω ) ℒ n π ℒ π n = - - - k R e ( ω ) - - k R e ( ω ) - - k R e ( ω ) FIG. 8: Real and imaginary parts of the transverse modes for three negative values of the product of the coupling terms,ˆ L πn ˆ L nπ = − . , − , −
4, for a static background, V = 0. We now perform the same analysis considering that the unperturbed fluid is moving with a velocity V in the, e.g. x –direction. Once again, we assume that the four-velocity and the perturbations are taken in the same direction, sothat u µ = γ (1 , V, ,
0) and k µ = ( ω, k, , ω and k . For the sake of illustration,the solutions of the resulting dispersion relation are plotted for two different cases: in Fig. 11, for negative valuesof ˆ L πn ˆ L nπ , and in Fig. 12, for positive values of ˆ L πn ˆ L nπ considering the following values of the background fluidvelocity, V = 0 . V = 0 . V = 0 . k = 0. This simplifies considerably the calculations and allows us to provide basicnecessary conditions for linear stability. In this limit, the dispersion relation reduces to − A ( γ ˆ ω ) + i B ( γ ˆ ω ) + (cid:104) C ( γ ˆ ωV ) (cid:105) ( γ ˆ ω ) − i ( γ ˆ ωV ) = 0 . Thus, it is possible to express the solutions in a simple analytical formˆ ω = 0 , γ ˆ ω ± = i (cid:0) B − V (cid:1) ± (cid:113) ( B − V ) − A − C V ) A − C V . (103)Naturally, the vanishing solution corresponds to the hydrodynamic mode, while the remaining solution correspondsto the non-hydrodynamic mode. The stability of the perturbations is usually governed by the non-hydrodynamicmode, which will either decay or increase exponentially depending on the choices of transport coefficients. If V = 0,we recover the previous solutions, ˆ ω ± = i/ ˆ τ π , i/ ˆ τ n .6 k I m ( ω ) ℒ n π ℒ π n = k I m ( ω ) ℒ n π ℒ π n = k I m ( ω ) ℒ n π ℒ π n = - - k R e ( ω ) - - k R e ( ω ) - - - k R e ( ω ) FIG. 9: Real and imaginary parts of the transverse modes for three positive values of the product of the coupling terms,ˆ L πn ˆ L nπ = 0 . , ,
6, for a static background, V = 0. k I m ( ω ) ℒ n π ℒ π n = - - k R e ( ω ) FIG. 10: Real and imaginary parts of the transverse modes for an unstable value for the product of the coupling terms,ˆ L πn ˆ L nπ = 10, for a static background, V = 0. Next, we investigate the necessary conditions the transport coefficients must satisfy in order for the non-hydrodynamic modes ˆ ω ± in Eq. (103) to be both stable in the homogeneous limit. Naturally, in order for themodes to have a positive imaginary part at k = 0, the numerator and denominator must carry the same sign. Thedenominator should not change sign as the velocity of the background fluid increases, otherwise resulting in an in-stability. Therefore, stable fluid-dynamical formulations must satisfy A > C V for all values of 0 ≤ V ≤
1, which isguaranteed by the following condition A > C . (104)The imaginary part of the numerator must also be positive. This condition requires, at least, that B − V > ≤ V ≤
1, leading to B > . (105)Finally, we study the term inside the square root in Eq. (103). One can show that for negative values of the productof the coupling terms, ˆ L nπ ˆ L πn ≤
0, which is the case considered here, such term is positive and smaller than
B − V ,7and does not affect the sign of the modes. Furthermore, considering also positive values for the product between thecoupling terms, ˆ L nπ ˆ L πn ≥
0, either the term inside the square root remains positive and smaller than
B − V or itbecomes negative, resulting in an imaginary value for the square root. Either way, the conditions given by Eqs. (104)and (105) are sufficient to maintain the theory linearly stable at k = 0. The conditions (104) and (105) can beexpressed as the following constraints for the transport coefficients B > ⇒ ˆ τ n + ˆ τ π > , (106) A > C = ⇒ ˆ L πn ˆ L nπ ≥ − τ n (ˆ τ π − . (107)The latter is equivalent to the causality condition obtained for perturbations around a background at rest, seeEq. (102). Combining the stability conditions above with those derived for a background at rest, see Eqs. (100) and(101), we then obtain ˆ τ n + ˆ τ π > , (108) − τ n (ˆ τ π − ≤ ˆ L πn ˆ L nπ < τ n ˆ τ n + ˆ τ π . (109)These inequalities can be further simplified by imposing that the product between the coupling terms is negative,ˆ L πn ˆ L nπ ≤
0. In this case, the last relation reduces to (cid:12)(cid:12)(cid:12) ˆ L πn ˆ L nπ (cid:12)(cid:12)(cid:12) ≤ τ n (ˆ τ π − , ˆ τ π ≥ , ˆ τ n ≥ L πn ˆ L nπ can also lead to stable theories. In this case, itwould even be possible to violate the condition ˆ τ π ≥
1. As a matter of fact, if ˆ τ π ≤
1, we see that only positive valuesof ˆ L πn ˆ L nπ are allowed.Then, removing the scaling factors from the transport coefficients, we obtain the following condition | (cid:96) πn (cid:96) nπ | ≤ τ n ( τ π − τ η ) . (111)This condition is satisfied in calculations from the Boltzmann equation [29, 30]. So far, we are not aware of anymicroscopic calculations that does not satisfy this condition. Longitudinal modes
In the presence of coupling terms, the linearized equations for the longitudinal modes can be expressed in thefollowing form
Ω 0 0 − κ
00 Ω − κ − κ Ω 0 − κ − i ˆ τ κ ˆ κ i ˆ τ n ˆΩ + 1 i ˆ L nπ ˆ κ − i ˆ κ − i ˆ L πn ˆ κ i ˆ τ π ˆΩ + 1 δ ˜ n B /n δ ˜ ε/w δ ˜ u (cid:107) δ ˜ ξ (cid:107) δ ˜ χ (cid:107) = 0 (112)We then obtain the following dispersion relation (cid:20)(cid:18) ˆΩ −
13 ˆ κ (cid:19) ( i ˆ τ π ˆΩ + 1) − i ˆ κ ˆΩ (cid:21) (cid:104) ˆΩ( i ˆ τ n ˆΩ + 1) − i ˆ τ κ ˆ κ (cid:105) −
23 ˆ L πn ˆ L nπ (cid:18) ˆΩ −
13 ˆ κ (cid:19) ˆΩˆ κ = 0 . (113)For the sake of convenience, we rewrite this expression as − A ˆΩ + i B ˆΩ + (cid:0) AS ˆ κ (cid:1) ˆΩ − i BD κ ˆΩ − (cid:0) E ˆ κ (cid:1) ˆΩˆ κ + i ˆ τ κ κ = 0 , (114)where A , B , and C were defined in the previous section, in Eqs. (90) and (91), and we further introduced the variables S ≡ A + 3ˆ τ π ˆ τ κ + 4 C A , (115) D ≡ B + 3ˆ τ κ + 4 B , (116) E ≡ τ κ + ˆ τ π ˆ τ κ + 43 ( C − ˆ τ n ) . (117)8 k I m ( ω ) V = k I m ( ω ) V = k I m ( ω ) V = - - k R e ( ω ) k R e ( ω ) k R e ( ω ) FIG. 11: Real and imaginary parts of the transverse modes considering a negative value for the product of the coupling terms,ˆ L πn ˆ L nπ = −
1, for three different values of the background velocity V = 0 . V = 0 .
4, and V = 0 . Other useful definitions that will be employed in the remaining of this paper are
M ≡ E A , R ≡ (cid:112) S − M . (118)Initially, we shall consider the case where the unperturbed system is at rest, i.e., u µ = (1 , , , − A ˆ ω + i B ˆ ω + (cid:16) AS ˆ k (cid:17) ˆ ω − i BD ˆ k ˆ ω − (cid:16) E ˆ k (cid:17) ˆ ω ˆ k + i ˆ τ κ k = 0 , (119)The solutions of this polynomial equation are rather complicated. As before, we shall look at such solutions in thelimits of small and large wavenumber in order to extract some information about the behavior of the modes. In thiscase, for small wavenumbers, ˆ k (cid:28)
1, we have that ω B L, + = i ˆ τ n + O (cid:16) ˆ k (cid:17) , (120) ω shear L = i ˆ τ π + O (cid:16) ˆ k (cid:17) , (121) ω B L, − = i ˆ τ κ ˆ k + O (cid:16) ˆ k (cid:17) , (122) ω sound ± = ± √ k + 23 i ˆ k + O (cid:16) ˆ k (cid:17) . (123)We have two non-hydrodynamic modes and three hydrodynamic modes. At small wavenumbers, the hydrodynamicmodes behave like the modes found in Navier-Stokes theory. Furthermore, for large wavenumber, the modes can bewritten as ˆ ω = ± ˆ k √S ± R + i A B ( S ± R ) − BD ( S ± R ) + ˆ τ κ
15 (
S ± R ) − S ( S ± R ) + 3 M + O (cid:18) k (cid:19) , (124)ˆ ω = i ˆ τ κ E + O (cid:18) k (cid:19) . (125)9 k I m ( ω ) V = k I m ( ω ) V = k I m ( ω ) V = - k R e ( ω ) k R e ( ω ) k R e ( ω ) FIG. 12: Real and imaginary parts of the transverse modes considering a positive value for the product of the coupling terms,ˆ L πn ˆ L nπ = 2, for three different values of the background velocity V = 0 . V = 0 .
4, and V = 0 . In order for these modes to be stable, it is required that their imaginary parts are positive and thus we have thefollowing necessary conditions:(i) ˆ τ κ / E > √S ± R is real;(iii) (cid:104) B ( S ± R ) − BD ( S ± R ) + ˆ τ κ (cid:105) / (cid:104)
15 (
S ± R ) − S ( S ± R ) + 3 M (cid:105) is positive.The causality condition further imposes the following constraints for the asymptotic group velocitylim k →∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ Re( ω ) ∂k (cid:12)(cid:12)(cid:12)(cid:12) = √S ± R ≤ . (126)In the following, we translate the conditions above into conditions for the transport coefficients.Condition (i) is guaranteed as long as E >
0, which leads to another constraint for the product of the couplingcoupling terms, given by ˆ L πn ˆ L nπ <
32 ˆ τ κ (ˆ τ π + 4) , (127)which is automatically satisfied if ˆ L πn ˆ L nπ <
0. Condition (ii) is satisfied if
S ≥ R (this condition is automaticallysatisfied if the product of the coupling terms is negative) and if R is real. The latter condition further implies that S > M . (128)This condition will be further developed when we discuss the modes for perturbations on top of a moving backgroundfluid. Finally, condition (iii) is satisfied as long as3 B ( S + R ) − BD ( S + R ) + ˆ τ κ > , (129)3 B ( S − R ) − BD ( S − R ) + ˆ τ κ < . (130)0If the product of the coupling terms is negative, these conditions are automatically satisfied and do not provide anynew constraints for the stability of the system. On the other hand, if the product of the coupling terms is positive,these conditions will lead to constraints for the transport coefficients, but we shall not discuss or explicitly derivethem in this work. Furthermore, all conditions listed above can be shown to be equivalent to the ones found by usingthe Routh-Hurwitz criterion [34–36]. In this sense, they are valid not only in the large wavenumber regime, but alsofor any values of wavenumber k . k I m ( ω ) ℒ n π ℒ π n = - - k R e ( ω ) k I m ( ω ) ℒ n π ℒ π n = - k R e ( ω ) FIG. 13: Real and imaginary parts of the longitudinal modes considering a negative and a positive values for the product ofthe coupling terms, ˆ L πn ˆ L nπ = − , V = 0. The solutions of Eq. (113) for a static background are displayed in Fig. 13 for a negative value of the coupling termˆ L πn ˆ L nπ = − L πn ˆ L nπ = 2 (bottom panels). The inclusion of the couplingproduces a similar behavior when compared to the transverse modes, see Figs. 8 and 9, where the sign of the productof the coupling terms dictates which modes merge at large values of wavenumber. Similarly to what was observed forthe transverse modes, the imaginary part of the longitudinal modes become constant at large wavenumber. The realparts of the longitudinal modes do not show any qualitative variation at large wavenumber as we change the sign ofthe coupling term.Now we discuss perturbations on top of a moving background fluid. Again, we assume that the four-velocityand the perturbations are taken in the same direction, e.g., the x axis. Therefore, we have u µ = γ (1 , V, ,
0) and k µ = ( ω, k, , κ given in Eqs. (61) and (62). In Figs. 14 and 15,we display the solutions of the dispersion relation in Eq. (113), substituting Eqs. (61) and (62), considering severalvalues of velocity and product of the coupling terms. In the following, we shall study the properties of the solutionsat vanishing wavenumber, ˆ k = 0, in order to extract a set of necessary linear stability conditions, as already done inthe previous sections. In this case, the dispersion relation reads( γ ˆ ω ) (cid:2) − A (cid:0) − S V + M V (cid:1) γ ˆ ω + i (cid:0) ˆ τ κ V + 3 B − B V (cid:1) γ ˆ ω + 3 − V (cid:3) = 0 , (131)1with 2 non-vanishing solutions, corresponding to the two longitudinal non-hydrodynamic modes of the theory, γ ˆ ω ± = i ˆ τ κ V + B (cid:0) − D V (cid:1) ± (cid:113) [ˆ τ κ V + B (3 − D V )] − A (3 − V ) (1 − S V + M V )3 A (1 − S V + M V ) . (132)As already discussed, in order to guarantee the stability of these modes, we must have that the imaginary part ofthe dispersion relation is positive for all possible values of the velocity. Evidently, this is achieved if both numeratorand denominator carry the same sign. The denominator in Eq. (132) is positive for V = 0, and thus it must remainpositive for all values of V – a change of sign would already imply in an instability. Therefore, we must have1 − S V + M V > , ∀ ≤ V ≤ . (133)Here, Eq. (133) is a polynomial function which is quadratic in V and positive at V = 0. The condition above issatisfied if the smallest root of this polynomial is larger than 1. This guarantees that the function is positive in thephysical interval of 0 ≤ V ≤ S − R ≥ M . (134)Note that this inequality leads to the causality condition derived in Eq. (126) for perturbations on top of a fluid atrest. This can be seen by taking the definition of M = S − R and then S − R ≥ ( S − R )( S + R ) = ⇒ S + R ≤ ⇒ −
32 ˆ τ π (ˆ τ n − ˆ τ κ ) − ˆ τ n (ˆ τ π − ≤ ˆ L πn ˆ L nπ , (135)where we used the previously derived stability conditions S − R > R ≥
0. Once again, we see that the stabilityconditions obtained for perturbations on top of a moving background are related to the causality conditions satisfiedby perturbations on top of a static background. The relations in Eq. (135) can be used to derive the additionalrelations M = ( S + R )( S − R ) ≤ ⇒
32 ˆ τ κ (ˆ τ π + 4) − A ≤ ˆ L πn ˆ L nπ , (136) S ≤ M ⇒ −
23 ˆ L πn ˆ L nπ ≤ (ˆ τ π −
2) (ˆ τ n − ˆ τ κ ) . (137)In order to have stable modes, the numerator in Eq. (132) must also be positive-definite. One can show that theterm inside the square root in Eq. (132) is positive-definite as long as we assume that the product of the coupling termsis negative. Therefore, in order for the modes ˆ ω + and ˆ ω − to have a positive-definite imaginary part it is sufficient toimpose that ˆ τ κ B V − D V + 1 ≥ , ∀ ≤ V ≤ . (138)In order for the inequality to be fulfilled, the smallest root of the polynomial above must be greater than 1, renderingthe function positive in the physical interval in which the background velocity does not exceed the speed of light. Theanalysis here is analogous to the one performed on Eq. (133), leading to a new constraint for the relaxation timesˆ τ π + ˆ τ n ≥ ˆ τ κ + 2 , (139)which is clearly satisfied for by the transport coefficients calculated from the Boltzmann equation [13, 29, 30] andused in the previous section of this paper.Assuming the coupling terms satisfy ˆ L πn ˆ L nπ <
0, the stability conditions derived for the longitudinal modes canbe summarized as ˆ τ π ≥ , ˆ τ n ≥ ˆ τ κ , | ˆ L πn ˆ L nπ | ≤
32 (ˆ τ π −
2) (ˆ τ n − ˆ τ κ ) . (140)From the transverse modes, we had the additional condition | ˆ L πn ˆ L nπ | ≤ τ n (ˆ τ π − . (141)2 k I m ( ω ) V = k I m ( ω ) V = k I m ( ω ) V = - k R e ( ω ) - - - k R e ( ω ) k R e ( ω ) FIG. 14: Real and imaginary parts of the longitudinal modes considering a negative value for the product of the coupling terms,ˆ L πn ˆ L nπ = −
1, for three different values of background velocity, V = 0 . V = 0 .
4, and V = 0 . However, one can demonstrate that this is contained in the inequalities in Eq. (140). Therefore, the conditions forˆ τ π and ˆ τ n are the same as the ones obtained without coupling terms and we just obtain an additional constraint forˆ L πn ˆ L nπ . If we use the relaxation times calculated from the Boltzmann equation, | ˆ L πn ˆ L nπ | ≤ ≈ . . (142)For the sake of completeness, the solution of Eq. (113) for perturbations of top of a moving background are displayedin Figs. 14 and 15 considering a negative and a positive value for the product of the coupling terms, ˆ L πn ˆ L nπ = − L πn ˆ L nπ = 2, respectively. We also show several examples of configurations that are driven unstable by thecoupling terms in Fig. 16, taking the following values for the product of the coupling terms that violate Eq. (142),ˆ L πn ˆ L nπ = − , − , − , − VI. CONCLUSIONS AND REMARKS
In this work we presented a linear stability analysis of Israel-Stewart theory including the effects of shear-stresstensor and net-baryon current. In particular, we investigate the effects that second order terms (that are linear inthe dissipative currents) that couple one dissipative current with the other, the diffusion-viscous couplings, can haveon the stability and causality of the theory. In order to achieve this goal, we considered small perturbations arounda global equilibrium state with energy density ε , vanishing net-baryon number, and a finite fluid 4-velocity u µ . Thestability condition is that the system returns to this global equilibrium state after being perturbed.We first considered the case where the coupling terms are zero. In this case, the modes related to fluctuations ofenergy, momentum and net-baryon number decouple. The dispersion relation for the modes related to fluctuations ofenergy and momentum obtained here are identical to the ones first derived in Ref. [11]. We thus obtained the samestability conditions for the shear relaxation time, as in Ref. [11], given by τ π ≥ ηε + P . Furthermore, since we considered the effects of net-baryon current, we also obtained a stability condition for the3 k I m ( ω ) V = k I m ( ω ) V = k I m ( ω ) V = - k R e ( ω ) - - k R e ( ω ) k R e ( ω ) FIG. 15: Real and imaginary parts of the longitudinal modes considering a positive value for the product of the coupling term,ˆ L πn ˆ L nπ = 2, for three different values of background velocity V = 0 . V = 0 .
4, and V = 0 . ℒ n π ℒ π n = - ℒ n π ℒ π n = - ℒ n π ℒ π n = - ℒ n π ℒ π n = - - - - - - - - k I m ( ω ) FIG. 16: Imaginary part of the unstable shear modes for different negative values of the coupling term, ˆ L πn ˆ L nπ = − , − , − , −
60, for V = 0 . net-baryon diffusion relaxation time, τ n ≥ κ n ¯ n B . We then investigated how the introduction of the aforementioned coupling terms affects these stability conditions.In this scenario, we proved that the theory is not stable for arbitrary values of such coupling terms. In order to be4consistent with kinetic theory calculations and the derivation of fluid-dynamics from second law of thermodynamics,we assumed that the product of the coupling terms is negative. We then showed that the stability conditions for therelaxation times derived in the presence of the coupling term are not modified by the inclusion of the coupling terms,and remain being the ones listed above. Furthermore, we obtained a stability condition that must be satisfied by thecoupling terms themselves, given by | (cid:96) πn (cid:96) nπ | ≤ (cid:18) τ π − ηε + P (cid:19) (cid:18) τ n − κ n ¯ n B (cid:19) . Once more, we emphasize that these conditions are obtained assuming (cid:96) πn (cid:96) nπ <
0. However, we showed some examplethat the system can be stable for (cid:96) πn (cid:96) nπ >
0, but these cases were not studied thoroughly.We confirmed in this paper that the causality conditions obtained for perturbations around a background fluid atrest lead to instabilities for perturbations on top of a moving fluid, as it was first derived in Ref. [11]. Finally, wenote that the stability conditions obtained by Olson in Ref. [26] also include the effects of diffusion-viscous couplingand even bulk viscous pressure, which was neglected here. However, in this case the constraints for the transportcoefficients were written in a more convoluted form (constraints for the diffusion-viscous coupling terms were notexplicitly derived) and the dispersion relations were not explicitly calculated.
VII. ACKNOWLEDGEMENTS
The authors thank J. Noronha for insightful discussions. GSD thanks Conselho Nacional de DesenvolvimentoCient´ıfico e Tecnol´ogico (CNPq) for support and Funda¸c˜ao Carlos Chagas Filho de Amparo `a Pesquisa do Estadodo Rio de Janeiro (FAPERJ), grant number E-26/202.747/2018. CVB thanks Coordena¸c˜ao de Aperfei¸coamento dePessoal de N´ıvel Superior (CAPES) and FAPERJ for support. [1] L. Rezzolla and O. Zanotti,
Relativistic hydrodynamics (Oxford University Press, 2013).[2] K. Yagi, T. Hatsuda, and Y. Miake,
Quark-gluon plasma: From big bang to little bang (Cambridge University Press, 2005).[3] W. Florkowski, M. P. Heller, and M. Spaliski, Reports on Progress in Physics , 046001 (2018).[4] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. , 123 (2013), arXiv:1301.2826 [nucl-th] .[5] R. Derradi de Souza, T. Koide, and T. Kodama, Prog. Part. Nucl. Phys. , 35 (2016), arXiv:1506.03863 [nucl-th] .[6] C. Gale, S. Jeon, and B. Schenke, International Journal of Modern Physics A , 1340011 (2013).[7] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Second Edition: Volume 6 (Course of Theoretical Physics) , Course oftheoretical physics / by L. D. Landau and E. M. Lifshitz, Vol. 6 (Butterworth-Heinemann, 1987).[8] C. Eckart, Phys. Rev. , 919 (1940).[9] W. A. Hiscock and L. Lindblom, Phys. Rev. D , 725 (1985).[10] G. S. Denicol, T. Kodama, T. Koide, and P. Mota, Journal of Physics G: Nuclear and Particle Physics , 115102 (2008).[11] S. Pu, T. Koide, and D. H. Rischke, Phys. Rev. D , 114039 (2010).[12] H. Grad, Communications on Pure and Applied Mathematics , 331 (1949).[13] W. Israel, Annals Phys. , 310 (1976).[14] W. Israel and J. Stewart, Annals Phys. , 341 (1979).[15] I. M¨uller and J. M. Mart´ı, Living Rev. Rel. , 1 (1999).[16] I. M¨uller, Z. Phys. , 329 (1967).[17] I. S. Liu, I. M¨uller, and T. Ruggeri, Ann. Phys. (N.Y.) , 191 (1986).[18] B. Carter, Proc. R. Soc. London, Ser A , 45 (1991).[19] M. Grmela and H. C. ¨Ottinger, Phys. Rev. E , 6620 (1997).[20] R. Baier, P. Romatschke, D. T. Son, A. O. Starinets, and M. A. Stephanov, JHEP , 100100 (2008).[21] G. Denicol, T. Koide, and D. Rischke, Phys. Rev. Lett. , 162501 (2010).[22] E. Molnr, H. Niemi, G. Denicol, and D. Rischke, Phys. Rev. D , 074010 (2014).[23] G. Denicol, H. Niemi, E. Molnar, and D. Rischke, Phys. Rev. D , 114047 (2012), [Erratum: Phys.Rev.D 91, 039902(2015)].[24] F. S. Bemfica, M. M. Disconzi, and J. Noronha, Phys. Rev. D , 104064 (2018).[25] W. A. Hiscock and L. Lindblom, Annals of Physics , 466 (1983).[26] T. S. Olson, Annals of Physics , 18 (1990).[27] F. S. Bemfica, M. M. Disconzi, and J. Noronha, Phys. Rev. Lett. (2019).[28] F. S. Bemfica, M. M. Disconzi, V. Hoang, J. Noronha, and M. Radosz, “Nonlinear Constraints on Relativistic Fluids FarFrom Equilibrium,” (2020), arXiv:2005.11632 [hep-th] .[29] G. S. Denicol, H. Niemi, E. Molnr, and D. H. Rischke, Phys. Rev. D (2012). [30] G. Denicol, E. Molnr, H. Niemi, and D. Rischke, Eur. Phys. J. A , 170 (2012).[31] A. Muronga, Phys. Rev. C (2004).[32] S. Pu, T. Koide, and Q. Wang, AIP Conf. Proc. , 186 (2010).[33] J. D. Jackson, Classical electrodynamics , 3rd ed. (Wiley, New York, NY, 1999).[34] E. J. Routh,