# Causality and stability in First-order Conformal Anisotropic Hydrodynamics

Fábio S. Bemfica  Escola de Ciências e Tecnologia, Universidade Federal do Rio
Grande do Norte, 59072-970, Natal, RN, Brazil  Mauricio Martinez  Department
of Physics, North Carolina State University, Raleigh, NC 27695, USA  Masoud
Shokri  Institute for Theoretical Physics, Max-von-Laue-Str. 1, D-60438
Frankfurt am Main, Germany

###### Abstract

We formulate the first-order dissipative anisotropic hydrodynamical theory for
a relativistic conformal uncharged fluid, which generalizes the Bemfica-
Disconzi-Noronha-Kovtun (BDNK) first-order viscous fluid framework. Our
approach maintains causal behavior in the nonlinear regime with or without
general relativity coupling, and we derive and analyze the constraints on
transport coefficients imposed by causality. We demonstrate the causal and
stable behavior of our theory in specific cases, including the discussion of
nonlinear causality as well as stability for linearized perturbations. We
apply our newly developed first order anisotropic theory to the Bjorken flow
and show how causality and stability impose constraints on the behavior of the
early-time attractor.

##  I Introduction

Causality is one of the most important guiding principles in physics. In the
bottom-up approach to constructing modern effective field theories, causality
mandates that effective quantum field theories meet certain requirements,
including microcausality, causal propagation limited to the speed of light,
stability of the vacuum, and possible constraints on commutation relations de
Rham  _et al._ (2022). Like other effective field theories, hydrodynamics is
also subject to the restrictions imposed by causality. In this context,
causality mandates that a solution to the relativistic fluid equations at a
specific space-time point $x$ is entirely defined by the past space-time
region that is causally connected to $x$ Hawking and Ellis (2023); Wald
(1984). While causality is an important physical requirement for relativistic
fluid dynamical equations of motion, it is not the only one. Two additional
requirements are necessary for the equations of motion: to be locally well-
posed and stable. The latter demands that small perturbations around the
thermal state decay over time, while the former ensures that the system
follows a well-defined space-time evolution for a given set of initial
conditions. The original first order dissipative hydrodynamical equations of
motion, known as the Navier-Stokes (NS) equations and derived by Landau and
Lifschitz Landau and Lifshitz (2013) as well as Eckart Eckart (1940) showed to
be acausal in the linear and nonlinear regimes Hiscock and Lindblom (1985a);
Pichon (1965). To address the acausality and stability concerns of the Navier-
Stokes equations, second-order hydrodynamics theories were introduced by
Israel and Stewart (IS)Israel (1976); Israel and Stewart (1979). Since their
seminal work, more recent formulations of second-order hydrodynamics have been
developed Baier  _et al._ (2008); Denicol  _et al._ (2012); Hiscock and
Lindblom (1983). However, it is unclear if second order theories are causal in
the full non-linear regime.

A recent development in relativistic hydrodynamics is the BDNK theory proposed
by Bemfica, Disconzi and Noronha Bemfica  _et al._ (2018, 2019, 2021a, 2021b,
2022) alongside with Kovtun Kovtun (2019); Hoult and Kovtun (2020). The BDNK
theory has offered a practical and straightforward approach to addressing the
longstanding issues of causality, stability, and local well-posedness in
relativistic fluids. Its development has created new opportunities for
advancing our understanding of the fundamental principles underlying
relativistic fluid dynamics. Basically BDNK theory is a generalization of the
first order NS relativistic theory which is causal in the linear and non-
linear regimes 111The rigorous mathematical demonstrations of these statements
are found in Refs. Bemfica  _et al._ (2021a, b, 2022). as well as locally
well-posed in Sobolev spaces in the presence and/or absence of gravity. BDNK
theory is also linearly stable around global equilibrium. An interesting
feature of BDNK theory is that the definition of the hydrodynamical variables,
such as the energy density, receive out-of-equilibrium corrections. In this
sense BDNK first order theory differs from the standard approaches based on
the Landau and Eckart frames.

On the other hand, anisotropic hydrodynamics has emerged as a very successful
phenomenological model for describing non-equilibrium fluids with large
spatial anisotropies in high-energy nuclear collisions Martinez and Strickland
(2010); Florkowski and Ryblewski (2011); Bazow  _et al._ (2014); Molnár  _et
al._ (2016); Molnar  _et al._ (2016). It has been used to accurately model the
space-time evolution of the fireball in these collisions and has shown good
agreement with experimental results (see Sect. 10 of Ref.Alqahtani  _et al._
(2018) and references therein). Additionally, anisotropic hydrodynamics has
passed rigorous numerical tests against exact kinetic theory models based on
the Boltzmann equation (See Sect. 7-8 of Ref.Alqahtani  _et al._ (2018) and
references therein), providing confidence in its efficacy for studying fluids
in far-from-equilibrium situations. Despite these successes, the foundational
aspects of anisotropic hydrodynamics remain incompletely understood,
particularly with respect to the role of causality and its constraints.
Motivated by this gap in understanding, we introduce a new approach in this
work to investigate causality in anisotropic hydrodynamics. We develop a novel
first order theory for a conformal, uncharged fluid near thermal equilibrium
that exhibits a spatial anisotropy characterized by a space-like vector
$l^{\mu}$.

The paper is organized as follows: in section II we discuss the most general
form of a first order anisotropic conformal fluid theory. We show in Sec. III
that physical requirements of the second law of thermodynamics up to the order
of validity of the theory leads unambiguously to the anisotropic BDNK theory
of a conformal fluid. Section IV examines the conditions that causality
imposes on the transport coefficients, both with and without gravity. Linear
stability of our novel first order anisotropic theory is analyzed in Sec. V.
The application of our novel anisotropic first order theory for a fluid
undergoing Bjorken flow is presented in Sect. VI. Our conclusions are
summarized in Sec. VII. Technical details are found in the appendices.

#### Notations and conventions

We use the natural unites in which $\hbar=c=k_{B}=1$, and adopt the Lorentzian
metric $g_{\mu\nu}$ with signature $-+++$. The standard covariant derivative
is denoted by $\nabla$, and the conformal covariant derivative by
$\mathcal{D}$. We use the standard symmetrization and antisymmetrization
notations, such that, for example, for a rank-2 tensor we obtain that
$A_{\mu\nu}$ is $A_{(\mu\nu)}=\tfrac{1}{2}\left(A_{\mu\nu}+A_{\nu\mu}\right)$
and $A_{[\mu\nu]}=\tfrac{1}{2}\left(A_{\mu\nu}-A_{\nu\mu}\right)$,
respectively. The Riemann tensor is defined as
$R^{\sigma}_{\rho\mu\nu}=2\left(\partial_{[\mu}\Gamma^{\sigma}_{\nu]\rho}+\Gamma^{\sigma}_{[\mu\beta}\Gamma^{\beta}_{\nu]\rho}\right)$,
and the Ricci tensor as $R_{\mu\nu}=R_{\rho\mu\sigma\nu}g^{\rho\sigma}$.

##  II First order anisotropic conformal BDNK theory

A common approach in hydrodynamics is to expand physical observables in terms
of gradients and truncate the series expansion at a given order in the
derivatives. For example, if one considers the energy-momentum tensor
$T^{\mu\nu}$, the gradient expansion takes the form
$T^{\mu\nu}=\mathcal{O}(1)+\mathcal{O}(\partial)+\cdots$, where
$\mathcal{O}(1)$ corresponds to the ideal fluid contributions while
$\mathcal{O}(\partial^{n})$ are the viscous corrections of order $n$ in
derivatives (time and space derivatives written in covariant form) of the
dynamical variables such as energy, density, etc. Therefore, first-order
theories contain only the first-order derivative corrections
$\mathcal{O}(\partial)$ 222The gradient expansion is equivalent to the Knudsen
expansion in kinetic theory Cercignani (1988); Groot (1980).. Since quantities
such as temperature and number density are only unambiguously defined in
equilibrium, their out-of-equilibrium corrections may differ due to different
choices of $\mathcal{O}(\partial)$. In fact, different first-order theories
are connected by transformations in these out-of-equilibrium variables 333For
example, consider a fluid’s theory defined by a set of $N$ out-of-equilibrium
thermodynamic variables ${\psi_{a}}$ with $a=1,...,N$, and suppose we perform
a transformation $\psi_{a}\to\psi_{a}+\delta\psi_{a}$, where $\delta\psi_{a}$
is first order in the derivative of the thermodynamic variables $\psi_{a}$.
This leads to the existence of two different first-order theories that are, in
fact, equivalent up to the first-order corrections.. The connection between
the Landau-Lifshitz and Eckart frames, for instance, can be understood through
such transformations Israel (1976). Recent works have provided a comprehensive
discussion on this subtle yet crucial aspect Kovtun (2019); Bemfica  _et al._
(2022); Dore  _et al._ (2022).

We shall develop the first order theory of a conformal uncharged fluid which
has a spatial anisotropy along an arbitrary direction determined by a space-
like vector $l^{\mu}$. We start by introducing the most general form of the
energy-momentum tensor that describes this particular fluid.

###  II.1 Anisotropic energy momentum tensor

Here, we turn to the problem of constructing a nonlinearly causal anisotropic
extension of the conformal BDNK hydrodynamics for an uncharged fluid Bemfica
_et al._ (2018). The most general form for the energy-momentum tensor
$T^{\mu\nu}$, decomposed in the directions parallel and orthogonal to the
time-like velocity flow $u^{\mu}$ ($u^{\mu}u_{\mu}=-1$) and the anisotropic
space-like vector $l^{\mu}$ ($l^{\mu}l_{\mu}=1$) which is invariant under the
little group $SO(2)$, reads as Molnár  _et al._ (2016); Molnar  _et al._
(2016)

$\displaystyle T^{\mu\nu}$ $\displaystyle=$
$\displaystyle\mathcal{E}\,u^{\mu}u^{\nu}+\mathcal{P}_{l}\,l^{\mu}l^{\nu}+\mathcal{P}_{\perp}\,\Xi^{\mu\nu}+2\,M\,u^{\left(\mu\right.}l^{\left.\nu\right)}+2\,W_{\perp
u}^{\left(\mu\right.}u^{\left.\nu\right)}$ (1) $\displaystyle+2\,W_{\perp
l}^{\left(\mu\right.}l^{\left.\nu\right)}+\pi_{\perp}^{\mu\nu}$
$\displaystyle=$
$\displaystyle\mathcal{E}\,u^{\mu}u^{\nu}+(\mathcal{P}_{l}-\mathcal{P}_{\perp})\,l^{\mu}l^{\nu}+\mathcal{P}_{\perp}\,\Delta^{\mu\nu}+2\,M\,u^{\left(\mu\right.}l^{\left.\nu\right)}+2\,W_{\perp
u}^{\left(\mu\right.}u^{\left.\nu\right)}$ $\displaystyle+2\,W_{\perp
l}^{\left(\mu\right.}l^{\left.\nu\right)}+\pi_{\perp}^{\mu\nu}\;.$

Here, $\Delta^{\mu\nu}=g^{\mu\nu}+u^{\mu}u^{\nu}$ and
$\Xi^{\mu\nu}=\Delta^{\mu\nu}-l^{\mu}l^{\nu}$ are the projectors orthogonal to
$u^{\mu}$, and to $u^{\mu}$ and $l^{\mu}$ together, respectively. Quantities
$\mathcal{E}=\mathcal{P}_{l}+2\mathcal{P}_{\perp},\,M,\,W_{\perp
u}^{\mu},\,W_{\perp l}^{\mu}$, and $\pi_{\perp}^{\mu\nu}$ contain four
scalars, two vectors, and one tensor build out of the leading-order
hydrodynamics fields and their first-order derivatives, with their concrete
forms to be determined subsequently. In particular,
$\mathcal{E}=\varepsilon+\mathcal{E}^{(1)}$, where $\mathcal{E}^{(1)}$ is the
first order correction to the energy and $\varepsilon$ is the equilibrium
energy density. Assuming $T$ to be an effective temperature, then, in the
conformal case, $\varepsilon=\varepsilon(T)\propto T^{4}$ and
$\varepsilon=P_{l}+2P_{\perp}$, where $P_{l}$ and $P_{\perp}$ are the leading
order contributations to $\mathcal{P}_{l}$ and $\mathcal{P}_{\perp}$,
respectively. Due to the conformal invariance, the derivatives of the scalar
fields are related through
$\nabla_{\mu}P_{l}=(P_{l}/\varepsilon)\nabla_{\mu}\varepsilon$ and
$\nabla_{\mu}P_{\perp}=(P_{\perp}/\varepsilon)\nabla_{\mu}\varepsilon$.

In Eq. (1) the vectors $W_{\perp\,l}^{\mu}$, $W_{\perp\,u}^{\mu}$, and the
symmetric traceless tensor $\pi^{\mu\nu}_{\perp}$ are orthogonal to $u$ and
$l$, i.e.

$V_{\mu}W_{\perp\,u}^{\mu}=V_{\mu}W_{\perp\,l}^{\mu}=V_{\mu}\pi^{\mu\nu}_{\perp}=0\;,$

with $V=\\{u,l\\}$. Note that the above expressions leads to
$\Xi_{\mu\nu}\pi_{\perp}^{\mu\nu}=0$. In the absence of other conserved
currents, the fluid’s evolution is governed by the energy-momentum
conservation $\nabla_{\nu}T^{\mu\nu}=0$, which will be referred to as the
equations of motion (EOM) and due to the conformal invariance can be
equivalently expressed as $\mathcal{D}_{\nu}T^{\mu\nu}=0$, where
$\mathcal{D}_{\mu}$ is the conformal covariant derivative compatible with
$u^{\mu}$ Loganayagam (2008). In particular

$\mathcal{D}_{\mu}u_{\nu}=\Delta^{\alpha}_{\mu}\nabla_{\nu}u_{\nu}-\Delta_{\mu\nu}\nabla_{\alpha}u^{\alpha}/3\,,$

and

$\mathcal{D}_{\mu}\varepsilon=\nabla_{\mu}\varepsilon+4\varepsilon(u^{\alpha}\nabla_{\alpha}u_{\mu}-u_{\mu}\nabla_{\alpha}u^{\alpha}/3)\;.$

The EOM may be decomposed into the directions parallel to $u$ and $l$ and the
directions perpendicular to both as

