Phenomenological Relativistic Second-Order Hydrodynamics for Multiflavor Fluids

In this work, we perform a phenomenological derivation of the first- and second-order relativistic hydrodynamics of dissipative fluids. To set the stage, we start with a review of the ideal relativistic hydrodynamics from energy-momentum and particle number conservation equations. We then go on to discuss the matching conditions to local thermodynamical equilibrium, symmetries of the energy-momentum tensor, decomposition of dissipative processes according to their Lorentz structure, and finally, the definition of the fluid velocity in the Landau and Eckart frames. With this preparatory work, we first formulate the first-order (Navier-Stokes) relativistic hydrodynamics from the entropy flow equation, keeping only the first-order gradients of thermodynamical forces. A generalized form of diffusion terms is found with a matrix of diffusion coefficients describing the relative diffusion between various flavors. The procedure of finding the dissipative terms is then extended to the second order to obtain the most general form of dissipative function for multiflavor systems up to the second order in dissipative fluxes. The dissipative function now includes in addition to the usual second-order transport coefficients of Israel-Stewart theory also second-order diffusion between different flavors. The relaxation-type equations of second-order hydrodynamics are found from the requirement of positivity of the dissipation function, which features the finite relaxation times of various dissipative processes that guarantee the causality and stability of the fluid dynamics. These equations contain a complete set of nonlinear terms in the thermodynamic gradients and dissipative fluxes arising from the entropy current, which are not present in the conventional Israel-Stewart theory.


