# Stability of multi-component relativistic viscous hydrodynamics from Israel-
Stewart and reproducing DNMR from maximizing the entropy

Dekrayat Almaalol  Illinois Center for Advanced Studies of the Universe,
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL
61801, USA  Travis Dore  Illinois Center for Advanced Studies of the Universe,
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL
61801, USA  Jacquelyn Noronha-Hostler  Illinois Center for Advanced Studies of
the Universe, Department of Physics, University of Illinois at Urbana-
Champaign, Urbana, IL 61801, USA

###### Abstract

The Quark Gluon Plasma produced in heavy-ion collisions has three relevant
conserved charges: baryon number (B), strangeness (S), and electric charge
(Q). Here we derive the Israel-Stewart framework for BSQ diffusion coupled to
shear and bulk viscosity and the thermodynamic derivatives needed to couple
this to a BSQ equation of state. We reproduce a subset of the ${\mathcal{J}}$
terms in the DNMR approach and propose a technique to derive other terms from
the maximum entropy principle. Finally, we preform a stability analysis that
can be used to constraint transport coefficients in BSQ hydrodynamics.

##  I Introduction

Within the standard model, baryon number (B) and electric charge (Q) are
always conserved. Additionally, the time scales of relativistic heavy-ion
collisions are significantly shorter than the strange quark weak decay such
that strangeness (S) is also conserved. Thus, dynamical simulations of heavy-
ion collisions must conserve BSQ, which is especially important at finite
baryon densities $\rho_{B}$. Already, initial studies using only one conserved
charged Denicol  _et al._ (2018); Du and Heinz (2020), ideal BSQ dynamics
Schäfer  _et al._ (2021), or lower dimensions Fotakis  _et al._ (2020); Du
_et al._ (2021) have begun to make model-to-data comparisons. From these
initial studies it was shown that the inclusion of strangeness neutrality can
play a strong role Monnai  _et al._ (2021), fluctuations of BSQ conserved
charges may affect strangeness flow Martinez  _et al._ (2019a); Carzon  _et
al._ (2022), and that BSQ diffusion affects the rapidity dependence of charge
distribution Fotakis  _et al._ (2020). Because baryon diffusion has very non-
trivial behavior at the critical point (aka baryon opalescence) wherein it
drops to 0 Rougemont  _et al._ (2015); Monnai  _et al._ (2017); Rougemont  _et
al._ (2017); Grefa  _et al._ (2022), it is that much more urgent that we
understand the interplay between BSQ diffusion and how that couples to other
transport coefficients (such as the bulk viscosity that diverges at the
critical point Monnai  _et al._ (2017); Martinez  _et al._ (2019b); Rajagopal
_et al._ (2020); Dore  _et al._ (2020, 2022a)).

At this point different approaches exist for deriving the relativistic viscous
hydrodynamic equations of motion i.) Phenomenological Israel-Stewart approach
wherein the equations of motion are derived directly from enforcing that the
entropy current must be positive Israel and Stewart (1979) ii.) the DNMR
approach Denicol  _et al._ (2012a) which derives the equations of motion from
the Boltzmann equation using the method of moments by applying a systematic
power-counting scheme in Knudsen and inverse Reynolds number iii.) the BDNK
equations of motion wherein a generic first-order expansion of the energy-
momentum tensor and charge current is employed, leading to stable and causal
equation of motion in non-traditional hydrodynamic frames Bemfica  _et al._
(2018, 2022); Kovtun (2019); Hoult and Kovtun (2020), iv.) anisotropic
hydrodynamics Martinez and Strickland (2010); Martinez  _et al._ (2012);
Alqahtani  _et al._ (2017) which implements a re-summation scheme using the
moment method of kinetic theory. It has been shown in a number of different
approaches for toy model scenarios that DNMR approaches an attractor more
quickly than Israel-Stewart Heinz  _et al._ (2016); Strickland  _et al._
(2018); Behtash  _et al._ (2018). Since BDNK is such a new theory, studies of
its attracting like behavior are not as developed, although there has been
some work Bemfica  _et al._ (2018). In this work, we will only consider
second-order theories.

Because heavy-ion collisions probe the QCD phase diagram dynamically, the
search for the QCD critical point depends not only on universality in the EoS
but also the dynamical universality class. Thus, the effects of critical
scaling of transport coefficients may probe far-from-equilibrium dynamics
close to the critical point. In fact, it is known that the hydrodynamic
picture breaks down exactly at the critical point Stephanov and Yin (2018);
Nahrgang  _et al._ (2019); Rajagopal  _et al._ (2020); Dore  _et al._ (2020),
although uncertainties remain in terms of how close one can approach the
critical point and still remain in a hydrodynamic regime. While a number of
advances have been made on critical scaling within hydrodynamics Stephanov and
Yin (2018); An  _et al._ (2020); Pradeep  _et al._ (2022), there is not yet
one standardized approach that is agreed on within the community that is also
known to preserve the properties of causality and stability in the equations
of motion.

Thus, our approach is to use hydrodynamic equations of motion that are the
most stable even far-from-equilibrium. Previously, it was confirmed Dore  _et
al._ (2020) that DNMR equations of motion appear to be the most stable (by
studying attractors) when it comes to critical scaling of transport
coefficients when compared to phenomenological Israel-Stewart. For the case of
BSQ charge conservation, new relativistic viscous fluid dynamics equations
must be derived taking into account the BSQ diffusion matrix Greif  _et al._
(2018). Initial studies have done this using a phenomenological Israel-Stewart
approach in Monnai (2012), and based on DNMR in Fotakis  _et al._ (2022)
(although not with all possible DNMR terms). However, there have not yet been
studies that can systematically connect the Israel-Stewart approach and DNMR,
nor have there been causality and stability analyses (a la Bemfica  _et al._
(2021); Plumberg  _et al._ (2021); Chiu and Shen (2021)) applied to BSQ
diffusion yet. Here we use a new approach that can reconcile nearly all terms
in the equation of motion between phenomenological Israel-Stewart and DNMR.
Then, we preform a linear stability analysis to find constraints between both
the transport coefficients and the EoS. We discuss the implications of these
constraints in specific limits. For instance, while the exploration of meta-
stable states requires other physics beyond what is discussed here, we can at
least determine if certain features of a meta-stable region like $c_{s}^{2}<0$
or a negative enthalpy are consistent with our stability analysis. Indeed, we
find that nothing prevents such features within our stability constraints.

This paper is organized as follows. In Sec. II we derive the Israel-Stewart-
based BSQ relativistic viscous hydrodynamic equations of motion and connect
them to DNMR. In Sec. III we derive our stability constraints for BSQ
hydrodynamics using two different approaches. Finally, in Sec. IV we provide
our conclusions and outlook.

##  II multi-component transient fluid dynamics with BSQ charges

Hydrodynamics is an effective field theory based on the local conservation
laws of energy-momentum tensor $T^{\mu\nu}(x)$ and charge currents
$N^{\mu}_{q}(x)$,

$\displaystyle D_{\mu}T^{\mu\nu}=0\,;\qquad\qquad
D_{\mu}N^{\mu}_{q}=0\,,\qquad q\in\\{B,S,Q\\}\,,$ (1)

where $D_{\mu}$ is the covariant derivative which gives for a general
coordinate system

$\displaystyle
D_{\mu}T^{\mu\nu}=\frac{1}{\sqrt{g}}\partial_{\mu}(\sqrt{g}T^{\mu\nu})+\Gamma^{\nu}_{\mu\alpha}T^{\mu\alpha}=0\,,$
(2) $\displaystyle
D_{\mu}N_{q}^{\mu}=\frac{1}{\sqrt{g}}\partial_{\mu}(\sqrt{g}T^{\mu\nu})=0\,,$
(3)

where the metric tensor $g^{\mu\nu}$ is defined with negative signature
$(+,-,-,-$). By extending to multiple conserved charges, the set of equations
of motion which evolve the system increases to 4+$N_{q}$ equations, which
describe the space-time evolution of $10+4\,N_{q}$ independent components of
the energy momentum tensor and the charges current(s), respectively.

The energy-momentum tensor and the charge currents consist of an equilibrium
contribution and a dissipative part,

$\displaystyle T^{\mu\nu}$
$\displaystyle=T^{\mu\nu}_{0}+\Pi^{\mu\nu}\,,\qquad\qquad
N^{\mu}_{q}=N^{\mu}_{q,0}+n^{\mu}_{q}\,;\qquad q\in\\{B,S,Q\\}$ (4)

where the “ 0 ” subscript refer to the equilibrium contribution. The non-
equilibrium contributions of the energy moment tensor is
$\Pi^{\mu\nu}=\Pi\Delta^{\mu\nu}+\pi^{\mu\nu}$ where
$\Pi=\frac{1}{3}\operatorname{Tr}\left[\Pi^{\mu\nu}\right]$ is the bulk
pressure and the symmetric, traceless component of $\Pi^{\mu\nu}$ is
$\pi^{\mu\nu}$, known as the shear stress tensor. The non-equilibrium
contribution of the charge current for a given conserved charge $q$ is
$n^{\mu}_{q}$. The decomposition of the energy momentum tensor and charge
current in the Landau frame is

$\displaystyle T^{\mu\nu}$ $\displaystyle=\varepsilon
u^{\mu}u^{\nu}-(P+\Pi)\Delta^{\mu\nu}+\pi^{\mu\nu}\,;$ (5) $\displaystyle
N^{\mu}_{q}$ $\displaystyle=\rho_{q}u^{\mu}+n^{\mu}_{q}\;,\qquad\qquad
q\in\\{B,S,Q\\}$ (6)

provided the baryon (B), strange (S), and electric charge (Q) are conserved
simultaneously. The time-like future directed flow vector $u^{\mu}$ is
normalized such that $u^{\mu}u_{\mu}=1$, and $u^{\mu}u^{\nu}$ is defined as
the temporal projection operator. Similarly,
$\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ is a spatial projection operator.
The thermodynamic variables are then energy density, $\varepsilon$, the
pressure, $P$, and the individual charge densities $\rho_{q}$.

The 4-flow vector $u^{\mu}$ is defined as the eigenvector for the energy-
momentum tensor $u_{\mu}T^{\mu\nu}=\varepsilon\,T^{\mu\nu}$ with the energy
density $\varepsilon$ as its eigenvalue. The energy and charge densities are
hence defined in the Landau frame by

$\varepsilon\
u^{\mu}=u_{\nu}T^{\mu\nu}\;;\qquad\rho_{q}=u_{\mu}N^{\mu}_{q}\;,$ (7)

where, in the local rest frame of the fluid, $u^{\mu}_{LRF}=(1,0,0,0)$. With
this particular choice of “ hydrodynamic frame ”, we have defined
$\varepsilon$ as the equilibrium energy density , and $\rho_{q}$ as the
equilibrium charge density. Eq. (7) is sometimes referred to as “ Landau
matching ” since taken together to hold generally, they imply the use of the
Landau frame. Where we then have $u_{\mu}\pi^{\mu\nu}=0;\qquad
u_{\mu}n^{\mu}_{q}=0\,\,\forall\,\,q\in\\{B,S,Q\\}$ by definition.

The best constraints on relativistic viscous hydrodynamics come from heavy-ion
collisions at the LHC/top RHIC energies that corresponds to a vanishing baryon
density. Thus, Landau matching has been a natural choice in heavy-ion
collisions because it eliminates extra transport coefficients and because the
Eckart frame is not well defined for zero net charge. However, in the case of
diffusion, there will be now transport coefficients but one can still use the
previously determined constraints in the limit of vanishing baryon densities
as a starting place.

Applying the temporal and spatial projections $u^{\mu}u^{\nu}$ and
$\Delta^{\mu\nu}$ to the conservation laws in Eq. (1) and using the
decomposition of $T^{\mu\nu}$ and $N^{\mu}$ Eqs. (5) and (6) gives

$\displaystyle D\varepsilon$
$\displaystyle=-(\varepsilon+P)\theta-\Pi\theta+\pi_{\mu\nu}\sigma^{\mu\nu}\,;\qquad\qquad\qquad\qquad\qquad\quad\rm``energy\,\,conservation"$
(8) $\displaystyle(\varepsilon+P)Du^{\mu}+\Pi Du^{\mu}$
$\displaystyle=\nabla^{\mu}(P{+}\Pi)-\Delta^{\mu\nu}\nabla^{\lambda}\pi_{\nu\lambda}+\pi^{\mu\nu}Du_{\nu}\,;\quad\qquad\qquad\rm``momentum\,\,conservation"$
(9) $\displaystyle D\rho_{q}$
$\displaystyle=-\rho_{q}\theta-\nabla_{\mu}n^{\mu}_{q}\;.\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\quad\rm``charge\,\,conservation"$
(10)

where $D\equiv u_{\mu}D^{\mu}$ is the covariant time derivative, $\theta\equiv
D_{\mu}u^{\mu}$ is the expansion scalar, and
$\nabla^{\mu}\equiv\Delta_{\mu\nu}\partial^{\nu}$ is the spatial gradient in
the local rest frame of the fluid. From Eq. (9) it becomes clearer why we
refer to $\Pi$ as the bulk pressure because it enters as $P+\Pi$ such that one
obtains an effective pressure.

For an ideal BSQ-charged fluid Eqs. (8 \- 10) give a set of 7 first-order PDEs
for the equilibrium fields $(\varepsilon,P,u^{\mu},\rho_{q})$, which represent
8 independent variables; therefore, one needs an additional input to the
equations. The local equilibrium assumption of hydrodynamics means that the
thermodynamic equation of state, $P=P(\varepsilon,\rho_{q})$, should be used
to reduce the space of variables. With that, the system is completely
determined provided the initial values of
$(\varepsilon_{0},P_{0},u^{\mu}_{0},\rho_{0}^{q})$ are known.

###  II.1 Phenomenological approach

When including viscous effects, one also needs to be able to describe the
dynamical evolution of the bulk, shear stress tensor, and the diffusive BSQ
charge-currents $(\Pi$, $\pi_{\mu\nu}$, $n^{\mu}_{q})$ extending the dynamical
space of variables to $(8+1+5+3N_{q})$ which requires an additional $6+3N_{q}$
equations for the dissipative currents to evolve the 23 - 1 independent
variables where the EoS closes the equations through the relation
$P=P(\varepsilon,\rho_{q})$.

Now that the equations of motions from conservation laws are established in
Eqs. (8 \- 10), we need the time evolution for the out of equilibrium
contributions: $\pi^{\mu\nu}$, $\Pi$, and $n_{q}^{\mu}$. We derive the
dynamical equations for these currents using the phenomenological approach of
extended thermodynamics originally formulated by Israel-Stewart Israel and
Stewart (1979). In this approach, The second law of thermodynamics

$\displaystyle
D_{\mu}S^{\mu}_{0}=\beta_{\nu}D_{\mu}T^{\mu\nu}_{0}-\sum_{q}^{B,S,Q}\alpha_{q}D_{\mu}N^{\mu}_{q,0}\geq
0,$ (11)

is extended out of equilibrium using the local conservation laws (1) in the
form (4) while constraining the entropy production by the second law of
thermodynamics which then gives

$\displaystyle D_{\mu}S^{\mu}_{0}=-\beta_{\nu}D_{\mu}T^{\mu\nu}_{non-
eq}+\sum_{q}^{B,S,Q}\alpha_{q}D_{\mu}N^{\mu}_{q,non-eq}\geq 0,$ (12)

where $S^{\mu}_{non-eq}$ is the entropy current out of equilibrium. The
thermodynamic quantities are defined by $\beta_{\nu}=\frac{u^{\nu}}{T}$, with
T being the equilibrium temperature. The thermodynamic potential associated
with the conserved charge $q$ is $\alpha_{q}=\frac{\mu_{q}}{T}$ where
$\mu_{q}$ is the chemical potential. Finally, $T^{\mu\nu}_{non-eq}$ and
$N^{\mu}_{q,non-eq}$ are the non-equilibrium contributions to the full energy
momentum tensor and charge currents from Eq. (4).

The starting point is the proposition that the second-order entropy current
for a multi-component fluid can be written as the following

$\displaystyle{S}^{\mu}_{non-eq}$
$\displaystyle=su^{\mu}-\sum_{q}^{B,S,Q}\alpha_{q}{n}^{\mu}_{q}-\frac{1}{2}u^{\mu}\left(\beta_{\Pi}{\Pi}^{2}+\beta_{\pi}{\pi}^{\mu\nu}{\pi}_{\mu\nu}+\sum_{q}^{B,S,Q}\beta_{n}^{qq^{\prime}}n_{q}^{\mu}{n}_{\mu}^{q^{\prime}}\right)+\sum_{q}^{B,S,Q}\left(\delta_{n\Pi}^{q}n^{\mu}_{q}{\Pi}+\delta_{n\pi}^{q}n^{\nu}_{q}{\pi}^{\mu}_{\nu}\right)$
$\displaystyle+u_{\nu}(\delta_{\Pi\pi}\Pi\pi^{\mu\nu})+higher\,\,order\,\,terms\,\,,$
(13)

where the equilibrium entropy density, s, is given by Gibbs–Duhem equation

$\varepsilon+p=sT+\sum_{q}^{B,S,Q}n_{q}\mu_{q}\,,$ (14)

and the equilibrium entropy current is given by $s^{\mu}_{0}=su^{\mu}$. The
second order transport coefficients are $\\{\beta_{\Pi},\beta_{\pi}\\}$ for
bulk and shear stress tensor, and $\\{\delta_{\Pi\pi}\\}$ for coupling between
($\Pi,\pi_{\mu\nu}$). The matrix $\beta_{n}^{qq^{\prime}}$ is a symmetric
matrix for diffusion of charge $q\leftrightarrow q^{\prime}$, and
$\\{\delta_{n\Pi}^{q},\delta_{n\pi}^{q}\\}$ are the second order coupling
coefficients between any fluid component q and bulk pressure and shear stress-
tensor, respectively.

Extending the standard IS treatment to multiple charges by including all three
charge currents $\\{n_{B},n_{S},n_{Q}\\}$ adds six additional independent
contributions from the $\beta_{qq^{\prime}}$ term which account for the
diffusive coupling among these charges. It will also extend the terms
$\\{\delta_{n\Pi}^{q},\delta_{n\pi}^{q}\\}$ similarly.

The entropy production rate at a space-time point $x$ is obtained from the
four divergence of the entropy current vector in Eq. (13), giving