$\displaystyle u^{\alpha}\mathcal{D}_{\alpha}\mathcal{E}$ $\displaystyle=$
$\displaystyle-(\mathcal{P}_{l}-\mathcal{P}_{\perp})l^{\mu}l^{\nu}\sigma_{\mu\nu}-M\mathcal{D}_{\alpha}l^{\alpha}-l^{\nu}\mathcal{D}_{\nu}M-\mathcal{D}_{\nu}W^{\nu}_{\perp
u}-2W_{\perp l}^{\nu}l^{\alpha}\sigma_{\mu\nu}$ (2a)
$\displaystyle-\pi_{\perp}^{\mu\nu}\sigma_{\perp\mu\nu},$ $\displaystyle
l^{\alpha}\mathcal{D}_{\alpha}\mathcal{P}_{l}$ $\displaystyle=$
$\displaystyle-(\mathcal{P}_{l}-\mathcal{P}_{\perp})\mathcal{D}_{\alpha}l^{\alpha}-M\,l^{\mu}l^{\nu}\sigma_{\mu\nu}-u^{\alpha}\mathcal{D}_{\alpha}M-l_{\nu}W_{\perp
u}^{\alpha}\mathcal{D}_{\alpha}u^{\nu}-\,l_{\nu}u^{\alpha}\mathcal{D}_{\alpha}W_{\perp
u}^{\nu}$ (2b) $\displaystyle-l_{\mu}l^{\nu}\mathcal{D}_{\nu}W^{\mu}_{\perp
l}-\mathcal{D}_{\alpha}W_{\perp
l}^{\alpha}+\pi_{\perp}^{\nu\alpha}\mathcal{D}_{\alpha}l_{\nu},$
$\displaystyle\Xi^{\mu\alpha}\mathcal{D}_{\alpha}\mathcal{P_{\perp}}$
$\displaystyle=$
$\displaystyle-(\mathcal{P}_{l}-\mathcal{P}_{\perp})\,\Xi^{\mu}_{\nu}l^{\alpha}\mathcal{D}_{\alpha}l^{\nu}-M\,\Xi^{\mu}_{\nu}l^{\alpha}\mathcal{D}_{\alpha}u^{\nu}-\,M\Xi^{\mu}_{\nu}u^{\alpha}\mathcal{D}_{\alpha}l^{\nu}-\Xi^{\mu}_{\nu}W_{\perp
u}^{\alpha}\mathcal{D}_{\alpha}u^{\nu}$ (2c)
$\displaystyle-\Xi^{\mu}_{\nu}u^{\alpha}\mathcal{D}_{\alpha}W_{\perp
u}^{\nu}-W_{\perp
l}^{\mu}\mathcal{D}_{\alpha}l^{\alpha}-\Xi^{\mu}_{\nu}W_{\perp
l}^{\alpha}\mathcal{D}_{\alpha}l^{\nu}-\Xi^{\mu}_{\nu}l^{\alpha}\mathcal{D}_{\alpha}W_{\perp
l}^{\nu}-\Xi^{\mu}_{\nu}\mathcal{D}_{\alpha}\pi_{\perp}^{\nu\alpha},$

where
$\sigma_{\mu\nu}=\mathcal{D}_{(\mu}u_{\nu)}=\Delta^{\alpha\beta}_{\mu\nu}\nabla_{\mu}u_{\nu}$
is the shear tensor, with
$\Delta^{\alpha\beta}_{\mu\nu}=[\Delta^{\alpha}_{\mu}\Delta^{\beta}_{\nu}+\Delta^{\beta}_{\mu}\Delta^{\alpha}_{\nu}-(2/3)\Delta^{\alpha\beta}\Delta_{\mu\nu}]/2$.
We have also introduced a transverse shear tensor as

$\sigma_{\perp\mu\nu}=\Xi^{\alpha\beta}_{\mu\nu}\mathcal{D}_{\alpha}u_{\beta}=\Xi^{\alpha\beta}_{\mu\nu}\nabla_{\alpha}u_{\beta}\;,\mbox{\quad
with\quad}\Xi^{\alpha\beta}_{\mu\nu}=(\Xi^{\alpha}_{\mu}\Xi^{\beta}_{\nu}+\Xi^{\beta}_{\mu}\Xi^{\alpha}_{\nu}-\Xi^{\alpha\beta}\Xi_{\mu\nu})/2\;.$

It is important to keep in mind that an expression that only includes
derivatives up to a certain order may be approximated at that order when the
equations of motion are considered. As an example, one can rewrite (2a) by
separating the leading and first order contributions from $\mathcal{E}$ and
$\mathcal{P}$ to obtain

$\displaystyle
u^{\alpha}\mathcal{D}_{\alpha}\varepsilon+(P_{l}-P_{\perp})l^{\mu}l^{\nu}\sigma_{\mu\nu}$
$\displaystyle=$
$\displaystyle-u^{\alpha}\mathcal{D}_{\alpha}\mathcal{E}^{(1)}-(\mathcal{P}_{l}^{(1)}-\mathcal{P}_{\perp}^{(1)})l^{\mu}l^{\nu}\sigma_{\mu\nu}-M\mathcal{D}_{\alpha}l^{\alpha}-l^{\alpha}\mathcal{D}_{\alpha}M-\mathcal{D}_{\nu}W^{\nu}_{\perp
u}$ (3) $\displaystyle-\,2W_{\perp
l}^{\alpha}l^{\nu}\sigma_{\mu\nu}-\pi_{\perp}^{\mu\nu}\sigma_{\perp\mu\nu}.$

Upon comparing the left and right hand sides of the equation above, it becomes
evident that the terms $u^{\alpha}\mathcal{D}_{\alpha}\varepsilon$ and
$l^{\mu}l^{\nu}\sigma_{\mu\nu}$ only contain first-order derivatives of the
zeroth-order quantities. However, the combination
$u^{\alpha}\mathcal{D}_{\alpha}\varepsilon+(P_{l}-P_{\perp})l^{\mu}l^{\nu}\sigma_{\mu\nu}$
is second-order on-shell, meaning that it becomes second order when the
equations of motion given by (3) are taken into account. It can be said that
the same combination is of the first-order off-shell. Similar arguments may be
applied in the remaining equations of motion in (2) above. In the next
section, we make use of generic arguments based on the second law of
thermodynamics together with the above on-shell analysis in order to determine
the remaining first order terms $M$, $E^{(1)}$, etc, by assuming the physical
applicability of a first-order theory.

##  III Entropy and entropy production

The second law of thermodynamics mandates that entropy should be at a maximum
in a state of equilibrium, and the entropy production should not have a
negative value. However, for a first-order theory one must only consider the
entropy up to the first-order on-shell derivative, and hence the entropy
production up to the second-order on-shell Kovtun (2019). This ensures that
the theory is being applied outside of its intended physical applicability.

The entropy current of an anisotropic conformal uncharged fluid reads as (see
Appendix A)

$S^{\mu}=\frac{P_{\perp}u^{\mu}-u_{\nu}T^{\mu\nu}}{T}\;,$ (4)

where $P_{\perp}$ is the leading order of $\mathcal{P}_{\perp}$, i.e.,
$\mathcal{P}_{\perp}=P_{\perp}+\mathcal{P}_{\perp}^{(1)}$. In this case, the
entropy density reads

$-u_{\mu}S^{\mu}=s+\frac{\mathcal{E}^{(1)}}{T}\;,$ (5)

where the leading order contribution to the entropy density is
$s=(\varepsilon+P_{\perp})/T$. From Eq. (4) we calculate the entropy
production

$\displaystyle\nabla_{\mu}S^{\mu}$ $\displaystyle=$
$\displaystyle-\frac{\pi_{\perp}^{\mu\nu}\sigma_{\perp\mu\nu}}{T}-\frac{2W_{\perp
l}^{\mu}l^{\nu}\sigma_{\mu\nu}}{T}-\frac{\varepsilon-3P_{\perp}}{4T}\frac{u^{\alpha}\mathcal{D}_{\alpha}\varepsilon}{\varepsilon}-\frac{(\mathcal{P}_{l}-\mathcal{P}_{\perp})l^{\mu}l^{\nu}\sigma_{\mu\nu}}{T}$
(6) $\displaystyle-\frac{\mathcal{E}^{(1)}u^{\nu}+W_{\perp
u}^{\nu}+M\,l^{\nu}}{4T}\frac{\mathcal{D}_{\nu}\varepsilon}{\varepsilon}\;.$

The choices $\pi_{\perp\,\mu\nu}\propto-\sigma_{\perp\,\mu\nu}$ and
$W_{\perp\,l}^{\mu}\propto-\Xi^{\mu\nu}l^{\alpha}\sigma_{\alpha\nu}$ give
positive contributions to the first two terms in (6). On the other hand, the
third term

$-\frac{\varepsilon-3P_{\perp}}{4T}\frac{u^{\alpha}\mathcal{D}_{\alpha}\varepsilon}{\varepsilon}$
(7)

cannot be positive definite since $\mathcal{D}_{\alpha}\varepsilon$ has no
definite sign. Furthermore, the order of derivatives of this term, as given by
Eq. (3), depends on the leading order of pressures and could be either first
or second order. Thus, positive entropy production in Eq. (6) demands the
elimination of the term in Eq. (7) by means of the choice

$P_{\perp}=P_{l}=P=\frac{\varepsilon}{3}\;.$ (8)

It is worth mentioning that this choice for the pressures applied to Eqs. (2)
makes $u^{\alpha}\mathcal{D}_{\alpha}\varepsilon$,
$l^{\alpha}\mathcal{D}_{\alpha}\varepsilon$, and
$\Xi^{\alpha\beta}\mathcal{D}_{\beta}\varepsilon$ of order $\partial^{3}$ on-
shell. This means that
$\mathcal{D}_{\alpha}\varepsilon=\mathcal{O}(\partial^{3})$ when Eqs. (2)
apply. After imposing (8), the term

$\frac{(\mathcal{P}_{l}-\mathcal{P}_{\perp})l^{\mu}l^{\nu}\sigma_{\mu\nu}}{T}=\frac{(\mathcal{P}_{l}^{(1)}-\mathcal{P}_{\perp}^{(1)})l^{\mu}l^{\nu}\sigma_{\mu\nu}}{T}$
(9)

becomes second order off-shell since $\pi_{\mu\nu}$ is of first order on- and
off-shell. Hence, the first order viscous corrections $\mathcal{P}_{l}^{(1)}$
and $\mathcal{P}_{\perp}^{(1)}$ must be chosen in such way that their
difference is positive up to second order when Eqs. (2) apply. Subsequently,
and by following the prescription outlined in the previous section, the
simplest anisotropic extension of the conformal BDNK hydrodynamics of an
uncharged fluid is formulated by the energy-momentum tensor

$\displaystyle T^{\mu\nu}$ $\displaystyle=$
$\displaystyle(\varepsilon+\mathcal{E}^{(1)})\,u^{\mu}u^{\nu}+(\mathcal{P}_{l}^{(1)}-\mathcal{P}_{\perp}^{(1)})\,l^{\mu}l^{\nu}+(P+\mathcal{P}_{\perp}^{(1)})\,\Delta^{\mu\nu}+2\,M\,u^{\left(\mu\right.}l^{\left.\nu\right)}+2\,W_{\perp
u}^{\left(\mu\right.}u^{\left.\nu\right)}$ (10) $\displaystyle+2\,W_{\perp
l}^{\left(\mu\right.}l^{\left.\nu\right)}+\pi_{\perp}^{\mu\nu}\;,$

where its components are written as

$\displaystyle\mathcal{E}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{3\chi}{4}\frac{u^{\mu}\mathcal{D}_{\mu}\varepsilon}{\varepsilon}\;,$
(11a) $\displaystyle\mathcal{P}_{l}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{\chi_{l}}{4}\frac{u^{\mu}\mathcal{D}_{\mu}\varepsilon}{\varepsilon}-2\eta_{ll}l^{\alpha}l^{\beta}\sigma_{\alpha\beta}\;,$
(11b) $\displaystyle\mathcal{P}_{\perp}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{\chi_{\perp}}{4}\frac{u^{\mu}\mathcal{D}_{\mu}\varepsilon}{\varepsilon}+\eta_{ll}l^{\alpha}l^{\beta}\sigma_{\alpha\beta},$
(11c) $\displaystyle\pi_{\perp}^{\mu\nu}$ $\displaystyle=$
$\displaystyle-2\eta_{\perp}\sigma_{\perp}^{\mu\nu}\;,$ (11d) $\displaystyle
W_{\perp l}^{\mu}$ $\displaystyle=$
$\displaystyle-2\eta_{l}\Xi^{\mu}_{\lambda}l_{\nu}\sigma^{\lambda\nu}\;,$
(11e) $\displaystyle W_{\perp u}^{\mu}$ $\displaystyle=$
$\displaystyle\frac{\lambda_{\perp}}{4}\Xi^{\mu\nu}\frac{\mathcal{D}_{\nu}\varepsilon}{\varepsilon}\;,$
(11f) $\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{\lambda_{l}}{4}\frac{l^{\nu}\mathcal{D}_{\nu}\varepsilon}{\varepsilon}\;.$
(11g)

In the previous expression we have 8 transport coefficients
$\\{\chi,\chi_{l},\chi_{\perp},\eta_{l},\eta_{ll},\eta_{\perp},\lambda_{l},\lambda_{\perp}\\}$
which are proportional to $\varepsilon^{3/4}$, due to the conformal
invariance. The conformal invariance also requires
$\chi_{l}+2\chi_{\perp}=3\chi$ to ensure $T^{\mu}_{\mu}=0$.

Plugging Eqs. (11) into the EOMs (2) and assuming condition (8), the entropy
density is given by
$-u_{\mu}S^{\mu}=(\varepsilon+P)/T+\mathcal{O}(\partial^{2})$ and the on-shell
entropy production is simply

$\nabla_{\mu}S^{\mu}=\frac{\pi_{\perp}^{\mu\nu}\pi_{\perp\mu\nu}}{2\eta_{\perp}T}+\frac{W_{\perp
l}^{\mu}W_{\perp
l\,\mu}}{\eta_{l}T}+\frac{\left(\mathcal{P}_{\perp}^{(1)}-\mathcal{P}_{l}^{(1)}\right)^{2}}{3\eta_{ll}T}+\mathcal{O}(\partial^{3})\;.$
(12)

We conclude that the on-shell entropy production given by the previous
expression is positive up to second order in derivatives if

$\eta_{\perp}\geq 0\;,\qquad\eta_{ll}\geq 0\;,\qquad\eta_{l}\geq 0\;.$ (13)

###  III.1 The isotropic conformal limit of BDNK

Once the space-like anisotropic vector $l$ is chosen, the isotropic limit for
the conformal uncharged fluid is reproduced by setting

$\eta_{\perp}=\eta_{l}=\eta_{ll}=\eta\;,\qquad\chi_{l}=\chi_{\perp}=\chi\;,\qquad\lambda_{\perp}=\lambda_{l}=\lambda\;.$
(14)

If we plug the above condition into (10), the energy-momentum tensor reduces
to the one derived in BDNK theory (Bemfica  _et al._ , 2018), i.e.,

$\displaystyle T^{\mu\nu}$ $\displaystyle=$
$\displaystyle\left[\varepsilon+\frac{3\chi}{4}\frac{u^{\mu}\mathcal{D}_{\mu}\varepsilon}{\varepsilon}\right]u^{\mu}u^{\nu}+\left[P+\frac{\chi}{4}\frac{u^{\mu}\mathcal{D}_{\mu}\varepsilon}{\varepsilon}\right]\Delta^{\mu\nu}-2\eta\sigma^{\mu\nu}+\frac{\lambda
u^{(\mu}\Delta^{\nu)\alpha}}{2}\frac{\mathcal{D}_{\alpha}\varepsilon}{\varepsilon}\;,$
(15)

where $\chi,\lambda,\eta\propto\varepsilon^{3/4}$. The equivalence between the
expressions derived above and the corresponding one in BDNK theory Bemfica
_et al._ (2018) becomes clearer when considering the limit established by
Eqs.(14) while writing $\Delta^{\mu\nu}=\Xi^{\mu\nu}+l^{\mu}l^{\nu}$ in
Eq.(15) together with the following identity Molnár  _et al._ (2016); Molnar
_et al._ (2016)

$\displaystyle\sigma_{\mu\nu}$ $\displaystyle=$
$\displaystyle\sigma_{\perp\mu\nu}-\frac{1}{2}\Xi_{\mu\nu}l^{\alpha}l^{\beta}\sigma_{\alpha\beta}+2\Xi^{\alpha}_{(\mu}l_{\nu)}l^{\beta}\sigma_{\alpha\beta}+l_{\mu}l_{\nu}l^{\alpha}l^{\beta}\sigma_{\alpha\beta}.$
(16)

Note that the traceless components of the energy momentum tensor (15) are
$\sigma_{\perp\mu\nu}$,
$\Xi^{\alpha}_{(\mu}l_{\nu)}l^{\beta}\sigma_{\alpha\beta}$, and the
combination
$\Xi_{\mu\nu}l^{\alpha}l^{\beta}\sigma_{\alpha\beta}-2l_{\mu}l_{\nu}l^{\alpha}l^{\beta}\sigma_{\alpha\beta}$.
These terms are needed and explain why these multiply the same coefficient
$\eta_{ll}$ in Eqs. (11b) and (11c).