Introduction
Relativistic hydrodynamics has been widely applied in recent years in heavy-ion physics [1][2][3][4] and astrophysics [5][6][7]. For example, the experiments at Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) motivated significantly the development of the relativistic hydrodynamics with the data being well-described in terms of a fluid with low shear viscosity; for reviews, see [8,9]. Relativistic hydrodynamics is also used in the general relativistic simulations of compact stars in isolation or binaries, particularly in the context of gravitational wave emission by neutron star mergers observed in 2017 [10]. A further astrophysical area of applications of relativistic hydrodynamics is the physics of compact star rotational dynamics, specifically glitches and their relaxations; for reviews, see [11,12].
The hydrodynamic description of fluids is valid close to the local thermal equilibrium. The hydrodynamic state of a relativistic fluid is described through its energy-momentum tensor and currents of conserved charges, which in the low-frequency and long-wavelength can be Taylor-expanded around their equilibrium values in thermodynamic gradients (socalled thermodynamic forces). The validity of such gradient expansion is guaranteed due to the clear separation between the typical microscopic and macroscopic scales of the system. The zeroth-order term in this expansion corresponds to the limit of the ideal hydrodynamics.
The truncation of the gradient expansion at the first order leads to the relativistic Navier-Stokes (NS) theory, which was worked out by Eckart [13] and Landau-Lifshitz [14]. It is known that the solutions of the relativistic NS equations are acausal and unstable [15][16][17][18]. The reason for the acausality is the parabolic structure of NS equations, which originated from the linear constitutive relations between the dissipative fluxes and the thermodynamic forces. Recent work demonstrated that the acausalities and instabilities in relativistic hydrodynamics are a consequence of the matching procedure to the local equilibrium reference state. More general matching conditions were used to render the theory causal and stable at first order [19][20][21].
The problem of acausality can be solved in the second-order theory, where additional terms appear that contain higher (second)-order derivatives in thermodynamic quantities. For non-relativistic fluids, the second-order theory was proposed by Müller [22], and then rediscovered and extended to relativistic systems by Israel and Stewart [23,24]. In these theories, the dissipative fluxes are treated as independent state variables that satisfy relaxation-type equations derivable from the entropy principle. The relaxation terms which appear in these equations recover the causality of the theory [18,25]. The relaxation equations for dissipative fluxes in the second-order theories and their numerical studies in simulations have been discussed extensively in the literature; see [26][27][28][29][30][31][32][33][34][35][36][37].
The relativistic second-order hydrodynamics can be obtained from moments of the Boltzmann equation for the distribution function [38][39][40]. This theory provides a systematic way of evaluating the new coefficients describing the relaxation of dissipative quantities at weak coupling and in the quasiparticle limit. An alternative approach valid in the strong coupling limit is Zubarev's non-equilibrium statistical operator formalism, which was recently applied to obtain the second-order hydrodynamics and/or Kubo-type formulae for transport coefficients [41][42][43]. Related approaches based on field theory methods are given in Refs. [44,45].
Recent applications of relativistic hydrodynamics focused on the systems featuring spin, polarization, and vorticity. A systematic calculation of the corrections of the stressenergy tensor and currents up to the second order in thermal vorticity were given in Refs. [46,47] and the relevant Kubo formulae were derived. Similarly, the problem of relativistic hydrodynamics under strong spin was studied on the basis of entropy-current analysis in Ref. [48] where the seventeen transport coefficients of highly anisotropic relativistic hydrodynamics were identified. Ref. [49] formulated relativistic spin hydrodynamics of Dirac fermions in a torsionful curved background and derived the relevant Kubo formulae associated with the correlation functions. The spin hydrodynamics was also derived from the Boltzmann equation using the method of moments in Ref. [50] up to the second order. Plasmas with vorticity were studied in second-order dissipative hydrodynamics in the presence of a chiral imbalance in Ref. [51] and the dispersion relations of chiral vortical waves were obtained.
Specific applications of second-order hydrodynamics include, for example, the computation of the transport coefficients from a Chapman-Enskog-like expansion for a system of massless quarks and gluons [34]. The second-order hydrodynamics was recently applied to multicomponent systems with hard-sphere interactions in Ref. [52], where the multicomponent nature of the system was taken into account. The second-order hydrodynamics was applied to quark-gluon plasma in heavy ion collisions, for example, in Ref. [53], showing the coupling between diffusion and shear and bulk viscosities in a multi-component system. A discussion of the use of relativistic anisotropic hydrodynamics to study the physics of ultrarelativistic heavy-ion collisions was given in Ref. [54], using as the main ingredient the quasiparticle anisotropic hydrodynamics model for quark-gluon plasma. Anomalies in the field theory were included in the description of heavy-ion collisions in Ref. [55]. The strong gravity regime of second-order hydrodynamics, relevant for astrophysical applications, was studied in Ref. [56]. The first-order theory was applied in the context of the out-of-equilibrium dynamics of viscous fluids in a spatially flat cosmology [57].
Recent work also extended the current formulations of relativistic hydrodynamics to include the effect of magnetic fields. The first-order dissipative effects and relevant transport coefficients were derived, and a numerical method was developed using the method of moments in Ref. [58]. The effects of anomalies in the magnetohydrodynamics were included in Ref. [59]. Charge diffusion in the second-order theories was recently studied in Ref. [60]. Applications to neutron stars and relevant transport coefficients were derived, for example, in Refs. [61][62][63]. A comprehensive review of the new developments in this field is given in Ref. [64].
More formal recent developments of relativistic hydrodynamics include formulations that avoid specific frames. (The choice of frames in relativistic heavy ion collisions was discussed in Ref. [65].) First-and second-order viscous relativistic hydrodynamics were formulated without specific frame conditions in Refs. [21,66], respectively, with the causality and stability conditions explicitly formulated in the conformal regime in the second work. Earlier, the stability of first-order hydrodynamics was established in Ref. [20].
The purpose of this article is twofold. First, we review the phenomenological derivation of the second-order dissipative hydrodynamics from general principles of conservation laws and the second law of thermodynamics. The dissipative equations are obtained from expansions of the energy-momentum tensor and currents around their equilibrium values in thermodynamic gradients (so-called thermodynamic forces) order by order. Second, we consider the theory in the general case of a system with l independent flavors of conserved charges, keeping all possible second-order terms which arise from the entropy current and the second law of thermodynamics. In the case where the nonlinear terms in the thermodynamic gradients and dissipative fluxes are dropped, these equations are reduced to the original Israel and Stewart hydrodynamics [23,24].
This work is organized as follows. In Section 2, we present the ideal hydrodynamics of non-dissipative fluids. Dissipation in the fluids is introduced in Section 3. We discuss the first-order hydrodynamics and obtain the Navier-Stokes equations in Section 4. The second-order Israel-Stewart theory in the case of a system with l independent flavors of conserved charges is given in Section 5. A summary is provided in Section 6.