$\displaystyle\partial_{\mu}{S}^{\mu}$
$\displaystyle=\left(-\beta\theta\Pi-\beta_{\Pi}\dot{{\Pi}}\Pi-\frac{\dot{\beta}_{\Pi}}{2}\Pi^{2}-\frac{\beta_{\Pi}}{2}{\Pi^{2}}\theta\right)+\left(\beta\pi_{\mu\nu}\sigma^{\mu\nu}-\beta_{\pi}\dot{{\pi}}^{\mu\nu}{\pi}_{\mu\nu}-\frac{\dot{\beta}_{\pi}}{2}{\pi}^{\mu\nu}{\pi}_{\mu\nu}+\frac{\beta_{\pi}}{2}{\pi}^{\mu\nu}{\pi}_{\mu\nu}\theta\right)$
$\displaystyle-\sum_{q}\left(n^{\mu}_{q}\nabla_{\mu}\alpha_{q}+\beta_{qq^{\prime}}\dot{{n}}^{\mu}_{q}n_{\mu}^{q^{\prime}}+\frac{\dot{\beta}_{qq^{\prime}}}{2}{n}^{\mu}_{q}n_{\mu}^{q^{\prime}}+\frac{\beta_{qq^{\prime}}}{2}{n}^{\mu}_{q}{n}_{\mu}^{q^{\prime}}\theta\right)-\bigg{(}(\partial_{\mu}u_{\nu})\delta_{\Pi\pi}\Pi\pi_{\mu\nu}+\delta_{\Pi\pi}\Pi(u_{\nu}\partial_{\mu}{\pi}^{\mu\nu})\bigg{)}$
$\displaystyle-\sum_{q}\bigg{(}\delta_{n\Pi}^{q}n_{q}^{\mu}\nabla_{\mu}{\Pi}+\delta_{n\Pi}^{q}\Pi\partial_{\mu}{n}^{\mu}_{q}+\lambda_{\Pi}^{q}{n}_{q}^{\mu}\Pi\nabla_{\mu}\delta_{n\Pi}^{q}-\delta_{n\pi}^{q}n^{\mu}_{q}\nabla^{\nu}{\pi}_{\mu\nu}+\delta^{q}_{n\pi}\pi_{\mu\nu}\nabla^{\langle\mu}n_{q}^{\nu\rangle}+\lambda_{\pi}^{q}\pi_{\mu\nu}n_{q}^{\langle\mu}\nabla^{\nu\rangle}\delta^{q}_{n\pi}\bigg{)}\,,$
(15)

where the subscript “ $\dot{}\,$ ” denotes the covariant time derivative
$u^{\mu}D_{\mu}$ of a particular term. The angular brackets across indices is
to emphasize that after contracting with $\pi^{\mu\nu}$, only the trace-less
symmetric part survives.

The standard treatment of Israel-Stewart is to define the tensors
$\\{\Pi,\pi_{\mu\nu},n^{\mu}_{q}\\}$ in such a way which introduces an
effective description of these dissipative currents. This definition
reproduces the generic structural form of second law of thermodynamics in
order to guarantee the positivity of the entropy production. We follow the
same procedure here and start by regrouping the dissipative terms in the
entropy production to get

$\displaystyle\partial_{\mu}{S}^{\mu}$
$\displaystyle=\beta{\Pi}\left\\{-\theta-\frac{\beta_{\Pi}}{\beta}\dot{{\Pi}}-\frac{\dot{\beta}_{\Pi}}{2\beta}{\Pi}-\frac{\beta_{\Pi}}{2\beta}{\Pi}\theta-\frac{1}{\beta}\sum_{q}\left(\delta_{n\Pi}^{q}\partial_{\mu}{n}^{\mu}_{q}+\lambda_{\Pi}^{q}{n}_{q}^{\mu}\nabla_{\mu}\delta_{n\Pi}^{q}\right)-\frac{\delta_{\Pi\pi}}{2\beta}(\partial_{\mu}u_{\nu})\pi^{\mu\nu}\right\\}$
$\displaystyle+\beta\pi_{\mu\nu}\left\\{\sigma^{\mu\nu}-\frac{\beta_{\pi}}{\beta}\dot{{\pi}}^{\mu\nu}-\frac{\dot{\beta}_{\pi}}{2\beta}{\pi}^{\mu\nu}+\frac{\beta_{\pi}}{2\beta}{\pi}^{\mu\nu}\theta-\frac{1}{\beta}\sum_{q}\left(\delta^{q}_{n\pi}\nabla^{\langle\mu}n_{q}^{\nu\rangle}+\lambda_{\pi}^{q}n_{q}^{\langle\mu}\nabla^{\nu\rangle}\delta^{q}_{1})\right)-\frac{\delta_{\Pi\pi}}{2\beta}(\partial_{\mu}u_{\nu})\Pi\right\\}$
$\displaystyle+\sum_{q}{n}^{\mu}_{q}\sum_{q^{\prime}}\left\\{-\nabla_{\mu}\alpha_{q^{\prime}}-\beta_{qq^{\prime}}\dot{{n}}_{\mu}^{q^{\prime}}-\frac{\dot{\beta}_{qq^{\prime}}}{2}{n}_{\mu}^{q^{\prime}}-\frac{\beta_{qq^{\prime}}}{2}{n}_{\mu}^{q^{\prime}}\theta-\left(\delta_{n\Pi}^{q}\nabla_{\mu}{\Pi}-(1-\lambda_{\Pi}^{q}){\Pi}\nabla_{\mu}\delta_{n\Pi}^{q}\right)-\delta_{n\pi}^{q}\nabla^{\nu}{\pi}_{\mu\nu}\right.$
$\displaystyle-(1-\lambda_{\pi}^{q}){\pi}_{\mu\nu}\nabla^{\nu}\delta_{n\pi}^{q}\left.\right\\}\,,$
(16)

where while grouping we redistributed the terms
$\lambda_{\Pi}^{q}{n^{\mu}_{q}\Pi}\nabla_{\mu}\delta_{n\Pi}^{q}$ and
$(1-\lambda_{\Pi}^{q}){n^{\mu}_{q}\Pi}\nabla_{\mu}\delta_{n\Pi}^{q}$ between
the bulk and the charge currents relaxation equations are to account for the
arbitrary contribution of the second order coupling terms $n^{\mu}_{q}\Pi$,
and similarly for $n^{\mu}_{q}\pi_{\mu\nu}$.

Both $\lambda_{\Pi}^{q},\lambda_{\pi}^{q}$ can vary in the range between
$\\{0,1\\}$ to set their relevant weight. Note the the coupling terms
$n^{\mu}_{q}\Pi$ and $n^{\mu}_{q}\pi_{\mu\nu}$ have been ignored previously
for studies of baryon diffusion e.g . Denicol  _et al._ (2018); however, there
has been no check made to determine their role or significance. We include
these terms here for their relevance especially when diffusion exists for
multiple charges.

Imposing the second law of thermodynamics guarantee positive definiteness of
the entropy production which, therefore, must be quadratic in dissipative
currents leading to a non-negative divergence of the entropy current
$S^{\mu}$. This requires making everything in ”big brackets” in Eq. (16) equal
to its corresponding dissipative current such that

$\displaystyle\nabla_{\mu}S^{\mu}=\frac{\Pi\Pi}{\zeta
T}+\frac{\pi^{\mu\nu}\pi_{\mu\nu}}{2\eta
T}+\frac{n^{\mu}_{q}n_{\mu}^{q^{\prime}}}{\kappa_{qq^{\prime}}T^{2}}\,\geq
0\,$ (17)

Thus, to ensure positive entropy produce, we define

$\displaystyle\Pi$
$\displaystyle=\zeta\left\\{-\theta-\frac{\beta_{\Pi}}{\beta}\dot{{\Pi}}-\frac{\dot{\beta}_{\Pi}}{2\beta}{\Pi}-\frac{\beta_{\Pi}}{2\beta}{\Pi}\theta-\frac{1}{\beta}\sum_{q}\left(\delta_{n\Pi}^{q}\partial_{\mu}{n}^{\mu}_{q}+\lambda_{\Pi}^{q}{n}_{q}^{\mu}\nabla_{\mu}\delta_{n\Pi}^{q}\right)-\frac{\delta_{\Pi\pi}}{2\beta}\pi_{\mu\nu}\sigma_{\mu\nu}\right\\}\,,$
(18) $\displaystyle{\pi}^{\mu\nu}$
$\displaystyle=2\eta\left\\{\sigma^{\mu\nu}-\frac{\beta_{\pi}}{\beta}\dot{{\pi}}^{\mu\nu}-\frac{\dot{\beta}_{\pi}}{2\beta}{\pi}^{\mu\nu}+\frac{\beta_{\pi}}{2\beta}{\pi}^{\mu\nu}\theta-\frac{1}{\beta}\sum_{q}\left(\delta^{q}_{n\pi}\nabla^{\langle\mu}n_{q}^{\nu\rangle}+\lambda_{\pi}^{q}n_{q}^{\langle\mu}\nabla^{\nu\rangle}\delta^{q}_{n\pi})\right)-\frac{\delta_{\Pi\pi}}{2\beta}\Pi\sigma_{\mu\nu}\right\\}\,,$
(19) $\displaystyle{n}_{\mu}^{q}$
$\displaystyle=\sum_{q^{\prime}}\kappa_{qq^{\prime}}\left\\{-\nabla_{\mu}\alpha_{q^{\prime}}-\beta_{lq^{\prime}}\dot{{n}}_{\mu}^{l}-\frac{\dot{\beta}_{lq^{\prime}}}{2}{n}_{\mu}^{l}-\frac{\beta_{lq^{\prime}}}{2}{n}_{\mu}^{l}\theta-\delta_{n\Pi}^{q^{\prime}}\nabla_{\mu}{\Pi}-\delta_{n\pi}^{q^{\prime}}\nabla^{\nu}{\pi}_{\mu\nu}-\tilde{\lambda}_{\Pi}^{q^{\prime}}{\Pi}\nabla_{\mu}\delta_{0}^{q^{\prime}}-\tilde{\lambda}_{\pi}^{q^{\prime}}{\pi}_{\mu\nu}\nabla^{\nu}\delta_{1}^{q^{\prime}}\right\\}\,,$
(20)

By construction, Eqs. (18-20) should be constrained to reduce to the Navier-
Stokes (NS) equations when the higher order terms in the entropy vanish
leading to

$\displaystyle\Pi^{NS}\equiv-\zeta\theta\,,\qquad\qquad\pi_{\mu\nu}^{NS}\equiv
2\eta\sigma_{\mu\nu}\,,\qquad\qquad\left(n^{\mu}_{q}\right)^{NS}\equiv\kappa_{qq^{\prime}}\nabla^{\mu}\alpha_{q^{\prime}},$
(21)

which leads to the definition of the relaxation times for the shear, bulk, and
diffusion and transport coefficients,

$\displaystyle\tau_{\pi}=2\beta_{\pi}\eta/\beta,\,\qquad\qquad\tau_{\Pi}=\beta_{\Pi}\zeta/\beta,\qquad\qquad\,\,\tau_{qq^{\prime}}=\beta_{ql}\kappa_{q^{\prime}l}/\beta\,$
(22)

where now both $\tau_{i}$ and the transport coefficients go to zero as the
second order terms vanish and the system approaches its local equilibrium
state. Note that $\tau_{qq^{\prime}}$ is now a $3\times 3$ matrix defined via
Eq. (22). We had made the assumption that $\kappa_{qq^{\prime}}$ is an
invertible matrix in Eq.(17) which combined with the condition that
$\beta_{qq^{\prime}}$ is also an invertible matrix defines the relaxation time
matrix for diffusion to be invertible as well. Note that, physically, all of
these matrices should be real and symmetric and therefore the assumption of
invertability makes sense.

Finally, we obtain the equations of motion for the bulk pressure $\Pi$,

$\displaystyle\tau_{\Pi}\dot{\Pi}+\Pi$
$\displaystyle=-(\zeta+\frac{\tau_{\Pi}}{2}\Pi)\theta-\frac{\tau_{\Pi}}{2\beta_{\Pi}}\dot{\beta}_{\Pi}\Pi-\frac{\zeta\delta_{n\Pi}^{q}}{\beta}\partial_{\mu}n^{\mu}_{q}+\frac{\zeta\lambda_{\Pi}^{q}}{\beta}n_{q}^{\mu}\nabla_{\mu}\delta_{n\Pi}^{q}-\frac{\zeta\delta_{\Pi\pi}}{2\beta}\pi^{\mu\nu}\sigma_{\mu\nu}\,,$
(23)

the shear stress tensor $\pi^{\mu\nu}$,

$\displaystyle\tau_{\pi}\dot{\pi}^{\mu\nu}+\pi^{\mu\nu}$
$\displaystyle=2\eta\sigma^{\mu\nu}+\frac{\tau_{\pi}}{2}\pi^{\mu\nu}\theta-\frac{\tau_{\pi}\dot{\beta}_{\pi}}{2\beta_{\pi}}\pi^{\mu\nu}-\frac{2\eta\delta^{q}_{n\pi}}{\beta}\nabla^{\langle\mu}n_{q}^{\nu\rangle}+\frac{2\eta\lambda_{\pi}^{q}}{\beta}n_{q}^{\langle\mu}\nabla^{\nu\rangle}\delta^{q}_{n\pi}-\frac{\eta\delta_{\Pi\pi}}{\beta}\Pi\sigma_{\mu\nu}\,,$
(24)

and the BSQ diffusive charge currents $n^{\mu}_{q}$,

$\displaystyle\tau_{qq^{\prime}}\dot{n}^{\mu}_{q^{\prime}}+n^{\mu}_{q}$
$\displaystyle=-\kappa_{qq^{\prime}}\nabla^{\mu}\alpha_{q^{\prime}}+\frac{\tau_{qq^{\prime}}n^{\mu}_{q^{\prime}}}{2\beta}\theta-\frac{\tau_{qq^{\prime}}}{2\beta}\dot{\beta}_{q^{\prime}l}n^{\mu}_{l}-\frac{\kappa_{qq^{\prime}}\delta_{n\Pi}^{q^{\prime}}}{\beta}\nabla^{\mu}\Pi-\frac{\kappa_{qq^{\prime}}\delta_{n\pi}^{q^{\prime}}}{\beta}\nabla_{\nu}\pi^{\mu\nu}$
$\displaystyle-\frac{\kappa_{qq^{\prime}}\tilde{\lambda}_{\Pi}^{q^{\prime}}}{\beta}\Pi\nabla^{\mu}\delta_{n\Pi}^{q^{\prime}}-\frac{\kappa_{qq^{\prime}}\tilde{\lambda}_{\pi}^{q^{\prime}}}{\beta}\pi^{\mu\nu}\nabla_{\nu}\delta_{n\pi}^{q^{\prime}}\,.$
(25)

where $\\{\tau_{\Pi},\tau_{\pi},\tau_{qq^{\prime}}\\}$ are the microscopic
time scale of the dissipative fluxes; therefore, Eqs. (23-25) are relaxation
type equations of the viscous and diffusive currents, where basically, the
effective fluid description of the system emerges as the these currents relax
towards their equilibrium counterparts.

These equations, in addition to the conservation laws Eqs. (8-10) provide the
space-time evolution for the transient Israel-Stewart theory multi-component
fluid with BSQ charges. An addition of the formalism we introduced here is
introduction of the $\pi_{\mu\nu}\Pi$ term which we discuss further in Sec.
(II.4). In this formalism, we focus on the proper treatment of the second
order coupling terms between shear-bulk, shear-diffusion, and bulk-diffusion.
This will allow us to systematically scrutinize and understand the interplay
between these dynamics. Finally, having included the full BSQ charges mixes
the time derivatives of the charge currents through the term
$\tau_{qq^{\prime}}\dot{n}^{\mu}_{q^{\prime}}$ which renders these dynamics
sensitive to the chemical content of the fluid.

###  II.2 Thermodynamics input: Equation of State

The equations of motion up to this point are generic for any charged fluid and
carry no information on the type of matter (beyond the number of conserved
charges within the system) nor the type of microscopic interactions the
constituents of the system will have in order to evolve back to equilibrium.
The relevant information for the thermodynamics of the QCD matter and the
transport properties are provided through the Equation of State (EoS) and
transport coefficients, respectively. Ideally, the transport coefficients
would be linked directly to the EoS such that one could calculate the 4D EoS
directly from lattice QCD and all corresponding transport coefficients. The
reality is that not all EoS can directly provide the transport coefficients
(for instance, the 4D lattice QCD reconstructed EoS is only a Taylor series
Noronha-Hostler  _et al._ (2019); Monnai  _et al._ (2019) and does not provide
any information about transport coefficients).

For 3 conserved charges, one requires a 4-dimensional (4D) EoS since one then
must include temperature and BSQ chemical potentials
$\left(T,\mu_{B},\mu_{S},\mu_{Q}\right)$. For most EoS models, it is more
natural to calculate thermodynamic observables along grids of
$\left(T,\mu_{B},\mu_{S},\mu_{Q}\right)$. Thus, it is also significantly
easier to take derivatives along trajectories across
$\left(T,\mu_{B},\mu_{S},\mu_{Q}\right)$, rather than along lines of constant
entropy, $s$, for instance. However, relativistic hydrodynamics does not
directly obtain $\left(T,\mu_{B},\mu_{S},\mu_{Q}\right)$, instead it is more
natural to use energy density $\varepsilon$ or entropy density $s$ combined
with the 3 BSQ densities i.e. $\left(\varepsilon,n_{B},n_{S},n_{Q}\right)$ or
$\left(s,n_{B},n_{S},n_{Q}\right)$. Thus, there is always a mapping issue
wherein one must translate from the more natural variables of the EoS
$\left(T,\mu_{B},\mu_{S},\mu_{Q}\right)$ into the more natural variables for
hydrodynamics $\left(s,n_{B},n_{S},n_{Q}\right)$.

###  II.3 Transport coefficients