##  IV Causality

Now, we shall show that there exist conditions for the transport parameters
defined in Eqs. (11) such that the resulting hydrodynamic theory exhibits
nonlinear causality. It is worth noting that the EOM
$\nabla_{\nu}T^{\mu\nu}=0$ give rise to a system of quasi-linear partial
differential equations (PDE), i.e., PDEs that do not contain products of their
highest-order derivative terms. The independent variables in the energy-
momentum tensor include the energy density $\varepsilon$ and the components of
the fluid’s velocity $u^{\mu}$, with the anisotropic vector $l^{\mu}$ 444The
components of the vector $l^{\mu}$ are constrained by being a unitary space-
like vector $l^{\mu}l_{\mu}=1$ and being orthogonal to $u^{\mu}$. For
instance, at the earliest stages of a heavy-ion collision, it is common to
choose $u^{\mu}=\gamma(v_{z})(1,0,0,v_{z})$ with
$l^{\mu}=\gamma(v_{z})(v_{z},0,0,1)$ Molnar  _et al._ (2016). being chosen
accordingly. To investigate causality in quasi-linear systems, we analyze the
principal part of the system of equations. This part contains only the terms
of the highest order in each variable and determines the order of the partial
differential equations in each variable. As an example, a quasi-linear system
that is of order one in $\varepsilon$ and two in $u^{\mu}$ does not contain
terms like $(\partial\varepsilon)^{2}$, $(\partial^{2}u)^{2}$, or
$\partial\varepsilon\,\partial^{2}u$, and its principal part can contain only
terms of form $\partial\varepsilon$ or $\partial^{2}u$. As explained in
Bemfica  _et al._ (2018), we assume $u^{\mu}$ to have four independent
components, with the constraint $u^{\mu}u_{\mu}=-1$ being imposed at an
initial time and being preserved through the fluid’s evolution. The EOM is
then decomposed in directions parallel and perpendicular to $u^{\mu}$. This
gives rise to 5 equations of the 5 independent variables

$-u_{\mu}\nabla_{\nu}T^{\mu\nu}=0,\quad\Delta^{\mu}_{\alpha}\nabla_{\nu}T^{\alpha\nu}=0\;,$
(17)

with the constraint $u^{\mu}u_{\mu}=-1$ being used when required 555This is
equivalent of considering $\nabla_{\nu}T^{\mu\nu}=0$ together with the
evolution equation for the constraint
$u^{\beta}\nabla_{\beta}[u^{\alpha}\nabla_{\alpha}(u^{\mu}u_{\mu})]=0$. Let us
now consider Eqs. (17) together with the constitutive relation (10) for the
energy-momentum tensor, and the dissipative fluxes given in (11), coupled with
gravity through Einstein’s equation

$R^{\mu\nu}-\frac{1}{2}g^{\mu\nu}R=8\pi GT^{\mu\nu}\;,$ (18)

and the gauge freedom fixed by assuming the harmonic gauge
$g^{\mu\nu}\Gamma_{\mu\nu}^{\alpha}=0$. As explained in Appendix B, the
causality of the system is determined by the vectors
$\xi_{\mu}=\nabla_{\mu}\Phi$ normal to the characteristic hypersurface
$\Phi(x)=0$, which are the roots of the characteristic equation (66). With the
roots of the characteristic equation obtained as $\xi_{0}=\xi_{0}(\xi_{i})$,
the system is causal if (C1) they are real,

$\xi_{0}\in\mathbb{R}\;,$ (19)

and (C2) $\xi_{\alpha}=(\xi_{0},\xi_{i})$ is not timelike, i.e.,

$\xi_{\alpha}\xi^{\alpha}\geq 0\;.$ (20)

In our case, the characteristic equation has 30 roots, of which 20 are light-
like and thus causal, arising from pure gravity. The remaining pieces of the
characteristic equation are spacelike roots, which we refer to as the matter
sector, can be further decomposed into two parts. The first part, which
contains four roots, is

$A^{2}=\left(\lambda_{\perp}a^{2}-\eta_{\perp}v^{2}+\delta\eta_{\perp
l}b^{2}\right)^{2}=0\;,$ (21)

where

$\delta\eta_{\perp l}=\eta_{\perp}-\eta_{l}\;,\qquad
a=u^{\mu}\xi_{\mu}\;,\qquad b=l^{\mu}\xi_{\mu}\;,\qquad
v^{\mu}=\Delta^{\mu\nu}\xi_{\mu}\;.$ (22)

The second part which contains the six remaining roots reads as

$\displaystyle
H^{\parallel}_{\parallel}(V,\xi)\left[A^{2}+A\left(U_{1}^{\mu}l_{\mu}+U_{2}^{\mu}\xi_{\mu}\right)+U_{1}^{\mu}l_{\mu}U_{2}^{\nu}\xi_{\nu}-U_{1}^{\mu}\xi_{\mu}U_{2}^{\nu}l_{\nu}\right]=0\;.$
(23)

The explicit forms of the terms $H^{\parallel}_{\parallel}(V,\xi)$,
$U_{1}^{\mu}$, and $U_{2}^{\mu}$ are calculated explicitly in Appendix B,
which for convenience we listed below

$\displaystyle
U_{1}^{\mu}l_{\mu}=b^{2}\delta\eta_{lll\perp}+a^{2}\delta\lambda+\delta\eta_{\perp
l}v^{2}+b^{2}\delta\eta_{lll}-\frac{a^{2}b^{2}\delta\lambda(\chi_{\perp}+\lambda_{\perp}+\delta\lambda+\delta\chi)}{4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)}\;,$ (24a) $\displaystyle
U_{1}^{\mu}\xi_{\mu}=b^{3}\delta\eta_{lll\perp}+a^{2}b\delta\lambda+(\delta\eta_{\perp
l}+\delta\eta_{lll})bv^{2}-\frac{a^{2}b\delta\lambda[(\chi_{\perp}+\lambda_{\perp})v^{2}+b^{2}(\delta\lambda+\delta\chi)]}{4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)}\;,$ (24b) $\displaystyle
U_{2}^{\mu}l_{\mu}=b\left(\frac{\delta\chi}{3}+\delta\eta_{lll}\right)+\frac{(\chi_{\perp}-\eta_{ll})b}{3}-\frac{a^{2}b(\chi+\lambda_{\perp})(\chi_{\perp}+\lambda_{\perp}+\delta\lambda+\delta\chi)}{4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)}\;,$ (24c) $\displaystyle
U_{2}^{\mu}\xi_{\mu}=b^{2}\left(\frac{\delta\chi}{3}+\delta\eta_{lll}\right)+\frac{(\chi_{\perp}-\eta_{ll})v^{2}}{3}-\frac{a^{2}(\chi+\lambda_{\perp})[(\chi_{\perp}+\lambda_{\perp})v^{2}+b^{2}(\delta\lambda+\delta\chi)]}{4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)}\;,$ (24d) $\displaystyle
H^{\parallel}_{\parallel}(V,\xi)=\frac{3\chi
a^{2}+\lambda_{\perp}v^{2}+\delta\lambda\,b^{2}}{4\varepsilon}\;.$ (24e)

Here, $\delta\eta_{lll\perp}=4\eta_{l}-3\eta_{ll}-\eta_{\perp}$,
$\delta\eta_{lll}=\eta_{ll}-\eta_{l}$, $\delta\chi=\chi_{l}-\chi_{\perp}$, and
$\delta\lambda=\lambda_{l}-\lambda_{\perp}$. Because $\xi_{\mu}$ needs to be
spacelike according with Eq. (20), no root with $v=0$ is allowed. Note that
$v=0$ yields to the following result

$\xi_{\mu}\xi^{\mu}=(-u^{\mu}u^{\nu}+\Delta^{\mu\nu})\xi_{\mu}\xi_{\nu}=-a^{2}+v^{2}=-a^{2}<0\;,$

which clearly violates causality according with Eq. (20). Upon inspection of
Eqs. (21) and (23), it is apparent that if either $\lambda_{\perp}$ or $\chi$,
which give rise to the leading order power in $a$, are zero, causality
violating roots with $b=v=0$ arise. Building upon this observation and the
known isotropic conformal case Bemfica  _et al._ (2018), we can justify the
following choice to ensure causality, i.e.

$\lambda_{\perp}>0\;,\qquad\chi>0\;,\qquad\varepsilon>0\;.$ (25)

We can rewrite $A$ by writing $b=l_{\mu}v^{\mu}=v\cos\theta$
($0\leq\theta\leq\pi$) using the Cauchy-Schwarz inequality, as both vectors
are orthogonal to $u^{\mu}$. In this case, $\Delta^{\mu\nu}$ defines a real
inner product among these vectors. This results in the rewriting of $A$, Eq.
(21), as follows

$A=\lambda_{\perp}(a^{2}-\tau v^{2})\;,$ (26)

where

$\tau=\frac{\eta_{\perp}-\delta\eta_{\perp
l}\cos^{2}\theta}{\lambda_{\perp}}\;.$ (27)

When considering the constraints (25), the roots of $A=0$ obey (19) and (20)
if, and only if 666See for instance Ref. (Bemfica  _et al._ , 2022).,
$0\leq\tau<1$ 777The equality $\tau=1$ may be used if the particles are
massless.. Thus, if the conditions (25) hold, the roots of $A=0$ are causal
if, and only if,

$\displaystyle\lambda_{\perp}>\max(\eta_{l},\,\eta_{\perp})\geq 0\;.$ (28)

However, the previous inequality does not lead to a sufficient condition for
$\lambda_{\perp}$ when comparing the isotropic limit of this constraint, i.e.,
$\lambda>\eta$, and comparing with the corresponding one derived in the
conformal isotropic case Bemfica  _et al._ (2018). Therefore, it is needed to
analyze the roots of (23) as well. In general this analysis becomes
cumbersome, and as is shown in Appendix C, it requires to finding the roots of
the following polynomial in $\varrho\equiv a^{2}/v^{2}$

$p(\varrho)=\sum_{i=0}^{3}\alpha_{i}\varrho^{i}\;.$ (29)

This polynomial is found by writing Eq. (23) in the generic form given by Eq.
(73), with the coefficients $\alpha_{i}$ to be determined from the
aforementioned equality. The causality condition is then stated in the
statement C.1: Assuming $\alpha_{3}$ to be positive, causality then requires
$p(\varrho)$ to be positive for $\varrho\geq 1$ and negative for $\varrho<0$,
while satisfying the following inequality

$18\alpha_{0}\alpha_{1}\alpha_{2}\alpha_{3}-4\alpha_{2}^{3}\alpha_{0}+\alpha_{2}^{2}\alpha_{1}^{2}-4\alpha_{3}\alpha_{1}^{3}-27\alpha_{3}^{2}\alpha_{0}^{2}\geq
0\;.$ (30)

We have explained the application of the generic conditions of statement C.1.
To simplify the analysis of the causality conditions, it is more instructive
to consider specific cases that provide clearer results. From the constitutive
relations (11), we realize that the spatial anisotropy appears explicitly in
different dissipative fluxes of the energy-momentum tensor: in the scalar
sector with $\delta\chi\neq 0$, in the vector sector with $\delta\lambda\neq
0$, or in the tensor sector through the differences between 3 different
transport parameters, i.e, $\eta_{\perp}$, $\eta_{l}$, and $\eta_{ll}$. If the
only source of anisotropy is the scalar sector, then, as stated in statement
C.2, causality requires the 5 inequalities given by Eqs. (77), together with
Eqs. (13) and (28), to be satisfied. The aforementioned conditions (77) are,
for example, satisfied by the following choices

$\lambda=\chi=10\eta\;,\qquad\chi_{\perp}=\delta\chi=15\eta/2\;,$ (31)

with $\eta_{\perp}=\eta_{l}=\eta_{ll}$, and $\delta\lambda=0$. Taking the
isotropic limit (14), the conditions (77) reduce to the ones of the conformal
isotropic BDNK Bemfica  _et al._ (2018), i.e.,

$\displaystyle\chi>4\eta\;,\qquad\lambda>\frac{3\chi\eta}{\chi-\eta}\;.$ (32a)

On the other hand, if the anisotropy is only in the vector sector, and the
transport parameters of the tensor sector vanish, i.e., the system is
shearless, as stated in statement C.3, the condition (25) is sufficient for
causality, without any further constraint on $\lambda_{l}$.

To check for causality in more general cases, one can use statement C.1. For
instance, in the example worked out entirely in Appendix D one finds the
following conditions that respect causality

$\eta_{\perp}=\eta\,,\quad\eta_{l}=\frac{2\eta}{3}\,,\quad\eta_{ll}=\frac{5\eta}{6}\,,\quad\lambda_{\perp}=\frac{13\eta}{2}\,,\quad\lambda_{l}=6\eta\,,\quad\chi=5\eta\,,\quad\chi_{\perp}=\frac{11\eta}{2}\,,\quad\chi_{l}=\frac{16\eta}{3}\,,\mbox{\quad
with\quad}\eta>0\;.$ (33)

The proof of linear stability for this choice of parameters may be found in
Appendix E.

##  V Linear stability

In this section, we study the stability of the linearized EOM for the
hydrodynamic theory developed in Secs. II and III. We adopt standard methods
Hiscock and Lindblom (1985b); Bemfica  _et al._ (2019) and consider small
perturbations of the hydrodynamic fields $\varepsilon$ and $u^{\mu}$ around a
homogeneous background, which corresponds to a global equilibrium with
constant hydrodynamic fields in flat space-time. In particular, we assume the
energy density to be $\varepsilon_{0}+\delta\varepsilon(t,x^{i})$, with
$\varepsilon_{0}$ being constant, and the fluid’s velocity to be
$u^{\mu}_{0}+\delta u^{\mu}(t,x^{i})$, with $u^{\mu}_{0}\delta
u_{\mu}=\order{\delta^{2}}$. The fluid velocity in equilibrium is
$u^{\mu}_{0}=\gamma\left(1,v^{i}\right)$, with $v^{i}$ being the components of
the 3-velocity $\mathbf{v}$, and $\gamma=1/\sqrt{1-\mathbf{v}^{2}}$ the
Lorentz factor. We then expand the equations of motion up to the first order
in perturbations, which are assumed to be in the form of plane waves, i.e.,
$\delta\varepsilon(t,x^{i})\to
e^{-iT_{0}x^{\mu}k_{\mu}}\delta\varepsilon(k^{\mu})$ and $\delta u(t,x^{i})\to
e^{-iT_{0}x^{\mu}k_{\mu}}\delta u(k^{\mu})$ , where $k^{\mu}=(i\Gamma,k^{i})$
and the presence of the equilibrium temperature $T_{0}$ in the exponent makes
the modes $k^{\mu}$ dimensionless. Non-trivial solutions to the EOM may lead
to imaginary solutions $\Gamma=\Gamma(k^{i})$, where linear stability is
verified if, and only if, $\real(\Gamma)\leq 0$.

The roots $\Gamma=\Gamma(k^{i})$ that come from the EOM are usually
cumbersome, so the analysis becomes very difficult. However, new results show
that causality + linear stability in the local rest frame (LRF) leads to
linear stability in any boosted frame. This relation has been shown to be true
for strongly hyperbolic systems of equations, as demonstrated in Ref.(Bemfica
_et al._ , 2022). A recent study by Gavassino Gavassino (2022) demonstrated
the validity of this relation for the general case. The author showed that if
a mode grows for an observer A, it can be thought of as a parametrization of
the time coordinate in that frame. If the theory is causal, this growth should
be preserved between frames, otherwise there would be an inversion in the time
direction. Accordingly, and since we already have the conditions for
causality, we shall study linear stability in the LRF and, in some specific
cases, linear stability in a homogeneous boosted frame (where $k^{i}=0$ but
$v^{i}\neq 0$).

