Stability analysis for a thermodynamically consistent model of relativistic fluid dynamics

The search for thermodynamic admissibility moreover reveals a fundamental difference between liquids and gases in relativistic fluid dynamics, as the reversible convection mechanism is much simpler for liquids than for gases. In relativistic fluid mechanics, positive entropy production is known to be insufficient for guaranteeing stability. Much more restrictive criteria for thermodynamic admissibility have become available in nonequilibrium thermodynamics. We here perform a linear stability analysis for a model of relativistic hydrodynamics that is based on the GENERIC (general equation for the nonequilibrium reversible-irreversible coupling) framework of nonequilibrium thermodynamics. Assuming a simple form of the entropy function, we find stability for the entire range of physically meaningful model parameters. For relativistic fluid dynamics, full thermodynamic admissibility indeed leads to stability.


I. INTRODUCTION
Relativistic fluid dynamics is important in astrophysics and cosmology as, for example, it allows us to describe the collapse of stars into neutron stars, flows around black holes, jets with relativistic speeds originating from the core of active galactic nuclei, the formation of galaxies or the expansion of the entire universe. The first use of relativistic fluid dynamics in the field of high-energy physics appeared in the 1950s, in the works of Landau [1] and Khalatnikov [2]. In recent years, the topic came to a renewed attention due to its potential to describe the experiments on quark-gluon plasmas, produced in Relativistic Heavy Ion Colliders (RHIC) [3][4][5][6].
The first attempts to generalize the laws of fluid mechanics and thermodynamics to the relativistic context were based on an extended Fourier law for heat conduction [7,8]. These formulations suffered from two fundamental flaws: they are unstable [33,34] and the parabolic nature of the classical differential equation for heat flow leads to instantaneous propagation of heat, hence violating causality [9]. To overcome this issue, some authors added an ad hoc relaxation term to Fourier's law, turning the parabolic equation into a hyperbolic equation with a finite propagation speed [10,11].
A big step forward in formulating relativistic fluid dynamics resulted from the insight that the entropy should depend on additional structural variables (related to the momentum and energy fluxes) and that the heat flux should not be proportional to the entropy flux. In 1967, Müller demonstrated in the context of heat flow that such a more complete description of nonequilibrium states solves the problem of infinite propagation speeds, and he moreover showed that such a description emerges naturally from kinetic gas theory [12]. It was then natural to extend Grad's moment method for the Boltzmann equation [13] to the relativistic context [14][15][16][17]. The derivation of moment equations from the relativistic Boltzmann equation culminated in what is now known as the Israel-Stewart model [18,19]. All work along these lines suffers from the infamous closure problem for the moment hierarchy.
The problem of relativistic fluid dynamics has also been approached with the tools of nonequilibrium thermodynamics.
The basic idea of EIT (extended irreversible thermodynamics) is to consider the dissipative fluxes (e.g., the heat flux vector and the viscous pressure tensor) as independent variables and to formulate convection-relaxation equations for such fluxes [20]. In the limit of slow phenomena, these equations reduce to the classical constitutive laws, but they are also suitable to describe fast phenomena, since they lead to hyperbolic equations with finite speeds of propagation for thermal and viscous perturbations. Such models are compatible with the kinetic theory of Grad's 14-moment method [21]; they are applicable to non-stationary processes, compatible with causality and stable under specific constraints [22]. The GENERIC (general equation for the nonequilibrium reversible-irreversible coupling) framework of nonequilibrium thermodynamics is more versatile in the choice of variables than EIT, but more restrictive in the structure of the equations, in particular, by imposing a Hamiltonian structure on reversible dynamics as the basis of a clear and general reversible-irreversible separation. A GENERIC model of relativistic fluid dynamics has been developed in [23,24], and gravity has been included in [25,26]. Further possible approaches include Carter's theory [27], Van's model [28] and the conformal field theory formulation [29].
In the present paper we study the linear stability of the GENERIC model developed in [23,24]. We first present this model and offer a sketch of the heuristic motivation of these equations for relativistic fluid dynamics offered in the recent textbook [30].