Starting with the bulk $\Pi$ equations of motion Eq. (23) , we have 10
independent transport coefficients for 3 conserved charges (note we count
$\tau_{\Pi}$ and $\zeta$ separately even though they can be related because
there is a coefficient that depends on the microscopic framework). If we have
X number of conserved charges, then the number of transport coefficients for
$\Pi$ is $N^{\Pi}_{trans\;coef.}=4+2\cdot X$ such that for $X=3$, we obtain
$N^{\Pi}_{trans\;coef.}=10$. For the shear stress tensor $\pi^{\mu\nu}$
relaxation equation Eq. (24), we again have 10 transport coefficients,
however, one also appears in the $\Pi$ equation of motion so we only have 9
new independent transport coefficients. We can then write the formula for the
number of $\pi^{\mu\nu}$ transport coefficients as
$N^{\pi^{\mu\nu}}_{trans\;coef.}=4+2\cdot X$ and if $\textrm{$\zeta>0$
then}\;N^{\pi^{\mu\nu}}_{trans\;coef.}=N^{\pi^{\mu\nu}}_{trans\;coef.}-1$ such
that if the bulk viscosity is nonzero then we subtract one transport
coefficient because $\delta_{\pi\Pi}$ appears in both equations of motion.
Finally, for the diffusion equations 25, if we have X number of conserved
charges, we gain new diffusion transport coefficients according to
$N^{diff}_{trans\;coef.}=2\cdot\left(\sum_{i=1}^{X}i\right)+X\cdot 4$ such
that for $X=3$ conserved charges we end up 24 tranport coefficients. Our final
number of transport coefficients is then

$N_{trans\;coef.}=N^{\Pi}_{trans\;coef.}+\left(N^{\pi^{\mu\nu}}_{trans\;coef.}-1\right)+N^{diff}_{trans\;coef.}$
(26)

where for the BSQ conserved charges we have 10+9+24=43 independent transport
coefficients.

The transport coefficients
$\\{\beta_{i},\tau_{i},\delta_{i}\\}\,{\forall}\,i\in\\{\Pi,\pi^{\mu\nu},n^{\mu}\\}$
are dynamical quantities and it is not possible to know the functional form of
an associated transport coefficients based on the Israel-Stewart theory.
Instead, we require input from a microscopic theory. Calculations of transport
coefficients has been done in a weakly-coupled approach from kinetic theory
Denicol  _et al._ (2012b, a); Molnár  _et al._ (2014); Fotakis  _et al._
(2020); Greif  _et al._ (2018); Rose  _et al._ (2020), hadron resonance gas
Noronha-Hostler  _et al._ (2009); Pal (2010); Khvorostukhin  _et al._ (2010);
Tawfik and Wahba (2010); Alba  _et al._ (2015); Ratti  _et al._ (2011); Tiwari
_et al._ (2012); Noronha-Hostler  _et al._ (2012); Kadam and Mishra (2014,
2015, 2016); Kadam and Pawar (2019); Mohapatra  _et al._ (2019), or transport
theory Wesp  _et al._ (2011); Ozvenchuk  _et al._ (2013); Rose  _et al._
(2018); Rais  _et al._ (2020) as well strongly coupled approach from
holography Kovtun  _et al._ (2005); Casalderrey-Solana  _et al._ (2014);
Finazzo and Noronha (2014); Finazzo  _et al._ (2015); Rougemont  _et al._
(2015, 2017); Grefa  _et al._ (2022). Current holographic calculations include
one conserved charge such as baryon conservation but can at least consider
small perturbations in the strangeness or electric charge sector. On the
kinetic theory side, recent effort have been done to understand and calculate
the dependence of the second order transport coefficients on chemical
potential. The final form for the relaxation times will rely on the approach
we take to extract the transport coefficients which in turn will directly
affect the dynamics and trajectories. Here we list the functional form of the
transport coefficients and arrange the proper matching with our equations.

In the case of $\mu_{q}=0$, the shear and bulk relaxation times are provided
by Denicol  _et al._ (2014)

$\displaystyle\frac{\zeta}{\tau_{\Pi}}$
$\displaystyle=15\left(\frac{1}{3}-c^{2}_{s}\right)^{2}(\varepsilon+p)\,,\qquad\qquad\qquad\frac{\eta}{\tau_{\pi}}=\frac{\varepsilon+p}{5}\,,$
(27)

And for the second order transport coefficients

$\displaystyle\frac{\delta_{\Pi\Pi}}{\tau_{\Pi}}=\frac{2}{3};\qquad\qquad\qquad\frac{\delta_{\pi\pi}}{\tau_{\pi}}$
$\displaystyle=\frac{4}{3}\,;\qquad\qquad\qquad\quad\frac{\lambda_{\pi\Pi}}{\tau_{\pi}}=\frac{6}{5};$
(28) $\displaystyle\frac{\tau_{\pi\pi}}{\tau_{\pi}}$
$\displaystyle=\frac{10}{7};\qquad\qquad\quad\quad\frac{\lambda_{\Pi\pi}}{\tau_{\Pi}}=\frac{8}{5}\left(\frac{1}{3}-c^{2}_{s}\right)\,.$

For a one component massless gas in the Boltzmann limit Denicol  _et al._
(2018)

$\displaystyle\frac{\kappa_{B}}{\tau_{B}}$
$\displaystyle=\,n_{B}\left(\frac{1}{3}\coth\left(\alpha_{B}\right)-\frac{n_{B}T}{\varepsilon+p}\right).$
(29)

where at this point, only diagonal terms of the diffusion matrix available
from kinetic theory calculations Fotakis  _et al._ (2022). A particular
feature of the Israel-Stewart formalism is the inclusion of time and spatial
derivative of thermodynamic potentials $\\{\alpha_{q}\\}$ as well as
derivatives of the second order transport coefficients
$\\{\beta_{\Pi},\beta_{\pi},\beta_{qq^{\prime}},\delta_{0}^{q},\delta_{1}^{q}\\}$.
These transport coefficients turns out to be complex combinations of
thermodynamics quantities for the case of finite charges. While these terms
have been ignored mostly in the literature, Dore  _et al._ (2020) showed that
they could impact the trajectories of the dynamical evolution. This might end
up an important effect when it comes to the search for critical point;
therefore, we take a step forward and calculate these thermodynamic derivative
quantities. To do that an additional thermodynamics input from the equation of
state is needed, but in principle calculate therms like ${\dot{\beta_{\Pi}}}$
with the aid of Eq. (8). Doing so we can see that it incorporates some higher
order terms beyond second order in dissipative currents. In addition, these
terms could be in particular important at the transition region, for instance,
at the QCD critical point it is expected to have $c_{s}^{2}\rightarrow 0$
(e.g. see Parotto  _et al._ (2020); Karthein  _et al._ (2021)). Since the
$\zeta$ scales with $\left(\frac{1}{3}-c^{2}_{s}\right)^{2}$ then a dip in
$c_{s}^{2}$ leads to a peak in $\zeta$. As we will see in Sec. III, our choice
of the underlying microscopic interactions will definitely impact the regime
of applicability of these equations through their dependence and connection to
the EoS. Whether the analysis will be biased to a particular set or not is
something we will explore in the future.

###  II.4 Connections with DNMR

To make connections to the DNMR framework Denicol  _et al._ (2012a), we start
by writing the general form of the relaxation equations as a systematic
expansion in Knudsen and inverse Reynolds numbers such as

$\displaystyle\dot{Re}_{i}^{-1}+Re_{i}^{-1}$ $\displaystyle=\textrm{Navier-
Stokes}+{\mathcal{J}}+{\mathcal{K}}+{\mathcal{R}}\,,$ (30)

where $Re_{i},Kn_{i}$ are the Reynolds and Knudsen numbers for the particular
dissipative current i$\in\\{\Pi,\pi^{\mu\nu},n^{\mu}_{q}\\}$. The
${\mathcal{J}}$, ${\mathcal{K}}$, and ${\mathcal{R}}$ terms are proportional
to $Kn_{i}\,Re_{i}^{-1}$, $Kn_{i}Kn_{j}$, and $Re_{i}^{-1}Re_{j}^{-1}$,
respectively. The general form of the relaxation-type equations in IS theories
has particular features which differ from DNMR, and we would like to discuss
some of them in this section.

Consider for example the bulk relaxation equation Eq. (23),

$\displaystyle\tau_{\Pi}\dot{\Pi}+\Pi$
$\displaystyle=-(\zeta+\frac{\tau_{\Pi}}{2}\Pi)\theta-\frac{\tau_{\Pi}}{2\beta_{\Pi}}\dot{\beta}_{\Pi}\Pi-\frac{\zeta\delta_{n\Pi}^{q}}{\beta}\partial_{\mu}n^{\mu}_{q}+\frac{\zeta\lambda_{\Pi}^{q}}{\beta}n_{q}^{\mu}\nabla_{\mu}\delta_{n\Pi}^{q}-\frac{\zeta\delta_{\Pi\pi}}{2\beta}\pi^{\mu\nu}\sigma_{\mu\nu}\,,$
(31)

For convenience, we write the equation in the following form,

$\displaystyle\tau_{\Pi}\dot{Re}^{-1}_{\Pi}+Re^{-1}_{\Pi}$
$\displaystyle=K_{\Pi}-\frac{\tau_{\Pi}\dot{\beta_{\Pi}}}{2\beta_{\Pi}}Re^{-1}_{\Pi}+\frac{\tau_{\Pi}}{2}K_{\Pi}Re^{-1}_{\Pi}+\frac{\zeta\lambda_{\Pi}^{q}\nabla_{\mu}\delta_{n\Pi}^{q}}{\beta}Re^{-1}_{n_{q}}-\frac{\zeta\delta_{n\Pi}}{\beta}\partial_{\mu}{Re}^{-1}_{n_{q}}-\frac{\zeta\delta_{\Pi\pi}}{2\beta}K_{\pi}Re^{-1}_{\pi}\,,$
(32)

with $K_{\Pi}\equiv\frac{\tau_{\Pi}}{\theta}$ and $Re_{\Pi}^{-1}\propto\Pi$.
It has been well known previously that IS theories fail to reproduce the
systematic expansion of DNMR theories. DNMR theories on the other hand does
not have those terms with time derivatives of thermodynamics potentials
$\\{\dot{\beta}_{\Pi},\dot{\beta}_{\pi},\dot{\beta}_{qq^{\prime}}\\}$.

In addition, the DNMR formalism implements an irreducible set of basis to
build the moments of the distribution function out of equilibrium. This
guarantees that dissipative fluxes created from the moments of the Boltzmann
equation are independent and orthogonal set of dynamical degrees of freedom.
This feature has not been enforced in IS theories where no symmetry arguments
was imposed on the dissipative fields.

The authors of DNMR have addressed in Denicol  _et al._ (2012a); Molnár  _et
al._ (2014) the relative importance of the above ${\mathcal{J}}$,
${\mathcal{K}}$, and ${\mathcal{R}}$ terms in their framework. The
${\mathcal{K}}$ terms were shown to render the equations of motion non-
hyperbolic leading to instability of the evolution, and this issue would be
recovered only by including higher order terms. The ${\mathcal{R}}$ terms on
the other hand arise from the non-linearities of the collision kernel in the
Boltzmann equation which signals their potential relevance in the far from
equilibrium regime; otherwise, it was shown in Molnár  _et al._ (2014) that,
in the Boltzmann limit of mass-less gas with a constant cross section, these
contribution can be ignored. This leaves us with the ${\mathcal{J}}$ $\propto$
$Kn_{i}\,Re_{j}^{-1}$ which has been used in most codes that implement DNMR
theory. Interestingly, in the Israel-Stewart approach, one can actually still
reproduce at least a subset of the DNMR terms by including contributions to
the entropy production which has a zero contribution to the entropy current
itself. An example is the term $\propto\pi_{\mu\nu}\sigma^{\mu\nu}$ from Eq.
(32) which does not exist in the standard IS formalism. It is not clear to us
why these important terms in DNMR comes up from zero net-entropy current terms
in the IS formalism.

One interpretation of the shear-bulk coupling term in the equations of motion
is that it originates from a thermodynamically reversible process i.e. this
process should be symmetric under reversing the flow velocity. For this
reason, one would just associate a zero transport coefficient for this term at
the level of the entropy current. Provided the above argument, one could look
at the level of entropy production where the 2nd law does not uniquely specify
the definition of entropy production and it only constrain it to be positive
definite. In this sense, we have some freedom in adding and subtracting terms
as long as we respect the order of the entropy current. Other zero entropy
terms can be treated on the same footing to generate the remaining
${\mathcal{J}}$ terms. For the term $\propto\pi_{\mu\nu}\sigma^{\mu\nu}$, it
was interesting to see that it is actually the original first order NS term
where as $\tau_{\pi}\rightarrow\infty$ and as the shear stress approaches it
NS limit $\pi_{\mu\nu}=2\eta\sigma_{\mu\nu}$ the above contribution becomes
$\propto\sigma_{\mu\nu}\sigma^{\mu\nu}$.

##  III Stability analysis of the BSQ hydrodynamics

A critical benchmark for any hydrodynamics theory is to pass the stability
test, which warrant the existence of a steady stable equilibrium state. Any “
small ” perturbations around this equilibrium background should decay quickly
and the system should asymptotically converge back to its equilibrium state.
Growth of perturbations, on the other hand, signals an instability of the
equilibrium state which then means that the configuration space cannot be
physically bounded. For Israel-Stewart theory, a linear stability analysis has
been done in both Landau and Eckart frames for the case of a single charge
Hiscock and Lindblom (1983); Olson (1990). In these works, the authors showed
that in the linear regime Israel-Stewart theories are causal if they are
stable. In a more recent work, Gavassino et al Gavassino  _et al._ (2022a),
further studied the relationship between causality and stability and
demonstrated that while causality alone does not guarantee stability, a theory
can be proven to be causal if it is thermodynamically stable.

In this section, we follow the method developed in Gavassino (2021) to (1)
derive the Lyapunov functional for the entropy current in the Landau frame and
generalize it to a fluid with multiple charges, and (2) systematically extract
the conditions for the positivity of the energy functional. We finally proceed
with the potential importance of the results by (a) making connections with
another fundamental concept of the physical system which is the causality of
signal propagation in this theory, and (b) check if the analysis is biased to
a particular choice of transport coefficients. The entropy current in Eq. (13)
accounts for all dissipative tensors up to second order as well as all
possible coupling terms with non-zero entropy production. The total entropy of
a space-like hyper-surface $\Sigma$ is given by

$\displaystyle S[\Sigma]$ $\displaystyle=\int_{\Sigma}{S}^{\mu}d\Sigma_{\mu},$
(33)

where for a perfect fluid this entropy current is conserved; however,
dissipation processes and diffusion lead to an increase in the entropy. For an
isolated system, and starting from a non-equilibrium state, the total entropy
is a Lyapounov function and must be positive definite function of time. If we
evaluate Eq. (33) on any later surface $\Sigma^{\prime}$ its value must
increase

$\displaystyle S(\Sigma^{\prime})-S(\Sigma)$
$\displaystyle=\int_{\Sigma}\nabla_{\mu}{S}^{\mu}dV\,.$ (34)

which means as the system evolve, the entropy increases until it reaches its
maximum value. This maximum is unique and is achieved only for the global
thermal equilibrium state. In this sense, the equilibrium state is an
attractor for all thermodynamically consistent configurations of the non-
equilibrium states. If we now “ linearly ” perturb the global thermodynamic
state with a small perturbation, we can quantify the entropy difference
between the perturbed and non-perturbed state as

$\displaystyle\Phi_{E}\equiv S(\Sigma)-\tilde{S}(\Sigma)$
$\displaystyle=\int_{\Sigma}{E}^{\mu}d\Sigma_{\mu}\,,$ (35)

where the “$\,\,\tilde{}\,\,$” refers to the entropy of the perturbed state,
$\Phi_{E}$ is the flux of the current $E^{\mu}$ which carries information on
the magnitude and propagation of perturbations. Our goal then is to quantify
the conditions that bound the fluctuations from growing indefinitely and also
to drive the system back to its equilibrium state as $\,t\rightarrow\infty$.

We start by first linearly perturbing our hydrodynamical fields about the
state of equilibrium up to second order in fluctuations. Then, defining the
equilibrium state by setting all disspative currents to zero, see App. B,
gives the following

$\displaystyle\delta S^{\mu}_{non-eq}$
$\displaystyle=\left(\varepsilon+p-\mu_{q}\,\rho_{q}\right)\delta
u^{\mu}+\left(\delta\varepsilon-\mu_{q}\,\delta\rho_{q}-\delta\mu_{q}\,\delta\rho_{q}+\delta
s\delta
T\right)u^{\mu}+\left(\delta\varepsilon-\mu_{q}\,\delta\rho_{q}\right)\delta
u^{\mu}$ $\displaystyle-\frac{\mu_{q}}{T}\delta
n_{q}^{\mu}-\frac{1}{2}u^{\mu}\left(\beta_{\Pi}\delta\Pi\delta\Pi+\beta_{\pi}\delta{\pi}^{\mu\nu}\delta{\pi}_{\mu\nu}+\sum_{q}\beta_{n}^{qq^{\prime}}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}\right)$
$\displaystyle-\sum_{q}\left(\delta_{n\Pi}^{q}\delta
n^{\mu}\delta\Pi+\delta_{n\pi}^{q}\delta
n_{\nu}^{q}\delta\pi^{\mu\nu}\right)\,,$ (36)

here the $\delta$ refer to the perturbed field. To drive the system to its
equilibrium state, these perturbations must strictly conserve both energy-
momentum and the BSQ charge currents during the full evolution. Therefore, we
impose a sort of “ landau matching” on the fluctuations. Starting from the
perturbed form of the landau condition $u_{\mu}T^{\mu\nu}=\varepsilon u^{\nu}$
one gets

$\displaystyle u_{\nu}\delta T^{\mu\nu}$ $\displaystyle=\varepsilon\delta
u^{\mu}+\delta\varepsilon u^{\mu}+\delta\varepsilon\delta u^{\mu}-\delta
u_{\nu}T^{\mu\nu}-\delta u_{\nu}\delta T^{\mu\nu}\,,$ (37)

which provides the zero energy flux contributions to the entropy fluctuations.
Likewise, the charge current $N^{\mu}_{q}=\rho u^{\mu}+n^{\mu}_{q}$ can be
expanded in perturbed form to give

$\displaystyle\delta N^{\mu}_{q}$ $\displaystyle=\rho_{q}\,\delta
u_{\mu}+\delta\rho_{q}\,u_{\mu}+\delta\rho_{q}\,\delta u_{\mu}+\delta
n^{\mu}_{q}\,,$ (38)

which provides the zero charge flux contributions to the entropy fluctuations.
Note that this condition will impose conditions on each of the charges
separately Subtracting out the zero flux contributions using the killing
vector and conservation laws Eq. (7), we can finally define the information
current $E^{\mu}$ such that

$\displaystyle E^{\mu}$ $\displaystyle=\delta s^{\mu}-\textrm{zero energy
flux}-\textrm{zero number flux }\,,$ (39)

where now using the perturbed entropy form Eq. (36) then gives the final form
of the energy functional