We begin by studying linear stability in the LRF, i.e., where
$u^{\mu}_{0}=(1,0,0,0)$ and, as a consequence, $\delta u^{0}=0$. Furthermore,
since $l$ is orthogonal to $u$, it can only have spatial components in the
LRF, i.e., $l^{\mu}_{0}=\left(0,l^{i}\right)$, with keeping $l^{\mu}l_{\mu}=1$
in mind. We note that although $l^{\mu}$ must also be perturbed in order to
preserve the aforementioned orthogonality, its perturbation does not
contribute to the dissipative fluxes of (11) up to the first order in
perturbations.

Let us define $\delta\bar{\varepsilon}=\delta\varepsilon/(\varepsilon+P)$ and
$\delta\bar{T}^{\mu\nu}=\delta T^{\mu\nu}/(\varepsilon+P)$, together with the
dimensionless quantities $\bar{\eta}_{l}=\eta_{l}/s$,
$\bar{\eta}_{ll}=\eta_{ll}/s$, $\bar{\eta}_{\perp}=\eta_{\perp}/s$,
$\bar{\lambda}_{l}=\lambda_{l}/s$, $\bar{\lambda}_{\perp}=\lambda_{\perp}/s$,
$\bar{\chi}=\chi/s$, $\bar{\chi}_{l}=\chi_{l}/s$, and
$\bar{\chi}_{\perp}=\chi_{\perp}/s$, where
$s=(\varepsilon+P)/T=4\varepsilon/3T$ is the equilibrium entropy density. We
may also define $\bar{k}=l_{i}k^{i}$,
$\kappa^{i}=\Xi^{i}_{j}k^{j}=k^{i}-l^{i}\bar{k}$, and
$\kappa=\sqrt{\kappa_{i}\kappa^{i}}$ in order to perform the decomposition

$\delta u^{i}=l^{i}\delta u_{L}+\frac{\kappa^{i}}{\kappa}\delta
u_{\parallel}+\delta u^{i}_{\perp},$

where $\delta u_{L}=l_{i}\delta u^{i}$, $\delta u_{\parallel}=\kappa_{i}\delta
u^{i}/\kappa$, and $\delta u_{\perp}^{i}=\Xi^{i}_{j}\delta
u^{j}-\kappa^{i}\delta u_{\parallel}$. With that in mind, we obtain the
following equations for the modes:

$\displaystyle\partial_{\nu}\delta\bar{T}^{0\nu}=M_{11}\delta\bar{\varepsilon}+M_{12}\delta
u_{L}+M_{13}\delta u_{\parallel}=0\;,$ (34a) $\displaystyle
l_{i}\partial_{\nu}\delta\bar{T}^{i\nu}=M_{21}\delta\bar{\varepsilon}+M_{22}\delta
u_{L}+M_{23}\delta u_{\parallel}=0\;,$ (34b)
$\displaystyle\frac{\kappa_{i}}{\kappa}\partial_{\nu}\delta\bar{T}^{i\nu}=M_{31}\delta\bar{\varepsilon}+M_{32}\delta
u_{L}+M_{33}\delta u_{\parallel}=0\;,$ (34c)
$\displaystyle\omega^{i}_{j}\partial_{\nu}\delta\bar{T}^{j\nu}=\left[\bar{\lambda}_{\perp}\Gamma^{2}+\Gamma+\bar{\eta}_{\perp}k^{2}+(\bar{\eta}_{l}-\bar{\eta}_{\perp})\bar{k}^{2}\right]\delta
u^{i}_{\perp}=0\;,$ (34d)

where $\omega^{ij}=\Xi^{ij}-\kappa^{i}\kappa^{j}/\kappa^{2}$ is the projector
orthogonal to $\kappa^{i}$ and $l^{i}$, $k^{2}=k_{i}k^{i}$, and

$\displaystyle M_{11}$ $\displaystyle=$
$\displaystyle\bar{\chi}\Gamma^{2}+\Gamma-\frac{k^{2}\bar{\lambda}_{\perp}}{3}+\frac{\bar{k}^{2}}{3}(\bar{\lambda}_{\perp}-\bar{\lambda}_{l})\;,$
(35a) $\displaystyle M_{12}$ $\displaystyle=$ $\displaystyle
i\bar{k}\left[\left(\bar{\lambda}_{l}+\bar{\chi}\right)\Gamma+1\right]\;,$
(35b) $\displaystyle M_{13}$ $\displaystyle=$ $\displaystyle
i\kappa\left[\left(\bar{\lambda}_{\perp}+\bar{\chi}\right)\Gamma+1\right]\;,$
(35c) $\displaystyle M_{21}$ $\displaystyle=$
$\displaystyle\frac{i}{3}\bar{k}\left[\left(\bar{\lambda}_{l}+\bar{\chi}_{l}\right)\Gamma+1\right]\;,$
(35d) $\displaystyle M_{22}$ $\displaystyle=$
$\displaystyle\bar{\lambda}_{l}\Gamma^{2}+\Gamma+\kappa^{2}\bar{\eta}_{l}+\frac{(4\bar{\eta}_{ll}-\bar{\chi}_{l})\bar{k}^{2}}{3}\;,$
(35e) $\displaystyle M_{23}$ $\displaystyle=$
$\displaystyle\frac{\bar{k}\kappa}{3}\left(3\bar{\eta}_{l}-\bar{\chi}_{l}-2\bar{\eta}_{ll}\right)\;,$
(35f) $\displaystyle M_{31}$ $\displaystyle=$
$\displaystyle\frac{i\kappa}{3}\left[\left(\bar{\chi}_{l}+\bar{\lambda}_{\perp}\right)\Gamma+1\right]\;,$
(35g) $\displaystyle M_{32}$ $\displaystyle=$
$\displaystyle\frac{\bar{k}\kappa}{3}\left(3\bar{\eta}_{l}-\bar{\chi}_{l}-2\bar{\eta}_{ll}\right)\;,$
(35h) $\displaystyle M_{33}$ $\displaystyle=$
$\displaystyle\Gamma^{2}\bar{\lambda}_{\perp}+\Gamma+\frac{k^{2}}{3}\left(-\bar{\chi}_{l}+\bar{\eta}_{ll}+3\bar{\eta}_{\perp}\right)+\frac{\bar{k}}{3}\left(3\bar{\eta}_{l}+\bar{\chi}_{l}-\bar{\eta}_{ll}-3\bar{\eta}_{\perp}\right)\;.$
(35i)

Note that the transverse modes $\delta u_{\perp}^{i}$ decouple from the rest
and give the shear channel polynomial equation

$\bar{\lambda}_{\perp}\Gamma^{2}+\Gamma+\bar{\eta}_{\perp}k^{2}+(\bar{\eta}_{l}-\bar{\eta}_{\perp})\bar{k}^{2}=0\;.$
(36)

The equations for the longitudinal modes (34a)–(34c) are nontrivial, i.e.,
give nonzero solutions for the perturbations, when

$\det(M)=0\;,$ (37)

where $M=[M_{ij}]_{3\times 3}$. The roots of Eq. (37) give the wave modes of
the sound channel. It is worth mentioning that the mode equations in a Lorentz
boosted frame are obtained by performing a boost which yields to the following
change

$\displaystyle\bar{k}\to l_{\mu}k^{\mu}\,,\quad
k^{2}\to-\gamma^{2}(\Gamma+ik_{i}v^{i})^{2}+\Gamma^{2}+k^{i}k_{i}\,,\quad\kappa^{2}\to
k^{2}-(l_{\mu}k^{\mu})^{2}\;,\quad\Gamma\to\gamma(\Gamma+ik_{i}v^{i})\;.$ (38)

In the boosted frame, $u^{\mu}=\gamma(1,v^{i})$ and
$l^{\mu}=(\tilde{l}^{k}v_{k},\tilde{l}^{i})/\sqrt{1-v^{l}\tilde{l}_{l}}$, with
$\tilde{l}^{i}$ being a unitary 3-vector that coincides with the anisotropic
unitary 3-vector $l^{i}$ in the LRF ($v=0$).

In the LRF case, by means of the Cauchy-Schwarz inequality, one may set
$l_{i}k^{i}=kx$, where $-1\leq x\leq 1$ covers all possible directions of
$k^{i}$. Then, Eq. (36) can be written as

$\bar{\lambda}_{\perp}\Gamma^{2}+\Gamma+\left[\bar{\eta}_{\perp}\left(1-x^{2}\right)+\bar{\eta}_{l}\right]k^{2}=0.$
(39)

The roots of the above equation have zero or negative real parts if, and only
if, $\bar{\lambda}>0$ and
$\bar{\eta}_{\perp}+x^{2}(\bar{\eta}_{l}-\bar{\eta}_{\perp})\geq 0$. Both
conditions are guaranteed by the constraints (25) and (28) because
$1-x^{2}\geq 0$, and we have that either $\bar{\eta}_{\perp}>\bar{\eta}\geq 0$
or $\bar{\eta}>\bar{\eta}_{\perp}\geq 0$. Hence, Eqs. (25) and (28) lead to
the linear stability of the shear channel in the LRF.

Now, let us consider (39) in the boosted homogeneous case, i.e., with
vanishing $k^{i}$. By employing again the Cauchy-Schwarz inequality, we can
write $(l_{\mu}k^{\mu})^{2}=(\Delta^{\mu\nu}k_{\mu}k_{\nu})y^{2}$, with
$-1\leq y\leq 1$, and then apply (38) to obtain

$\Gamma\left[\left(\bar{\lambda}_{\perp}-(1-y^{2})v^{2}\bar{\eta}_{\perp}-y^{2}v^{2}\bar{\eta}_{l}\right)\gamma\Gamma+1\right]=0\;,$
(40)

which has roots for all $v^{2}\in[0,1)$ with nonpositive real parts if
$\bar{\lambda}_{\perp}\geq\bar{\eta}_{\perp}\left(1-y^{2}\right)+\bar{\eta}_{l}y^{2}$
and $\bar{\lambda}_{\perp}>0$. Both inequalities are guaranteed by the
causality condition (28) and Eq. (25), according to which
$\bar{\eta}_{\perp}\left(1-y^{2}\right)+\bar{\eta}_{l}y^{2}\leq\max(\bar{\eta}_{l},\bar{\eta}_{\perp})<\bar{\lambda}_{\perp}$.
In this condition we clearly see the connection between stability in any frame
and causality+stability in the LRF Gavassino (2022).

The polynomial for the sound channel in (37) is of power six, with complicated
coefficients that depend on all transport parameters as well as $\bar{k}$ and
$k^{2}$. The analysis of this type of polynomials is extremely complex, so it
is more convenient and better to examine each set of parameters separately. As
an example, the stability of the causal set of parameters (33) is demonstrated
in Appendix E.

##  VI Bjorken flow

In this section, we apply the formalism developed in this work to the case of
Bjorken flow, i.e., when the fluid’s velocity is
$u^{\mu}=\left(1,0,0,0,\right)$ in the so-called Milne coordinates
$\left(\tau,x,y,\xi\right)$, which are related to the usual Minkowski
coordinates via $\tau=\sqrt{t^{2}-z^{2}}$ and
$\xi=\frac{1}{2}\,\log[(t+z)/(t-z)]$. From a mathematical perspective, one can
choose the spacelike anisotropic vector to be in any of the $x$, $y$, and
$\eta$ directions. However, the physical picture of the heavy ion collisions
suggests

$l^{\mu}=\frac{1}{\tau}\,\left(0,0,0,1\right)\,,$ (41)

to be the right choice. Alternatively, other choices may be possible, but one
should be cautious of potential non-physical issues when making such a choice,
as illustrated in Appendix F. By equating $u$ and $l$ into Eqs. (11) and using
the Bjorken symmetries, the dissipative fluxes reduce to

$\displaystyle\mathcal{E}^{(1)}$ $\displaystyle=$
$\displaystyle\tilde{\chi}^{3}T^{3}\left(\frac{1}{\tau}\,+\frac{3\dot{T}}{T}\right)\,,$
(42a) $\displaystyle\mathcal{P}_{l}^{(1)}$ $\displaystyle=$ $\displaystyle
T^{3}\left(\frac{\tilde{\chi}_{l}-4\tilde{\eta}_{ll}}{3\tau}+\frac{3\tilde{\chi}_{l}\dot{T}}{T}\right)\,,$
(42b) $\displaystyle\mathcal{P}_{\perp}^{(1)}$ $\displaystyle=$ $\displaystyle
T^{3}\left(\frac{\tilde{\chi}_{\perp}+2\tilde{\eta}_{ll}}{3\tau}+\frac{3\tilde{\chi}_{\perp}\dot{T}}{T}\right)\,,$
(42c)

and

$\pi_{\perp}^{\mu\nu}=0\,,\qquad W_{\perp l}^{\mu}=0\,,\qquad W_{\perp
u}^{\mu}=0\,,\qquad M=0\,.$ (43)

Namely, anisotropy only appears in the scalar sector with three independent
relevant transport parameters, $\chi=\tilde{\chi}T^{3}$,
$\chi_{\perp}=\tilde{\chi}_{\perp}T^{3}$, and
$\eta_{ll}=\tilde{\eta}_{ll}T^{3}$. Recall that in the isotropic case, the
relevant transport parameters for the Bjorken flow are shear viscosity $\eta$,
and $\chi$ Bemfica  _et al._ (2018); Shokri and Taghinavaz (2020). Similar to
the cases of IS Heller and Spalinski (2015) and isotropic conformal BDNK, the
fluid’s evolution is governed by only one equation,

$\displaystyle
9\tilde{\chi}\frac{\tau^{2}\ddot{T}}{T}+18\tilde{\chi}\frac{\tau^{2}\dot{T}^{2}}{T^{2}}+\left(\frac{3\tau(9\tilde{\chi}-\tilde{\chi}_{\perp})}{T}+12\tau^{2}\right)\dot{T}+4\tau
T+3\tilde{\chi}-2\tilde{\chi}_{\perp}-4\tilde{\eta}_{ll}=0\,.$ (44)

At late times, the solution to the previous equation can be written as a power
series,

$T=\frac{\Lambda}{(\Lambda\tau)^{1/3}}\left(1-\frac{\tilde{\eta}_{ll}}{2(\Lambda\tau)^{2/3}}-\frac{\tilde{\eta}_{ll}(\tilde{\chi}_{l}+5\tilde{\chi}_{\perp})}{24(\Lambda\tau)^{2/3}}+\cdots\right)\,,$
(45)

where $\Lambda$ is a constant with energy dimensions. We clarify that when we
refer to late times, we mean $\Lambda\tau\gg 1$. One may assume that up to the
first-order, the power series solution is equal to the one of the isotropic
conformal case, to reduce the number of free parameters. Such an assumption
gives rise to $\eta_{ll}=\eta$.

To gain more insight into the physical implications of (44), especially in the
far from equilibrium regime, we assume the dimensionless parameters $w=T\tau$
888The variable $w$ is proportional to the inverse Knudsen number for the
conformal case. and $f(w)=\frac{\tau}{w}\derivative{w}{\tau}$ Heller and
Spalinski (2015). The EOM (44) reduces to a first-order nonlinear differential
equation,

$\displaystyle\frac{9\tilde{\chi}}{4}f(w)^{2}+wf(w)\left(1+\frac{3}{4}\tilde{\chi}f^{\prime}(w)\right)-\frac{6\tilde{\chi}+\tilde{\chi}_{\perp}}{2}f(w)+\frac{3\tilde{\chi}+\tilde{\chi}_{\perp}-\tilde{\eta}_{ll}}{3}-\frac{2w}{3}=0\,.$
(46)