II. MODEL
We restrict ourselves to special relativity. For the metric, we adopt the Minkowski tensor with the convention (η µν ) = diag(−1, 1, 1, 1). For the indexes of four-vectors and tensors, we use Greek letters to denote components from 0 to 3 (i.e. time and space components) and Latin letters to denote components from 1 to 3 (i.e. spatial components only). We define a dimensionless velocity four-vector by means of its where the v i are the components of the local velocity of the fluid v, c is the speed of light and γ is the Lorentz factor, chosen in such a way that that is The relativistic continuity equation reads with ρ the mass density and the subscript f denoting quantities evaluated in the co-moving local reference frame of the fluid. To define such a frame, different choices are possible: in the Eckart framework, the particle flow vanishes in the local rest frame of the fluid [7], while in the Landau framework, the momentum density vanishes in the local rest frame [8].
We decide to operate within the Eckart framework, i.e. we consider as the local reference frame the frame in which the spatial components of the velocity vanish u 0 f = 1, u i f = 0. The tensorη µν = η µν + u µ u ν denotes the spatial projection operator. In the local rest frame such an operator selects the spatial components of vectors and tensors it is applied to.
The energy and momentum balances in the relativistic setting are combined and can be expressed as where T µν is the energy-momentum tensor. Such a tensor has to incorporate energy and momentum densities as well as energy and momentum fluxes. To avoid the presence of spatial derivatives in the formulation of the fluxes, in the spirit of the description of nonrelativistic complex fluids [30], we introduce the relaxing structural variables α µν and ω µ , namely the momentum and energy fluxes respectively. At this point the exact physical meaning of α µν and ω µ is left open, but it will emerge in the following. We require Thus, in the co-moving reference frame, the temporal compo- The four-tensor α µν is dimensionless and symmetric, therefore it has only six independent components. By making use of the upper convected Maxwell model for a structural variable [30], we express the time relaxation of α µν as The left hand side represents the upper convected derivative of α µν and structurally reminds of a Lie algebra. In the case of a space-like vector, this would be equivalent to the Lie derivative of α µν in a rigorous sense [31]. The upper convected derivative appearing at the left hand side of Eq. (7) describes the variation in time of the tensor α µν , when following the trajectory of the fluid particle it refers to (translation) as well as the rotation and deformation of α µν along such a trajectory. In particular, u λ ∂ λ α µν is the so-called substantial derivative of α µν accounting for the time derivative and the convection term, while the remaining two terms on the left-hand side account for rotations and deformations. The parameters λ 0 and λ 2 are two constants, representing the characteristic relaxation times of the trace-partᾱ µν and the trace-free partα µν of α µν , The underlying assumption of Eq. (7) is that the trace and traceless part of α µν relax independently from each other. Similarly, the time relaxation equation for the four-vector ω µ is formulated as where the left-hand side is again an upper convected derivative and λ 1 is a characteristic relaxation time. From now on, we renameη µν ω ν =ω µ . The main reason why we introduced η µν is related to the non-relaxing nature of ω 0 in the localrest frame, since (ω 0 ) f corresponds to the local equilibrium temperature, as we will show below. The introduction ofη µν guarantees that both sides of Eq. (9) vanish, upon contraction with u µ , sinceω The appropriate constitutive equation for the energymomentum tensor T µν , is derived from the entropy balance equation, with the assumption that the entropy density in the local rest frame s f depends on the mass density ρ f and internal energy density f , as well as on the structural variables, i.e. it can be expressed as s f (ρ f , f , α µν ,ω µ ). The spatial projection ω µ is used here instead of ω µ because the temporal part of the latter includes thermodynamic information and would be redundant with ρ f and f . From the entropy balance equation, a constraint on the ω µ arises [30], where T f is the nonequilibrium temperature in the local rest frame, defined through Eq. (12) is formally similar to the classical definition of thermodynamic temperature in the nonrelativistic setting, but entails some additional ambiguity. In the equilibrium setting, the derivative of s f with respect to f is performed keeping ρ f constant. Here, we additionally have to keep the structural variables constant. Therefore the definition of the temperature depends on the particular choice of the structural variables.
The final constitutive equation for the energy-momentum tensor obtained from the entropy balance reads where p f is the pressure in the local rest frame. From a physical point of view, u µ T µν u ν = ρ f c 2 + f can be interpreted as the total internal energy density, with the first term on the right-hand side associated with the mass density and the second term associated with the thermodynamic state of the system. The dependence of the energy-momentum tensor on the form of the entropy function is a characteristic feature of nonequilibrium thermodynamics. The associated entropy equation has the balance form, with the time derivative and flux (i.e. flux four-vector) of entropy at the left-hand side, and the entropy production at the righthand side and it reads The entropy-flux four-vector has a contribution proportional toω µ , showing thatω µ is related to a heat flux. Thus, Eq. (11) and Eq. (14) help to clarify the physical meaning of the fourvector ω ν , whose temporal component in the local rest frame represents the temperature, while the spatial components represent the heat flux. For our stability analysis, we need a particular functional form of the entropy function. We choose the quadratic function where H α and H ω are positive constants, to guarantee the concavity of the entropy with respect to the variables it depends on. A quadratic entropy function is not only consistent with the second-order Israel-Stewart model and the standard formulation of EIT, but should also be sufficient for producing the linearized equations for our stability analysis. For the entropy function (15), the energy-momentum tensor (13) becomes Note that this tensor T µν is manifestly symmetric.
For the linear stability analysis, we choose to take ρ f and f as independent thermodynamic variables. All thermodynamic properties are expressed in terms of three quantities: adiabatic thermal expansivity α s , isothermal speed of sound c T and heat capacity ratio γ. The following expressions hold [32] ρ Summarizing, our model consists of 14 equations: the continuity equation Eq. (4), the energy-momentum tensor conservation equation Eq. (5), the relaxation equations of the structural variables α µν and ω µ , respectively Eq. (7) and Eq. (9). The energy-momentum tensor is expressed as Eq. (13). A thermodynamically admissible entropy balance equation Eq. (14) automatically emerges. To close the model, two state equations are required, relating the independent thermodynamic variables f , ρ f to T f and p f and providing the relations in Eqs. (17), (18) and (19). The unknown variables are 14 in total: ρ f , f , u µ , α µν and ω µ (1 + 1 + 3 + 6 + 3 = 14), given the symmetry of the α µν and the constraints in Eqs. (2), (6), (11).