$\displaystyle E^{\mu}$ $\displaystyle=-(\varepsilon+{p})\delta u^{\nu}\delta
u_{\nu}u^{\mu}+\delta{p}\,\delta
u^{\mu}+\frac{u^{\mu}}{2(\varepsilon+{p})}\left.\frac{\partial\varepsilon}{\partial{p}}\right|_{s}(\delta{p})^{2}+\frac{u^{\mu}}{2(\varepsilon+{p})}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial{p}}{\partial s}\right|_{\alpha_{q}}(\delta
s)^{2}+\left.\frac{\partial\alpha_{q}}{\partial{p}}\right|_{s}\delta{p}\,\delta
n_{q}^{\mu}$
$\displaystyle+\left.\frac{\partial}{\partial\alpha_{q}}\right|_{p}\delta
s\,\delta n_{q}^{\mu}+\delta\Pi\,\delta u^{\mu}+\delta\pi^{\mu\nu}\delta
u_{\nu}+\frac{1}{2}u^{\mu}\left(\beta_{\Pi}\delta\Pi\delta\Pi+\beta_{\pi}\delta{\pi}^{\mu\nu}\delta{\pi}_{\mu\nu}+\sum_{q}\beta_{n}^{qq^{\prime}}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}\right)$
$\displaystyle+\sum_{q}\left(\delta_{n\Pi}^{q}\delta
n^{\mu}\delta\Pi+\delta_{n\pi}^{q}\delta
n_{\nu}^{q}\delta\pi^{\mu\nu}\right)\,.$ (40)

The information current $E^{\mu}$ carries all information about the
perturbation and can be written as

$\displaystyle\nabla_{\mu}E^{\mu}=-\left[\frac{\delta\Pi\delta\Pi}{\zeta
T}+\frac{\delta\pi^{\mu\nu}\delta\pi_{\mu\nu}}{2\eta T}+\frac{\delta
n^{\mu}_{q}\delta n_{\mu}^{q}}{\kappa_{qq^{\prime}}T^{2}}\right]\,,$ (41)

where now one can state that

$\displaystyle\nabla_{\mu}E^{\mu}\leq 0;\qquad\iff\rm
provided\,if\,and\,only\,if\qquad\nabla_{\mu}S^{\mu}\geq 0\,,$ (42)

which is a direct consequence of the second law of thermodynamics. The
relation Eq. (42) is therefore interpreted as the loss of information about
the perturbations as the system evolves towards its equilibrium state which
minimizes its energy and maximizes its entropy. The total energy associated
with $E^{\mu}$ can then be found by integration of the above expression over
the space-like surface $\Sigma$

$\displaystyle E(\Sigma)$
$\displaystyle=\int_{\Sigma}E^{\mu}d\Sigma_{\mu}\qquad=\int
e\beta^{\mu}dx_{\mu}\,,$

here $\beta^{\mu}=u^{\mu}/T$ and e is defined as the energy density of
$E^{\mu}$ given by $e=\frac{TE^{\mu}n_{\mu}}{u^{\nu}n_{\nu}}$. The vector
$n_{\mu}$ is the future directed unit normal to $\Sigma$. Since $E$ can not
increase in time, we have that

$\displaystyle E^{\mu}n_{\mu}\geq 0\,,$ (43)

The stability of the theory is now contingent on the proof that the energy
density of the energy current Eq. (40) is positive definite.

###  III.1 Stability conditions

We include the detailed method of derivation and the steps for extracting the
conditions in App. B. After contracting the energy current $E^{\mu}$ with the
future directed time-like vector $n_{\mu}$, the energy density is given by

$\displaystyle e$ $\displaystyle=-(\varepsilon+{p})\delta u^{\nu}\delta
u_{\nu}+\delta{p}\,\delta
u^{\mu}\lambda_{\mu}+\frac{1}{2(\varepsilon+{p})}\left.\frac{\partial\varepsilon}{\partial{p}}\right|_{s}(\delta{p})^{2}+\frac{1}{2(\varepsilon+{p})}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial{p}}{\partial s}\right|_{\alpha_{q}}(\delta
s)^{2}+\left.\frac{\partial\alpha_{q}}{\partial{p}}\right|_{s}\delta{p}\,\delta\,n_{q}^{\mu}\lambda_{\mu}$
$\displaystyle+\left.\frac{\partial\alpha_{q}}{\partial s}\right|_{p}\delta
s\,\delta n_{q}^{\mu}\lambda_{\mu}+\delta\Pi\,\delta
u^{\mu}\lambda_{\mu}+\delta\pi^{\mu\nu}\delta
u_{\nu}\lambda_{\mu}+\frac{1}{2}\left(\beta_{\Pi}\delta\Pi\delta\Pi+\beta_{\pi}\delta{\pi}^{\mu\nu}\delta{\pi}_{\mu\nu}+\sum_{q}\beta_{n}^{qq^{\prime}}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}\right)$
$\displaystyle+\sum_{q}\left(\delta_{n\Pi}^{q}\lambda_{\mu}\delta
n^{\mu}\delta\Pi+\delta_{n\pi}^{q}\lambda_{\mu}\delta
n_{\nu}^{q}\delta\pi^{\mu\nu}\right)\,,$ (44)

here $\lambda^{\mu}=\frac{\Delta^{\mu\nu}n_{\nu}}{n_{c}u^{c}}$. In order to
prove the positivity of the energy functional, we first write it in terms of
linearly independent combinations. Then, we assure these combinations are
positive definite by completing the squares of all coupled terms.

To do so, we first write the energy functional Eq. (44) in the form of
$k\times k$ dimensional matrix where k is the dimensionality of the phase
space or state configurations. Thus, proving positivity of this matrix is a
direct proof of stability of the dynamical system.

$\left(\begin{array}[]{c}\delta u^{\mu}\\\ \delta p\\\ \delta s\\\ \delta
n_{q}^{\mu}\\\ \delta\Pi_{\mu\nu}\\\ \delta\pi\\\
\end{array}\right)^{T}\left(\begin{array}[]{ccc|ccccc}\frac{(\varepsilon+p)}{2}&\frac{\lambda_{\mu}}{2}&0&0&\frac{\lambda_{\mu}}{2}&\frac{\lambda_{\nu}}{2}\\\
\frac{\lambda_{\mu}}{2}&\frac{1}{(\varepsilon+p)}\left.\frac{\partial\varepsilon}{\partial
p}\right|_{s}&0&\frac{\lambda_{\mu}}{2}\left.\frac{\partial\alpha_{q}}{\partial
p}\right|_{s}&0&0\\\
0&0&\left.\frac{1}{(\varepsilon+p)}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\frac{\partial p}{\partial
s}\right|_{\alpha_{q}}&\frac{\lambda_{\mu}}{2}\left.\frac{\partial\alpha_{q}}{\partial
s}\right|_{p}&0&0&\\\ \\\ \hline\cr\\\
0&\frac{\lambda_{\mu}}{2}\left.\frac{\partial\alpha_{q}}{\partial
p}\right|_{s}&\frac{\lambda_{\mu}}{2}\left.\frac{\partial\alpha_{q}}{\partial
s}\right|_{p}&\beta_{n}^{qq^{\prime}}&\frac{\lambda_{\mu}}{2}\delta_{n\Pi}^{q}&\frac{\lambda_{\mu}}{2}\delta_{n\pi}^{q}\\\
\frac{\lambda_{\mu}}{2}&0&0&\delta_{n\Pi}^{q}\frac{\lambda_{\mu}}{2}&\beta_{\Pi}&0\\\
\frac{\lambda_{\nu}}{2}&0&0&\delta_{n\pi}^{q}\frac{\lambda_{\mu}}{2}&0&\beta_{\pi}\\\
\end{array}\right)\left(\begin{array}[]{c}\delta u^{\mu}\\\ \delta p\\\ \delta
s\\\ \delta n_{q}^{\mu}\\\ \delta\Pi_{\mu\nu}\\\ \delta\pi\\\
\end{array}\right)\vspace{0.5cm}$

For a symmetric matrix, the proof that the determinants of all upper left
matrices are positive suffice to prove its positivity.

The conditions relate in a non-trivial way and can be sorted into: expressions
combining thermodynamics quantities such as $(\varepsilon,P)$, thermodynamics
potentials and derivatives $(\alpha_{q},c_{s}^{2})$ as well as dissipative
corrections $(\Pi,\pi_{\mu\nu})$, and transport coefficients
$(\beta_{\Pi},\beta_{\pi},...)$. Here we only list the final conditions and
discuss their physical implications. a detailed derivation is done in Sec. B.
We classify the conditions as follow:

  * •

Thermodynamics constraints:

$\displaystyle
c1:\;\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
p}\right|_{s}\geq 0\,;\qquad\qquad\qquad
c2:\;\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial s}\right|_{\alpha_{q}}\geq 0\,.$
(45)

where Eq. (45) is the thermodynamic constraint on the speed of sound since
$c_{s}^{2}=\left.\frac{\partial p}{\partial\varepsilon}\right|_{s}$. For most
regions of the QCD phase diagram, it is reasonable to assume that the enthalpy
$w=\varepsilon+p\geq 0$ and $c_{s}^{2}\geq 0$ are positive such that these
relations in Eq. (45) are rather trivially satisfied. However, if one passes
through a metstable phase by spinodal decomposition then it is possible to
reach a regime where $c_{s}^{2}<0$ or $w<0$. The only configuration where
stability is maintained while $c_{s}^{2}$ is allowed to be negative is if the
enthalpy is negative i.e $\varepsilon+p<0$, which implies that $p<0$ and
$|\varepsilon|<|p|$ since the energy density must be positive. In this case,
the speed of sound is imaginary i.e. $c_{s}\in\Im$ and $c_{s}^{2}<0$. The
consequence of a negative pressure and $c_{s}^{2}$, is that it would imply
that $\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial s}\right|_{\alpha_{q}}<0$ must
be negative. Conversely, if we know that $c_{s}^{2}>0$, then $w>0$ such that
it implies $\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial s}\right|_{\alpha_{q}}>0$ must
be positive. Then, $\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial s}\right|_{\alpha_{q}}$ is a
thermodynamic condition on heat conduction.

  * •

Transport constraints:

$\displaystyle c3:\;\beta_{\Pi}\geq 0\,,\qquad\qquad\qquad
c4:\;\beta_{\pi}\geq 0\,;\qquad\qquad$ (46)

where $\lambda$ is the norm of $\lambda^{\mu}$. The first two conditions
constrain the shear and bulk viscosities to be positive.

  * •

Constraints connecting thermodynamics with transport:

These condition have a more complicated form; however, one can use the
previous inequalities to simplify them. For example,

$\displaystyle
c5:\;\frac{\varepsilon+p}{\lambda^{2}}\geq\frac{(\delta_{n\pi}^{q}+1)^{2}}{2\beta_{\pi}}-\frac{\rho_{q}^{2}\beta_{n}^{qq^{\prime}}}{\lambda^{2}}$
(47)

that nicely ties together a relation between the enthalpy $w=\varepsilon+p$,
diffusion matrix $\beta_{n}^{qq^{\prime}}$, and the coupling of diffusion to
the shear stress $\delta_{n\pi}^{q}$. Since the charge density and the
coupling coefficients are squared and $\beta_{\pi}$ is constrained to be
positive by Eq. (46), their sign does not directly impact the matrix elements
but their magnitude does. For example for $w>0$, the constraint c5 allows for
a negative $\beta_{n}^{qq^{\prime}}$ as long as the total
$\frac{(\delta_{n\pi}^{q}+1)^{2}}{2\beta_{\pi}}-\frac{\rho_{q}^{2}\beta_{n}^{qq^{\prime}}}{\lambda^{2}}$
is less than the entropy $w$.

In contrast, if $w<0$, then we have:

$\displaystyle\textrm{If
}\varepsilon+p<0,\;\textrm{then}\quad\frac{(\delta_{n\pi}^{q}+1)^{2}}{2\beta_{\pi}}<\rho_{q}^{2}\beta_{n}^{qq^{\prime}}\,\textrm{and}\;\beta_{n}^{qq^{\prime}}>0$
(48)

since the term $\frac{(\delta_{n\pi}^{q}+1)^{2}}{2\beta_{\pi}}$ is positive
definite due to c4. Then, Eq. (48) implies that for a meta-stable region we
must have $\beta_{n}^{qq^{\prime}}>0$. Indeed, $\beta_{n}^{qq^{\prime}}$ will
have a dependence on the choice of the hyper-surface through the norm
$\lambda^{2}$; however, in the above relation, we used the most stringent
constraint by setting $\lambda^{2}\rightarrow 1$.

The next condition states that the relation between transport and diffusion
coefficients is constrained by the thermodynamics such that if we look at

$\displaystyle
c6:\;\frac{\beta_{n}^{qq^{\prime}}}{2\lambda^{2}}-\frac{({\delta_{n\pi}^{q}})^{2}}{\beta_{\pi}}-\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}\geq{\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}}-T^{2}(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}\,.$ (49)

One finds can find additional constraints from c6 in very specific scenarios.

Scenario 1 the speed of sound is $c_{s}^{2}>0$ and the right-hand side of Eq.
(49) is positive:

$\displaystyle\textrm{if}\;c_{s}^{2}>0\quad\textrm{and}\quad$
$\displaystyle\frac{1}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}\geq\underbrace{T^{2}\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}}_{>0\textrm{ from c2}}$ then
$\displaystyle\frac{\beta_{n}^{qq^{\prime}}}{2\lambda^{2}}\geq\frac{({\delta_{n\pi}^{q}})^{2}}{\beta_{\pi}}+\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}\,,$
(50)

that ensures the positivity of the right hand side of the inequality in Eq.
(49).

Scenario 2 the speed of sound is $c_{s}^{2}>0$ and the the right-hand side of
Eq. (49) is negative:

$\displaystyle\textrm{if}\;c_{s}^{2}>0\quad\textrm{and}\quad$
$\displaystyle\frac{1}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}<T^{2}\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}$ then
$\displaystyle\frac{\beta_{n}^{qq^{\prime}}}{2\lambda^{2}}\;\textrm{may be
positive or negative}$ (51)

that ensures the positivity of the right hand side of the inequality in Eq.
(49).

Scenario 3 for $c_{s}^{2}<0$ we know from Eq. (48) that
$\beta_{n}^{qq^{\prime}}>0$ such that there is always a possibility for the
left-hand side of Eq. (49) to be positive or negative, depends on if
$\beta_{n}^{qq^{\prime}}$ is larger or smaller than
$\frac{({\delta_{n\pi}^{q}})^{2}}{\beta_{\pi}}+\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}$.
If the right-hand side of Eq. (49) is negative, we do not have any additional
constraints. However, if the right-hand side of Eq. (49) is positive, then we
find that

$\displaystyle\textrm{if}\;c_{s}^{2}<0;\quad\textrm{and}$
$\displaystyle\frac{1}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}\geq T^{2}\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}$ (52) then
$\displaystyle\frac{\beta_{n}^{qq^{\prime}}}{2\lambda^{2}}\geq\frac{({\delta_{n\pi}^{q}})^{2}}{\beta_{\pi}}+\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}\,,$
(53)

where we have then find that regardless of the sign of $c_{s}^{2}$, we can
assume $\beta_{n}^{qq^{\prime}}>0$ if Scenario 1 or Scenario 3 are fulfilled.
However, it is possible for $\beta_{n}^{qq^{\prime}}<0$ in Scenario 2.

Our last condition is the most complex and comes from the coupling between the
charge current and the flow velocity which states the following inequality

$\displaystyle c7:\;$
$\displaystyle(\varepsilon+p)\bigg{(}\frac{1}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}-\frac{1}{2\lambda^{2}}\bigg{)}-\frac{1}{\beta_{\Pi}}-{\frac{2}{3\beta_{\pi}}}$
$\displaystyle\geq\frac{\\{-\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}-\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}-\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\\}^{2}}{\lambda^{2}\,\\{\frac{\beta_{n}^{qq^{\prime}}}{2\lambda^{2}}-\frac{({\delta_{n\pi}^{q}})^{2}}{\beta_{\pi}}-\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}-{\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}}+T^{2}(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}\\}}$ (54)

Note that we have shown the constraints on the positivity of the denominator
on the right hand side in the previous condition.

The constraints $c4,c5,c6,c7$ depend on the space-like vector
$\lambda^{2}=\lambda_{\mu}\lambda^{\mu}$ where different choices of the hyper-
surface will affect the particular value of this vector. As we show in the
text, one can enforce the strongest constraint by setting $\lambda^{2}=1$. We
finally note that the different structure of the conditions in the Landau
frame compared to the Ekart frame Hiscock and Lindblom (1983); Olson (1990)
comes from the fact that charges could diffuse independently and are not
restricted to co-move with the flow velocity. This creates a linear dependence
between the fluctuations of $\delta u^{\mu}$ and $\delta n^{\mu}_{q}$, which
we considered when deriving the constraints. While the general structure of
the conditions may seem different, it has been recently shown in Gavassino
_et al._ (2022b) that one can map from one set of conditions to the other by a
change of frame argument. A consequence of this effect for example is that it
relaxes the constraint on the diffusion matrix c5 to allow for negative
diffusion in some scenarios. One might still recover the positive definiteness
of $\beta_{qq^{\prime}}$ in the case of ideal BSQ hydrodynamics where the
charge fluctuations co-move with the energy flow.

##  IV Discussion and Conclusions

In this paper, we derived the relativistic viscous hydrodynamics equations of
motion with multiple conserved charges using the phenomenological approach
from Israel-Stewart. For the first time, we were able to recover the DNMR
$(K_{n}^{(1)}R_{n}^{(1)})$ terms from the phenomenological approach from
Israel-Stewart approach by adding new terms to the entropy production that
exactly cancel but that enter different dissipative equations of motion e.g.
the bulk pressure and the charge current equations. This theoretical
description is relevant for low energy heavy ion collisions and the search for
critical point/first-order phase transition and may be of interest for neutron
star mergers as well Most  _et al._ (2021, 2022).

After our equations of motion were derived, we performed a stability analysis
where the conditions are sufficient for hydrodynamic stability and are more
general than a mode analysis. As the work in Gavassino  _et al._ (2022a) show,
using the energy functional method, the conditions of thermodynamic stability
at linear order directly imply the linear causality. This means that the
conditions we derive here stand as penchmarks for stability as well as
causality of the BSQ Israel-Stewart theory in the Landau frame.