Relativistic Ideal Hydrodynamics
Relativistic hydrodynamics describes the state of the fluid employing its energymomentum tensor T µν and currents of conserved charges N µ a , such as the baryonic, electric, etc. Here, we consider the general case of a system with l independent flavors of conserved charges, which are labeled by the index a. The equations of relativistic hydrodynamics are contained in the conservation laws for the energy-momentum tensor and the charge currents In dissipative hydrodynamics, also the entropy principle ∂ µ S µ ≥ 0 should be applied to close the system (1), where S µ is the entropy 4-current.
Ideal hydrodynamics corresponds to the zeroth-order expansion of the energymomentum tensor and charge currents with respect to the thermodynamic forces. In this case, each fluid element maintains the local thermal equilibrium during its evolution [14,67]. The macroscopic state of the fluid is therefore fully described through the fields of the energy density (x) and the charge densities n a (x). The fluid 4-velocity is defined as where x µ is the displacement 4-vector of a fluid element, dτ = g µν dx µ dx ν is the infinitesimal change in the proper time, v(x) is the fluid 3-velocity and γ(x) = (1 − v 2 ) −1/2 is the relevant Lorentz factor. The 4-velocity given by Equation (2) is normalized by the condition u µ (x)u µ (x) = 1 and has only three independent components. Because the thermal equilibrium is maintained locally, each fluid element can be assigned well-defined local values of the temperature T(x), chemical potentials µ a (x) (conjugate to charge densities n a (x)), entropy density s(x) and pressure p(x). These quantities are related to the local energy and charge densities via an equation of state p = p( , n a ) and the standard thermodynamic relations (first law) (3) dp = sdT + ∑ a n a dµ a , (Gibbs − Duhem relation) where h is the enthalpy density. The fluid 4-velocity is defined in such a way that in the fluid rest frame, the energy and charge flows vanish: N i = 0 and T 0i = 0, if u i = 0. These conditions together with the spatial isotropy imply the following form for the energy-momentum tensor and charge currents [14,67]: where the index 0 labels the quantities in ideal hydrodynamics; ∆ µν = g µν − u µ u ν is the projection operator onto the 3-space orthogonal to u µ and has the properties From Equations (6) and (7), we obtain In the fluid rest frame, u µ = (1, 0, 0, 0); therefore, ∆ µν = diag(0, −1, −1, −1) and ∆ ij = ∆ ij = −δ ij . In this case, Equation (8) simplifies to We introduce also the energy 4-current or momentum density by the formula which is the relativistic generalization of the 3-momentum density (mass current). As seen from Equations (6) and (10), all charge currents are parallel to each other and to the energy flow, which is due to the possibility of a unique definition of the velocity field for ideal fluids. The equations of ideal hydrodynamics are obtained by substituting the expressions (6) into the conservation laws (1) Contracting the second equation in (11) once with u ν and then with the projector ∆ αν and taking into account Equation (7) and the condition u ν ∂ µ u ν = 0, we obtain Dn a + n a θ = 0, where we introduce the covariant time derivative, covariant spatial derivative, and velocity 4-divergence via D = u µ ∂ µ , ∇ α = ∆ αβ ∂ β and θ = ∂ µ u µ , respectively. The velocity 4divergence θ quantifies how fast the fluid is expanding (θ > 0) or contracting (θ < 0); it is called the fluid expansion rate. In the case of an incompressible flow of the fluid, we have θ = 0. It is not difficult to recognize in the first two relations in (12) the covariant expressions for the charge conservation law and the energy conservation law, respectively. The third equation is nothing more than the relativistic generalization of the ordinary Euler equation familiar from nonrelativistic hydrodynamics. From these equations, one can deduce that in relativistic hydrodynamics, the role of rest mass density is taken over by the enthalpy density h, which thus provides the correct inertia measure for relativistic fluids. This fact illustrates the importance of the quantity (10) as the relativistic analog of the momentum flux.
The system (12) contains l + 4 equations for l + 5 variables , p, n a and u µ . To close the system, one still needs to specify an equation of state p = p( , n a ), which relates the pressure to the conserved thermodynamic variables.
It is easy to show that the equations of ideal hydrodynamics lead automatically to entropy conservation. The entropy flux can be written as Using the thermodynamic relations (3)-(5) and the equations of motion (12), we obtain which is the second law of thermodynamics for a non-dissipative system. Using Equations (5), (6), (10) and (13), we can rewrite the entropy current as where we defined The expression (15) is the covariant form of the relation (5). To proceed further, it is convenient to modify Equations (3) and (4) using the definitions (16). We obtain for Equation (4) where we used Equation (5) in the second step. Now the first law of thermodynamics and the Gibbs-Duhem relation can be written in an alternative form: With the aid of Equations (6) and (15), these relations can be cast into a covariant form: where we used the second relation in Equation (18). One should note that, despite their vector form, these equations do not contain more information than the scalar thermodynamic relations (contraction of Equations (19) and (20) with the projector ∆ µν leads to identities).