The late-time expansion (45) is written in terms of the variable $w$ as

$f(w)=\frac{2}{3}+\frac{\tilde{\eta}_{ll}}{3w}+\frac{\tilde{\eta}_{ll(\tilde{\chi}_{l}+5\tilde{\chi}_{\perp})}}{18w^{2}}+\order{\frac{1}{w^{3}}\,}\,,$
(47)

which is valid when $w\gg 1$. The pressure anisotropy,
${\mathcal{A}}=\left({\mathcal{P}}_{\perp}-{\mathcal{P}}_{l}\right)/P$, in the
isotropic conformal BDNK, and in contrast to IS theory, is purely determined
by shear viscosity and does not depend on $f$. In the anisotropic case, the
situation is different because ${\mathcal{A}}$ receives contribution from
$f(w)$,

${\mathcal{A}}=\frac{2}{3}\frac{\tilde{\chi}_{\perp}-\tilde{\chi}_{l}}{w}\left(f-\frac{2}{3}\right)+\frac{6\tilde{\eta}_{ll}}{w}\,.$
(48)

Note that at late times $f(w)\sim 2/3+\order{w^{-1}}$ so both, isotropic and
anisotropic first order BDNK theories share the same forward attractor. This
is expected since asymptotically the system relaxes towards the thermal
equilibrium. Following Ref. Heller and Spalinski (2015), we assume a
correction to the first-order on-shell terms in (47),

$f(w)=\frac{2}{3}+\frac{\tilde{\eta}_{ll}}{3w}+\delta f(w)\,.$ (49)

By equating the previous expression into Eq. (46) while assuming $\delta f\ll
f$, and expanding in $1/w$ around $w\to\infty$, we obtain

$\delta
f(w)\sim\exp(-\frac{2w}{\tilde{\chi}})w^{\frac{\tilde{\eta}_{ll}+\tilde{\chi}_{\perp}}{\tilde{\chi}}}\,.$
(50)

Since $\chi>0$ due to causality, at late times the perturbation (50) decays
faster than the perturbative terms of the late-time expansion. The relation
between the coefficients in the above form and the analytical structure of the
Borel resummation might be of interest, but it will not be discussed here.

Equation (46) can also be studied for the existence of pullback attractors.
Following the transasymptotic and dynamical systems methods outlined in Refs.
Behtash  _et al._ (2018, 2019a, 2019b, 2021) one can show that the initial
value which rise to the attractor solution is found by expanding (46) around
$w=0$,

$f(w\ll
1)=\frac{7}{9}-\frac{\tilde{\chi}-\tilde{\chi}_{\perp}}{9\tilde{\chi}}+\frac{\sqrt{12\tilde{\eta}_{ll}\tilde{\chi}+\tilde{\chi}_{\perp}^{2}}}{9\tilde{\chi}}\,.$
(51)

In the following we use the slow roll approximation Heller and Spalinski
(2015). Namely, we assume $\absolutevalue{f^{\prime}}$ to be much smaller than
$\absolutevalue{f}$ and expand (46) in terms of $f^{\prime}/f$ to obtain an
algebraic equation. This leads to two solutions for $f(w)$ that we expand
around $w\to\infty$ and compare with (47), to recognize the stable one,

$f(w)_{\rm
slowroll}=\frac{7}{9}-\frac{\tilde{\chi}-\tilde{\chi}_{\perp}}{9\tilde{\chi}}-\frac{2w}{9\tilde{\chi}}+\frac{\sqrt{\left(2w-\tilde{\chi}_{\perp}\right)^{2}+12\tilde{\eta}_{ll}\tilde{\chi}}}{9\tilde{\chi}}\,.$
(52)

If at early times, $f$ gets larger than 1 then the fluid experiences
reheating. To prevent reheating we must have

$\chi_{l}>4\eta_{ll}>0\,,$ (53)

which can be shown that is equivalent to the condition for stability and
causality and reduces to $\chi>4\eta$ in the isotropic limit. Therefore,
causality forbids reheating for the Bjorken model discussed in this section.

Finally, we might assume an off-shell, or nonphysical, entropy current
$S^{\mu}_{\rm off}$, that is, a current determined from (4) without
considering the power-counting discussed in Sec. III. If causal and stable
parameters are chosen for the attractor solution (52), $\nabla\cdot S_{\rm
off}$ will be negative at early times. However, if causality and stability are
violated, the divergence will always be positive. $\nabla\cdot S_{\rm off}$
approaches zero rapidly and then changes sign in the stable and causal case.

##  VII Conclusions

In this work, we explore the constraints imposed by causality in anisotropic
hydrodynamics. Our approach involves deriving the simplest and most general
form of the first-order anisotropic energy-momentum tensor, based on the
invariance under the little group $SO(2)$ and the second law of
thermodynamics. We demonstrate that our theory is both nonlinearly causal
whether or not it is coupled to gravity and that it is stable in the
linearized regime. It is worth mentioning that causality is not only in the
heart of any relativistic theory but also is essential for well-posedness of a
system of relativistic covariant equations since it must be valid in any
Lorentz frame. Furthermore, we show that the standard isotropic BDNK theory
can be recovered as a limit of our novel approach.

We verify the linear stability in the local rest frame (LRF) and use causality
to apply the recent result obtained by Gavassino Gavassino (2022), which
guarantees that linear stability in the LRF together with causality results in
linear stability in any boosted frame. We give an illustration of such
stability-causality connection by verifying these results in the homogeneous
boosted frame. For the specific case of Bjorken flow, we investigate the
causality conditions necessary to ensure the existence of forward and pullback
attractors. Our findings reveal that the behavior of these attractors at early
and late times is constrained by causality conditions. Violation of these
conditions results in a reheating effect, which is the increase of temperature
from its initial value at very early times.

There are several potential avenues for future research that can build on our
work. Firstly, it would be valuable to investigate how our novel anisotropic
theory can be derived from a coarse-graining approach, as is typically done in
relativistic kinetic theory. Additionally, extending our analysis to the most
general non-conformal case for both charged and uncharged fluids, with and
without gravity, could have significant implications in other fields such as
cosmology Bemfica  _et al._ (2023). Our work may also impact the development
of a BDNK-type theory for resistive dissipative magnetohydrodynamics Denicol
_et al._ (2019); Armas and Camilloni (2022), where the magnetic field serves a
similar purpose to the anisotropy vector $l^{\mu}$.

Finally, while our approach has been successful in analyzing the near-
equilibrium regime, it would be valuable to extend our techniques to the
extreme far-from-equilibrium dynamics. This could provide valuable insights
into how causality affects the space-time evolution of these systems. These
fascinating research problems are important and require further investigation
in the future.

###### Acknowledgements.

M. S. was supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through the Collaborative Research Center CRC-TR 211
“Strong-interaction matter under extreme conditions” - project number
315477589 - TRR 211 and by the State of Hesse within the Research Cluster
ELEMENTS (Project ID 500/10.006). M. M. was supported in part by the US
Department of Energy Grant No. DE-FG02-03ER41260 and BEST (Beam Energy Scan
Theory) DOE Topical Collaboration.

##  Appendix A Anisotropic free energy density

In this appendix, we present the thermodynamic relations of an anisotropic
fluid. To this end, let us assume the leading-order generating functional of
the fluid to be Jensen  _et al._ (2012); Banerjee  _et al._ (2012)

$W_{0}=\int\differential[4]{x}\sqrt{-g}{\mathcal{F}}(T,\ell)\,.$ (54)

Here, ${\mathcal{F}}$ is the free energy that is to be determined and
$\ell=\sqrt{g_{\mu\nu}l^{\mu}l^{\nu}}$, which is set to unity at the end
999$\ell$ is similar to the parameter $\xi$ for the superfluid case in Jensen
_et al._ (2012). Note that, $\ell$ varies only because of the metric
variation.. The energy-momentum tensor can then be derived from

$T^{\mu\nu}=\frac{2}{\sqrt{-g}}\functionalderivative{W}{g_{\mu\nu}}=\Bigg{(}{\mathcal{F}}\delta\sqrt{-g}+\sqrt{-g}\left(\frac{\partial{{\mathcal{F}}}}{\partial{T}}\right)_{\ell}\delta
T+\sqrt{-g}\left(\frac{\partial{{\mathcal{F}}}}{\partial{\ell}}\right)_{T}\delta\ell\Bigg{)}\,.$
(55)

The variation of temperature and $\ell$ are

$\delta T=\frac{1}{2}\,Tu^{\mu}u^{\nu}\delta
g_{\mu\nu}\,,\qquad\delta\ell=\frac{1}{2}\,l^{\mu}l^{\nu}\delta g_{\mu\nu}\,.$

Plugging the above relations into (55), we obtain

$T^{\mu\nu}_{0}=\left(-{\mathcal{F}}+T\left(\frac{\partial{{\mathcal{F}}}}{\partial{T}}\right)\Big{|}_{\ell}\right)u^{\mu}u^{\nu}+\left(\frac{\partial{{\mathcal{F}}}}{\partial{\ell}}\right)\Big{|}_{T}l^{\mu}l^{\nu}+{\mathcal{F}}\Delta^{\mu\nu}\,.$
(56)

The pressure and longitudinal pressure are defined as

$\displaystyle
P\equiv\frac{1}{3}\,\Delta_{\mu\nu}T^{\mu\nu}={\mathcal{F}}+\frac{1}{3}\,\left(\frac{\partial{{\mathcal{F}}}}{\partial{\ell}}\right)_{T}\,,\qquad
P_{l}\equiv
l_{\mu}l_{\nu}T^{\mu\nu}={\mathcal{F}}+\left(\frac{\partial{{\mathcal{F}}}}{\partial{\ell}}\right)_{T}\,.$

Plugging (56) into the above, we recognize the free energy,

${\mathcal{F}}=\frac{3P-P_{l}}{2}\equiv P_{\perp}\,.$ (57)

After identifying the free energy with the transverse pressure, one can
observe the relations between thermodynamic quantities,

$\displaystyle\epsilon\equiv
u_{\mu}u_{\nu}T^{\mu\nu}=-P_{\perp}+T\left(\frac{\partial{P_{\perp}}}{\partial{T}}\right)_{\ell}\,,\qquad\left(\frac{\partial{P_{\perp}}}{\partial{\ell}}\right)_{T}=P_{l}-P_{\perp}\,.$
(58)

The first equation above is Euler’s equation for the anisotropic fluid that
determines its entropy density,

$s=\frac{\epsilon+{\mathcal{F}}}{T}=\left(\frac{\partial{P_{\perp}}}{\partial{T}}\right)_{\ell}\,.$
(59)

From the above definition of the entropy density and (58), we find

$\displaystyle\differential{\epsilon}$ $\displaystyle=$
$\displaystyle-\left(P_{l}-P_{\perp}\right)\differential{\ell}+T\differential{s}\,\qquad\differential{P_{\perp}}=\left(P_{l}-P_{\perp}\right)\differential{\ell}+s\differential{T}\,.$

At this point, we set $\ell=1$ to find the Gibbs-Duhem relation and the first
law of thermodynamics for the anisotropic fluid

$\displaystyle\differential{\epsilon}$ $\displaystyle=$ $\displaystyle
T\differential{s}\,\qquad\differential{P_{\perp}}=s\differential{T}\,.$ (60)

Finally, the entropy current is found from the covariant form of (59)

$Ts_{\mu}=P_{\perp}u_{\mu}-u^{\nu}T_{\mu\nu}\,.$ (61)

##  Appendix B Details of causality analysis

In this appendix we give the mathematical details of causality analysis that
are not presented in Sec. IV. The system of equations arising from (17) and
(18) can be written as 101010One may notice that causality is blind to the
leading order anisotropy of pressure since it does not appear in the principal
part of the system of PDEs.

$\displaystyle\frac{3\chi
u^{\alpha}u^{\beta}+\lambda_{\perp}\Delta^{\alpha\beta}+\delta\lambda\,l^{\alpha}l^{\beta}}{4\varepsilon}\partial_{\alpha}\partial_{\beta}\varepsilon+\left[(\chi+\lambda_{\perp})u^{(\alpha}\delta^{\beta)}_{\nu}+\delta\lambda\,l^{(\alpha}u^{\beta)}l_{\nu}\right]\partial_{\alpha}\partial_{\beta}u^{\nu}$
$\displaystyle\hskip
113.81102pt+h_{a}^{\parallel,\,\alpha\beta}(\varepsilon,u,g)\partial_{\alpha}\partial_{\beta}g_{a}=b^{\parallel}(\partial\varepsilon,\partial
u,\partial g)\;,$ (62a)
$\displaystyle\frac{(\chi_{\perp}+\lambda_{\perp})\Delta^{\mu(\alpha}u^{\beta)}+(\delta\lambda+\delta\chi)\,l^{\mu}l^{(\alpha}u^{\beta)}}{4\varepsilon}\partial_{\alpha}\partial_{\beta}\varepsilon+C^{\mu\alpha\beta}_{\nu}\partial_{\alpha}\partial_{\beta}u^{\nu}+h_{a}^{\mu,\,\alpha\beta}(\varepsilon,u,g)\partial_{\alpha}\partial_{\beta}g_{a}$
$\displaystyle\hskip 113.81102pt=b^{\mu}(\partial\varepsilon,\partial
u,\partial g)\;,$ (62b) $\displaystyle
g^{\alpha\beta}\partial_{\alpha}\partial_{\beta}g_{a}=b_{a}(\partial\varepsilon,\partial
u,\partial g)\;,$ (62c)

where there is an implicit sum over repeated $a,b=\mu\nu$, with $\mu\leq\nu$,
over the values $00,\,01,\,02,\,03,\,11,\,12,\,13,\,22,\,23,$ and $33$.
Furthermore, $\delta\lambda=\lambda_{l}-\lambda_{\perp}$,
$\delta\chi=\chi_{l}-\chi_{\perp}$, and we have defined

$\displaystyle
C^{\mu\alpha\beta}_{\nu}\equiv\delta\eta_{lll\perp}\,l^{\mu}l^{\beta}l^{\alpha}l_{\nu}+\left(\frac{\delta\chi}{3}+\delta\eta_{lll}\right)\,l^{\mu}l^{(\beta}\delta^{\alpha)}_{\nu}+\delta\eta_{lll}\,\Delta^{\mu(\beta}l^{\alpha)}l_{\nu}+\frac{(\chi_{\perp}-\eta_{ll})\Delta^{\mu(\beta}\delta^{\alpha)}_{\nu}}{3}$
$\displaystyle\hskip
56.9055pt+\delta\lambda\,l^{\mu}l_{\nu}u^{\alpha}u^{\beta}+\delta\eta_{\perp
l}\Delta^{\alpha\beta}l^{\mu}l_{\nu}+\left[\lambda_{\perp}u^{\alpha}u^{\beta}-\eta_{\perp}\Delta^{\alpha\beta}+\delta\eta_{\perp
l}l^{\alpha}l^{\beta}\right]\delta^{\mu}_{\nu}\;,$ (63)

wherein $\delta\eta_{lll\perp}=4\eta_{l}-3\eta_{ll}-\eta_{\perp}$,
$\delta\eta_{lll}=\eta_{ll}-\eta_{l}$ and $\delta\eta_{\perp
l}=\eta_{\perp}-\eta_{l}$. We may rewrite (62) in the matrix form

$H^{I}_{J}(V,\partial)V^{J}+b^{I}=0,$ (64)

where $I,J$ take the values $\parallel$, $\mu=0,1,2,3$, and
$a=00,01,\cdots,33$, $V^{I}=(\varepsilon,u^{\nu},g_{b})$,

$H(V,\partial)^{I}_{J}=h(V)^{I,\,\alpha\beta}_{J}\partial_{\alpha}\partial_{\beta}$