In this analysis, we find 7 stability constraints that either affect the
thermodynamic variables, the transport coefficients, or a mixature of the
two.We explore the possibility of the the meta-stable region across a first-
order phase transition where $c_{s}^{2}<0$ is possible and find that we should
still maintain hydrodynamic stability within that region as long as the
enthalpy is also negative and the Navier Stokes transport coefficients for
shear viscosity, bulk viscosity, and charge diffusion are positive definite.
Interestingly enough, under very specific conditions outside of the meta-
stable region (i.e. $c_{s}^{2}>0$) we can allow for a negative charge
diffusion.

At large baryon densities the role of transport coefficients are potentially
even more crucial than at $\mu_{B}=0$. Unfortunately, due to the fermion sign
problem we have no guidance from first principle lattice QCD calculations on
their behavior. That being said studies of critical scaling found that at the
critical point the bulk viscosity scales with the cubic power of the
correlation Moore and Saremi (2008); Monnai  _et al._ (2017); Martinez  _et
al._ (2019b) and is already large because of the dip in $c_{s}^{2}\rightarrow
0$ Dore  _et al._ (2020, 2022b). Across a first-order phase transition, one
expects a jump in the transport coefficients Soloveva  _et al._ (2021); Grefa
_et al._ (2022) (although, for a meta stable region things may be more
complicated and we are not aware of papers exploring that possibility yet in
the context of heavy-ion collisions). Thus, any constraints that are possible
on their behavior are extremely important for simulations of heavy-ion
collisions.

These conditions imply the causality of the solution at least in the linear
regime. As shown in Sec. III.1, the stability criteria creates a direct
correspondence between the EoS and the microscopic non-equilibrium physics.
The distinction from the non-linear causality constraints derived in Bemfica
_et al._ (2021) is that our stability constraints are expressed in terms of
thermodynamics and transport coefficients only while the conditions for
causality include the magnitude of the dissipative fluxes. In the IS
formalism, the transport coefficients set the magnitude of those
contributions, so a constraint on the transport coefficient already
constraints the relevant flux term.

One surprising outcome from this analysis is that we find that the terms that
only appear in DNMR but not previous works on Israel-Stewart have zero
contributions to the entropy current such that they do not affect
thermodynamic stability. From the second law of thermodynamics, irreversible
processes increase the entropy production. Since those terms do not increase
the entropy, they might not transfer any energy and momentum on the hyper-
surface; and therefore, they carry no information. The question then is what
role they play, and how do they arise from the Boltzmann equation in the DNMR
framework.

For the reasons mentioned above, the extra DNMR terms would only be
constrained from non-linear causality analysis. Previous work have found that
these extra terms better handle far-from-equilibrium effects that can occur at
the critical point Dore  _et al._ (2020) (i.e. the inclusion of these terms
more quickly reached an attractor) but that work also never found that the
lack of these terms lead to any thermodynamic instabilities. Thus, it appears
that these extra terms shift out-of-equilibriumness between different
dissipitative currents. For instance, imagine a large peak in the bulk
viscosity. Then these terms could shift some the enhancement over to the shear
stress tensor or charge currents to avoid issues with a large bulk viscosity
such that the out-of-equilibriumness is more balanced between different
disspitive currents. We leave a more concrete study of these terms to an
upcoming future work.

##  V Acknowledgements:

We thank G. S. Denicol and J. Noronha for helpful comments. J.N.H, T.D., D.A.
are supported by the US-DOE Nuclear Science Grant No. DE-SC0020633.

## References

  * Denicol  _et al._ (2018) G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C 98, 034916 (2018), arXiv:1804.10557 [nucl-th] .
  * Du and Heinz (2020) L. Du and U. Heinz, Comput. Phys. Commun. 251, 107090 (2020), arXiv:1906.11181 [nucl-th] .
  * Schäfer  _et al._ (2021) A. Schäfer, I. Karpenko, X.-Y. Wu, J. Hammelmann, and H. Elfner, (2021), arXiv:2112.08724 [hep-ph] .
  * Fotakis  _et al._ (2020) J. A. Fotakis, M. Greif, C. Greiner, G. S. Denicol, and H. Niemi, Phys. Rev. D 101, 076007 (2020), arXiv:1912.09103 [hep-ph] .
  * Du  _et al._ (2021) L. Du, X. An, and U. Heinz, Phys. Rev. C 104, 064904 (2021), arXiv:2107.02302 [hep-ph] .
  * Monnai  _et al._ (2021) A. Monnai, B. Schenke, and C. Shen, Int. J. Mod. Phys. A 36, 2130007 (2021), arXiv:2101.11591 [nucl-th] .
  * Martinez  _et al._ (2019a) M. Martinez, M. D. Sievert, D. E. Wertepny, and J. Noronha-Hostler, (2019a), arXiv:1911.10272 [nucl-th] .
  * Carzon  _et al._ (2022) P. Carzon, M. Martinez, M. D. Sievert, D. E. Wertepny, and J. Noronha-Hostler, Phys. Rev. C 105, 034908 (2022), arXiv:1911.12454 [nucl-th] .
  * Rougemont  _et al._ (2015) R. Rougemont, J. Noronha, and J. Noronha-Hostler, Phys. Rev. Lett. 115, 202301 (2015), arXiv:1507.06972 [hep-ph] .
  * Monnai  _et al._ (2017) A. Monnai, S. Mukherjee, and Y. Yin, Phys. Rev. C 95, 034902 (2017), arXiv:1606.00771 [nucl-th] .
  * Rougemont  _et al._ (2017) R. Rougemont, R. Critelli, J. Noronha-Hostler, J. Noronha, and C. Ratti, Phys. Rev. D 96, 014032 (2017), arXiv:1704.05558 [hep-ph] .
  * Grefa  _et al._ (2022) J. Grefa, M. Hippert, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, Phys. Rev. D 106, 034024 (2022), arXiv:2203.00139 [nucl-th] .
  * Martinez  _et al._ (2019b) M. Martinez, T. Schäfer, and V. Skokov, Phys. Rev. D 100, 074017 (2019b), arXiv:1906.11306 [hep-ph] .
  * Rajagopal  _et al._ (2020) K. Rajagopal, G. Ridgway, R. Weller, and Y. Yin, Phys. Rev. D 102, 094025 (2020), arXiv:1908.08539 [hep-ph] .
  * Dore  _et al._ (2020) T. Dore, J. Noronha-Hostler, and E. McLaughlin, Phys. Rev. D 102, 074017 (2020), arXiv:2007.15083 [nucl-th] .
  * Dore  _et al._ (2022a) T. Dore, J. M. Karthein, I. Long, D. Mroczek, J. Noronha-Hostler, P. Parotto, C. Ratti, and Y. Yamauchi, (2022a), arXiv:2207.04086 [nucl-th] .
  * Israel and Stewart (1979) W. Israel and J. M. Stewart, Annals Phys. 118, 341 (1979).
  * Denicol  _et al._ (2012a) G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, Phys. Rev. D 85, 114047 (2012a), [Erratum: Phys.Rev.D 91, 039902 (2015)], arXiv:1202.4551 [nucl-th] .
  * Bemfica  _et al._ (2018) F. S. Bemfica, M. M. Disconzi, and J. Noronha, Phys. Rev. D 98, 104064 (2018), arXiv:1708.06255 [gr-qc] .
  * Bemfica  _et al._ (2022) F. S. Bemfica, M. M. Disconzi, and J. Noronha, Phys. Rev. X 12, 021044 (2022), arXiv:2009.11388 [gr-qc] .
  * Kovtun (2019) P. Kovtun, JHEP 10, 034 (2019), arXiv:1907.08191 [hep-th] .
  * Hoult and Kovtun (2020) R. E. Hoult and P. Kovtun, JHEP 06, 067 (2020), arXiv:2004.04102 [hep-th] .
  * Martinez and Strickland (2010) M. Martinez and M. Strickland, Nucl. Phys. A 848, 183 (2010), arXiv:1007.0889 [nucl-th] .
  * Martinez  _et al._ (2012) M. Martinez, R. Ryblewski, and M. Strickland, Phys. Rev. C 85, 064913 (2012), arXiv:1204.1473 [nucl-th] .
  * Alqahtani  _et al._ (2017) M. Alqahtani, M. Nopoush, R. Ryblewski, and M. Strickland, Phys. Rev. Lett. 119, 042301 (2017), arXiv:1703.05808 [nucl-th] .
  * Heinz  _et al._ (2016) U. Heinz, D. Bazow, G. S. Denicol, M. Martinez, M. Nopoush, J. Noronha, R. Ryblewski, and M. Strickland, Nucl. Part. Phys. Proc. 276-278, 193 (2016), arXiv:1509.05818 [nucl-th] .
  * Strickland  _et al._ (2018) M. Strickland, J. Noronha, and G. Denicol, Phys. Rev. D 97, 036020 (2018), arXiv:1709.06644 [nucl-th] .
  * Behtash  _et al._ (2018) A. Behtash, C. N. Cruz-Camacho, and M. Martinez, Phys. Rev. D 97, 044041 (2018), arXiv:1711.01745 [hep-th] .
  * Stephanov and Yin (2018) M. Stephanov and Y. Yin, Phys. Rev. D 98, 036006 (2018), arXiv:1712.10305 [nucl-th] .
  * Nahrgang  _et al._ (2019) M. Nahrgang, M. Bluhm, T. Schaefer, and S. A. Bass, Phys. Rev. D 99, 116015 (2019), arXiv:1804.05728 [nucl-th] .
  * An  _et al._ (2020) X. An, G. Başar, M. Stephanov, and H.-U. Yee, Phys. Rev. C 102, 034901 (2020), arXiv:1912.13456 [hep-th] .
  * Pradeep  _et al._ (2022) M. Pradeep, K. Rajagopal, M. Stephanov, and Y. Yin, Phys. Rev. D 106, 036017 (2022), arXiv:2204.00639 [hep-ph] .
  * Greif  _et al._ (2018) M. Greif, J. A. Fotakis, G. S. Denicol, and C. Greiner, Phys. Rev. Lett. 120, 242301 (2018), arXiv:1711.08680 [hep-ph] .
  * Monnai (2012) A. Monnai, Phys. Rev. C 86, 014908 (2012), arXiv:1204.4713 [nucl-th] .
  * Fotakis  _et al._ (2022) J. A. Fotakis, E. Molnár, H. Niemi, C. Greiner, and D. H. Rischke, Phys. Rev. D 106, 036009 (2022), arXiv:2203.11549 [nucl-th] .
  * Bemfica  _et al._ (2021) F. S. Bemfica, M. M. Disconzi, V. Hoang, J. Noronha, and M. Radosz, Phys. Rev. Lett. 126, 222301 (2021), arXiv:2005.11632 [hep-th] .
  * Plumberg  _et al._ (2021) C. Plumberg, D. Almaalol, T. Dore, J. Noronha, and J. Noronha-Hostler,  (2021), 10.1103/PhysRevC.105.L061901, arXiv:2103.15889 [nucl-th] .
  * Chiu and Shen (2021) C. Chiu and C. Shen, Phys. Rev. C 103, 064901 (2021), arXiv:2103.09848 [nucl-th] .
  * Noronha-Hostler  _et al._ (2019) J. Noronha-Hostler, P. Parotto, C. Ratti, and J. M. Stafford, Phys. Rev. C 100, 064910 (2019), arXiv:1902.06723 [hep-ph] .
  * Monnai  _et al._ (2019) A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C 100, 024907 (2019), arXiv:1902.05095 [nucl-th] .
  * Denicol  _et al._ (2012b) G. S. Denicol, E. Molnár, H. Niemi, and D. H. Rischke, Eur. Phys. J. A 48, 170 (2012b), arXiv:1206.1554 [nucl-th] .
  * Molnár  _et al._ (2014) E. Molnár, H. Niemi, G. S. Denicol, and D. H. Rischke, Phys. Rev. D 89, 074010 (2014), arXiv:1308.0785 [nucl-th] .
  * Rose  _et al._ (2020) J.-B. Rose, M. Greif, J. Hammelmann, J. A. Fotakis, G. S. Denicol, H. Elfner, and C. Greiner, Phys. Rev. D 101, 114028 (2020), arXiv:2001.10606 [nucl-th] .
  * Noronha-Hostler  _et al._ (2009) J. Noronha-Hostler, J. Noronha, and C. Greiner, Phys. Rev. Lett. 103, 172302 (2009), arXiv:0811.1571 [nucl-th] .
  * Pal (2010) S. Pal, Phys. Lett. B 684, 211 (2010), arXiv:1001.1585 [nucl-th] .
  * Khvorostukhin  _et al._ (2010) A. S. Khvorostukhin, V. D. Toneev, and D. N. Voskresensky, Nucl. Phys. A 845, 106 (2010), arXiv:1003.3531 [nucl-th] .
  * Tawfik and Wahba (2010) A. Tawfik and M. Wahba, Annalen Phys. 522, 849 (2010), arXiv:1005.3946 [hep-ph] .
  * Alba  _et al._ (2015) P. Alba, R. Bellwied, M. Bluhm, V. Mantovani Sarti, M. Nahrgang, and C. Ratti, Phys. Rev. C 92, 064910 (2015), arXiv:1504.03262 [hep-ph] .
  * Ratti  _et al._ (2011) C. Ratti, S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo (Wuppertal-Budapest), Nucl. Phys. A 855, 253 (2011), arXiv:1012.5215 [hep-lat] .
  * Tiwari  _et al._ (2012) S. K. Tiwari, P. K. Srivastava, and C. P. Singh, Phys. Rev. C 85, 014908 (2012), arXiv:1111.2406 [hep-ph] .
  * Noronha-Hostler  _et al._ (2012) J. Noronha-Hostler, J. Noronha, and C. Greiner, Phys. Rev. C 86, 024913 (2012), arXiv:1206.5138 [nucl-th] .
  * Kadam and Mishra (2014) G. P. Kadam and H. Mishra, Nucl. Phys. A 934, 133 (2014), arXiv:1408.6329 [hep-ph] .
  * Kadam and Mishra (2015) G. P. Kadam and H. Mishra, Phys. Rev. C 92, 035203 (2015), arXiv:1506.04613 [hep-ph] .
  * Kadam and Mishra (2016) G. P. Kadam and H. Mishra, Phys. Rev. C 93, 025205 (2016), arXiv:1509.06998 [hep-ph] .
  * Kadam and Pawar (2019) G. Kadam and S. Pawar, Adv. High Energy Phys. 2019, 6795041 (2019), arXiv:1802.01942 [hep-ph] .
  * Mohapatra  _et al._ (2019) R. K. Mohapatra, H. Mishra, S. Dash, and B. K. Nandi, (2019), arXiv:1901.07238 [hep-ph] .
  * Wesp  _et al._ (2011) C. Wesp, A. El, F. Reining, Z. Xu, I. Bouras, and C. Greiner, Phys. Rev. C 84, 054911 (2011), arXiv:1106.4306 [hep-ph] .
  * Ozvenchuk  _et al._ (2013) V. Ozvenchuk, O. Linnyk, M. I. Gorenstein, E. L. Bratkovskaya, and W. Cassing, Phys. Rev. C 87, 064903 (2013), arXiv:1212.5393 [hep-ph] .
  * Rose  _et al._ (2018) J. B. Rose, J. M. Torres-Rincon, A. Schäfer, D. R. Oliinychenko, and H. Petersen, Phys. Rev. C 97, 055204 (2018), arXiv:1709.03826 [nucl-th] .
  * Rais  _et al._ (2020) J. Rais, K. Gallmeister, and C. Greiner, Phys. Rev. D 102, 036009 (2020), arXiv:1909.04522 [hep-ph] .
  * Kovtun  _et al._ (2005) P. Kovtun, D. T. Son, and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005), arXiv:hep-th/0405231 .
  * Casalderrey-Solana  _et al._ (2014) J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal, and U. A. Wiedemann,  _Gauge/String Duality, Hot QCD and Heavy Ion Collisions_ (Cambridge University Press, 2014) arXiv:1101.0618 [hep-th] .
  * Finazzo and Noronha (2014) S. I. Finazzo and J. Noronha, Phys. Rev. D 89, 106008 (2014), arXiv:1311.6675 [hep-th] .
  * Finazzo  _et al._ (2015) S. I. Finazzo, R. Rougemont, H. Marrochio, and J. Noronha, JHEP 02, 051 (2015), arXiv:1412.2968 [hep-ph] .
  * Denicol  _et al._ (2014) G. S. Denicol, S. Jeon, and C. Gale, Phys. Rev. C 90, 024912 (2014), arXiv:1403.0962 [nucl-th] .
  * Parotto  _et al._ (2020) P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang, J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer, and M. Stephanov, Phys. Rev. C 101, 034901 (2020), arXiv:1805.05249 [hep-ph] .
  * Karthein  _et al._ (2021) J. M. Karthein, D. Mroczek, A. R. Nava Acuna, J. Noronha-Hostler, P. Parotto, D. R. P. Price, and C. Ratti, Eur. Phys. J. Plus 136, 621 (2021), arXiv:2103.08146 [hep-ph] .
  * Hiscock and Lindblom (1983) W. A. Hiscock and L. Lindblom, Annals Phys. 151, 466 (1983).
  * Olson (1990) T. S. Olson, Annals Phys. 199, 18 (1990).
  * Gavassino  _et al._ (2022a) L. Gavassino, M. Antonelli, and B. Haskell, Phys. Rev. Lett. 128, 010606 (2022a), arXiv:2105.14621 [gr-qc] .
  * Gavassino (2021) L. Gavassino, Class. Quant. Grav. 38, 21LT02 (2021), arXiv:2104.09142 [gr-qc] .
  * Gavassino  _et al._ (2022b) L. Gavassino, M. Antonelli, and B. Haskell, Phys. Rev. D 106, 056010 (2022b).
  * Most  _et al._ (2021) E. R. Most, S. P. Harris, C. Plumberg, M. G. Alford, J. Noronha, J. Noronha-Hostler, F. Pretorius, H. Witek, and N. Yunes, Mon. Not. Roy. Astron. Soc. 509, 1096 (2021), arXiv:2107.05094 [astro-ph.HE] .
  * Most  _et al._ (2022) E. R. Most, A. Haber, S. P. Harris, Z. Zhang, M. G. Alford, and J. Noronha, (2022), arXiv:2207.00442 [astro-ph.HE] .
  * Moore and Saremi (2008) G. D. Moore and O. Saremi, JHEP 09, 015 (2008), arXiv:0805.4201 [hep-ph] .
  * Dore  _et al._ (2022b) T. Dore, J. Karthein, D. Mroczek, P. Parotto, J. Noronha-Hostler, and C. Ratti, EPJ Web Conf. 259, 10001 (2022b), arXiv:2109.05098 [nucl-th] .
  * Soloveva  _et al._ (2021) O. Soloveva, D. Fuseau, J. Aichelin, and E. Bratkovskaya, Phys. Rev. C 103, 054901 (2021), arXiv:2011.03505 [nucl-th] .