Matching Conditions
Ideal hydrodynamics relies on the strong assumption of the local thermodynamic equilibrium of each fluid element. However, the local equilibrium cannot be maintained permanently since all real systems react to the non-uniformities of the fluid on finite time scales by generating irreversible (dissipative) fluxes. These fluxes tend to eliminate the local gradients to drive the system toward global thermal equilibrium, leading thereby to energy dissipation and entropy increase.
As in the case of ideal hydrodynamics, the equations of relativistic dissipative hydrodynamics can be obtained from the conservation laws given by Equation (1). The expressions (6) then need to be generalized by taking into account all possible effects of dissipation. The energy-momentum tensor obtains an anisotropic contribution because of irreversible momentum exchange between different fluid elements, as well as heat transfer due to the relative motion of energy and charge currents. The diffusion processes between different charge species in their turn introduce dissipative terms in the charge currents N µ a . As a result, the energy-momentum tensor and the charge currents for a dissipative fluid can be written in the following form: where j µ a and τ µν are dissipative terms, and the tensor τ µν is symmetric. Note that since the system is out of thermodynamic equilibrium, the thermodynamic parameters are not well-defined anymore, and one needs to impose additional conditions to specify what should be exactly understood under , n a and p in Equations (21) and (22). The thermodynamic variables in non-equilibrium states can be defined only by means of a fictitious equilibrium state, which should be constructed in such a way as to satisfy the thermodynamic relations (3)-(5) [24]. In order to construct such an equilibrium state for the given values of T µν and N µ a , we define first the energy and charge densities via the so-called matching (fitting) conditions: Equations (23) imply simply that and n a are the time-like eigenvalues of the energymomentum tensor and the charge currents, respectively, measured by a local observer comoving with the fluid element. We remark that the quantities and n a depend in general on the choice of u µ , which is the consequence of the ambiguity of the definition of the velocity field for dissipative fluids (see Section 3.3 for details).
We define in the next step an equilibrium entropy density via an equation of state s = s( , n a ), which is chosen to be the same function of the parameters and n a , as it would be in full thermodynamic equilibrium [15,23,24,68]. The rest of the local thermodynamic quantities can be defined through standard thermodynamic relations (3)-(5), i.e., which implies that all thermodynamic parameters are defined in non-equilibrium states via the same functions of and n a as their equilibrium counterparts. The thermodynamic parameters defined above form a fictitious equilibrium state for which the reversible thermodynamic relations are formally valid. However, it is worthwhile to stress that only the energy and charge densities can be ascribed a definite physical meaning, whereas the quantities s, p, T, and µ a do not retain their usual physical meaning when the system is out of thermal equilibrium. These quantities are mathematically convenient to use since they have the approximate physical meaning of their equilibrium counterparts for small departures from the local equilibrium [23,24].
For instance, the quantity p is not the actual thermodynamic pressure (i.e., the work done by a unit change of the fluid volume), but differs from the latter by an additional non-equilibrium term, which is of the first order in velocity gradients [14,24]. Similarly, the quantity s is not the actual non-equilibrium entropy density since it is defined through reversible thermodynamic relations. In other words, s is not the proper quantity which should increase in non-equilibrium processes according to the second law of thermodynamics. However, the first non-equilibrium correction, in this case, appears only at the second order in gradients, and, therefore, can be ignored in the first-order theory [14,23,24]. Indeed, the correction to s should be negatively defined since the entropy attains its maximum in equilibrium. The only scalar quantity which might contribute to the entropy density at the first order in gradients is the velocity 4-divergence θ; therefore, the non-equilibrium entropy density should be s = s( , n a ) − aθ, with some thermodynamic coefficient a. Since θ can be of both signs, we conclude that a = 0.

Decomposition in Different Dissipative Processes
To identify different dissipative processes, it is convenient to separate scalar, vector, and traceless parts of the tensor τ µν . We first note that the matching conditions (23) together with Equations (21) and (22) impose the following orthogonality conditions on the dissipative terms: The tensor τ µν can be further decomposed into its irreducible components parallel and orthogonal to the fluid 4-velocity u µ . For that purpose, it is useful to introduce a fourth-rank traceless projector orthogonal to u µ via which has the properties The most general tensor decomposition of τ µν consistent with the orthogonality condition (25) can be now written as where we defined new dissipative quantities via and used the properties (7). From Equations (21), (22) and (28), we obtain the most general decompositions of the energy-momentum tensor and the charge currents into their irreducible components: The energy flow defined in Equation (10) can be generalized for dissipative fluids as The dissipative terms j µ a , q µ , π µν and Π in Equations (30) and (31) are called charge diffusion fluxes, energy diffusion flux, shear stress tensor, and bulk viscous pressure, respectively. The shear stress tensor is the traceless spatial part of the energy-momentum tensor and describes the dissipation by anisotropic momentum flow, whereas the bulk viscous pressure is the non-equilibrium part of the pressure and is responsible for dissipation during isotropic expansion or compression. The dissipative currents satisfy the following conditions: the first three of which reflect the fact that the dissipation in the fluid should be spatial. Finally, all quantities on the right-hand sides of Equations (21) and (22) can be obtained by the relevant projections of T µν and N µ a : as follows from Equations (21)-(23), (29) and the properties (7) and (27). In the fluid rest frame, we have From these expressions, it can be seen that all dissipative quantities are, as expected, purely spatial in the rest state of the fluid. Each of the vectors q µ and j µ a thus has 3 independent components, and the shear stress tensor π µν has 5 independent components since its trace vanishes. Thus, the total number of independent quantities in Equations (36) and (37) is 4l + 10, which corresponds to the degrees of freedom of the energy-momentum tensor and the charge currents (remembering that p is determined by an equation of state). As for the velocity of the fluid, u µ should not be treated as an independent variable but should be associated with one of the physical currents. We devote the next subsection to the discussion of the possible definitions of the velocity field.