III. LINEAR STABILITY ANALYSIS
We here formulate a linear stability analysis for the abovementioned model of relativistic fluid dynamics. Following standard procedures [22,28,33,34], we express a generic perturbation of any variable q from the equilibrium rest state q 0 in the form δq = δqe ωt+ikx , with t the time and x one of the space directions in the rest frame.

A. Equilibrium rest-frame
We consider an equilibrium state corresponding to the rest state. The values of the variables in such a state are ρ eq = ρ f,0 , p eq = p f,0 , T eq = T f,0 , with constant values ρ f,0 , T f,0 and p f,0 for a particular equilibrium state. The first order variations of the variables, namely the unknowns of the linearized problem, are δρ, δT, δu x , δu y , δu z , δω x , δω y , δω z , δα xx , δα xy , δα xz , δα yy , δα yz , δα zz . The following relationships hold

B. Linearized equations
To prepare the stability analysis, we write the linearized evolution equations of the model as where the subscript 0 has been dropped in T f,0 , ρ f,0 and f,0 for simplicity of notation. We define the following parameters and, from Eq. (24), the equation of motion where Natural ranges of parameters are: 1 < γ < 2 , X > 0, 0 < Y < 1, Q > 0, λ j > 0. The constraints on Y arise from Eq. (28): the temperature should increase when the system is compressed (∂ i u i < 0), hence Y < 1. From the physical point of view, Q gives an indication on the relativistic character of a fluid. The limit Q → 0 corresponds to a nonrelativistic fluid. As the speed of sound c T approaches the speed of light c, the parameter Q grows. The time scales λ 0 and λ 2 can be chosen independently, without the need to impose any further inequality, as it can be seen by splitting the tensorial relaxation equation (26) for the α jk into its tracefree and trace part For small relaxation times λ j , the quantitiesα jk , α ii and ω j relax quickly to the quasi-stationary solutions; hence, from Eqs. (31), (32) and (25) respectively, one obtains These equations lead to the proper interpretation of the relaxing structural variables α jk and ω j : Eqs. (33) and (34) contain the proper combinations of the velocity gradient tensor for the Newtonian shear and dilatational stresses expected in the nonrelativistic limit; Eq. (35) contains the heat flux (multiplied by a factor). In particular, −ω j consists of the convective and conductive flux contributions in the temperature equation.
C. Block structure The linearized system can be rewritten in the matrix form M · δq = 0, where M is a 14 × 14 matrix with block-structure (6 + 3 + 3 + 1 + 1) and δq = (δρ, δT, δu x , δω x , δα xx , δα yy +δα zz 2 , δu y , δω y , δα xy , δu z , δω z , δα xz , δα yy − δα zz , δα yz ) T . Without loss of generality, we can choose ρ f = T f = c = λ 1 = 1, which is equivalent to choosing the units for density, temperature, velocity and time. The matrix M is 3] and N = [6 × 6]. In particular In the above equations, we defined Φ = ω/k as well as the additional dimensionless parameters To assess the stability, we look for the roots Φ n of the characteristic polynomial detM = 0, with n = 1, ..., 14, i.e. we look for the values of Φ for which the system has nontrivial solutions for any initial perturbation of the variables δq n respect to the equilibrium rest state q 0 . We assume k R + , therefore, if Re(Φ n ) < 0 also ω n < 0, with ω n = kΦ n ; hence any initial perturbations δq n will decay in time and the system will be stable.
Given the block structure of the matrix M the characteristic polynomial can be factorized as We consider the different blocks separately. The solutions of the characteristic polynomial Eq. (41) are found for Eq. (42) implies Φ = −3F , which is always real and negative, thus corresponding to decaying perturbations. In particular, it corresponds to ω = −1/λ 2 , i.e. a constant decaying rate related to the characteristic relaxation constant ofα µν and independent of k, the spatial wavelength of the perturbation.
This condition characterizes the relaxation of shear and second normal stress difference.
In the physically acceptable range of parameters, Eq. (43) has only solutions featuring a negative real part. To proof this, we apply the Routh-Hurwiz criterion [35] for a 3 rd order polynomial a 0 + a 1 Φ + a 2 Φ 2 + Φ 3 . According to this criterion, such a polynomial has only solutions with a negative real part if a i > 0 for ∀i and a 1 a 2 > a 0 . The first condition can be verified from Eq. (43), given the physically acceptable ranges of the parameters, which must all be real and positive. The second condition can be rewritten as which is true for real valued, positive F , G, X, Y , Z. In the limits of Z = 0 or Y = 0, corresponding to the limit of T f = 0 (since γ = 1 and α s = 0 would provide unphysical values of Y , for nonzero values of H ω ), one solution is the purely real negative Φ = −3F , corresponding to an additional shear mode ω = −1/λ 2 . The same would happen in the limit case Q = 0 (nonrelativistic limit) or in the case that F = G, i.e. λ 1 = λ 2 . In particular in the nonrelativistic limit case of Q = 0, additional roots of Eq. (43) will be Φ = −3G(1+Y Z) and Φ = 0. The latter corresponds to a case of marginal stability, while the former corresponds to a perturbation decaying in time with ω = −(1 + T 3 f H ω /c T )/λ 1 , independent of k. This can be interpreted as a non-hydrodynamic mode, where the decaying occurs faster at higher temperatures. Physically, one can imagine a higher thermal agitation of the fluid particles to correspond to a higher dissipation rate.
We now consider Eq. (44). To determine the sign of the roots, we rearrange Eq. (44) as which can be compactly rewritten as X = −N/D with N the numerator and D the denominator of the fraction. In the case of Φ purely real, positive values of Φ would lead to both N > 0 and D > 0, hence X < 0. Since these are physically unacceptable values of X, purely real positive roots are not possible and purely real roots Φ can only be negative.
Purely imaginary values of Φ correspond to an oscillatory non-growing solution of the linear system and do not concern its stability. Therefore further investigation is required only for the case of a complex-valued root, with R = Re(Φ) = 0 and I = Im(Φ) = 0. We substitute Eq. (47) inside Eq. (46) and we impose Im(X) = 0, since we know from Eq. (27) that X has to be a real number. The condition for such a requirement to be met is expressed in terms of Q = Q(R, I, E, F, G, Y, Z) and inserted into Eq. (46), thus finding with We compactly rewrite Eq. (48) as X = −N /D , where N and D are the numerator and the denominator of the fraction, respectively. From the equations above, one can see that T 1 > 0 and T 2 > 0, in the case of R > 0. Thus N > 0, D > 0 and X < 0, which is against the assumptions on the acceptable values of X. Therefore, it has to be R = Re(Φ) ≤ 0 and the proof of the unconditional stability of the system is complete, within the physically acceptable range of values of the parameters. In Figs In the case of Q = 0.001, close to the nonrelativistic limit, we find several modes that have been previously identified by other authors [22]: the line labeled 1 (green online) is characterized by a proportionality to −k 2 near k = 0 and identifies a diffusive mode for k < k c , where k c is a characteristic wavelength. The line 2 (brown online), represents a relaxation mode related to the heat fluxes ω j : this can be seen from the fact that, for k = 0, it takes the value ω = −1/λ 1 (see Eq. (25)). The real parts of 1 and 2 merge for k > k c , while the complex parts becomes nonzero and conjugate, turning both 1 and 2 into two non-hydrodynamic modes. The modes 3 and 4 (orange online) are complex conjugate with negative real part, hence they represent two sound modes. The modes 5 (magenta online) and 6 (yellow online) have negative real part, reaching an asymptote as k increases and null complex part, thus they can be identified as two non-propagating relaxation modes. In particular, at k = 0, the mode 5 takes the value ω = −1/λ 0 while the mode 6 takes the value ω = −1/λ 2 . Hence we conclude that the mode 5 is related to the relaxation of α ii (see Eq. (32)) and the mode 6 is related to the relaxation ofα jk (see Eq. (31)). From the 3×3 block R, three modes appear, represented by the dashed lines: two shear modes, which turn into non-hydrodynamic modes for k larger than a characteristic value k c (dark-grey and light-grey, blue and red online) and one non-hydrodynamic mode (black). For Q = 1, the lines 1 and 2 represent two shear modes for k < k c , since they are shifted downwards respect to the origin of the axes, while the other modes retain a similar meaning. In the case of Q = 0.1, an intermediate situation between the two previous cases appears and the different modes are combined in a nontrivial way.