##  Appendix A shear-bulk coupling

In this appendix, we calculate the divergence of the zero entropy current term

$\displaystyle\partial_{\mu}(\delta_{\Pi\pi}u_{\nu}\Pi\pi^{\mu\nu})$
$\displaystyle=\delta_{\Pi\pi}\Pi\pi^{\mu\nu}(\partial_{\mu}u_{\nu})+u_{\nu}\Pi{\pi^{\mu\nu}}(\partial_{\mu}{\delta}_{\Pi\pi})+\delta_{\Pi\pi}u_{\nu}\pi^{\mu\nu}(\partial_{\mu}{\Pi})+\delta_{\Pi\pi}\Pi
u_{\nu}(\partial_{\mu}{\pi}^{\mu\nu})$
$\displaystyle=\delta_{\Pi\pi}\Pi\pi^{\mu\nu}\partial_{\mu}u_{\nu}+\delta_{\Pi\pi}\Pi
u_{\nu}\partial_{\mu}{\pi}^{\mu\nu}$
$\displaystyle=\delta_{\Pi\pi}\Pi\pi^{\mu\nu}\nabla_{\mu}u_{\nu}-\delta_{\Pi\pi}\Pi\pi^{\mu\nu}\nabla_{\mu}u_{\nu}$
(55)

Where in going from second to third line in Eq. (55) we used the relation

$\displaystyle\partial_{\mu}\underbrace{\left(u_{\nu}\pi^{\mu\nu}\right)}_{\textrm{=0}}$
$\displaystyle=$ $\displaystyle
u_{\nu}\partial_{\mu}\pi^{\mu\nu}+\pi^{\mu\nu}\partial_{\mu}u_{\nu}$

Which trivially gives

$\displaystyle u_{\nu}\partial_{\mu}\pi^{\mu\nu}$ $\displaystyle=$
$\displaystyle-\pi^{\mu\nu}\partial_{\mu}u_{\nu}$

And manipulated the relation

$\displaystyle\partial^{\mu}-\nabla^{\mu}$ $\displaystyle=$
$\displaystyle\partial_{\alpha}g^{\mu\alpha}-\Delta^{\mu\alpha}\partial_{\alpha}$
(56) $\displaystyle=$
$\displaystyle\left[g^{\mu\alpha}-\Delta^{\mu\alpha}\right]\partial_{\alpha}$
$\displaystyle=$ $\displaystyle u^{\mu}u^{\alpha}\partial_{\alpha}$
$\displaystyle=$ $\displaystyle u^{\mu}D$

such that

$\displaystyle\pi^{\mu\nu}\partial_{\mu}u_{\nu}$ $\displaystyle=$
$\displaystyle\pi^{\mu\nu}(u_{\mu}D+\nabla_{\mu})u_{\nu}$ (57)
$\displaystyle=$ $\displaystyle\pi^{\mu\nu}\nabla^{\mu}u^{\nu}.$

Back to Eq. (55), we see that therm $\nabla_{\mu}u_{\nu}$ can be divided into
its symmetric $(\dots)$ and anti-symmetric $\\{\dots\\}$ part

$\displaystyle\nabla_{\mu}u_{\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\nabla_{\mu}u_{\nu}+\nabla_{\nu}u_{\mu}\right)+\frac{1}{2}\left(\nabla_{\mu}u_{\nu}-\nabla_{\nu}u_{\mu}\right)$
(58) $\displaystyle=$
$\displaystyle\sigma_{\mu\nu}+\frac{1}{3}\Delta_{\mu\nu}\theta+\Omega_{\mu\nu}$

where

$\displaystyle\sigma_{\mu\nu}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\left(\nabla_{\mu}u_{\nu}+\nabla_{\nu}u_{\mu}-\frac{2}{3}\Delta_{\mu\nu}\theta\right)$
(59) $\displaystyle=$
$\displaystyle\nabla_{(\mu}u_{\nu)}-\frac{1}{3}\Delta_{\mu\nu}\theta$
$\displaystyle\theta$ $\displaystyle\equiv$ $\displaystyle\nabla_{\mu}u^{\mu}$
(60) $\displaystyle\Omega_{\mu\nu}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\left(\nabla_{\mu}u_{\nu}-\nabla_{\nu}u_{\mu}\right)$
(61)

where the shear stress tensor $\sigma_{\mu\nu}$ is symmetric, the expansion
rate $\Theta$ is symmetric, and the vorticity $\Omega_{\mu\nu}$ is anti-
symmetric. Inserting this symmetric form back into the term
$\delta_{\Pi\pi}\Pi\pi^{\mu\nu}\nabla_{\mu}u_{\nu}$ will get us

$\displaystyle\delta_{\Pi\pi}\Pi\pi^{\mu\nu}\nabla_{\mu}u_{\nu}$
$\displaystyle=\delta_{\Pi\pi}\Pi\pi^{\mu\nu}\big{(}\sigma_{\mu\nu}+\frac{1}{3}\Delta_{\mu\nu}\theta+\Omega_{\mu\nu}\big{)}$
$\displaystyle=\delta_{\Pi\pi}\Pi\pi^{\mu\nu}\sigma_{\mu\nu}$ (62)

where $\pi^{\mu\nu}\Omega_{\mu\nu}=0$ and
$\frac{1}{3}\pi^{\mu\nu}\Delta_{\mu\nu}\theta=0$ by symmetry. This term is the
new addition to the bulk and shear dynamical relaxation equations which arises
from the relevant zero entropy current contribution.

##  Appendix B Lyapunov functional

To construct an energy functional which is quadratic in the perturbation
fields, we start first by expressing the entropy current Eq. (13) as
${S}^{\mu}={S}^{\mu}_{0}+{S}^{\mu}_{1}+{S}^{\mu}_{2}$.

Then, one can express the fluctuations around the equilibrium state as

$\delta{S}^{\mu}=\delta{S}^{\mu}_{0}+\delta{S}^{\mu}_{1}+\delta{S}^{\mu}_{2}.$
(63)

The perturbed fields are as follows

$\displaystyle\tilde{\varepsilon}$
$\displaystyle=\varepsilon+\delta\varepsilon;\qquad\qquad\qquad\qquad\qquad\qquad\tilde{T}=T+\delta
T;\qquad\qquad\qquad\qquad\qquad\qquad\tilde{\Pi}=\Pi+\delta\Pi$
$\displaystyle\tilde{P}$ $\displaystyle=p+\delta
p;\qquad\qquad\qquad\qquad\qquad\qquad\tilde{s}=s+\delta
s;\qquad\qquad\qquad\qquad\qquad\qquad\tilde{\pi}^{\mu\nu}=\pi^{\mu\nu}+\delta\pi^{\mu\nu}$
$\displaystyle\tilde{\rho_{q}}$
$\displaystyle=\rho_{q}+\delta\rho_{q};\qquad\qquad\qquad\qquad\qquad\tilde{\mu_{q}}=\mu_{q}+\delta\mu_{q};\qquad\qquad\qquad\qquad\qquad\qquad\tilde{n}^{\mu}_{q}={n}^{\mu}_{q}+\delta{n}^{\mu}_{q}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\tilde{u}^{\mu}={u}^{\mu}+\delta{u}^{\mu}$
(64)

for perturbations at $0^{th}$ order of the equilibrium entropy current

$S^{\mu}_{0}=su^{\mu},$ (65)

with

$sT=\varepsilon+p-\sum_{q}\mu_{q}\,\rho_{q},\qquad\qquad q\in\\{B,S,Q\\}\,.$
(66)

Therefore, the perturbed expression for this entropy becomes

$\displaystyle S^{\mu}_{0}+\delta S^{\mu}_{0}$ $\displaystyle=(s+\delta
s)(u^{\mu}+\delta u^{\mu})$ $\displaystyle S^{\mu}_{0}+\delta S^{\mu}_{0}$
$\displaystyle=su^{\mu}+s\delta u^{\mu}+\delta su^{\mu}+\delta s\delta
u^{\mu}$ $\displaystyle\delta S^{\mu}_{0}$ $\displaystyle=s\delta
u^{\mu}+\delta su^{\mu}+\delta s\delta u^{\mu}.$ (67)

while the perturbed form of the thermodynamic relation 66 gives

$\displaystyle\delta s\,T$
$\displaystyle=(\varepsilon+p-\sum_{q}\mu_{q}\,\rho_{q}-sT)+\delta\varepsilon+\delta
p-\sum_{q}(\mu_{q}\,\delta\rho_{q}+\delta\mu_{q}\,\rho_{q}+\delta\mu_{q}\,\delta\rho_{q})-s\delta
T-\delta s\delta T\,,$ $\displaystyle=\delta\varepsilon+\delta
p-\mu_{q}\,\delta\rho_{q}-\delta\mu_{q}\,\rho_{q}-\delta\mu_{q}\,\delta\rho_{q}-s\delta
T-\delta s\delta T\,.$ (68)

The first term in parentheses go to zero by the identity Eq. (66). Inserting
the expression into Eq. (67) and keeping only up to second order terms gives

$\displaystyle T\delta S^{\mu}_{0}$
$\displaystyle=\left[\varepsilon+p-\mu_{q}\,\rho\right]\delta
u^{\mu}+\left[\delta\varepsilon-\mu_{q}\,\delta\rho+\delta
p-\delta\mu_{q}\,\rho-s\delta T-\delta\mu_{q}\,\delta\rho+\delta s\delta
T\right]u^{\mu}$
$\displaystyle+\left[\delta\varepsilon-\mu_{q}\,\delta\rho+\delta
p-\delta\mu_{q}\,\rho-s\delta T\right]\delta u^{\mu}\,.$ (69)

Finally eliminating the third term in the third square brackets using the
thermodynamics identity (Gibbs-Duhem) gives

$\displaystyle T\delta S^{\mu}_{0}$
$\displaystyle=\left[\varepsilon+p-\mu_{q}\,\rho\right]\delta
u^{\mu}+\left[\delta\varepsilon-\mu_{q}\,\delta\rho-\delta\mu_{q}\,\delta\rho+\delta
s\delta
T\right]u^{\mu}+\left[\delta\varepsilon-\mu_{q}\,\delta\rho\right]\delta
u^{\mu}\,.$ (70)

Now, for the 1st order perturbations

$\displaystyle{S}^{\mu}_{1}$
$\displaystyle=\alpha_{q}n^{\mu}_{q}\,,\qquad\qquad q\in{B,S,Q}$
$\displaystyle\delta{S}^{\mu}_{1}$
$\displaystyle=\alpha_{q}n^{\mu}_{q}+\alpha_{q}\delta
n^{\mu}_{q}+\delta\alpha_{q}\delta n^{\mu}_{q}+\delta\alpha_{q}\delta
n^{\mu}_{q}\,,$ (71)

Finally, for the 2nd order perturbations

${S}^{\mu}_{2}=-\frac{1}{2}u^{\mu}\left(\beta_{\Pi}{\Pi}^{2}+\beta_{\Pi}{\pi}^{\mu\nu}{\pi}_{\mu\nu}+\sum_{q}\beta_{n}^{qq^{\prime}}n_{q}^{\mu}{n}_{\mu}^{q^{\prime}}\right)-\sum_{q}\left(\delta_{n\pi}^{q}n^{\mu}_{q}{\Pi}+\delta_{n\pi}^{q}n_{\nu}^{q}{\pi}^{\mu\nu}\right)\,,$
(72)

which then can be expanded as

$\displaystyle\delta{S}^{\mu}_{2}$
$\displaystyle=-\frac{1}{2}\left(u^{\mu}+\delta
u^{\mu}\right)\left[\beta_{\Pi}\left(\Pi+\delta\Pi\right)^{2}+\beta_{\Pi}\left({\pi}^{\mu\nu}+\delta{\pi}^{\mu\nu})({\pi}_{\mu\nu}+\delta{\pi}_{\mu\nu}\right)+\sum_{q}\beta_{n}^{qq^{\prime}}(n_{q}^{\mu}+\delta
n_{q}^{\mu})({n}_{\mu}^{q^{\prime}}+\delta{n}_{\mu}^{q^{\prime}})\right]$
$\displaystyle-\sum_{q}\left[\delta_{n\Pi}^{q}(n^{\mu}_{q}+\delta
n_{q}^{\mu})(\Pi+\delta\Pi)+\delta_{n\pi}^{q}(n_{\nu}^{q}+\delta
n_{\nu}^{q})({\pi}^{\mu\nu}+\delta{\pi}^{\mu\nu})\right]\,.$ (73)

At this point, we set all terms proportional to the currents
$\\{\Pi,\pi^{\mu\nu},n^{\mu}_{q};q\in\\{B,S,Q\\}\\}$ to zero by definition of
the equilibrium state and keep only terms at second order in the perturbations
to arrive at

$\displaystyle\delta{S}^{\mu}_{2}$
$\displaystyle=-\frac{1}{2}u^{\mu}\left(\beta_{\Pi}\delta\Pi\delta\Pi+\beta_{\Pi}\delta{\pi}^{\mu\nu}\delta{\pi}_{\mu\nu}+\sum_{q}\beta_{n}^{qq^{\prime}}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}\right)-\sum_{q}\left(\delta_{n\Pi}^{q}\delta
n^{\mu}_{q}\delta\Pi+\delta_{n\pi}^{q}\delta
n_{\nu}^{q}\delta\pi^{\mu\nu}\right).$ (74)

To get the full perturbed non-equilibrium entropy, we add up all terms in Eqs.
(70, 71, and 74) to get

$\displaystyle T\delta S^{\mu}_{non-eq}$
$\displaystyle=\left[\varepsilon+p-\sum_{q}\mu_{q}\,\rho\right]\delta
u^{\mu}+\left[\delta\varepsilon-\sum_{q}(\mu_{q}\,\delta\rho-\delta\mu_{q}\,\delta\rho)+\delta
s\delta
T\right]u^{\mu}+\left[\delta\varepsilon-\sum_{q}\mu_{q}\,\delta\rho\right]\delta
u^{\mu}$ $\displaystyle-\alpha_{q}\delta
n_{q}^{\mu}-\frac{1}{2}u^{\mu}\left(\beta_{\Pi}\delta\Pi\delta\Pi+\beta_{\Pi}\delta{\pi}^{\mu\nu}\delta{\pi}_{\mu\nu}+\sum_{q}\beta_{n}^{qq^{\prime}}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}\right)-\sum_{q}\left(\delta_{n\Pi}^{q}\delta
n^{\mu}\delta\Pi+\delta_{n\pi}^{q}\delta n_{\nu}^{q}\delta\pi^{\mu\nu}\right)$
(75)

The idea is that the energy-momentum and number should be conserved for all
physical state configurations. These ”zero flux contributions” are calculated
using the next two expressions

####  B.0.1 zero energy-momentum flux contribution

We start from the perturbed form of the Landau condition,
$u_{\mu}T^{\mu\nu}=\varepsilon u^{\nu}$, imposing that no fluctuations can
propagate

$\displaystyle(u_{\nu}+\delta u_{\nu})(T^{\mu\nu}+\delta T^{\mu\nu})$
$\displaystyle=(\varepsilon+\delta\varepsilon)(u^{\mu}+\delta u^{\mu})\,,$
$\displaystyle u_{\nu}T^{\mu\nu}+u_{\nu}\delta T^{\mu\nu}+\delta
u_{\nu}T^{\mu\nu}+\delta u_{\nu}\delta T^{\mu\nu}$ $\displaystyle=\varepsilon
u^{\mu}+\varepsilon\delta u^{\mu}+\delta\varepsilon
u^{\mu}+\delta\varepsilon\delta u^{\mu}\,,$ $\displaystyle u_{\nu}\delta
T^{\mu\nu}$ $\displaystyle=\varepsilon\delta u^{\mu}+\delta\varepsilon
u^{\mu}+\delta\varepsilon\delta u^{\mu}-\delta u_{\nu}T^{\mu\nu}-\delta
u_{\nu}\delta T^{\mu\nu}\,,$ $\displaystyle=(\rm zfc)_{energy-momentum}\,.$
(76)

Now we turn to the charge current(s) zero flux contribution. We start again
from the definition of Landau frame of the currents expanded in perturbed form
$N^{\mu}_{q}=\rho u^{\mu}+n^{\mu}_{q}$. Note that this will impose conditions
on each of the charges separately

$\displaystyle N^{\mu}_{q}+\delta N^{\mu}_{q}$
$\displaystyle=(\rho_{q}+\delta\rho_{q})(u^{\mu}+\delta
u^{\mu})+(n^{\mu}_{q}+\delta n^{\mu}_{q})\,,$
$\displaystyle=\rho_{q}u_{\mu}+\rho_{q}\delta u_{\mu}+\delta\rho
u_{\mu}+\delta\rho\delta u_{\mu}+n^{\mu}_{q}+\delta n^{\mu}_{q}\,,$
$\displaystyle\delta N^{\mu}_{q}$ $\displaystyle=\rho_{q}\,\delta
u_{\mu}+\delta\rho_{q}\,u_{\mu}+\delta\rho_{q}\,\delta u_{\mu}+\delta
n^{\mu}_{q}\,.$ $\displaystyle=(\rm zfc)_{charges}$ (77)

From energy flux contributions equation, we get

$\displaystyle{(\rm zfc)}_{T^{\mu\nu}}$ $\displaystyle=\varepsilon\delta
u^{\mu}+\delta\varepsilon u^{\mu}+\delta\varepsilon\delta u^{\mu}-\delta
u_{\nu}T^{\mu\nu}-\delta u_{\nu}\delta T^{\mu\nu}\,,$ (78)

Similarly, from the number flux contributions equation one finds

$\displaystyle{(\rm zfc)}_{N_{q}^{\mu}}$ $\displaystyle=\rho_{q}\,\delta
u_{\mu}+\delta\rho_{q}\,u_{\mu}+\delta\rho_{q}\,\delta u_{\mu}+\delta
n^{\mu}_{q}\,,$

and we obtain