Definition of Flow Velocity
Another important question in dissipative fluid dynamics is the proper definition of fluid velocity. The choice of the frame in the case of ideal hydrodynamics is simple, as it has the energy and charge current flowing parallel to each other. Then, the fluid rest frame is defined via the requirement that these currents vanish identically. The situation is different in the case of dissipative fluids because one is faced with simultaneous flows of energy and particle currents. There are two simple and also natural ways to define the fluid rest state in dissipative hydrodynamics, which we describe in this subsection in turn.
In the Landau frame (or L-frame) [14] the net energy flow vanishes, the definition of fluid velocity u µ in this frame is then given through the time-like eigenvector of T µν : which in combination with Equations (23), (30), and (33), leads to the following relations: where the subscript L indicates that the quantities are evaluated according to the Landau definition of u µ . With such a choice of the velocity field, the energy diffusion flux is zero, whereas the heat transport is accommodated in the particle diffusion fluxes j µ La . A different form of Equation (38) is obtained by noting that the energy flux (32) in the L-frame is parallel to the flow velocity: Next, consider the velocity of fluid u µ in a generic frame. It can be related to u µ L upon noticing that in a generic fluid rest frame the current (32) is given by P µ = (h, q i ) and, therefore, the boost velocity from an arbitrary rest state (u i = 0) to the Landau rest frame (u i L = 0) is given by v i L = q i /h = O 1 . We follow the convention of Refs. [23,24] to denote the n-th order quantities in deviation from equilibrium by the symbol O n . The transformation of the charge currents into the L-frame is then given by Setting here N i a = j i a , N i a = j i La and N 0 a = n a , we find the charge diffusion fluxes in the Landau rest state: It is worthwhile to note that the current on the left-hand side j i La is evaluated at a transformed coordinate denoted by x . However, the difference j i La can be neglected, as it is of the third order in gradients. Because in the fluid rest frame j 0 La = j 0 a = q 0 = 0, we can express Equation (42) in a covariant form as which is valid already in an arbitrary frame, i.e., not exclusively in the fluid rest frame. The 4-currents are the charge diffusion fluxes with respect to the energy flow, i.e., these are the charge currents in the absence of energy flow. Note that the combinations (44) are invariant under first-order changes in u µ despite the fact that the energy-diffusion flux q µ and the chargediffusion fluxes j µ a depend on the chosen velocity field. It will be demonstrated below (see in Section 4), that the energy dissipation in irreversible processes is associated precisely with these currents.
The remaining thermodynamic variables which appear in Equations (30) and (31) remain unchanged if we neglect the second-order changes in thermodynamic gradients induced by u µ [23,24]. Thus, we can summarize the relation between u and u L as which follows straightforwardly upon comparing Equations (32) and (40). The Eckart frame (E-frame) is defined as the frame in which the velocity field is parallel to one of the conserved currents N µ a . For a fluid with a single conserved charge (for example, the net particle number), one has N µ = nu µ + j µ and, therefore, the 4-velocity is defined as [13] which in combination with Equations (23) and (31) gives i.e., the particle diffusion flux vanishes. The subscript E in Equations (46) and (47) indicates that the quantities are evaluated in the E-frame. The boost velocity from an arbitrary rest frame to the Eckart rest frame is v i E = j i /n. This implies that the velocities u µ and u µ E can be related by We are now in a position to transform the energy flux into the Eckart rest frame; we find It follows then that the energy diffusion flux in the E-frame reads The quantity can be seen as the energy flow with respect to the particle flow; therefore, it can be identified with heat flux. The relation (51) demonstrates that heat conduction and particle diffusion correspond to the same phenomenon, which is observed from different reference frames. This is true only when the higher than the first-order deviations from equilibrium can be neglected.
To order O 1 , the relation between the L-frame and the E-frame can be established from Equations (45) and (48) as Finally, let us note that the extension of the discussion above to the case where multiple conserved charges are present is straightforward. This is achieved by attaching a reference frame to each of them, which formally then leads to the definition which, in turn, enforces the vanishing of the corresponding diffusion flux, i.e., j µ a = 0. In what follows below, we will keep the fluid velocity generic (if not stated otherwise) without specifying any particular reference frame.