IV. SUMMARY AND DISCUSSION
In the present paper we considered the thermodynamically consistent model for relativistic fluid dynamics developed within the GENERIC framework [23,24], where the treatment of nonequilibrium phenomena (momentum and heat transport) is described by means of additional structural variables to prevent the problem of the infinite propagation speeds. We performed a linear stability analysis of the equilibrium rest state and we analytically proved that, in the entire range of the physically meaningful parameters, the model is unconditionally stable. We verified that, in the nonrelativistic limit, the standard shear, sound, diffusive and nonpropagating modes are retrieved while, as relativistic effects become more relevant, the modes start to mix with each other in a nontrivial way.
In nonrelativistic fluid dynamics, the equations for liquids and gases possess the same Navier-Stokes-Fourier form and differ only by the characteristic values of the transport coefficients. This is a consequence of the fact that the basic equations express the balance laws for the conserved quantities mass, momentum and energy in a straightforward form. In relativistic fluid dynamics, we are faced with additional relaxation equations for non-conserved structural variables so that there is no reason to expect universal equations that encompass both liquids and gases.
The GENERIC equations of relativistic fluid dynamics include a Poisson-bracket formulation of the convection mechanism. A similar Poisson bracket has not been found for the 13-or 14-moment equations derived from the nonrelativistic or relativistic Boltzmann equation. The problem is that the convection terms for increasingly higher moments in momentum space are strongly coupled and hence very complicated so that they become difficult, if not impossible, to truncate. The situation is very different if the structural variables describe features related to forces and differences in position space. This is the reason why we expect a serious difference between fluid dynamics of relativistic liquids and gases. The heuristic derivation of the GENERIC equations based on the knowledge of complex fluids in [30] underlines their liquidlike background, whereas the Israel-Stewart model and the EIT equations are based on moment expansions for the Boltzmann equation and hence are appropriate for gas-like fluids. Note that, in the case of quark-gluon plasmas, recent experimental findings showed that a liquid-like description may be more accurate [36]. For gas-like fluids, GENERIC structure has so far been found only in the Boltzmann equation or the equivalent infinite moment hierarchy.
With the present paper we performed one step forward in showing that the thermodynamically consistent model developed within the GENERIC framework [23,24] should be considered as a valuable option, a worthwhile alternative to the existing ones, when addressing the study of relativistic fluids. To gain full confidence, linear stability should also be investigated in a boosted system [22], which is considerably more complicated and hence beyond the scope of the present analytical study (more terms in evolution equations, angle between boost and wave vector of perturbation). The elegance associated with the Poissonian structure of the convection terms in the relaxation equations (7) and (9) for the structural variables is an important argument in favor of the equations derived within the GENERIC framework. More physical entropy functions, presumably containing logarithmic rather than quadratic contributions, open the door to more realistic modeling, where the form (13) of the energy-momentum tensor is prescribed by thermodynamic admissibility.