$\displaystyle\mu_{q}\,\,\rho_{q}\delta u_{\mu}$ $\displaystyle=({\rm
zfc})-\mu_{q}\,\,\delta\rho_{q}u_{\mu}-\mu_{q}\,\,\delta\rho_{q}\delta
u_{\mu}-\mu_{q}\,\delta n^{\mu}_{q}\,.$

The final expression is now

$\displaystyle T\delta S^{\mu}_{non-eq}$ $\displaystyle=({\rm
zfc})-\delta\varepsilon u^{\mu}-\delta\varepsilon\delta u^{\mu}+\delta
u_{\nu}T^{\mu\nu}+\delta u_{\nu}\delta T^{\mu\nu}+p\delta u^{\mu}-({\rm
zfc})+{\mu_{q}\,\delta\rho}\,u_{\mu}$
$\displaystyle+{\mu_{q}\,\delta\rho}\,\delta u_{\mu}+{\mu_{q}\,\delta
n^{\mu}_{q}}+\left[\delta\varepsilon-\mu_{q}\,\delta\rho_{q}-\delta\mu_{q}\,\delta\rho+\delta
s\delta
T\right]u^{\mu}+\left[\delta\varepsilon-\mu_{q}\,\delta\rho\right]\delta
u^{\mu}$ $\displaystyle-\mu_{q}\delta n_{q}^{\mu}-\delta\mu_{q}\delta
n_{q}^{\mu}+\frac{\mu_{q}\delta n_{q}^{\mu}\delta
T}{T}-\frac{1}{2}u^{\mu}\left(\beta_{\Pi}\delta\Pi\delta\Pi+\beta_{\Pi}\delta{\pi}^{\mu\nu}\delta{\pi}_{\mu\nu}+\sum_{q}\beta_{n}^{qq^{\prime}}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}\right)$
$\displaystyle-\sum_{q}\left(\gamma_{n\Pi}^{q}\delta
n^{\mu}\delta\Pi+\gamma_{n\pi}^{q}\delta n_{\nu}^{q}\delta\pi^{\mu\nu}\right)$
(80)

To further simplify the above expression, we note first that the energy
momentum tensor (in the mostly negative convention) is decomposed as

$\displaystyle T^{\mu\nu}=\varepsilon
u^{\mu}u^{\nu}-(p+\pi)\delta^{\mu\nu}+\pi^{\mu\nu}$ (81)

Then by projecting into $\delta u^{\mu}$ and manipulating the relations
$u^{\nu}\delta u_{\nu}=\delta u^{\nu}\delta u_{\nu}/2$, one finds

$\displaystyle\delta u_{\nu}T^{\mu\nu}$ $\displaystyle=\varepsilon\delta
u_{\nu}u^{\mu}u^{\nu}-(p+\Pi)\delta u_{\nu}\delta^{\mu\nu}+\delta
u_{\mu}\pi^{\mu\nu}$ $\displaystyle\delta u_{\nu}T^{\mu\nu}$
$\displaystyle=\varepsilon\delta u_{\nu}u^{\mu}u^{\nu}-p\delta
u_{\nu}\delta^{\mu\nu}-\Pi\delta u_{\nu}\delta^{\mu\nu}+\delta
u_{\nu}\pi^{\mu\nu}$ $\displaystyle=\varepsilon\delta
u_{\nu}u^{\mu}u^{\nu}-p\delta u_{\nu}g^{\mu\nu}+u_{\mu}u_{\nu}p\delta
u_{\nu}-\Pi\delta u^{\mu}+\delta u_{\nu}\pi^{\mu\nu}$
$\displaystyle=-(\varepsilon+p)\frac{\delta u_{\nu}\delta
u^{\nu}}{2}u^{\mu}-p\delta u^{\mu}-\Pi\delta u^{\mu}+\delta
u_{\nu}\pi^{\mu\nu}$ (82)

Likewise

$\displaystyle T^{\mu\nu}+\delta T^{\mu\nu}$
$\displaystyle=(\varepsilon+\delta\varepsilon)u^{\mu}u^{\nu}-(p+\delta
p)\delta^{\mu\nu}-(\Pi+\delta\pi)\delta^{\mu\nu}+(\pi^{\mu\nu}+\delta\pi^{\mu\nu})$
$\displaystyle=\varepsilon u^{\mu}u^{\nu}+\delta\varepsilon
u^{\mu}u^{\nu}-p\delta^{\mu\nu}-\delta
p\delta^{\mu\nu}-\Pi\delta^{\mu\nu}-\delta\Pi\delta^{\mu\nu}+\pi^{\mu\nu}+\delta\pi^{\mu\nu}$

Expanding out all terms

$\displaystyle T^{\mu\nu}+\delta T^{\mu\nu}$
$\displaystyle=\varepsilon(u^{\mu}+\delta u^{\mu})(u^{\nu}+\delta
u^{\nu})+\delta\varepsilon(u^{\mu}+\delta u^{\mu})(u^{\nu}+\delta
u^{\nu})-p(g^{\mu\nu}-(u^{\mu}+\delta u^{\mu})(u^{\nu}+\delta u^{\nu}))$
$\displaystyle-\Pi(g^{\mu\nu}-(u^{\mu}+\delta u^{\mu})(u^{\nu}+\delta
u^{\nu}))+\pi^{\mu\nu}+\delta\pi^{\mu\nu}$
$\displaystyle=\varepsilon(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}+\delta
u^{\mu}u^{\nu}+\delta u^{\mu}\delta
u^{\nu})+\delta\varepsilon(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}+\delta
u^{\mu}u^{\nu}+\delta u^{\mu}\delta u^{\nu})$ $\displaystyle-
pg^{\mu\nu}-p(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}+\delta
u^{\mu}u^{\nu}+\delta u^{\mu}\delta u^{\nu})-\delta pg^{\mu\nu}-\delta
p(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}$ $\displaystyle+\delta
u^{\mu}u^{\nu}+\delta u^{\mu}\delta u^{\nu})-\Pi
g^{\mu\nu}-\Pi(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}+\delta
u^{\mu}u^{\nu}+\delta u^{\mu}\delta u^{\nu})-\delta\Pi g^{\mu\nu}$
$\displaystyle-\delta\Pi(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}+\delta
u^{\mu}u^{\nu}+\delta u^{\mu}\delta u^{\nu})+\pi^{\mu\nu}+\delta\pi^{\mu\nu}.$
(84)

Canceling the unperturbed terms on both sides and discarding terms at third
order in fluctuations gives

$\displaystyle\delta T^{\mu\nu}$ $\displaystyle=\varepsilon(u^{\mu}\delta
u^{\nu}+\delta u^{\mu}u^{\nu}+\delta u^{\mu}\delta
u^{\nu})+\delta\varepsilon(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}+\delta
u^{\mu}u^{\nu})$ $\displaystyle-p(u^{\mu}\delta u^{\nu}+\delta
u^{\mu}u^{\nu}+\delta u^{\mu}\delta u^{\nu})-\delta pg^{\mu\nu}-\delta
p(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}$ $\displaystyle+\delta
u^{\mu}u^{\nu})-\Pi(u^{\mu}\delta u^{\nu}+\delta u^{\mu}u^{\nu}+\delta
u^{\mu}\delta u^{\nu})-\delta\Pi g^{\mu\nu}$
$\displaystyle-\delta\Pi(u^{\mu}u^{\nu}+u^{\mu}\delta u^{\nu}+\delta
u^{\mu}u^{\nu})+\delta\pi^{\mu\nu}$ (85)

Finally, we project with $\delta u_{\nu}$ keeping in mind that it preserves
all symmetry arguments as the flow field $u_{\nu}$ and cancelling all third
order terms

$\displaystyle\delta u_{\nu}\delta T^{\mu\nu}$
$\displaystyle=\varepsilon(u^{\mu}\delta u^{\nu}\delta u_{\nu}+\delta
u^{\mu}\delta u_{\nu}u^{\nu})+\delta\varepsilon(u^{\mu}u^{\nu}\delta
u_{\nu})+p(u^{\mu}\delta u^{\nu}\delta u_{\nu}+\delta u^{\mu}u^{\nu}\delta
u_{\nu})$ $\displaystyle-\delta pg^{\mu\nu}\delta u_{\nu}+\delta
p(u^{\mu}u^{\nu}\delta u_{\nu})-\Pi(u^{\mu}\delta u^{\nu}\delta u_{\nu}+\delta
u^{\mu}u^{\nu}\delta u_{\nu})-\delta\Pi\delta u^{\mu}$
$\displaystyle-\delta\Pi(u^{\mu}u^{\nu}\delta u_{\nu}+u^{\mu}\delta
u^{\nu}\delta u_{\nu}+\delta u^{\mu}u^{\nu}\delta u_{\nu})+\delta
u_{\nu}\,\delta\pi^{\mu\nu}$ $\displaystyle=\varepsilon u^{\mu}\delta
u^{\nu}\delta u_{\nu}+pu^{\mu}\delta u^{\nu}\delta u_{\nu}-\delta p\delta
u^{\mu}-\Pi u^{\mu}\delta u^{\nu}\delta u_{\nu}-\delta\Pi\delta u^{\mu}+\delta
u_{\nu}\,\delta\pi^{\mu\nu}$ (86)

to arrive at

$\displaystyle\delta u_{\nu}T^{\mu\nu}+\delta u_{\nu}\delta T^{\mu\nu}$
$\displaystyle=-(\varepsilon+p)\frac{\delta u_{\nu}\delta
u^{\nu}}{2}u^{\mu}-p\delta u^{\mu}+\varepsilon u^{\mu}\delta u^{\nu}\delta
u_{\nu}+pu^{\mu}\delta u^{\nu}\delta u_{\nu}-\delta p\delta u^{\mu}$
$\displaystyle-\Pi\delta u^{\mu}+\delta u_{\nu}\pi^{\mu\nu}-\Pi u^{\mu}\delta
u^{\nu}\delta u_{\nu}-\delta\Pi\delta u^{\mu}+\delta
u_{\nu}\,\delta\pi^{\mu\nu}$ $\displaystyle=(\varepsilon+p-\Pi)\frac{\delta
u_{\nu}\delta u^{\nu}}{2}u^{\mu}-(p+\delta p+\Pi)\delta
u^{\mu}+(\pi^{\mu\nu}+\delta\pi^{\mu\nu})\delta u_{\nu}$ (87)

The final step in the derivation is to subtract those terms which will provide
us with the Lyabunov functional

$\displaystyle E^{\mu}$ $\displaystyle=({\rm zfc})-\delta S^{\mu}_{non-eq}$
(88)

Before we go further we note that some of our fluctuations are not independent
from the others. So first we rewrite all quantities in linear combinations of
independent ones to make sure our basis is complete and independent. First we
need to simplify

$\displaystyle\delta\mu_{q}\,\delta\rho_{q}+\delta s\delta T$ (89)
$\displaystyle ds=\frac{1}{T}d\varepsilon-\frac{\mu_{q}}{T}d\rho_{q}$ (90)

giving

$\displaystyle ds=\frac{1}{T}d\varepsilon-\frac{\mu_{q}}{T}d\rho_{q}$ (91)
$\displaystyle\delta\mu_{q}\,\delta\rho_{q}+\delta s\delta
T=\frac{1}{\varepsilon+p}\left[\left.\frac{\partial\varepsilon}{\partial
p}\right|_{s}(\delta p)^{2}+\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial s}\right|_{\mu_{q}/T}(\delta
s)^{2}\right]$ (92)
$\displaystyle\delta\frac{\mu_{q}}{T}=\left.\frac{\partial}{\partial
p}\frac{\mu_{q}}{T}\right|_{s}\delta p+\left.\frac{\partial}{\partial
s}\frac{\mu_{q}}{T}\right|_{p}\delta s$ (93)

##  Appendix C Stability constraints

The energy functional in Eq. (40) is written as function of fluctuations
including couplings which introduce nonlinear dependence between different
fluctuations. In order to prove the positivity of the energy functional, we
first write it in linearly independent combinations. Then, we assure these
combinations are positive definite by completing the squares of all coupled
terms.

I will start first by hydrodynamics fluctuations terms which are coupled to
thermodynamics

$\displaystyle e$ $\displaystyle=-(\varepsilon+p)\frac{\delta u_{\nu}\delta
u^{\nu}}{2}+{\delta p\,\delta u^{\mu}}\lambda_{\mu}\pm\delta\Pi\,\delta
u^{\mu}\lambda_{\mu}+\delta\pi^{\mu\nu}\delta
u_{\nu}\lambda_{\mu}+{\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
p}\right|_{s}(\delta p)^{2}}$
$\displaystyle+{\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial s}\right|_{\mu_{q}/T}(\delta
s)^{2}}+{\sum_{q}\left.\frac{\partial\alpha_{q}}{\partial
p}\right|_{s}\lambda_{\mu}\delta n_{q}^{\mu}\delta
p}+{\sum_{q}\left.\frac{\partial\alpha_{q}}{\partial
s}\right|_{p}\lambda_{\mu}\delta n_{q}^{\mu}\delta
s}+\frac{\beta_{\Pi}}{2}\delta\Pi\delta\Pi$
$\displaystyle+\frac{\beta_{\pi}}{2}\delta{\pi}^{\mu\nu}\delta{\pi}_{\mu\nu}+\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}+\sum_{q}\delta_{n\pi}^{q}\lambda_{\mu}\delta
n^{\mu}\delta\Pi+\sum_{q}\delta_{n\pi}^{q}\lambda_{\mu}\delta
n_{\nu}^{q}\delta\pi^{\mu\nu}$ (94)

Since we cannot constrain the magnitude of perturbations or how they
propagate, we need to rewrite them in a quadratic form to be able to extract
the conditions. Starting from the coupling between $(\delta u^{\mu},\delta
n_{q}^{\mu},\delta p)$ perturbations, one gets

$\displaystyle e_{u^{\mu},n_{q}^{\mu},p}$ $\displaystyle=\delta p\,\delta
u^{\mu}\lambda_{\mu}+\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
p}\right|_{s}(\delta p)^{2}+\sum_{q}\left.\frac{\partial\alpha_{q}}{\partial
p}\right||_{s}\lambda_{\mu}\delta n_{q}^{\mu}\delta p$
$\displaystyle=\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
p}\right|_{s}\left[\delta p+\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left(\delta
u^{\mu}\lambda_{\mu}+\left.\frac{\partial\alpha_{q}}{\partial
p}\right|_{s}\delta
n_{q}^{\mu}\lambda_{\mu}\right)\right]^{2}-\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}(\delta u^{\mu}\lambda_{\mu})^{2}$
$\displaystyle-\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}(\delta n_{q}^{\mu}\lambda_{\mu})(\delta
n_{q^{\prime}}^{\mu}\lambda_{\mu})-\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}(\delta u^{\mu}\lambda_{\mu})(\delta n_{q}^{\mu}\lambda_{\mu}).$
(95)

We follow the same procedure for the couplings between $({\delta
u^{\mu},\delta n_{q}^{\mu},\delta s})$ perturbations to get

$\displaystyle e_{n_{q}^{\mu},s}$
$\displaystyle=T\left.\frac{\partial{\alpha_{q}}}{\partial s}\right|_{p}\delta
s\,\delta
n_{q}^{\mu}\lambda_{\mu}+\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial s}\right|_{\alpha_{q}}(\delta
s)^{2}$
$\displaystyle=\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial
s}\right|_{\alpha_{q}}\left[\delta s+T(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}(\delta
n_{q}^{\mu}\lambda_{\mu})\right]^{2}-T^{2}(\varepsilon+p)$
$\displaystyle\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}(\delta n_{q}^{\mu}\lambda_{\mu})(\delta
n_{q^{\prime}}^{\mu}\lambda_{\mu})$ (96)

At this point we get a few linear independent combinations of the
thermodynamics quantities with the flow and charge currents fields. We
continue with coupled fluctuations among the pure hydrodynamics fields
starting with the bulk fluctuations $({\delta\Pi,\delta u^{\mu},\delta
n^{\mu}})$

$\displaystyle e_{\Pi,u^{\mu},n^{\mu}}$
$\displaystyle=\frac{\beta_{\Pi}}{2}\delta\Pi\delta\Pi+\delta\Pi\,\delta
u^{\mu}\lambda_{\mu}+\sum_{q}\delta_{n\Pi}^{q}\lambda_{\mu}\delta
n^{\mu}\delta\Pi$
$\displaystyle=\frac{\beta_{\Pi}}{2}\left[\delta\Pi+\frac{1}{\beta_{\Pi}}\delta
u^{\mu}\lambda_{\mu}\pm\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}\delta
n^{\mu}_{q}\lambda_{\mu}\right]^{2}-\frac{1}{\beta_{\Pi}}(\delta
u^{\mu}\lambda_{\mu})^{2}-\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}(\delta
n^{\mu}_{q}\lambda_{\mu})^{2}-\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}(\delta
n^{\mu}_{q}\lambda_{\mu})(\delta u^{\mu}\lambda_{\mu})$ (97)

as well as the couplings with shear stress fluctuations

$\displaystyle e_{\pi^{\mu\nu},u^{\mu},n^{\mu}}$
$\displaystyle=\frac{\beta_{\pi}}{2}\delta\pi_{\mu\nu}\delta\pi^{\mu\nu}\pm\delta\pi^{\mu\nu}\,\delta
u_{\nu}\lambda_{\mu}+\sum_{q}\delta_{n\pi}^{q}\lambda_{\mu}\delta
n_{\nu}\delta\pi^{\mu\nu}\,,$
$\displaystyle=\frac{\beta_{\pi}}{2}\left[\delta\pi^{\mu\nu}+\frac{1}{\beta_{\pi}}\lambda^{\mu}\delta
u^{\nu}-\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{\mu}\delta
n^{\nu}_{q}\right]^{2}-\frac{1}{\beta_{\pi}}\lambda^{2}\delta u^{\nu}\delta
u_{\nu}\pm\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}\lambda^{2}\delta
n^{\nu}_{q}\delta
n_{\nu}^{q}-\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{2}\delta
n_{\nu}^{q}\delta u^{\nu}\,.$ (98)

Inserting Eqs. (95-98) back into the general form Eq (44) then re-writing the
functional by grouping relevant terms give us