Equations of Relativistic Dissipative Hydrodynamics
Equations of relativistic dissipative hydrodynamics are obtained by substituting the decompositions (30) and (31) into the conservation laws (1). Using the same technique as in Section 2 and recalling the properties (33), we obtain where we introduced the velocity stress tensor as σ µν = ∆ αβ µν ∂ α u β . The system of Equations (54)-(56) contains l + 4 equations, as in the case of ideal hydrodynamics, while the number of independent variables is now 4l + 10. The additional unknown variables here are 3l components of the diffusion fluxes, and 5 components of the shear stress tensor and the bulk viscous pressure (recall that the equilibrium pressure is given by the equation of state, and one of the diffusion fluxes can always be eliminated). Thus, in order to solve the system of Equations (54)-(56), we need additional 3l + 6 equations for these dissipative quantities. These relations are described in the phenomenological theory via the second law of thermodynamics (the entropy principle).

Relativistic Navier-Stokes (First-Order) Theory
The entropy conservation law (14) derived in the framework of ideal fluid dynamics is no longer valid for dissipative fluids. In this case, it should be replaced by the second law of thermodynamics, which implies that the entropy production rate of an isolated system must be always non-negative: where equality holds only for reversible processes. By analogy with the decomposition of the charge currents (31), we can decompose S µ into contributions parallel and orthogonal to u µ where s is identified with the equilibrium entropy density (which is sufficient for the first-order accuracy, as explained in Section 3.1), and the vector j µ s satisfies the condition u µ j µ s = 0. The index 1 denotes that (58) is only the first-order approximation to the entropy flux. For small departures from equilibrium, it is natural to assume that j µ s is a linear combination of the energy and charge diffusion fluxes whereβ andα a are functions of thermodynamic variables and should be determined from the condition (57). This formulation of relativistic dissipative fluid dynamics was proposed by Eckart [13] and Landau-Lifshitz [14], and leads to the relativistic version of the NS theory. Inserting Equation (59) into Equation (58) and introducing the dissipation function via R ≡ T∂ µ S µ , we obtain where we employed the relations (3) and (5) and eliminated the terms D and Dn a using Equations (54) and (55). Requiring R ≥ 0, we identify from Equation (60) In the third and the fourth terms, we can replace ∂ µ → ∇ µ due to the orthogonality conditions (33). Then, we have We can further simplify the third term by approximating Du µ = h −1 ∇ µ p from Equation (12) since q µ is already of the order O 1 . Using Equation (18), we obtain Substituting this result into Equation (61) and recalling the definition (44), we obtain finally The second law of thermodynamics implies R ≥ 0, which requires the right-hand side of Equation (63) to be a quadratic form of thermodynamic forces θ, π µν and ∇ µ α a . Assuming linear dependence of the dissipative fluxes on the thermodynamic forces, we obtain the constitutive relations where η and ζ are called the shear and the bulk viscosities, respectively, and χ ab is the matrix of the diffusion coefficients. These relations together with Equations (54)-(56) constitute a closed system of equations, which are known as relativistic NS equations. From Equations (63) and (64), we obtain for the dissipation function The positivity of this expression is guaranteed by the positivity of the viscosity coefficients and the eigenvalues of the matrix χ ab (recall that the diffusion fluxes J µ a are spatial). We see from Equation (65) that the contribution of diffusion processes to the function R depends only on the relative diffusion currents J µ a , but not on the currents j µ a and q µ separately. This result, which was obtained in the framework of the relativistic NS (first-order) theory, is the direct consequence of the Lorentz invariance and indicates that the dissipation in the fluid is independent of the choice of the fluid velocity field, as expected.
The entropy current given by Equations (58) and (59) can now be written as where we used Equations (5), (16), (32) and (44) to obtain the second relation. One can give a simple interpretation of the second expression in Equation (66). Recalling that P µ /h = u µ L is the fluid velocity measured in the L-frame (see Equation (32)), we observe that the first term on the right-hand side of Equation (66) is the entropy current that is convected together with the energy. The second term arises as a result of the relative flow between the energy and the charges and, therefore, should be identified with the irreversible part of the entropy flow. Using Equations (5), (16), (32) and (44), we can rewrite Equation (66) also in the following form: which formally coincides with its counterpart of ideal hydrodynamics (15).
If we have only one sort of conserved charge (l = 1), then instead of the third relation in Equation (64), we have simply where we used the definition of the heat current given by Equation (51) and introduced the coefficient of thermal conductivity via Equation (69) establishes the relation between the thermal conductivity and the diffusion coefficient χ. Thus, there are only three independent transport coefficients in the first-order theory which relate the irreversible fluxes to the corresponding thermodynamic forces. Employing Equation (62), we can write the heat flux (68) in the following way which in the fluid rest frame reads Equation (71) is the relativistic generalization of the well-known Fourier law h i = −κ∂ i T of the non-relativistic hydrodynamics [14].
The expressions (65)- (67) in the case of l = 1 reduce to and The last relation in Equation (73) is the decomposition of the entropy flow into its reversible and irreversible components observed from the E-frame.
Note that in the case where there are no conserved charges, i.e., when l = 0, the heat conduction and/or diffusion phenomena are absent [69].