is a $15\times 15$ matrix linear operator, and

$\displaystyle h^{\parallel,\,\alpha\beta}_{\parallel}$ $\displaystyle=$
$\displaystyle\frac{3\chi
u^{\alpha}u^{\beta}+\lambda_{\perp}\Delta^{\alpha\beta}+\delta\lambda\,l^{\alpha}l^{\beta}}{4\varepsilon},$
(65a) $\displaystyle h^{\parallel,\,\alpha\beta}_{\nu}$ $\displaystyle=$
$\displaystyle(\chi+\lambda_{\perp})u^{(\alpha}\delta^{\beta)}_{\nu}+\delta\lambda\,l^{(\alpha}u^{\beta)}l_{\nu},$
(65b) $\displaystyle h^{\mu,\,\alpha\beta}_{\parallel}$ $\displaystyle=$
$\displaystyle\frac{(\chi_{\perp}+\lambda_{\perp})\Delta^{\mu(\alpha}u^{\beta)}+(\delta\lambda+\delta\chi)\,l^{\mu}l^{(\alpha}u^{\beta)}}{4\varepsilon},$
(65c) $\displaystyle h^{\mu,\,\alpha\beta}_{\nu}$ $\displaystyle=$
$\displaystyle C^{\mu\alpha\beta}_{\nu},$ (65d) $\displaystyle
h^{a,\,\alpha\beta}_{b}$ $\displaystyle=$
$\displaystyle\delta^{a}_{b}g^{\alpha\beta},\text{ where
}\delta^{a}_{b}=\delta^{\mu}_{\lambda}\delta^{\nu}_{\sigma}\text{ when
}a=\mu\nu\text{ and }b=\lambda\sigma,$ (65e) $\displaystyle
h^{a,\,\alpha\beta}_{\parallel}$ $\displaystyle=$ $\displaystyle
h^{a,\,\alpha\beta}_{\mu}=0.$ (65f)

The remaining expressions $h^{\mu,\,\alpha\beta}_{b}(V)$ and $b^{I}(\partial
V)$ are irrelevant to what follows. The principal part of each equation $I$ is
contained in $H(V,\partial)^{I}_{J}V^{J}$. It is worth mentioning that all
terms $b^{I}(\partial V)$ are functions of at most the first-order derivative
of the variables $V^{I}=(\varepsilon,u^{\nu},g_{a})$, with products containing
first-order derivatives such as $(\partial V^{K})^{n}$, $(\partial
V^{K})^{n}(\partial V^{L})^{m}$, among others, being allowed. On the other
hand, $h^{I,\,\alpha\beta}_{J}(V)$ are functions of the variables $V^{I}$ only
and not their derivatives. Therefore, this is a quasi-linear PDE system and
the usual tools to compute causality apply. The characteristic surfaces
$\\{\Phi(x)=0\\}$ are determined by the principal part of the equations by
solving the characteristic equation $\det[H(V^{k},\xi)]=0$, with
$\xi_{\mu}=\nabla_{\mu}\Phi$ Courant and Hilbert (1991); Choquet-Bruhat
(2009). Note that the components of the matrix $H(V,\xi)$ are
$H^{I}_{J}(V,\xi)=h(V)_{J}^{I,\alpha\beta}\xi_{\alpha}\xi_{\beta}$. The system
is causal if, for any given real $\xi_{i}$, (C1) the roots
$\xi_{0}=\xi_{0}(\xi_{i})$ of the characteristic equation are real and (C2)
$\xi_{\alpha}=(\xi_{0}(\xi_{i}),\xi_{i})$ is spacelike or lightlike, i.e.,
$\xi_{\mu}\xi^{\mu}\geq 0$. Condition (C2) guarantees that the hypersurfaces
$\\{\Phi(x)=0\\}$ are timelike or lightlike, ensuring that there is no
superluminal information. 111111Note that the matrix $H(V^{k},\xi)$ is
invertible only when $\xi$ is timelike, i.e., solutions are only possible over
spacelike or lightlike hypersurfaces $\Phi$.

We may now compute the characteristic equation, for which we must compute the
determinant of the matrix $H(V,\xi)$ that reads

$\displaystyle\det[H(V,\xi)]$ $\displaystyle=$
$\displaystyle\det\begin{bmatrix}H^{\parallel}_{\parallel}(V,\xi)&H^{\parallel}_{\nu}(V,\xi)&H^{\parallel}_{b}(V,\xi)\\\
H^{\mu}_{\parallel}(V,\xi)&H_{\nu}^{\mu}(V,\xi)&H_{b}^{\mu}(V,\xi)\\\
0_{10\times 1}&0_{10\times 4}&\xi_{\mu}\xi^{\mu}I_{10}\end{bmatrix}$ (66)
$\displaystyle=$ $\displaystyle(\xi_{\mu}\xi^{\mu})^{10}M\;,$

where

$M=\det\begin{bmatrix}H^{\parallel}_{\parallel}(V,\xi)&H^{\parallel}_{\nu}(V,\xi)\\\
H^{\mu}_{\parallel}(V,\xi)&H_{\nu}^{\mu}(V,\xi)\end{bmatrix},$ (67)

wherein $I_{10}$ is the $10\times 10$ identity matrix and $0_{m\times n}$ is
the $m\times n$ null matrix. Out of 30 roots of the characteristic equation,
20 are the real roots $\xi_{\mu}\xi^{\mu}=0$ coming from the pure gravity
sector and, thus, are causal. These are expected to be lightlike roots since
pure gravity fields have no mass. If matter instead of pure radiation, which
is permitted for a conformal fluid, is treated, the remaining roots are
expected to be spacelike. As for the matter sector, let us define
$a=u^{\mu}\xi_{\mu}$, $b=l^{\mu}\xi_{\mu}$,
$v^{\mu}=\Delta^{\mu\nu}\xi_{\mu}$, and
$v=\sqrt{\xi_{\mu}\xi_{\nu}\Delta^{\mu\nu}}$. Then, from (67) one obtains that

$\displaystyle
M=H^{\parallel}_{\parallel}(V,\xi)\det\begin{bmatrix}H_{\nu}^{\mu}(V,\xi)-\frac{H^{\mu}_{\parallel}(V,\xi)H^{\parallel}_{\nu}(V,\xi)}{H^{\parallel}_{\parallel}(V,\xi)}\\\
\end{bmatrix}$
$\displaystyle=H^{\parallel}_{\parallel}(V,\xi)\det\left[A\delta^{\mu}_{\nu}+U_{1}^{\mu}l_{\nu}+U_{2}^{\mu}\xi_{\nu}\right]$
$\displaystyle=A^{2}H^{\parallel}_{\parallel}(V,\xi)\Big{[}A^{2}+A\left(U_{1}^{\mu}l_{\mu}+U_{2}^{\mu}\xi_{\mu}\right)+U_{1}^{\mu}l_{\mu}U_{2}^{\nu}\xi_{\nu}-U_{1}^{\mu}\xi_{\mu}U_{2}^{\nu}l_{\nu}\Big{]}\;,$
(68)

where

$\displaystyle A$ $\displaystyle=$
$\displaystyle\lambda_{\perp}a^{2}-\eta_{\perp}v^{2}+\delta\eta_{\perp
l}b^{2}\;,$ (69a) $\displaystyle U_{1}^{\mu}$ $\displaystyle=$ $\displaystyle
b^{2}\delta\eta_{lll\perp}\,l^{\mu}+a^{2}\delta\lambda\,l^{\mu}+\delta\eta_{\perp
l}v^{2}l^{\mu}+b\delta\eta_{lll}\,v^{\mu}-\frac{ab\delta\lambda\,H^{\mu}_{\parallel}(V,\xi)}{H^{\parallel}_{\parallel}(V,\xi)}\;,$
(69b) $\displaystyle U_{2}^{\mu}$ $\displaystyle=$ $\displaystyle
b\left(\frac{\delta\chi}{3}+\delta\eta_{lll}\right)\,l^{\mu}+\frac{(\chi_{\perp}-\eta_{ll})v^{\mu}}{3}-\frac{a(\chi+\lambda_{\perp})H^{\mu}_{\parallel}(V,\xi)}{H^{\parallel}_{\parallel}(V,\xi)}\;,$
(69c)

In particular, let us write explicitly

$\displaystyle
U_{1}^{\mu}l_{\mu}=b^{2}\delta\eta_{lll\perp}+a^{2}\delta\lambda+\delta\eta_{\perp
l}v^{2}+b^{2}\delta\eta_{lll}-\frac{a^{2}b^{2}\delta\lambda(\chi_{\perp}+\lambda_{\perp}+\delta\lambda+\delta\chi)}{4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)}\;,$ (70a) $\displaystyle
U_{1}^{\mu}\xi_{\mu}=b^{3}\delta\eta_{lll\perp}+a^{2}b\delta\lambda+(\delta\eta_{\perp
l}+\delta\eta_{lll})bv^{2}-\frac{a^{2}b\delta\lambda[(\chi_{\perp}+\lambda_{\perp})v^{2}+b^{2}(\delta\lambda+\delta\chi)]}{4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)}\;,$ (70b) $\displaystyle
U_{2}^{\mu}l_{\mu}=b\left(\frac{\delta\chi}{3}+\delta\eta_{lll}\right)+\frac{(\chi_{\perp}-\eta_{ll})b}{3}-\frac{a^{2}b(\chi+\lambda_{\perp})(\chi_{\perp}+\lambda_{\perp}+\delta\lambda+\delta\chi)}{4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)}\;,$ (70c) $\displaystyle
U_{2}^{\mu}\xi_{\mu}=b^{2}\left(\frac{\delta\chi}{3}+\delta\eta_{lll}\right)+\frac{(\chi_{\perp}-\eta_{ll})v^{2}}{3}-\frac{a^{2}(\chi+\lambda_{\perp})[(\chi_{\perp}+\lambda_{\perp})v^{2}+b^{2}(\delta\lambda+\delta\chi)]}{4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)}\;,$ (70d) $\displaystyle
H^{\parallel}_{\parallel}(V,\xi)=\frac{3\chi
a^{2}+\lambda_{\perp}v^{2}+\delta\lambda\,b^{2}}{4\varepsilon}\;.$ (70e)

Thus, the matter sector contains 10 overall roots that must obey (C1) and
(C2), with 4 coming from

$A^{2}=\left(\lambda_{\perp}a^{2}-\eta_{\perp}v^{2}+\delta\eta_{\perp
l}b^{2}\right)^{2}=0\;,$ (71)

and the remaining 6 from 121212Note that in the product
$H^{\parallel}_{\parallel}(V,\xi)\left(U_{1}^{\mu}l_{\mu}U_{2}^{\nu}\xi_{\nu}-U_{1}^{\mu}\xi_{\mu}U_{2}^{\nu}l_{\nu}\right)$,
the term with denominator $H^{\parallel}_{\parallel}(V,\xi)$ cancels as
expected.

$\displaystyle
H^{\parallel}_{\parallel}(V,\xi)\left[A^{2}+A\left(U_{1}^{\mu}l_{\mu}+U_{2}^{\mu}\xi_{\mu}\right)+U_{1}^{\mu}l_{\mu}U_{2}^{\nu}\xi_{\nu}-U_{1}^{\mu}\xi_{\mu}U_{2}^{\nu}l_{\nu}\right]=0\;.$
(72)

##  Appendix C Causility conditions in the general and specific cases

The polynomial in (72) is of power 3 in $a^{2}$. Since for causality we must
obtain the roots in the form of $\varrho=\frac{a^{2}}{v^{2}}$, then we may
rewrite (72) as

$4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)\left[A^{2}+A\left(U_{1}^{\mu}l_{\mu}+U_{2}^{\mu}\xi_{\mu}\right)+U_{1}^{\mu}l_{\mu}U_{2}^{\nu}\xi_{\nu}-U_{1}^{\mu}\xi_{\mu}U_{2}^{\nu}l_{\nu}\right]=p(\varrho)v^{6};.$
(73)

We have multiplied (72) by $4\varepsilon$ to eliminate it from the denominator
in $H^{\parallel}_{\parallel}(V,\xi)$. Note that $p(\varrho)$ is a cubic
polynomial in $\varrho$. Since $l^{\mu}=\Delta^{\mu\nu}l_{\nu}$ and
$v^{\mu}=\Delta^{\mu\nu}\xi_{\nu}$ are vectors orthogonal to $u^{\mu}$,
$\Delta^{\mu\nu}$ define an inner product between them and, thus, we can apply
the Cauchy-Schwarz inequality to write $b=l_{\mu}v^{\mu}=\kappa\,v$, where
$\kappa\in[-1,1]$ depending on the root $\xi$ and the vector $l$. Causality of
the 6 roots of (73) follows from the statement:

######  Statement C.1 (The general case).

Let $p(\varrho)$ be defined by means of Eq. (73) and let us write it as

$p(\varrho)=\alpha_{3}\varrho^{3}+\alpha_{2}\varrho^{2}+\alpha_{1}\varrho+\alpha_{0}.$
(74)

Assume that

$\alpha_{3}>0\;.$ (75)

Then, causality requires that

$\displaystyle p(\varrho)>0,\quad\forall\varrho\geq 1$ (76a) $\displaystyle
p(\varrho)<0,\quad\forall\varrho<0,$ (76b) $\displaystyle
18\alpha_{0}\alpha_{1}\alpha_{2}\alpha_{3}-4\alpha_{2}^{3}\alpha_{0}+\alpha_{2}^{2}\alpha_{1}^{2}-4\alpha_{3}\alpha_{1}^{3}-27\alpha_{3}^{2}\alpha_{0}^{2}\geq
0,$ (76c) for all $\kappa\in[-1,1]$.

###### Proof.