$\displaystyle e$
$\displaystyle={\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
p}\right|_{s}\left[\delta p+\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left(\delta
u^{\mu}\lambda_{\mu}+\left.\frac{\partial\alpha_{q}}{\partial
p}\right|_{s}\delta n_{q}^{\mu}\lambda_{\mu}\right)\right]^{2}}$
$\displaystyle+{\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial
s}\right|_{\alpha_{q}}\left[\delta s+T(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\frac{\partial{\theta_{q}}}{\partial s}|_{p}(\delta
n_{q}^{\mu}\lambda_{\mu})\right]^{2}}$
$\displaystyle-{\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}(\delta
u^{\mu}\lambda_{\mu})^{2}}-{\frac{1}{\beta_{\Pi}}(\delta
u^{\mu}\lambda_{\mu})^{2}}$
$\displaystyle-{\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}(\delta
n_{q}^{\mu}\lambda_{\mu})^{2}}-{T^{2}(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}(\delta
n_{q}^{\mu}\lambda_{\mu})^{2}}\pm{\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}(\delta
n^{\mu}_{q}\lambda_{\mu})^{2}}$
$\displaystyle-{\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}(\delta u^{\mu}\lambda_{\mu})(\delta
n_{q}^{\mu}\lambda_{\mu})}-{\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}(\delta
n^{\mu}_{q}\lambda_{\mu})(\delta u^{\mu}\lambda_{\mu})}$
$\displaystyle+{\frac{\beta_{\pi}}{2}\left[\delta\pi^{\mu\nu}+\frac{1}{\beta_{\pi}}\lambda^{\mu}\delta
u^{\nu}-\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{\mu}\delta
n^{\nu}_{q}\right]^{2}}-{\frac{1}{\beta_{\pi}}\lambda^{2}\delta u^{\nu}\delta
u_{\nu}}-{\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}\lambda^{2}\delta
n^{\nu}_{q}\delta
n_{\nu}^{q}}-{\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{2}\delta
n_{\nu}^{q}\delta u^{\nu}}$
$\displaystyle+{\frac{\beta_{\Pi}}{2}\left[\delta\Pi+\frac{1}{\beta_{\Pi}}\delta
u^{\mu}\lambda_{\mu}-\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}\delta
n^{\mu}_{q}\lambda_{\mu}\right]^{2}}$
$\displaystyle-(\varepsilon+p)\frac{\delta u_{\nu}\delta
u^{\nu}}{2}+\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}$ (99)

As we can see from Eq. (99), in order to proceed further we need to project
out time like and space like components of the flow velocity and charge
currents. To do so, I will now manipulate the following relation

$\displaystyle\gamma^{\mu}_{\nu}$
$\displaystyle=g^{\mu}_{\nu}+u^{\mu}u_{\nu}-\frac{\lambda^{\mu}\lambda_{\nu}}{\lambda^{2}}\,,$
$\displaystyle=\Delta^{\mu}_{\nu}-\frac{\lambda^{\mu}\lambda_{\nu}}{\lambda^{2}}\,,$
$\displaystyle\gamma_{\nu}^{\mu}\delta u_{\mu}\delta u^{\nu}$
$\displaystyle=\Delta^{\mu}_{\nu}-\frac{\lambda^{\mu}\lambda_{\nu}}{\lambda^{2}}(\delta
u_{\mu}\delta u^{\nu})$ $\displaystyle=g^{\mu}_{\nu}\delta u_{\mu}\delta
u^{\nu}+u^{\mu}u_{\nu}\delta u_{\mu}\delta
u^{\nu}-\frac{\lambda^{\mu}\lambda_{\nu}}{\lambda^{2}}(\delta u_{\mu}\delta
u^{\nu})$ $\displaystyle=\delta u_{\nu}\delta u^{\nu}+(u_{\alpha}\delta
u^{\alpha})^{2}-\frac{\lambda^{\mu}\delta u_{\mu}}{\lambda^{2}}$ (100)

Using the relation $u^{\mu}\delta u_{\mu}=-\delta u_{\mu}\delta u^{\mu}/2$ we
cancel the middle term to get

$\displaystyle-\frac{\varepsilon+p}{2}\delta u_{\nu}\delta u^{\nu}$
$\displaystyle=-\frac{\varepsilon+p}{2}\gamma_{\nu}^{\mu}\delta u_{\mu}\delta
u^{\nu}-\frac{\varepsilon+p}{2}\frac{(\lambda_{\mu}\delta
u^{\mu})^{2}}{\lambda^{2}}$ (101)

Likewise

$\displaystyle\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}\delta
n_{q}^{\mu}\delta{n}_{\mu}^{q^{\prime}}$
$\displaystyle=\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}\gamma_{\nu}^{\mu}\delta
n^{q}_{\mu}\delta{n}^{\nu}_{q^{\prime}}+\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}\frac{(\lambda_{\mu}\delta
n_{q}^{\mu})^{2}}{\lambda^{2}}$

And finally, once more, I manipulate the relation

$\gamma_{\mu\nu}=g_{\mu\nu}+u_{\mu}u_{\nu}-\frac{\lambda_{\mu}\lambda_{\nu}}{\lambda^{2}}\,,$
(103)

to get

$\displaystyle\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}\lambda^{2}\delta
n^{\nu}_{q}\delta n_{\nu}^{q}$
$\displaystyle=-\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}\lambda^{2}\gamma_{\nu}^{\mu}\delta
n_{\mu}^{q}\delta
n^{\nu}_{q^{\prime}}-\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}(\lambda_{\mu}\delta
n^{\mu}_{q})^{2}$ (104) $\displaystyle-\frac{\lambda^{2}}{\beta_{\pi}}\delta
u_{\nu}\delta u^{\nu}$
$\displaystyle=-\frac{\lambda^{2}}{\beta_{\pi}}\gamma_{\nu}^{\mu}\delta
u_{\mu}\delta u^{\nu}-\frac{1}{\beta_{\pi}}(\lambda_{\mu}\delta u^{\mu})^{2}$
(105) $\displaystyle\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{2}\delta
n^{\nu}_{q}\delta u_{\nu}$
$\displaystyle=-\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{2}\gamma_{\nu}^{\mu}\delta
n^{\nu}_{q}\delta
u_{\nu}-\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}(\lambda_{\mu}\delta
n^{\mu}_{q})(\lambda_{\mu}\delta u^{\mu})$ (106)

Inserting the above expressions into the original form

$\displaystyle e$
$\displaystyle=\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
p}\right|_{s}\left[\delta p+\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left(\delta
u^{\mu}\lambda_{\mu}+\left.\frac{\partial\alpha_{q}}{\partial
p}\right|_{s}\delta n_{q}^{\mu}\lambda_{\mu}\right)\right]^{2}$
$\displaystyle+\frac{1}{\varepsilon+p}\left.\frac{\partial\varepsilon}{\partial
s}\right|_{p}\left.\frac{\partial p}{\partial
s}\right|_{\alpha_{q}}\left[\delta s+T(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}(\delta n_{q}^{\mu}\lambda_{\mu})\right]^{2}$
$\displaystyle-\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}(\delta
u^{\mu}\lambda_{\mu})^{2}-{\frac{1}{\beta_{\Pi}}(\delta
u^{\mu}\lambda_{\mu})^{2}}-\frac{\varepsilon+p}{2}\frac{(\lambda_{\mu}\delta
u^{\mu})^{2}}{\lambda^{2}}-{\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}(\delta n_{q}^{\mu}\lambda_{\mu})^{2}}$
$\displaystyle-{T^{2}(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}(\delta
n_{q}^{\mu}\lambda_{\mu})^{2}}\pm{\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}(\delta
n^{\mu}_{q}\lambda_{\mu})^{2}}-\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}(\lambda_{\mu}\delta
n^{\mu}_{q})^{2}+\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}\frac{(\lambda_{\mu}\delta
n_{q}^{\mu})^{2}}{\lambda^{2}}$
$\displaystyle-{\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}(\delta u^{\mu}\lambda_{\mu})(\delta
n_{q}^{\mu}\lambda_{\mu})}-{\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}(\delta
n^{\mu}_{q}\lambda_{\mu})(\delta
u^{\mu}\lambda_{\mu})}-\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}(\lambda_{\mu}\delta
n^{\mu}_{q})(\lambda_{\mu}\delta u^{\mu})$
$\displaystyle+{\frac{\beta_{\pi}}{2}\left[\delta\pi^{\mu\nu}+\frac{1}{\beta_{\pi}}\lambda^{\mu}\delta
u^{\nu}\pm\frac{\gamma_{n\pi}^{q}}{\beta_{\pi}}\lambda^{\mu}\delta
n^{\nu}_{q}\right]^{2}}+{\frac{\beta_{\Pi}}{2}\left[\delta\Pi+\frac{1}{\beta_{\Pi}}\delta
u^{\mu}\lambda_{\mu}\pm\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}\delta
n^{\mu}_{q}\lambda_{\mu}\right]^{2}}-\frac{\varepsilon+p}{2}\gamma_{\nu}^{\mu}\delta
u_{\mu}\delta u^{\nu}$
$\displaystyle+\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}\gamma_{\nu}^{\mu}\delta
n^{q}_{\mu}\delta{n}^{\nu}_{q^{\prime}}-{\frac{\lambda^{2}}{\beta_{\pi}}\gamma_{\nu}^{\mu}\delta
u_{\mu}\delta u^{\nu}-\frac{1}{\beta_{\pi}}(\lambda_{\mu}\delta
u^{\mu})^{2}}{-\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}\lambda^{2}\gamma_{\nu}^{\mu}\delta
n_{\mu}^{q}\delta
n^{\nu}_{q^{\prime}}}-{\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{2}\gamma_{\nu}^{\mu}\delta
n^{\nu}_{q}\delta u_{\nu}}$ (107)

Now by looking at the terms

$\displaystyle e$
$\displaystyle=\bigg{(}\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}{-\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}\lambda^{2}\bigg{)}\gamma_{\nu}^{\mu}\delta
n_{\mu}^{q}\delta
n^{\nu}_{q^{\prime}}}-{\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{2}\gamma_{\nu}^{\mu}\delta
n^{\nu}_{q}\delta u_{\mu}}$
$\displaystyle=\bigg{(}\sum_{q}\frac{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}{2\beta_{\pi}}\bigg{)}\gamma_{\nu}^{\mu}\delta
n_{\mu}^{q}\delta
n^{\nu}_{q^{\prime}}-\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{2}\gamma_{\nu}^{\mu}\delta
n^{\nu}_{q}\delta u_{\mu}$
$\displaystyle=\bigg{(}\sum_{q}\frac{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}{2\beta_{\pi}}\bigg{)}\bigg{[}\gamma_{\nu}^{\mu}\delta
n_{\mu}^{q}\delta
n^{\nu}_{q^{\prime}}-\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}\gamma_{\nu}^{\mu}\delta
n^{\nu}_{q}\delta u_{\mu}\bigg{]}$ (108)

where completing the square the gives

$\displaystyle e$ $\displaystyle=\bigg{[}\gamma_{\nu}^{\mu}\delta
n_{\mu}^{q}-\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}\gamma_{\nu}^{\mu}\delta
u_{\mu}\bigg{]}\bigg{[}\gamma^{\nu}_{\mu}\delta
n^{\mu}_{q}-\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}\gamma^{\nu}_{\mu}\delta
u^{\mu}\bigg{]}$ $\displaystyle=\gamma_{\nu}^{\mu}\gamma^{\nu}_{\mu}\delta
n_{\mu}^{q}\delta
n^{\mu}_{q^{\prime}}+\bigg{(}\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}\bigg{)}^{2}\gamma_{\nu}^{\mu}\delta
u_{\mu}\gamma^{\nu}_{\mu}\delta
u^{\mu}-2\bigg{(}\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}\bigg{)}\gamma_{\nu}^{\mu}\delta
u_{\mu}\gamma^{\nu}_{\mu}\delta u^{\mu}$
$\displaystyle*=\bigg{(}\sum_{q}\frac{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}{2\beta_{\pi}}\bigg{)}\bigg{[}\gamma_{\nu}^{\mu}\gamma^{\nu}_{\mu}\delta
n_{\mu}^{q}\delta
n^{\mu}_{q^{\prime}}+\bigg{(}\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}\bigg{)}^{2}\gamma_{\nu}^{\mu}\delta
u_{\mu}\gamma^{\nu}_{\mu}\delta u^{\mu}$
$\displaystyle-2\bigg{(}\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}\bigg{)}\gamma_{\nu}^{\mu}\delta
u_{\mu}\gamma^{\nu}_{\mu}\delta n_{q}^{\mu}\bigg{]}$
$\displaystyle*=\bigg{(}\sum_{q}\frac{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}{2\beta_{\pi}}\bigg{)}\bigg{[}\gamma_{\nu}^{\mu}\gamma^{\nu}_{\mu}\delta
n_{\mu}^{q}\delta
n^{\mu}_{q^{\prime}}\bigg{]}+\bigg{(}\frac{(\delta_{n\pi}^{q}\lambda^{2})^{2}}{2\beta_{\pi}(2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2})}\bigg{)}\bigg{[}\gamma_{\nu}^{\mu}\delta
u_{\mu}\gamma^{\nu}_{\mu}\delta u^{\mu}\bigg{]}$
$\displaystyle-2\bigg{(}\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}}\bigg{)}\gamma_{\nu}^{\mu}\delta
u_{\mu}\gamma^{\nu}_{\mu}\delta n_{q}^{\mu}\bigg{]}$ (109)

We then arrive at

$\displaystyle\bigg{(}\sum_{q}\frac{\beta_{n}^{qq^{\prime}}}{2}{-\frac{(\delta_{n\pi}^{q})^{2}}{\beta_{\pi}}\lambda^{2}\bigg{)}\gamma_{\nu}^{\mu}\delta
n_{\mu}^{q}\delta
n^{\nu}_{q^{\prime}}}-{\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}\lambda^{2}\gamma_{\nu}^{\mu}\delta
n^{\nu}_{q}\delta u_{\mu}}$
$\displaystyle=\bigg{(}\sum_{q}\frac{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}{2\beta_{\pi}}\bigg{)}$
$\displaystyle\times\bigg{[}\gamma_{\nu}^{\mu}\delta
n_{\mu}^{q}-\frac{\delta_{n\pi}^{q}\lambda^{2}}{2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2}}\gamma_{\nu}^{\mu}\delta
u_{\mu}\bigg{]}^{2}$
$\displaystyle-\bigg{(}\frac{(\delta_{n\pi}^{q}\lambda^{2})^{2}}{2\beta_{\pi}(2\beta_{\pi}\beta_{n}^{qq^{\prime}}-(\delta_{n\pi}^{q})^{2}\lambda^{2})}\bigg{)}\bigg{[}\gamma_{\nu}^{\mu}\delta
u_{\mu}\gamma^{\nu}_{\mu}\delta u^{\mu}\bigg{]},$ (110)

which provides two additional constraints. Likewise, I collect the terms

$\displaystyle e$
$\displaystyle=\bigg{(}\frac{\beta_{n}^{qq^{\prime}}}{2\lambda^{2}}-\frac{({\delta_{n\pi}^{q}})^{2}}{\beta_{\pi}}-\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}-{\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}}-T^{2}(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}\bigg{)}(\delta n^{\mu}_{q}\lambda_{\mu})^{2}$
$\displaystyle-\bigg{(}\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}+\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}-\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\bigg{)}(\delta n^{\mu}_{q}\lambda_{\mu})(\delta
u^{\mu}\lambda_{\mu}).$

To make things a bit more compact, I will define

$\displaystyle e_{n,u}$
$\displaystyle=\frac{\beta_{n}^{qq^{\prime}}}{2\lambda^{2}}-\frac{({\delta_{n\pi}^{q}})^{2}}{\beta_{\pi}}-\frac{(\delta_{n\Pi}^{q})^{2}}{\beta_{\Pi}}-{\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}\Bigg{)}^{2}}-T^{2}(\varepsilon+p)\left.\frac{\partial
s}{\partial\varepsilon}\right|_{p}\left.\frac{\partial s}{\partial
p}\right|_{\alpha_{q}}\Bigg{(}\left.\frac{\partial{\alpha_{q}}}{\partial
s}\right|_{p}\Bigg{)}^{2}$ (111) $\displaystyle C_{n,u}$
$\displaystyle=-\frac{\delta_{n\Pi}^{q}}{\beta_{\Pi}}+\frac{\delta_{n\pi}^{q}}{\beta_{\pi}}-\frac{\varepsilon+p}{2}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}\left.\frac{\partial{\alpha_{q}}}{\partial
p}\right|_{s}$ $\displaystyle e$ $\displaystyle=e_{n,u}\,\big{[}\delta
n_{\mu}^{q}\delta n^{\mu}_{q^{\prime}}-\frac{C_{n,u}}{e_{n,u}}\delta
n^{\mu}_{q}\delta u_{\mu}\big{]}$ $\displaystyle\bigg{(}\delta
n_{\mu}^{q}-\frac{C_{n,u}}{e_{n,u}}\delta u_{\mu}\bigg{)}^{2}=\delta
n_{\mu}^{q}\delta
n^{\mu}_{q^{\prime}}+\bigg{(}\frac{C_{n,u}}{e_{n,u}}\bigg{)}^{2}\delta
u^{\mu}\delta u_{\mu}-2\,\frac{C_{n,u}}{e_{n,u}}\delta n^{\mu}_{q}\delta
u_{\mu}$

to get to

$\displaystyle\delta n_{\mu}^{q}\delta
n^{\mu}_{q^{\prime}}-2\,\frac{C_{n,u}}{e_{n,u}}\delta n^{\mu}_{q}\delta
u_{\mu}=\bigg{(}\delta n_{\mu}^{q}-\frac{C_{n,u}}{e_{n,u}}\delta
u_{\mu}\bigg{)}^{2}-\bigg{(}\frac{C_{n,u}}{e_{n,u}}\bigg{)}^{2}\delta
u^{\mu}\delta u_{\mu}.$ (112)

Finally, I group all terms with $\delta u^{\mu}\lambda_{\mu}$ from Eq. (107)
and the expression in Eq. (112) to obtain the final condition

$\displaystyle-{\frac{\varepsilon+p}{4}\left.\frac{\partial
p}{\partial\varepsilon}\right|_{s}(\delta
u^{\mu}\lambda_{\mu})^{2}}-{\frac{1}{\beta_{\Pi}}(\delta
u^{\mu}\lambda_{\mu})^{2}}-{\frac{2}{3\beta_{\pi}}(\delta
u^{\mu}\lambda_{\mu})^{2}}-\frac{\varepsilon+p}{2}\frac{(\lambda_{\mu}\delta
u^{\mu})^{2}}{\lambda^{2}}-\frac{C_{n,u}^{2}}{e_{n,u}}\delta u^{\mu}\delta
u_{\mu}$ (113)

Which provides a constraint associated with $\delta u^{\mu}\lambda_{\mu}$.