Israel-Stewart (Second-Order) Theory
The first-order theory described in the previous section turns out to be acausal, and, therefore, cannot be regarded as a consistent theory of relativistic dissipative fluids. The origin of acausality lies in the constitutive relations (64), which imply that the thermodynamic forces generate dissipative fluxes instantaneously [15,17,22,23]. In addition, this theory suffers also from instability, which is a consequence of acausality [15][16][17][18]25]. However, as pointed out in the introduction, acausalities, and instabilities are a consequence of the matching procedure to the local-equilibrium reference state which can be generalized to obtain causal and stable first-order dissipative hydrodynamics [19][20][21].
It turns out that to recover the causality, the entropy current S µ is required to be at least a quadratic function of the dissipative fluxes. This idea of an extension of the entropy current up to the second order was first proposed by Müller [22] for nonrelativistic fluids. For relativistic fluids, a similar second-order theory was developed by Israel and Stewart [23,24]. In this subsection, we review briefly the Israel-Stewart (IS) formulation of causal hydrodynamics, following mainly Ref. [23].
As a starting point, the entropy current given by Equations (66) and (67) is extended up to the second order in dissipative quantities j µ a , q µ , π µν and Π. It is worth stressing that, despite the frame dependence of these quantities, the entropy current S µ along with the energy-momentum tensor T µν and the charge currents N µ a should be regarded as a primary variable and should be therefore frame-independent up to the second order [23,24] (we note that the expressions (66) and (67) are frame-independent only at the first order in deviations from equilibrium).
We write now the entropy current in the following form: where S µ 1 is the first-order contribution given by Equations (66) and (67), and the termsR µ and R µ collect all possible second-order corrections.
The most general form of the vectors R µ andR µ in a generic frame is [23] R µ = βu µ 2 where the new coefficients β Π , β π , β ab J = β ba J , α a Π and α a π are unknown functions of and n a . The vector R µ is frame-independent up to the second order, andR µ collects the second-order contributions to S µ which are not frame-independent to the order O 2 (note that we use different metric signature from Ref. [23]). Because the contribution S µ 1 is frame-independent only to the first order, the termR µ is necessary to provide the frame-independence of the total entropy current (74) [23].
The terms in Equations (75) and (76) which are proportional to u µ are responsible for the second-order corrections to the equilibrium entropy density s. Indeed, the non-equilibrium entropy density is identified with S = S µ u µ , which can be found from Equations (74)-(76) where we used Equations (66). The last inequality in Equation (77) requires the entropy density in non-equilibrium states to be smaller than its equilibrium value s. From here, we conclude that β Π ≥ 0, β π ≥ 0, and the matrix β ab J is positive-semidefinite. The term in Equation (77), which is proportional to q µ q µ represents the shift in the entropy density because of the change of the reference frame and is automatically negative. Those terms in Equations (75) and (76) which are orthogonal to u µ represent the irreversible entropy flux arising from couplings between the diffusion and viscous fluxes.
The phenomenological equations for the dissipative fluxes should be found again from the positivity condition of the dissipative function The first term in Equation (78) was already computed in Equation (61): Now we use Equation (56) to eliminate the acceleration term in Equation (79): where we used the second relation of Equation (18) where we used the notation ∇ µ = ∆ µν ∂ ν . Using Equations (76) and (81), we obtain Here, we introduced three additional coefficientsγ Π ,γ Π andγ π because there is an ambiguity in "sharing" the terms involving Πq µ and π µ ν q ν between the diffusion and viscous fluxes [15]. Note that Ref. [23] assumedγ π =γ Π =γ Π = 0. Now we take the divergence of Equation (75): where we introduced again two additional coefficients γ Π and γ π to "share" the terms containing ΠJ µ a and π µ ν J ν a between the diffusion and viscous fluxes [15]. We kept also all terms that contain gradients of transport coefficients which were neglected in Ref. [23].
Combining Equations (82) and (83), we obtain for the dissipative function (78) where we introduced the short-hand notationṡ and used the symmetries of the dissipative fluxes. The expressions in the square brackets in Equation (84) are called generalized or extended thermodynamic forces. Requiring R ≥ 0, fixing for simplicity the L-frame q µ = 0 and assuming linear relations between these forces and the dissipative fluxes, we obtain the following evolution equations: Equations (87)-(89) in this form for single type of conserved charges were derived in Ref. [15] and were written in the E-frame. In the original papers of Israel and Stewart [23,24,70], the terms which are nonlinear in thermodynamic forces and dissipative fluxes were dropped, although the general case of chemically reacting multicomponent mixtures was discussed in Ref. [23]. The current form generalizes the hydrodynamics equations derived in Refs. [15,23,24,70] to the case of multiple independent flavors of conserved charges, where all second-order terms arising from the entropy current are kept.
With the aid of Equations (87)-(89), the dissipative function (84) obtains the form which formally coincides with Equation (65). Defining relaxation times according to we can write Equations (87)-(89) in the following form: The first terms on the right-hand sides of these equations represent the corresponding NS contributions to the dissipative fluxes; see Equation (64). The first terms on the lefthand sides incorporate the relaxation of the dissipative fluxes to their NS values on finite time scales given by Equation (91). Thus, these relaxation terms imply a delay in the response of the dissipative fluxes to thermodynamic forces and recover the causality of the theory [18,25]. The rest of the terms in Equations (92)-(94) are responsible for spatial inhomogeneities in the dissipative fluxes as well as nonlinear couplings between different dissipative processes. We note that the derivation of the second-order hydrodynamics from the kinetic theory produces additional terms which are not obtained within the phenomenological theory [24]. The derivation of complete IS equations from kinetic theory is discussed in Refs. [30,38,[71][72][73].

Summary
We provided a review of the phenomenological theory of second-order relativistic hydrodynamics for systems with multiple conserved charges with an extension to multiflavor fluids. The hydrodynamic state of the system is described through the energy-momentum tensor and the 4-currents of conserved charges. We reviewed the derivation and content of the equations at zeroth, first, and second order in gradient expansion of the energymomentum tensor and currents, which led us to the ideal Navier-Stokes and Israel-Stewart hydrodynamics, respectively. From the positivity condition of the dissipative function, the most general set of dissipative processes was identified which contains also the relative diffusions between different conserved currents. We kept all the second-order gradient terms arising from the second law of thermodynamics, as well as the nonlinear terms in the thermodynamic gradients and dissipative fluxes, which were omitted in the canonical Israel-Stewart theory, but were kept in Ref. [15] in the case of one conserved flavor.
The hydrodynamics theory exposed in this article is phenomenological in nature. More formal but also complex derivations are available in the literature which utilize different concepts and approaches of statistical mechanics, for example, quasiparticle Boltzmann equation [38][39][40] or non-equilibrium statistical operator [41,43]. Nevertheless, the phenomenological theory can be reasonably expected to remain in the arsenal of theoretical tools that can be deployed in studies of fluid in new settings and/or under various external fields.
A few closing remarks are in order concerning the numerical implementations of the second-order dissipative hydrodynamics. First of all, the choice between the frames (i.e., Eckart vs. Landau) may be significantly influenced by the required computational cost, which in turn depends on such factors as system size, boundary conditions and the level of accuracy required. In general, the Landau frame is known to be computationally more expensive compared to the Eckart frame, but it also provides more insight into the behavior of the fluid. Eventually, the choice between the Landau and Eckart frames will depend on the specific problem at hand and the trade-off between accuracy and computational efficiency. A number of methods are available for the solution of the generalized (secondorder) Navier-Stokes equations for multifluids, for example, finite difference method, finite volume method, spectral methods, and Lagrangian particle methods; see Ref. [9] Section 5 for a pedagogical discussion in the context of relativistic heavy ion collisions. Again, the choice of the method will be dictated by the specifics of the problem at hand, the computational cost, and required accuracy.