First, we must ensure that all real roots of $p(\varrho)$ lie on the range
$[0,1)$ as demanded by (20). From (75), we must ensure that $p(\varrho)$ is
positive for all $\varrho\geq 1$ [condition (76a)] and negative for all
$\varrho<0$ [condition (76b). This guarantees that the real roots can only
occur in this desired range, and thus (20) is satisfied. As for (19), the
roots are real if the discriminant of the cubic polynomial (74) is greater or
equal to zero, which leads to the condition (76c). ∎

In what follows, we present causality conditions for two more specific cases.
The first case is when the anisotropy appears only in ${\mathcal{E}}^{(1)}$,
${\mathcal{P}}^{(1)}_{l}$, and ${\mathcal{P}}^{(1)}$:

######  Statement C.2 (Anisotropy in the $\chi$’s).

Consider the anisotropic conformal fluid theory defined by the energy-momentum
tensor in (1) and supplemented with Eqs. (11) with the choices
$\lambda_{\perp}=\lambda_{l}=\lambda,$ and
$\eta_{\perp}=\eta_{l}=\eta_{ll}=\eta$. Then, the corresponding EOM are causal
under assumption (25) if, and only if, condition (28) applies together with

$\displaystyle\kappa^{2}\delta\chi\lambda+\chi(4\eta+\lambda-\chi)+(\lambda+\chi)\chi_{\perp}\geq
0$ (77a)
$\displaystyle\kappa^{4}\delta\chi^{2}\lambda^{2}-2\kappa^{2}\lambda\chi\delta\chi(-4\eta+\lambda+\chi)+2(\lambda+\chi)\chi_{\perp}\left[\kappa^{2}\delta\chi\lambda+\chi(4\eta+\lambda-\chi)\right]$
$\displaystyle+\chi\left[16\eta^{2}\chi+8\eta\left(2\lambda^{2}+\lambda\chi-\chi^{2}\right)+\chi\left(-3\lambda^{2}-2\lambda\chi+\chi^{2}\right)\right]+(\lambda+\chi)^{2}\chi_{\perp}^{2}\geq
0,$ (77b) $\displaystyle\chi\geq 4\eta-\kappa^{2}\delta\chi,$ (77c)
$\displaystyle-\kappa^{2}\lambda\delta\chi+\chi(-4\eta+5\lambda+\chi)-(\lambda+\chi)\chi_{\perp}>0,$
(77d)
$\displaystyle-2\kappa^{2}\delta\chi\lambda-4\eta\lambda-12\eta\chi+7\lambda\chi-3(\lambda+\chi)\chi_{\perp}+3\chi^{2}>0,$
(77e)

for all $\kappa^{2}\in[0,1]$.

###### Proof.

The determinant $M$ in (B) can be rewritten as

$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{3\chi\lambda^{4}}{4\varepsilon}\prod_{a=1,\pm}(a^{2}-\tau_{a}v^{2})^{n_{a}},$
(78)

where $n_{1}=3$, $n\pm=1$, and

$\displaystyle\tau_{1}$ $\displaystyle=$ $\displaystyle\frac{\eta}{\lambda},$
(79a) $\displaystyle\tau_{\pm}$ $\displaystyle=$
$\displaystyle\frac{\alpha\pm\sqrt{\beta}}{6\lambda\chi},$ (79b)
$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\kappa^{2}\delta\chi\lambda+\chi(4\eta+\lambda-\chi)+(\lambda+\chi)\chi_{\perp},$
(79c) $\displaystyle\beta$ $\displaystyle=$
$\displaystyle\kappa^{4}\delta\chi^{2}\lambda^{2}-2\kappa^{2}\lambda\chi\delta\chi(-4\eta+\lambda+\chi)+2(\lambda+\chi)\chi_{\perp}\left[\kappa^{2}\delta\chi\lambda+\chi(4\eta+\lambda-\chi)\right]$
(79d)
$\displaystyle+\chi\left[16\eta^{2}\chi+8\eta\left(2\lambda^{2}+\lambda\chi-\chi^{2}\right)+\chi\left(-3\lambda^{2}-2\lambda\chi+\chi^{2}\right)\right]+(\lambda+\chi)^{2}\chi_{\perp}^{2}.$

The matter sector has two roots for (79a) with multiplicity 3 each and four
roots for (79b) (two for each $\tau_{\pm}$), a total of 10 roots. Now, (19) ad
(20) are observed if, and only if, $0\leq\tau_{a}<1$. For $\tau_{1}$ it is
guaranteed by (25) together with (28). As for $\tau_{\pm}$, it needs to be
real, i.e., $\beta\geq 0$, what corresponds to (77a), $\tau_{-}\geq 0$, and
$\tau_{+}<1$. For $\tau_{-}\geq 0$, we need that $\alpha\geq 0$ [condition
(77b)] together with $\alpha^{2}-\beta\geq 0$ [condition (77c)]. As for
$\tau_{+}<1$, it corresponds to $6\lambda\chi-\alpha\geq 0$ [condition (77d)]
together with $(6\lambda\chi-\alpha)^{2}-\beta>0$, i.e., condition (77e). All
the above conditions must be valid for all values of $\kappa\in[-1,1]$, i.e.,
for all possible values of the product $b=l_{\mu}v^{\mu}=\kappa v$. ∎

Finally, a much simpler case is given below:

######  Statement C.3 (The shearless case).

Consider the shearless anisotropic fluid ($\eta_{\perp}=\eta_{l}=\eta_{ll}=0$)
with $\delta\chi=0$ described by the energy-momentum tensor (10) and
supplemented by (11). The theory is causal if (25) is satisfied.

###### Proof.

In this particular case, the determinant in (B) becomes

$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{3\chi\lambda_{\perp}^{4}a^{6}}{4\varepsilon}\left(a^{2}-\frac{1}{3}v^{2}\right)^{2}=0.$
(80)

Note that there is one causal root $a=\xi_{\mu}u^{\mu}=0$ because in this case
$\xi_{\mu}\xi^{\mu}=-a^{2}+v^{2}=v^{2}>0$, with multiplicity 6, and two roots
$a^{2}=v^{2}/3$, also causal since $0<\tau=1/3<1$, with multiplicity 2 each,
completing the total 10 roots from the matter sector. Assumption (A1)
guarantees that the determinant $M$ is not trivially zero, which would give
any $\xi_{\mu}$ as a possible solution. ∎

##  Appendix D A causal example

In this appendix, we use statement C.1 to show that the set of parameters
(33), which are reproduced below for convenience, are causal

$\eta_{\perp}=\eta\,,\quad\eta_{l}=\frac{2\eta}{3}\,,\quad\eta_{ll}=\frac{5\eta}{6}\,,\quad\lambda_{\perp}=\frac{13\eta}{2}\,,\quad\lambda_{l}=6\eta\,,\quad\chi=5\eta\,,\quad\chi_{\perp}=\frac{11\eta}{2}\,,\quad\chi_{l}=\frac{16\eta}{3}\;.$

The above parameters satisfy (25) and (28) by construction, and we show that
satisfies the conditions of statement C.1 as well. For simplicity, we take the
overall factor $\eta^{3}/216$ out of the definition of $p(\varrho)$ in (73) to
obtain

$4\varepsilon
H^{\parallel}_{\parallel}(V,\xi)\left[A^{2}+A\left(U_{1}^{\mu}l_{\mu}+U_{2}^{\mu}\xi_{\mu}\right)+U_{1}^{\mu}l_{\mu}U_{2}^{\nu}\xi_{\nu}-U_{1}^{\mu}\xi_{\mu}U_{2}^{\nu}l_{\nu}\right]=\frac{\eta^{3}}{216}p(\varrho)v^{6}\;.$
(81)

Then, the coefficients in $p(\varrho)$ are given by

$\displaystyle\alpha_{0}$ $\displaystyle=$
$\displaystyle-364-1623\kappa^{2}+1674\kappa^{4}-119\kappa^{6}<0,\quad\forall\;\kappa^{2}\in[0,1]$
(82a) $\displaystyle\alpha_{1}$ $\displaystyle=$ $\displaystyle
16224+19601\kappa^{2}-19925\kappa^{4}>0,\quad\forall\;\kappa\in[0,1]$ (82b)
$\displaystyle\alpha_{2}$ $\displaystyle=$
$\displaystyle-18(7254-203\kappa^{2})<0,\quad\forall\;\kappa\in[0,1]$ (82c)
$\displaystyle\alpha_{3}$ $\displaystyle=$ $\displaystyle 126360\;,$ (82d)

where $\kappa$ is defined through $l_{\mu}v^{\mu}=\kappa v$ (see Appendix C).
For $\varrho<0$, we have

$p(\varrho)=-\left(\absolutevalue{\alpha_{0}}+\alpha_{1}\absolutevalue{\varrho}+\absolutevalue{\alpha_{2}}\varrho^{2}+\alpha_{3}\absolutevalue{\varrho}^{3}\right)<0\;,$
(83)

and, thus, (76b) is verified. On the other hand, for $\varrho\geq 1$

$p(\varrho)=\alpha_{0}+\alpha_{1}\varrho+(\alpha_{2}+\alpha_{3}\varrho^{2})\varrho\geq\alpha_{0}+\alpha_{1}+\alpha_{2}+\alpha_{3}\;,$
(84)

since $\alpha_{1},\alpha_{3}>0$. However,

$\alpha_{0}+\alpha_{1}+\alpha_{2}+\alpha_{3}=11648+21632\kappa^{2}-18251\kappa^{4}-119\kappa^{6}>0,\quad\forall\,\kappa^{2}\in[0,1],$
(85)

and, therefore, (76a) is also verified. Finally,

$\displaystyle
18\alpha_{0}\alpha_{1}\alpha_{2}\alpha_{3}-4\alpha_{2}^{3}\alpha_{0}+\alpha_{2}^{2}\alpha_{1}^{2}-4\alpha_{3}\alpha_{1}^{3}-27\alpha_{3}^{2}\alpha_{0}^{2}=324\left(2421746461615104-6266597114164608\kappa^{2}\right.$
$\displaystyle+24310721163158820\kappa^{4}-50747526702105948\kappa^{6}+59336411451755437\kappa^{8}-40187158939087070\kappa^{10}$
$\displaystyle\hskip
85.35826pt\left.+12398536143066361\kappa^{12}\right)>0\;,\quad\forall\,\kappa^{2}\in[0,1]\;.$
(86)

Hence, (76c) is also verified and the set of transport parameters (33) is
causal.

##  Appendix E Stability of the causal example (33)

In this appendix, we examine the stability of the sound channel (37) for the
causal set of parameters (33). In this case, up to an overall constant and,
because $\bar{\eta}>0$, by performing the changes $\Gamma\to\Gamma/\bar{\eta}$
and $k^{i}\to k^{i}/\bar{\eta}$ we obtain that Eq. (37) becomes

$a_{0}\Gamma^{6}+a_{1}\Gamma^{5}+a_{2}\Gamma^{4}+a_{3}\Gamma^{3}+a_{4}\Gamma^{2}+a_{5}\Gamma+a_{6}=0\;,$
(87)

where

$\displaystyle a_{0}$ $\displaystyle=$ $\displaystyle 21060\;,$ (88a)
$\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle 10962\;,$ (88b)
$\displaystyle a_{2}$ $\displaystyle=$ $\displaystyle
6\left[315+f_{1}(x^{2})k^{2}\right]\;,$ (88c) $\displaystyle a_{3}$
$\displaystyle=$ $\displaystyle 108+f_{2}(x^{2})k^{2}\;,$ (88d) $\displaystyle
a_{4}$ $\displaystyle=$
$\displaystyle\frac{k^{2}}{2}\left[f_{3}(x^{2})+f_{4}(x^{2})k^{2}\right],$
(88e) $\displaystyle a_{5}$ $\displaystyle=$ $\displaystyle
36k^{2}+f_{5}(x^{2})k^{4}\;,$ (88f) $\displaystyle a_{6}$ $\displaystyle=$
$\displaystyle k^{4}\left[f_{6}(x^{2})+f_{7}(x^{2})k^{2}\right]\;.$ (88g)

We have defined the functions

$\displaystyle f_{1}(x^{2})$ $\displaystyle=$ $\displaystyle 3498-70x^{2}\;,$
(89a) $\displaystyle f_{2}(x^{2})$ $\displaystyle=$ $\displaystyle
7248-270x^{2}\;,$ (89b) $\displaystyle f_{3}(x^{2})$ $\displaystyle=$
$\displaystyle 1680-36x^{2}\;,$ (89c) $\displaystyle f_{4}(x^{2})$
$\displaystyle=$ $\displaystyle 5548+6564x^{2}-6464x^{4}\;,$ (89d)
$\displaystyle f_{5}(x^{2})$ $\displaystyle=$ $\displaystyle
485+556x^{2}-553x^{4}\;,$ (89e) $\displaystyle f_{6}(x^{2})$ $\displaystyle=$
$\displaystyle 24+30x^{2}-30x^{4}\;,$ (89f) $\displaystyle f_{7}(x^{2})$
$\displaystyle=$ $\displaystyle 78+306x^{2}-310x^{4}+22x^{6}\;.$ (89g)

One may verify that all functions in (89) are positive for all
$x^{2}\in[0,1]$, which makes all $a_{I}$’s in (88) ($I=0,\cdots,6$) to be also
positive. Thus, pure real roots are in fact negative. As for imaginary roots,
we apply the Routh-Hurwitz criterion (RHC) Gradshteyn and Ryzhik (2007), which
in this case requires us to compute the following table

$\begin{array}[]{|c|c|c|c|c|}\hline\cr a_{0}&a_{2}&a_{4}&a_{6}&0\\\ \hline\cr
a_{1}&a_{3}&a_{5}&0&0\\\ \hline\cr b_{1}&b_{2}&b_{3}&0&0\\\ \hline\cr
c_{1}&c_{2}&0&0&0\\\ \hline\cr d_{1}&d_{2}&0&0&0\\\ \hline\cr e_{1}&0&0&0&0\\\
\hline\cr\end{array}$ (90)

where $b_{i}=(a_{1}a_{2i}-a_{0}a_{2i+1})/a_{1}$,
$c_{i}=(b_{1}a_{2i+1}-a_{1}b_{i+1})/b_{1}$,
$d_{i}=(c_{1}b_{i+1}-b_{1}c_{i+1})/c_{1}$, and
$e_{1}=(d_{1}c_{2}-c_{1}d_{2})/d_{1}$. Since $a_{I}>0$ for $I=0,\cdots,6$,
then $\real(\Gamma)<0$ if, and only if, $b_{1}>0$, $c_{1}>0$, $d_{1}>0$, and
$e_{1}>0$. Since $a_{1}>0$, then it is enough to obtain

$\displaystyle a_{1}b_{1}$ $\displaystyle=$ $\displaystyle
324\left[56925+(238974+3340x^{2})k^{2}\right]>0\;,$ (91a) $\displaystyle
b_{1}c_{1}$ $\displaystyle=$
$\displaystyle\frac{36}{203}\left[1024650+18g_{1}(x^{2})k^{2}+g_{2}(x^{2})k^{4}\right]\;,$
(91b) $\displaystyle b_{1}c_{1}d_{1}$ $\displaystyle=$
$\displaystyle\frac{36k^{2}}{203}\left[6147900g_{3}(x^{2})+27g_{4}(x^{2})k^{2}+6g_{5}(x^{2})k^{4}+g_{6}(x^{2})k^{6}\right]\;,$
(91c) $\displaystyle b_{1}c_{1}d_{1}e_{1}$ $\displaystyle=$
$\displaystyle\frac{36k^{4}}{203}\left[221324400g_{7}(x^{2})+972g_{8}(x^{2})k^{2}+27g_{9}(x^{2})k^{4}+6g_{10}(x^{2})k^{6}+g_{11}(x^{2})k^{8}\right],$
(91d)

where

$\displaystyle g_{1}(x^{2})$ $\displaystyle=$ $\displaystyle
1412154-77159x^{2}\;,$ (92a) $\displaystyle g_{2}(x^{2})$ $\displaystyle=$
$\displaystyle 174806118-143563237x^{2}+133959417x^{4}\;,$ (92b)
$\displaystyle g_{3}(x^{2})$ $\displaystyle=$ $\displaystyle 35-3x^{2}\;,$
(92c) $\displaystyle g_{4}(x^{2})$ $\displaystyle=$ $\displaystyle
422612235-188238658x^{2}+129974883x^{4}\;,$ (92d) $\displaystyle g_{5}(x^{2})$
$\displaystyle=$ $\displaystyle
20839149333-15122827103x^{2}+12589851858x^{4}+623807058x^{6}\;,$ (92e)
$\displaystyle g_{6}(x^{2})$ $\displaystyle=$ $\displaystyle
219646045656+96361454524x^{2}-432604607986x^{4}+620066970164x^{6}$ (92f)
$\displaystyle-290471866054x^{8}\;,$ $\displaystyle g_{7}(x^{2})$
$\displaystyle=$ $\displaystyle 23-18x^{2}+15x^{4}\;,$ (92g) $\displaystyle
g_{8}(x^{2})$ $\displaystyle=$ $\displaystyle
362387218-297049991x^{2}+243855688x^{4}-5179815x^{6}\,$ (92h) $\displaystyle
g_{9}(x^{2})$ $\displaystyle=$ $\displaystyle
216917135055-156807618144x^{2}+105968706122x^{4}+35462079400x^{6}$ (92i)
$\displaystyle-13719895953x^{8}\;,$ $\displaystyle g_{10}(x^{2})$
$\displaystyle=$ $\displaystyle
5000654914056-3108012289024x^{2}+1415965400581x^{4}+1798665419999x^{6}$ (92j)
$\displaystyle-134955905301x^{8}-625170319671x^{10}\;,$ $\displaystyle
g_{11}(x^{2})$ $\displaystyle=$ $\displaystyle
20097012270600-82381572326488x^{2}+329458308450878x^{4}-682695134979400x^{6}$
(92k)
$\displaystyle+795041858933812x^{8}-532773739193376x^{10}+161024463542854x^{12}\;.$

Since all 11 functions $g$ are greater than zero for all $x^{2}\in[0,1]$, then
all equations in (91) are positive, the RHC are verified and the system is
stable in the LRF. Since we proved in Sec. D that the system under these
parameters is causal, then the result in Gavassino (2022) ensures linear
stability in any boosted frame.

Just for the sake of illustration, let us verify stability in the boosted
homogeneous frame by taking the changes (38) in (87) and then taking
$k^{i}=0$. This will lead us to the roots $\Gamma=0$ of multiplicity 3 and the
roots of X.

$\beta_{0}(\gamma\Gamma)^{3}+\beta_{1}(\gamma\Gamma)^{2}+\beta_{2}\gamma\Gamma+\beta_{3}=0\;,$
(93)

where

$\displaystyle\beta_{0}$ $\displaystyle=$ $\displaystyle
2(1384+1698x^{2}-1461x^{4}-11x^{6})>0,\quad\forall\;x^{2}\in[0,1]\;,$ (94a)
$\displaystyle\beta_{1}$ $\displaystyle=$ $\displaystyle
4199+826x^{2}-553x^{4}>0,\quad\forall\;x^{2}\in[0,1]\;,$ (94b)
$\displaystyle\beta_{2}$ $\displaystyle=$ $\displaystyle
6(179+8x^{2}-5x^{4})>0,\quad\forall\;x^{2}\in[0,1]\;,$ (94c)
$\displaystyle\beta_{3}$ $\displaystyle=$ $\displaystyle 72\;.$ (94d)

Since all coefficients $\beta_{0,1,2,3}$ are positive, then there is only one
negative pure real root for the polynomial, as desired. As for the complex
roots, the remaining RHC to be computed is

$\beta_{1}\beta_{2}-\beta_{0}\beta_{3}=6(718405+140694x^{2}-78310x^{4}-8290x^{6}+2765x^{8})>0\;,$
(95)

which is greater than zero for all $x^{2}\in[0,1]$. Thus, linear stability is
also verified in the homogeneous boosted frame.

##  Appendix F The unphysical choice of $l$ in the Bjorken flow

Here, we repeat the study of Bjorken flow in Sec. VI, with an alternative
choice of spacelike anisotropy vector, $l=\partialderivative{x}$. The fluxes
of (11) in this case are

$\displaystyle\mathcal{E}^{(1)}$ $\displaystyle=$
$\displaystyle\tilde{\chi}^{3}T^{3}\left(\frac{1}{\tau}\,+\frac{3\dot{T}}{T}\right)$
(96a) $\displaystyle\mathcal{P}_{l}^{(1)}$ $\displaystyle=$ $\displaystyle
T^{3}\left(\frac{2\tilde{\eta}_{ll}+\tilde{\chi}_{l}}{3\tau}+\frac{3\tilde{\chi}_{l}\dot{T}}{T}\right)$
(96b) $\displaystyle\mathcal{P}_{\perp}^{(1)}$ $\displaystyle=$ $\displaystyle
T^{3}\left(\frac{\tilde{\chi}_{\perp}-\tilde{\eta}_{ll}}{3\tau}+\frac{3\tilde{\chi}_{\perp}\dot{T}}{T}\right)$
(96c) $\displaystyle\pi_{\perp}^{\mu\nu}$ $\displaystyle=$
$\displaystyle-2\eta_{\perp}T^{3}\mathrm{diag}\left(0,0,1/\tau,-1/\tau^{3}\right),$
(96d) $\displaystyle W_{\perp l}^{\mu}$ $\displaystyle=$ $\displaystyle 0,$
(96e) $\displaystyle W_{\perp u}^{\mu}$ $\displaystyle=$ $\displaystyle 0,$
(96f) $\displaystyle M$ $\displaystyle=$ $\displaystyle 0,$ (96g)

which gives rise to the following EOM

$\displaystyle
9\tilde{\chi}\frac{\tau^{2}\ddot{T}}{T}+18\tilde{\chi}\frac{\tau^{2}\dot{T}^{2}}{T^{2}}+\left(\frac{3\tau(6\tilde{\chi}+\tilde{\chi}_{\perp})}{T}+12\tau^{2}\right)\dot{T}+4\tau
T+\tilde{\chi}_{\perp}-3\tilde{\eta}_{\perp}-\tilde{\eta}_{ll}=0\,.$ (97)

The above can be expressed in terms of $w=T\tau$ and
$f(w)=\frac{\tau}{w}\derivative{w}{\tau}$ as

$\displaystyle\frac{9\tilde{\chi}}{4}f(w)^{2}+wf(w)\left(1+\frac{3}{4}\tilde{\chi}f^{\prime}(w)\right)-\frac{\tilde{\chi}-15\tilde{\chi}_{\perp}}{4}f(w)+\frac{18\tilde{\chi}-2\tilde{\chi}_{\perp}-3\tilde{\eta}_{\perp}-\tilde{\eta}_{ll}}{12}-\frac{2w}{3}=0\,.$
(98)

The pressure anisotropy reads

${\mathcal{A}}=\frac{2}{3}\frac{\tilde{\chi}_{\perp}-\tilde{\chi}_{l}}{w}\left(f-\frac{2}{3}\right)-\frac{3\tilde{\eta}_{ll}}{w}\,.$
(99)

As in Sec. VI, ${\mathcal{A}}$ has an “offshell” contribution, which cannot be
found in the isotropic conformal BDNK. However, the first-order term is
negative, in contrast to the isotropic case. The late-time expansion for $T$
is

$T=\frac{\Lambda}{(\Lambda\tau)^{1/3}}\left(1-\frac{\tilde{\eta}_{ll}+3\tilde{\eta}_{\perp}}{8(\Lambda\tau)^{2/3}}-\frac{(\tilde{\eta}_{ll}+3\tilde{\eta}_{\perp})(5\tilde{\chi}-\tilde{\chi}_{\perp})}{64(\Lambda\tau)^{2/3}}+\cdots\right)\,,$
(100)

and for $f$ is

$f(w)=\frac{2}{3}+\frac{\tilde{\eta}_{ll}+3\tilde{\eta}_{\perp}}{12w}+\frac{(\tilde{\eta}_{ll}+3\tilde{\eta}_{\perp})(5\tilde{\chi}-\tilde{\chi}_{\perp})}{48w^{2}}+\order{\frac{1}{w^{3}}\,}\,.$
(101)

We may consider the following linear perturbation at late times

$f(w)=\frac{2}{3}+\frac{\tilde{\eta}_{ll}+3\tilde{\eta}_{\perp}}{12w}+\delta
f(w)\,,$ (102)

which up to the first order in perturbation in late times, is

$\delta
f(w)\sim\exp(-\frac{2w}{\tilde{\chi}})w^{\frac{\tilde{\eta}_{ll}+3\tilde{\eta}_{\perp}+2(\tilde{\chi}_{\perp}+\tilde{\chi}_{l})}{4\tilde{\chi}}}\,.$
(103)

The numerical attractor can be found from the following initial condition

$f(w\ll
1)=\frac{7}{9}+\frac{\tilde{\chi}_{l}-\tilde{\chi}_{\perp}}{36\tilde{\chi}}+\frac{\sqrt{12\left(\tilde{\eta}_{ll}+3\tilde{\eta}_{\perp}\right)\tilde{\chi}+\left(\tilde{\chi}_{\perp}+\tilde{\chi}_{l}\right)^{2}}}{18\tilde{\chi}}\,,$
(104)

and the slow-roll attractor is

$f(w)_{\rm
slowroll}=\frac{7}{9}+\frac{\tilde{\chi}_{l}-\tilde{\chi}_{\perp}}{36\tilde{\chi}}-\frac{2w}{9\tilde{\chi}}+\frac{\sqrt{\left(4w-(\tilde{\chi}_{l}+\tilde{\chi}_{\perp})\right)^{2}+12\left(3\tilde{\eta}_{\perp}+\tilde{\eta}_{ll}\right)\tilde{\chi}}}{18\tilde{\chi}}$
(105)

Neglecting the power counting argument, we find $\nabla\cdot S_{\rm off}$ is
initially negative if

$\chi_{\perp}>\eta_{ll}+3\eta_{\perp}>0\,.$ (106)

In the isotropic limit, the above reproduces the stability condition

$\chi>4\eta\,.$

The condition $\nabla\cdot S<0$ at $w=0$ is equivalent to

$\frac{2}{3}<f(0)<1\,,$ (107)

which prevents early reheating for the attractor solution.

## References

  * de Rham  _et al._ (2022) C. de Rham, S. Kundu, M. Reece, A. J. Tolley, and S.-Y. Zhou, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.06805 [hep-th] .
  * Hawking and Ellis (2023) S. W. Hawking and G. F. R. Ellis,  _The Large Scale Structure of Space-Time_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2023).
  * Wald (1984) R. M. Wald,  _General Relativity_ (Chicago Univ. Pr., Chicago, USA, 1984).
  * Landau and Lifshitz (2013) L. Landau and E. Lifshitz,  _Fluid Mechanics: Landau and Lifshitz: Course of Theoretical Physics, Volume 6_, v. 6 (Elsevier Science, 2013).
  * Eckart (1940) C. Eckart, Phys. Rev. 58, 919 (1940).
  * Hiscock and Lindblom (1985a) W. A. Hiscock and L. Lindblom, Phys. Rev. D 31, 725 (1985a).
  * Pichon (1965) G. Pichon,  _Annales de l’IHP Physique théorique_ , 2, 21 (1965).
  * Israel (1976) W. Israel, Annals Phys. 100, 310 (1976).
  * Israel and Stewart (1979) W. Israel and J. M. Stewart, Annals Phys. 118, 341 (1979).
  * Baier  _et al._ (2008) R. Baier, P. Romatschke, D. T. Son, A. O. Starinets, and M. A. Stephanov, JHEP 04, 100 (2008), arXiv:0712.2451 [hep-th] .
  * Denicol  _et al._ (2012) G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, Phys. Rev. D 85, 114047 (2012), [Erratum: Phys.Rev.D 91, 039902 (2015)], arXiv:1202.4551 [nucl-th] .
  * Hiscock and Lindblom (1983) W. A. Hiscock and L. Lindblom, Annals Phys. 151, 466 (1983).
  * 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._ (2019) F. S. Bemfica, M. M. Disconzi, and J. Noronha, Phys. Rev. D 100, 104020 (2019), arXiv:1907.12695 [gr-qc] .
  * Bemfica  _et al._ (2021a) F. S. Bemfica, M. M. Disconzi, C. Rodriguez, and Y. Shao, Commun. Pure Appl. Anal. 20, 2279 (2021a).
  * Bemfica  _et al._ (2021b) F. S. Bemfica, M. M. Disconzi, and P. J. Graber, Commun. Pure Appl. Anal. 20, 2885 (2021b).
  * 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] .
  * Florkowski and Ryblewski (2011) W. Florkowski and R. Ryblewski, Phys. Rev. C 83, 034907 (2011), arXiv:1007.0130 [nucl-th] .
  * Bazow  _et al._ (2014) D. Bazow, U. W. Heinz, and M. Strickland, Phys. Rev. C 90, 054910 (2014), arXiv:1311.6720 [nucl-th] .
  * Molnár  _et al._ (2016) E. Molnár, H. Niemi, and D. H. Rischke, Phys. Rev. D 94, 125003 (2016), arXiv:1606.09019 [nucl-th] .
  * Molnar  _et al._ (2016) E. Molnar, H. Niemi, and D. H. Rischke, Phys. Rev. D 93, 114025 (2016), arXiv:1602.00573 [nucl-th] .
  * Alqahtani  _et al._ (2018) M. Alqahtani, M. Nopoush, and M. Strickland, Prog. Part. Nucl. Phys. 101, 204 (2018), arXiv:1712.03282 [nucl-th] .
  * Cercignani (1988) C. Cercignani,  _The Boltzmann Equation and Its Applications_ (Springer-Verlag, New York, 1988).
  * Groot (1980) S. R. D. Groot,  _Relativistic Kinetic Theory. Principles and Applications_ (Amsterdam, Netherlands: North-holland ( 1980) 417p, 1980).
  * Dore  _et al._ (2022) T. Dore, L. Gavassino, D. Montenegro, M. Shokri, and G. Torrieri, Annals Phys. 442, 168902 (2022), arXiv:2109.06389 [hep-th] .
  * Loganayagam (2008) R. Loganayagam, JHEP 05, 087 (2008), arXiv:0801.3701 [hep-th] .
  * Hiscock and Lindblom (1985b) W. A. Hiscock and L. Lindblom, Phys. Rev. D 31, 725 (1985b).
  * Gavassino (2022) L. Gavassino, Phys. Rev. X 12, 041001 (2022), arXiv:2111.05254 [gr-qc] .
  * Shokri and Taghinavaz (2020) M. Shokri and F. Taghinavaz, Phys. Rev. D 102, 036022 (2020), arXiv:2002.04719 [hep-th] .
  * Heller and Spalinski (2015) M. P. Heller and M. Spalinski, Phys. Rev. Lett. 115, 072501 (2015), arXiv:1503.07514 [hep-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] .
  * Behtash  _et al._ (2019a) A. Behtash, C. N. Cruz-Camacho, S. Kamata, and M. Martinez, Phys. Lett. B 797, 134914 (2019a), arXiv:1805.07881 [hep-th] .
  * Behtash  _et al._ (2019b) A. Behtash, S. Kamata, M. Martinez, and H. Shi, Phys. Rev. D 99, 116012 (2019b), arXiv:1901.08632 [hep-th] .
  * Behtash  _et al._ (2021) A. Behtash, S. Kamata, M. Martinez, T. Schäfer, and V. Skokov, Phys. Rev. D 103, 056010 (2021), arXiv:2011.08235 [hep-ph] .
  * Bemfica  _et al._ (2023) F. S. Bemfica, M. M. Disconzi, J. Noronha, and R. J. Scherrer, Phys. Rev. D 107, 023512 (2023), arXiv:2210.13372 [gr-qc] .
  * Denicol  _et al._ (2019) G. S. Denicol, E. Molnár, H. Niemi, and D. H. Rischke, Phys. Rev. D 99, 056017 (2019), arXiv:1902.01699 [nucl-th] .
  * Armas and Camilloni (2022) J. Armas and F. Camilloni, JCAP 10, 039 (2022), arXiv:2201.06847 [hep-th] .
  * Jensen  _et al._ (2012) K. Jensen, M. Kaminski, P. Kovtun, R. Meyer, A. Ritz, and A. Yarom, Phys. Rev. Lett. 109, 101601 (2012), arXiv:1203.3556 [hep-th] .
  * Banerjee  _et al._ (2012) N. Banerjee, J. Bhattacharya, S. Bhattacharyya, S. Jain, S. Minwalla, and T. Sharma, JHEP 09, 046 (2012), arXiv:1203.3544 [hep-th] .
  * Courant and Hilbert (1991) C. Courant and D. Hilbert,  _Methods of Mathematical Physics_ , 1st ed., Vol. 2 (John Wiley & Sons, Inc., 1991) p. 852.
  * Choquet-Bruhat (2009) Y. Choquet-Bruhat,  _General Relativity and the Einstein Equations_ (Oxford University Press, New York, 2009).
  * Gradshteyn and Ryzhik (2007) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ , seventh ed. (Elsevier/Academic Press, Amsterdam, 2007) pp. xlviii+1171, translated from the Russian, Translation edited and with a preface by Alan Jeffrey and Daniel Zwillinger, With one CD-ROM (Windows, Macintosh and UNIX).