# Wave Phenomena In General Relativistic Magnetohydrodynamics

Ankit Kumar Panda<EMAIL_ADDRESS>Victor Roy<EMAIL_ADDRESS>School of Physical Sciences, National Institute of Science Education and
Research, An OCC of Homi Bhabha National Institute, Jatni-752050, India

(August 27, 2024)

###### Abstract

Here we study the wave propagation and stability of general relativistic non-
resistive dissipative second-order magnetohydrodynamic equations in curved
space-time. We solve the Boltzmann equation for a system of particles and
antiparticles using the relaxation time approximation and the Chapman-Enskog-
like gradient expansion for the off-equilibrium distribution function,
truncating beyond second-order in curved space-time in electromagnetic fields.
Unlike holographic calculation [1], we show that the viscous evolution
equations do not explicitly depend on the curvature of space-time. Also, we
have tested the causality and stability of the second-order theory in curved
space-time in the presence of linearised metric perturbation and derived
dispersion relations for various modes. Interestingly, we found the coupling
of gravitational modes with the usual magneto-sonic modes in the small wave-
number limit. Also, we show additional non-hydrodynamical modes arise due to
gravity for a bulk-viscous fluid.

††preprint: APS/123-QED

##  I Introduction

General Relativistic Magnetohydrodynamics (GRMHD) is an essential frame-work
to study astrophysical systems such as binary neutron star mergers, which is
one of the source of the gravitational waves [2, 3, 4, 5]. Usually, in GRMHD
numerical studies, the ideal magnetohydrodynamics limit ( large magnetic
Reynolds number) is used, and other dissipative effects (due to bulk and shear
viscosity, etc.) are usually neglected [4, 6, 7, 8, 9]. The ideal-MHD (note
that we use ”ideal” to denote large electrical conductivity, and zero or
vanishing viscosity is termed as non-viscous fluid to avoid confusion) might
be a good approximation for large celestial bodies where the length scale
involved is large, and gradients are small [10, 11]. On the other hand, in
heavy-ion collisions, two heavy nuclei collide and create a hot and dense
Quark-Gluon-Plasma (QGP) that evolves for a few femtoseconds like a strongly
coupled fluid, gradients are usually large, and dissipative effects are not
negligible. In this case the positively charged nuclei move with relativistic
speed and generate one of the strongest transient magnetic fields in the
universe [12, 13, 14, 15]. The QGP produced in the heavy-ion collisions
behaves like a fluid with negligible specific shear viscosity [16, 17, 18,
19]. Theoretical studies show that the QGP has finite electrical conductivity;
hence the evolution of QGP in a strong electromagnetic background field
naturally indicates that relativistic magnetohydrodynamics is the appropriate
framework for the space-time evolution. Unlike astrophysical systems, the
space-time metric considered in heavy-ion collisions is flat. In short, in
GRMHD, we generally neglect the viscous effects, and in special relativistic
MHD (as in heavy-ion collisions), we ignore the curvature of space-time.

Any relativistically consistent theory, in principle, must obey the causality
condition, i.e., the speed of propagation of perturbations cannot be
superluminal. For the case of relativistic viscous hydrodynamics, it is well
known that relativistic extension of the Navier-Stokes leads to acausality
[20] and one of the remedies were suggested by Mueller-Israel-Stewart
(MIS)[21, 22] through the inclusion of second-order dissipative corrections in
entropy four-current. As we know, the second-order relativistic viscous
hydrodynamics equations give rise to a causal propagation (in a restricted
sense) of the perturbation in fluids [23, 24, 25, 26] and the same is true for
a charged fluid in an external electromagnetic field [27]. The second-order
causal dissipative ideal [28, 29] and resistive [30, 31] MHD formulation were
recently derived from the relativistic Boltzmann kinetic equation using the
moment method and the relaxation time approximations. In these MIS-type causal
theories, the dissipative fluxes are treated as independent variables and are
governed by relaxation-type equations. The two formulations in [28, 29, 30,
31] produce almost identical equations for the evolution of these dissipative
quantities, albeit with different transport coefficients and some extra
coefficients in RTA approximation. In both formulations, the space-time metric
was assumed to be flat. In the present study, we lift this restriction, allow
the space-time to be curved, and formulate the second-order dissipative
resistive GRMHD. We show that this gives rise to causal theory when subjected
to linear perturbations. Throughout the paper, we assume that the
gravitational field varies very little over the fluid particles’ mean free
path and mean free time.

The manuscript is arranged as follows: in Section (II), we describe the
Boltzmann equation in the Riemannian space, where we use the relaxation time
approximation for simplifying the collision kernel. In the next Section (III),
we define various macroscopic fluid variables from the kinetic theory using a
one-particle phase-space distribution function; we also discuss the
conservation of energy-momentum and the Maxwell equations in the covariant
form. The following Section (IV) contains the first and second-order general
relativistic magnetohydrodynamics equations derived from the relativistic
kinetic theory. In Section (V) we investigate the causality and stability
using linear perturbation considering dynamic space-time metric. Section (IV)
and (V) contain the main results. Finally, we summarise our work in Section
(VI).

Throughout the work we use metric with mostly negative sign i.e.,
$g_{\mu\nu}=diag\left(+,-,-,-\right)$. We also use the natural units where
$\hbar=c=k_{B}=\mu_{0}=1$. The gravitational constant is denoted by G, and the
cosmological constant is denoted by $\Lambda$.

##  II The Boltzmann equation in a gravitational and electromagnetic field

Here we give a short textbook introduction to the Boltzmann equation in
Riemannian space in the presence of an external force due to electromagnetic
field. We further assume the collision integral in the Boltzmann equation is
given by the relaxation time approximation (RTA). Though it is one of the
simplest forms of the collision kernel in use, it essentially captures almost
all the system’s important features initially away from equilibrium and
relaxing back to the nearest equilibrium state through the inter-particle
collisions. For simplicity, we also assume the relaxation time is a Lorentz
scalar and independent of the momentum of the particles. However, this
momentum independence of the relaxation time is not a necessary condition for
the RTA to be applicable [32, 33].

Now, we recall that in the Minkowski space the phase space volume element in
two different inertial reference frames are related as
$d^{3}xd^{3}p=d^{3}x^{\prime}d^{3}p^{\prime}$. The one-particle distribution
function $f(x^{\mu},p^{\mu})$ is usually written as $f(\bf{x},\bf{p},t)$ due
to the on shell condition $p^{0}=\sqrt{|{\bf{p}}|^{2}+m^{2}}$. In a Riemannian
space the following relations hold:
$\sqrt{-g^{\prime}}\frac{d^{3}p^{\prime}}{p_{0}^{\prime}}=\sqrt{-g}\frac{d^{3}p}{p_{0}}$,
$\sqrt{-g^{\prime}}d^{4}x^{\prime}=\sqrt{-g}d^{4}x$, The four velocity is
defined as $u^{\mu}=\Gamma(1,\bf{v})$, where the Lorentz gamma
$\Gamma=\frac{1}{\sqrt{g_{00}\left(1+\frac{g_{0i}}{g_{00}}\beta^{i}\right)^{2}-\beta^{2}}}$.
The covariant components are given by $u_{\mu}=g_{\mu\nu}u^{\nu}$ with the
normalisation $u_{\mu}u^{\mu}=1$.

Considering the Boltzmann equation in a curved space-time under an
electromagnetic force, the evolution of one-particle distribution function
$f({\bf x},{\bf p},t)$ is given by [34]

$\displaystyle p^{\mu}\partial_{\mu}f+qF^{\mu\nu}p_{\nu}\frac{\partial
f}{\partial p^{\mu}}-\Gamma^{\mu}_{\alpha\rho}p^{\alpha}p^{\rho}\frac{\partial
f}{\partial p^{\mu}}=-\frac{u\cdot p}{\tau_{c}}\delta f,$ $\displaystyle
p^{\mu}\partial_{\mu}f+F^{\mu}\frac{\partial f}{\partial
p^{\mu}}=-\frac{u\cdot p}{\tau_{c}}\delta f$ (1)

where we define $F^{\mu}\equiv
qF^{\mu\nu}p_{\nu}-\Gamma^{\mu}_{\alpha\rho}p^{\alpha}p^{\rho}$ which includes
the effect of gravity (which is approximated here as space-time curvature) and
$F^{\mu\nu}=g^{\mu\alpha}F_{\alpha\beta}g^{\beta\nu}$ is the antisymmetric
electromagnetic field tensor. The RHS of the Boltzmann equation Eq.(1) is the
collision kernel, which we have approximated to be of the form of relaxation
time approximation, with $\tau_{c}$ being the relaxation time.
For a system slightly away from the local equilibrium, the one-particle
distribution function $f$ can be decomposed into an equilibrium $f_{0}$ and an
off-equilibrium part $\delta f$ as $f=f_{0}+\delta f$. The equilibrium
distribution function for particles is given
by:$f_{0}=\frac{1}{e^{\beta\left(g_{\mu\nu}u^{\mu}p^{\nu}\right)-\alpha}+r}$,
$r=0,\pm 1$ correspond the boltzmann, the fermi-dirac, and the bosonic
distribution respectively. For anti-particles $\alpha\rightarrow-\alpha$ in
the above expression.
At this point, we would like to take a detour and discuss about the form of
the equilibrium distribution function $f_{0}$ in presence of gravity. We know
that for the equilibrium distribution function the RHS of the Boltzmann
equation vanish by definition. The LHS of the Boltzmann equation should also
vanish; this imposes some restriction on the hydrodynamic variables. If we put
$f_{0}$ in the LHS of the Boltzmann equation we get

$-p^{\alpha}p^{\beta}D_{\alpha}\left(\beta
u_{\beta}\right)+p^{\mu}D_{\mu}\alpha=0,$ (2)

here
$D_{\alpha}u_{\beta}=\partial_{\alpha}u_{\beta}-\Gamma_{\alpha\beta}^{\mu}u_{\mu}$
is the covariant derivative of $u_{\beta}$. In Eq.(2) all the terms are
written in terms of a polynomial equation of $p^{\mu}$, hence all the
components should vanish independently. It follows that

$\displaystyle D_{\mu}\alpha$ $\displaystyle=$ $\displaystyle 0,$ (3)

$-D_{\alpha}(\beta u_{\beta})=\begin{cases}0,&\text{if }m\neq 0\\\
f(x^{\delta})g_{\alpha\beta},&\text{if }m=0\end{cases}$

where $f(x^{\delta})$ is some arbitrary function of the space-time coordinate.
If we consider a massless gas then $g_{\alpha\beta}p^{\alpha}p^{\beta}=0$, and
for $f_{0}$ the LHS of the Boltzmann equation becomes zero. But for the case
of $m\neq 0$ from the above equation we get :

$\displaystyle g_{\sigma\beta}\partial_{\alpha}(\beta
u^{\sigma})+g_{\sigma\alpha}\partial_{\beta}(\beta u^{\sigma})+\beta
u^{\sigma}\partial_{\sigma}g_{\alpha\beta}$ $\displaystyle=$ $\displaystyle
0.$

The above equation is a Killing equation with $\beta u^{\mu}$ being the
Killing vector. If we define the components of the killing vector in a
particular reference frame as $\beta u^{\mu}=(\beta u^{0},\vec{0})$ where
$\beta u^{0}$ is constant, the killing equation leads to

$\displaystyle\partial_{0}g_{\alpha\beta}$ $\displaystyle=$ $\displaystyle 0.$
(4)

This suggests that the metric tensor must be time-independent for gas of non-
vanishing rest mass to define an equilibrium distribution function. However,
for the massless case, there are no such restrictions.

Another point to ponder is the Liouville’s theorem, which along with the on-
shell condition gives rise to the following equation:

$\displaystyle\frac{\partial}{\partial
x^{a}}\left(\frac{g}{p_{0}}p^{a}\right)$ $\displaystyle=$ $\displaystyle 0,$
(5)

where

$\displaystyle p^{a}\equiv(p^{\mu},-\Gamma^{i}_{\mu\nu}p^{\mu}p^{\nu}),$
$\displaystyle\frac{\partial}{\partial
x^{a}}\equiv\left(\frac{\partial}{\partial x^{\mu}},\frac{\partial}{\partial
p^{i}}\right).$

Eq.(5) is useful in proving the conservation of energy-momentum and particle
four-current.

##  III Definitions and conservation equations

Now that we have defined the equilibrium distribution function and the
Boltzmann equation in curved spacetime, we use these definitions to define
hydrodynamic variables and conservation laws that we will use later to derive
the equation of motion of a fluid. The particle four flow and the energy-
momentum tensors are defined as the first and second moment of the
distribution function:

$\displaystyle N^{\mu}_{f}$ $\displaystyle=$
$\displaystyle\int{}{}\frac{\mathcal{G}_{s}d^{3}p}{(2\pi)^{3}p_{0}}\sqrt{-g}p^{\mu}f,$
(6) $\displaystyle T^{\mu\nu}_{f}$ $\displaystyle=$
$\displaystyle\int{}{}\frac{\mathcal{G}_{s}d^{3}p}{(2\pi)^{3}p_{0}}\sqrt{-g}p^{\mu}p^{\nu}f.$
(7)

$\mathcal{G}_{s}$ here is the spin-degeneracy factor. Usually, the particle
four-current $N^{\mu}_{f}$ corresponds to some conserved charges, and in that
case, the conservation of particle and the energy-momentum tensor are obtained
by replacing the partial derivatives with the covariant derivatives wherever
applicable:

$\displaystyle D_{\mu}N^{\mu}_{f}$ $\displaystyle=$
$\displaystyle\partial_{\mu}N^{\mu}_{f}+\Gamma^{\mu}_{\mu\rho}N^{\rho}_{f}=0,$
(8) $\displaystyle D_{\mu}T^{\mu\nu}_{f}$ $\displaystyle=$
$\displaystyle\partial_{\mu}T^{\mu\nu}_{f}+\Gamma^{\mu}_{\mu\rho}T^{\nu\rho}_{f}+\Gamma^{\nu}_{\mu\rho}T^{\mu\rho}_{f}=0,$
(9)

here $\Gamma^{\mu}_{\nu\rho}$ is the Christoffel symbol of second kind and is
given by :
$\Gamma^{\mu}_{\nu\rho}=\frac{g^{\mu\lambda}}{2}\left(\partial_{\nu}g_{\rho\lambda}+\partial_{\rho}g_{\nu\lambda}-\partial_{\lambda}g_{\nu\rho}\right)$
, $u\cdot p$=$g_{\mu\nu}u^{\mu}p^{\nu}$ , and
$g_{\alpha\beta}u^{\alpha}u^{\beta}=1$.
Also, we use $D=u^{\mu}D_{\mu}$ for the covariant time derivative, and the
covariant derivative is decomposed as $D_{\mu}=u_{\mu}D+\nabla_{\mu}$, where
$\nabla_{\mu}=\Delta^{\alpha}_{\mu}D_{\alpha}$. We decompose the gradient of
four velocity in terms of symmetric traceless shear, non-zero trace bulk, and
antisymmetric second rank vorticity tensor as
$\nabla_{\mu}u_{\nu}=\sigma_{\mu\nu}+\Delta_{\mu\nu}\frac{\theta}{3}+\omega_{\mu\nu}$
where $\Delta_{\mu\nu}=g_{\mu\nu}-u_{\mu}u_{\nu}$ ,$\theta=D_{\mu}u^{\mu}$,
$D_{\mu}g^{\mu\nu}=0$. The electromagnetic field tensor is given by
$F_{\mu\nu}=D_{\mu}A_{\nu}-D_{\nu}A_{\mu}=E_{\mu}u_{\nu}-E_{\nu}u_{\mu}+B_{\mu\nu}$,
where $B_{\mu\nu}=\epsilon_{\mu\nu\alpha\beta}u^{\alpha}B^{\beta}$,
$B^{\mu\nu}=-Bb^{\mu\nu}$, $\epsilon_{\mu\nu\alpha\beta}=\sqrt{-g}$ , and
$\epsilon^{\mu\nu\alpha\beta}=\frac{-1}{\sqrt{-g}}$. Since we consider a
combined system of a fluid , electromagnetic fields, and the presence of a
conserved charged current, the corresponding energy-momentum tensor of the
fluid, the EM fields, and the four current of the conserved quantity are given
by:

$\displaystyle T^{\mu\nu}_{f}$ $\displaystyle=$ $\displaystyle\epsilon
u^{\mu}u^{\nu}-\left(P+\Pi\right)\Delta^{\mu\nu}+\pi^{\mu\nu},$ (10)
$\displaystyle T^{\mu\nu}_{EM}$ $\displaystyle=$
$\displaystyle-F^{\mu}_{\alpha}F^{\nu\alpha}+\frac{1}{4}g^{\mu\nu}F_{\alpha\beta}F^{\alpha\beta},$
(11) $\displaystyle N^{\mu}_{f}$ $\displaystyle=$ $\displaystyle
n_{f}u^{\mu}+V_{f}^{\mu}.$ (12)

The fluid though under consideration is non-polarisable and non-magnetisable.
While deriving the conservation equations from the Boltzmann equation Eq.(1)
we want to express the left hand side in terms of thermodynamic forces
(spatial gradients of velocities, temperature etc.). For that we note
$\partial_{\lambda}f_{0}=-f_{0}\tilde{f}_{0}\left[\beta
p^{\mu}D_{\lambda}u_{\mu}+\beta
u_{\mu}\Gamma^{\mu}_{\lambda\rho}P^{\rho}-\partial_{\lambda}\alpha+(p\cdot
u)\partial_{\lambda}\beta\right]$ (where we have used
$\tilde{f}_{0}=(1-rf_{0})$, $r=\pm 1$ and $\partial_{\lambda}p^{\nu}=0$),
which gives rise to terms such as $\dot{\alpha}$, $\dot{\beta}$,
$\dot{u_{\mu}}$. Such terms do not contribute in the expression for the
thermodynamic force, hence we need to eliminate them. This we achieve with the
help of the conservation equations given below along with the relativistic
thermodynamic integrals $J^{(m)\pm}_{nq}$ defined later in appendix(A).

$\displaystyle D_{\mu}T^{\mu\nu}_{EM}$ $\displaystyle=$
$\displaystyle-F^{\nu\lambda}J_{\lambda f},$ $\displaystyle D_{\mu}J^{\mu}$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle T^{\mu\nu}$
$\displaystyle=$ $\displaystyle T^{\mu\nu}_{f}+T^{\mu\nu}_{EM},$
$\displaystyle D_{\mu}T^{\mu\nu}_{f}$ $\displaystyle=$ $\displaystyle
F^{\nu\lambda}J_{\lambda f}.$

Now using the fluid conservation equation, the expressions for $\theta$ , $D$,
$\sigma_{\mu\nu}$, and $\omega_{\mu\nu}$ we get the following expressions for
$\dot{\alpha}$, $\dot{\beta}$ and $\dot{u}^{\mu}$ :

$\displaystyle\dot{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{D_{20}}\left[J_{20}^{(0)-}\theta\left(\epsilon+P+\Pi\right)-J_{30}^{(0)+}\left(n_{f}\theta+D_{\mu}V_{f}^{\mu}\right)+J_{20}^{(0)-}\left(-\pi^{\mu\nu}\sigma_{\mu\nu}+qE^{\mu}V_{f\mu}\right)\right],$
$\displaystyle\dot{\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{D_{20}}\left[J_{10}^{(0)+}\theta\left(\epsilon+P+\Pi\right)-J_{20}^{(0)-}\left(n_{f}\theta+D_{\mu}V_{f}^{\mu}\right)+J_{10}^{(0)+}\left(-\pi^{\mu\nu}\sigma_{\mu\nu}+qE^{\mu}V_{f\mu}\right)\right],$
$\displaystyle\dot{u}^{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon+P}\left[\frac{n_{f}}{\beta}\left(\nabla^{\mu}\alpha-h\nabla^{\mu}\beta\right)-\Pi\dot{u}^{\mu}+\nabla^{\mu}\Pi-\Delta^{\mu}_{\nu}D_{\rho}\pi^{\rho\nu}\right]+\frac{1}{\epsilon+P}\left[qn_{f}E^{\mu}-qBb^{\mu\nu}V_{f\nu}\right],$

where $D_{20}=J^{(0)+}_{30}J^{(0)+}_{10}-J_{20}^{(0)-}J_{20}^{(0)-}$,
$h=\frac{\epsilon+P}{n_{f}}$ and
$\sigma^{\mu\nu}=\Delta^{\mu\nu}_{\alpha\beta}\nabla^{\alpha}u^{\beta}$. These
relations will be used later to derive the off-equilibrium distribution
function $\delta f$.

Now the definition of viscous stresses given in Eq.(10) in terms of $\delta f$
are given by:

$\displaystyle\pi^{\mu\nu}$ $\displaystyle=$
$\displaystyle\Delta^{\mu\nu}_{\alpha\beta}\int\frac{\mathcal{G}_{s}d^{3}p\sqrt{-g}}{(2\pi)^{3}p_{0}}p^{\alpha}p^{\beta}(\delta
f),$ (13) $\displaystyle\Pi$ $\displaystyle=$
$\displaystyle-\frac{\Delta_{\alpha\beta}}{3}\int\frac{\mathcal{G}_{s}d^{3}p\sqrt{-g}}{(2\pi)^{3}p_{0}}p^{\alpha}p^{\beta}(\delta
f),$ (14) $\displaystyle V^{\mu}$ $\displaystyle=$
$\displaystyle\Delta^{\mu}_{\alpha}\int\frac{\mathcal{G}_{s}d^{3}p\sqrt{-g}}{(2\pi)^{3}p_{0}}p^{\alpha}(\delta
f).$ (15)

##  IV Order-by-order magnetohydrodynamic equations

###  IV.1 Calculation of $\delta f$

In this section we discuss the details of $\delta f$ which is used in Eqs.
(13),(14),(15) to derive the evolution equation for dissipative stresses upto
second order. As mentioned earlier, we use RTA collision kernel in the
Boltzmann equation to derive $\delta f$ as a series expansion around $f_{0}$
in terms of three smallness paramaters Knudsen number
$\mathrm{Kn}=\tau_{c}\partial_{\mu}$ [35] , $\chi=qB\tau_{c}/T$
,$\xi=qE\tau_{c}/T$ [29, 31, 36]. We note that this is not the unique way to
derive the $\delta f$, for example in [28, 37] a relativistic moment method
was used, and some other variant of RTA can be found in [38, 39, 40]. In
addition due to non-zero Christoffel symbols we also need to make sure that
$\tau_{c}\Gamma^{\mu}_{\alpha\beta}<<1$ holds (corresponds to nearly flat
geometry between two collision)which corresponds to a weak gravitaional force
limit, $\mathcal{F}_{g}^{\mu}<<T^{3}$, where
$\mathcal{F}_{g}^{\mu}=\Gamma^{\mu}_{\alpha\beta}p^{\alpha}p^{\beta}$. With
the above assumptions we are ready to derive the first order correction
$\delta f^{(1)}$ by noting : $f=f^{(0)}+\delta f^{(1)}+\delta f^{(2)}+...$,
where $\delta f^{(n)}$ is obtained from the Boltzmann equation

$\displaystyle p^{\mu}\partial_{\mu}f+qF^{\mu\nu}p_{\nu}\frac{\partial
f}{\partial p^{\mu}}-\Gamma^{\mu}_{\alpha\rho}p^{\alpha}p^{\rho}\frac{\partial
f}{\partial p^{\mu}}=-\frac{u.p}{\tau_{c}}\delta f,$

which leads to

$\delta
f^{(n)}=\sum_{n=0}^{\infty}{\left[-\frac{\tau_{c}}{u.p}\left(p^{\mu}\partial_{\mu}+\mathcal{F}^{\mu}\frac{\partial}{\partial
p^{\mu}}-\mathcal{F}_{g}^{\mu}\frac{\partial}{\partial
p^{\mu}}\right)\right]}^{n}f_{0},$ (16)

where $\mathcal{F}^{\mu}\equiv qF^{\mu\nu}p_{\nu}$. The detailed derivation
and explicit expression for $\delta f^{(1)}$ and $\delta f^{(2)}$ are given in
appendix (B).

Using the expression for $D_{\mu}u_{\alpha}$( see Eq.(50)) in the expression
for $\delta f^{(1)}$ (Eq.(49)) we have

$\displaystyle\delta f^{(1)}$ $\displaystyle=$ $\displaystyle
f_{0}\tilde{f}_{0}\frac{\tau_{c}}{u.p}\left((u.p)p^{\mu}\partial_{\mu}\beta+\beta
p^{\mu}p^{\alpha}\left(u_{\mu}\dot{u}_{\alpha}+\omega_{\mu\alpha}+\sigma_{\mu\alpha}+\frac{\Delta_{\mu\alpha}\theta}{3}\right)-p^{\mu}\partial_{\mu}\alpha\right)$
(17) $\displaystyle+f_{0}\tilde{f}_{0}\frac{\tau_{c}}{u.p}{\beta
qF^{\mu\nu}\delta^{s}_{\mu}u_{s}p_{\nu}}.$

Similarly we can calculate $\delta f^{(2)}$:

$\displaystyle\delta f^{(2)}$ $\displaystyle=$
$\displaystyle-\frac{\tau_{c}}{u.p}\left(p^{\mu}\partial_{\mu}\delta
f^{(1)}+qF^{\mu\nu}p_{\nu}\frac{\partial\delta f^{(1)}}{\partial
p^{\mu}}-\Gamma^{\mu}_{\nu\rho}p^{\nu}p^{\rho}\frac{\partial\delta
f^{(1)}}{\partial p^{\mu}}\right).$ (18)

The explicit expression is quite long and hence given in appendix(B).

###  IV.2 First Order Constitutive relations

Using $\delta f^{(1)}$ we get the first-order constitutive relations for
dissipative stresses as,

$\displaystyle\Pi_{(1)}$ $\displaystyle=$
$\displaystyle-\tau_{c}\beta_{\Pi}\theta,$ (19) $\displaystyle V^{\mu}_{(1)}$
$\displaystyle=$ $\displaystyle\tau_{c}\beta_{V}\left(\nabla^{\mu}\alpha+\beta
qE^{\mu}\right),$ (20) $\displaystyle\pi^{\mu\nu}_{(1)}$ $\displaystyle=$
$\displaystyle 2\tau_{c}\beta_{\pi}\sigma^{\mu\nu}.$ (21)

where
$\beta_{\Pi}=\frac{5\beta}{3}J_{42}^{(1)+}+\mathcal{X}J_{31}^{(0)+}-\mathcal{Y}J_{21}^{(0)-}$
with

$\mathcal{X}=\frac{J_{10}^{(0)+}\left(\epsilon+P\right)-J_{20}^{(0)-}n_{f}}{D_{20}},$

and

$\mathcal{Y}=\frac{J_{20}^{(0)-}\left(\epsilon+P\right)-J_{30}^{(0)+}n_{f}}{D_{20}}.$

, $\beta_{V}=\frac{1}{h}J_{21}^{(0)-}-J_{21}^{(1)-}$ and $\beta_{\pi}=\beta
J_{42}^{(1)+}$. Not surprisingly, the results are similar to that of [31]
apart from the fact that the partial derivatives are now replaced by covariant
derivative.

###  IV.3 Second order evolution equations

In a similar manner using the expression of $\delta f^{(2)}$ we find the
second-order evolution equation,

$\displaystyle\frac{\Pi}{\tau_{c}}$ $\displaystyle=$
$\displaystyle-\dot{\Pi}-\delta_{\Pi\Pi}\Pi\theta+\lambda_{\Pi\pi}\pi^{\mu\nu}\sigma_{\mu\nu}-\tau_{\Pi
V}V\cdot\dot{u}-\lambda_{\Pi V}V\cdot\nabla\alpha-l_{\Pi V}D\cdot
V-\beta_{\Pi}\theta-qB\lambda_{\Pi VB}b^{\mu\beta}V_{\beta}V_{\mu}$ (22)
$\displaystyle+\tau_{c}\tau_{\Pi
VB}\dot{u}_{\alpha}qBb^{\alpha\beta}V_{\beta}-q\delta_{\Pi
VB}\nabla_{\mu}\left(\tau_{c}Bb^{\mu\beta}V_{\beta}\right)-q^{2}\tau_{c}\chi_{\Pi
EE}E^{\mu}E_{\mu},$ $\displaystyle\frac{V^{\mu}}{\tau_{c}}$ $\displaystyle=$
$\displaystyle-\dot{V}^{\langle\mu\rangle}-V_{\nu}\omega^{\nu\mu}+\lambda_{VV}V^{\nu}\sigma^{\mu}_{\nu}-\delta_{VV}V^{\mu}\theta+\lambda_{V\Pi}\Pi\nabla^{\mu}\alpha-\lambda_{V\pi}\pi^{\mu\nu}\nabla_{\nu}\alpha-\tau_{V\pi}\pi^{\mu}_{\nu}\dot{u^{\nu}}-qB\delta_{VB}b^{\mu\gamma}V_{\gamma}$
(23)
$\displaystyle+\tau_{V\Pi}\Pi\dot{u^{\mu}}+l_{V\pi}\Delta^{\mu\nu}D_{\gamma}\pi^{\gamma}_{\nu}-l_{V\Pi}\nabla^{\mu}\Pi+\beta_{V}\nabla^{\mu}\alpha+\tau_{c}qBl_{V\pi
B}b^{\sigma\mu}D^{\kappa}\pi_{\kappa\sigma}-q\tau_{c}\lambda_{VVB}Bb^{\gamma\nu}V_{\nu}\sigma^{\mu}_{\gamma}$
$\displaystyle+\tau_{c}qB\tau_{V\Pi
B}b^{\gamma\mu}\Pi\dot{u}_{\gamma}-\tau_{c}qBl_{V\Pi
B}b^{\gamma\mu}\nabla_{\gamma}\Pi-q\tau_{c}\delta_{VVB}Bb^{\mu\nu}V_{\nu}\theta-q\tau_{c}\mathbf{\rho}_{VVB}Bb^{\gamma\nu}V_{\nu}\omega^{\mu}_{\gamma}$
$\displaystyle+\chi_{VE}qE^{\mu}+q\Delta^{\mu}_{\alpha}\chi_{VE}D\left(\tau_{c}E^{\alpha}\right)-q\tau_{c}\rho_{VE}E^{\mu}\theta-q\tau_{VVB}\Delta^{\mu}_{\gamma}D\left(\tau_{c}Bb^{\gamma\nu}V_{\nu}\right),$
$\displaystyle\frac{\pi^{\mu\nu}}{\tau_{c}}$ $\displaystyle=$
$\displaystyle-\dot{\pi}^{\left<\mu\nu\right>}+2\beta_{\pi}\sigma^{\mu\nu}+2\pi^{\langle\mu}_{\gamma}\omega^{\nu\rangle\gamma}-\tau_{\pi\pi}\pi^{\langle\mu}_{\gamma}\sigma^{\nu\rangle\gamma}-\delta_{\pi\pi}\pi^{\mu\nu}\theta+\lambda_{\pi\Pi}\Pi\sigma^{\mu\nu}-\tau_{\pi
V}V^{\langle\mu}\dot{u}^{\nu\rangle}-\tau_{c}qB\tau_{\pi
VB}\dot{u}^{\langle\mu}b^{\nu\rangle\sigma}V_{\sigma}$ (24)
$\displaystyle+\lambda_{\pi V}V^{\langle\mu}\nabla^{\nu\rangle}\alpha-l_{\pi
V}\nabla^{\langle\mu}V^{\nu\rangle}+\delta_{\pi
B}\Delta^{\mu\nu}_{\eta\beta}qBb^{\gamma\eta}g^{\beta\rho}\pi_{\gamma\rho}-qB\lambda_{\pi
VB}V_{\gamma}b^{\gamma\langle\mu}V^{\nu\rangle}-q\delta_{\pi
VB}\nabla^{\langle\mu}\left(\tau_{c}B^{\nu\rangle\gamma}V_{\gamma}\right)$
$\displaystyle+q^{2}\tau_{c}\chi_{\pi
EE}\Delta^{\mu\nu}_{\sigma\rho}E^{\sigma}E^{\rho}.$

Here the evolution equations are similar to that found in [31] apart from the
fact that the partial derivatives are now replaced by covariant derivatives.
This shows no explicit dependence on the space-time curvature (if present) as
was also inferred in [1] and as was taken in [11] . However, the explicit
expression given in the appendix (B) might be helpful for numerical
simulations for a given geometry.

##  V Causality and stability for GRMHD

Here we study the causality and stability of non-resistive-MHD formulation in
presence of curved space-time and small metric perturbation. For convenience
we decompose the perturbation in the fluid variables and the space-time metric
in Fourier modes of the form ${\delta\tilde{X}}=\delta
Xe^{ig_{\mu\nu}k^{\mu}x^{\nu}}$, where $\delta X\equiv(\delta\epsilon,\delta
P,\delta u^{\mu}..)$. We decompose the metric as
$g_{\mu\nu}=\eta_{\mu\nu}+\delta\tilde{g}_{\mu\nu}$, where $\eta_{\mu\nu}$ is
the flat space-time Minkowski metric, and $\delta\tilde{g}_{\mu\nu}$ is the
small perturbation. Due to the fact that
$g^{\mu\alpha}g_{\alpha\nu}=\delta^{\mu}_{\nu}$ we have

$\displaystyle\delta[g^{\mu\alpha}g_{\alpha\nu}]$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle\delta\tilde{g}^{\mu\nu}$ $\displaystyle=$
$\displaystyle\delta\tilde{g}_{\alpha\beta}\eta^{\alpha\mu}\eta^{\beta\nu}.$

Before we linearise the conservation equations for fluid, let us discuss the
linearised form of the Einstein equation that is given by :

$\displaystyle\delta\tilde{R}_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}\eta^{\alpha\beta}\delta\tilde{R}_{\alpha\beta}+\Lambda\delta\tilde{g}_{\mu\nu}$
$\displaystyle=$ $\displaystyle-8\pi G\delta\tilde{T}_{\mu\nu},$
$\displaystyle-\partial_{\nu}\delta\tilde{\Gamma}^{\lambda}_{\mu\lambda}+\partial_{\lambda}\delta\tilde{\Gamma}^{\lambda}_{\mu\nu}+\Lambda\delta\tilde{g}_{\mu\nu}$
$\displaystyle=$ $\displaystyle-8\pi G\delta\tilde{T}_{\mu\nu}.$ (25)

We use the explicit form of $\delta\tilde{\Gamma}^{\mu}_{\nu\alpha}$ (the
details are given appendix(53)) in Eq.(25) to derive Eq.(29), Eq.(30).

For the sake of simplicity, without losing generalisation, one can use Lorentz
gauge $\delta\tilde{g}^{\mu\nu}k_{\nu}=0$ along with transverse
$\delta\tilde{g}^{\mu\nu}u_{\nu}=0$, and traceless gauge
$\delta\tilde{g}^{\mu}_{\mu}=0$ to reduce the number of independent variables.
With the following choice for
$k^{\mu}=\left(\omega,ksin\theta,0,kcos\theta\right)$ , and
$u^{\mu}=(1,0,0,0)$ we have the following form for the above mentioned three
gauges:

$\displaystyle\delta\tilde{g}^{\mu t}$ $\displaystyle=$ $\displaystyle 0,$
(26) $\displaystyle\delta\tilde{g}^{\mu x}k_{x}+\delta\tilde{g}^{\mu z}k_{z}$
$\displaystyle=$ $\displaystyle 0,$ (27) $\displaystyle\
\delta\tilde{g}^{tt}-\delta\tilde{g}^{xx}-\delta\tilde{g}^{yy}-\delta\tilde{g}^{zz}$
$\displaystyle=$ $\displaystyle 0.$ (28)

For this case we have $\delta\tilde{g}^{xx}$ and $\delta\tilde{g}^{xy}$ as two
independent components; hence from Eq.(25) we get the following equations for
$\delta\tilde{g}^{xx}$ and $\delta\tilde{g}^{xy}$:

$\displaystyle-\left(\Lambda+\frac{1}{2}\left(\omega^{2}-k^{2}\right)-8\pi
G\left(P_{0}+\frac{B^{2}_{0}}{2}\right)\right)\delta\tilde{g}^{xx}=$
$\displaystyle-8\pi
G\left(c_{s}^{2}\delta\tilde{\epsilon}+\delta\tilde{\Pi}+\delta\tilde{b}^{z}B^{2}_{0}-\delta\tilde{g}^{xx}\frac{(k_{x})^{2}}{(k_{z})^{2}}\frac{B^{2}_{0}}{2}\right)$
(29) $\displaystyle\left(\Lambda+\frac{\left(\omega^{2}-k^{2}\right)}{2}-8\pi
G\left(P_{0}+\frac{B^{2}_{0}}{2}\right)\right)\delta\tilde{g}^{xy}=0.$ (30)

Here we consider only bulk viscosity in the $\delta\tilde{T}^{\mu\nu}$. Next,
we linearise the Maxwell’s equation $D_{\mu}\tilde{F}^{\mu\nu}=0$ by taking
$b^{\mu}\equiv B^{\mu}/B=(0,0,0,1)$ . Also note that we have considered
$\Pi_{0}=0$.

The resulting linearised equations are :

$\displaystyle ik_{x}\delta\tilde{b}^{x}+ik_{z}\delta\tilde{b}^{z}$
$\displaystyle=$ $\displaystyle 0,$ (31) $\displaystyle
ik_{z}\delta\tilde{u}^{x}-i\omega\delta\tilde{b}^{x}$ $\displaystyle=$
$\displaystyle 0,$ (32) $\displaystyle
ik_{z}\delta\tilde{u}^{y}-i\omega\delta\tilde{b}^{y}$ $\displaystyle=$
$\displaystyle 0,$ (33)
$\displaystyle-i\omega\delta\tilde{b}^{z}-ik_{x}\delta\tilde{u}^{x}$
$\displaystyle=$ $\displaystyle 0.$ (34)

Next, we linearise the energy-momentum conservation equation for the fluid and
field:

$\displaystyle\partial_{\mu}\delta\tilde{T}^{\mu\nu}+\delta\tilde{\Gamma}^{\mu}_{\mu\alpha}T^{\alpha\nu}_{0}+\delta\tilde{\Gamma}^{\nu}_{\mu\alpha}T^{\mu\alpha}_{0}$
$\displaystyle=$ $\displaystyle 0$ (35)

The resulting linearised equations are :

$\displaystyle
i\omega\delta\tilde{\epsilon}-i\omega\delta\tilde{\Pi}+ik_{x}\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{x}+ik_{z}\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{z}=0,$
(36) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}+B^{2}_{0}-\frac{k_{x}^{2}}{\omega^{2}}B^{2}_{0}\right)\delta\tilde{u}^{x}+ik_{x}c_{s}^{2}\delta\tilde{\epsilon}+ik_{x}\delta\tilde{\Pi}-ik_{x}\frac{k_{x}^{2}+k_{z}^{2}}{k_{z}^{2}}\delta
B_{0}^{2}\tilde{g}^{xx}-ik_{z}B^{2}_{0}\delta\tilde{b}^{x}=0,$ (37)
$\displaystyle
i\omega\left(\epsilon_{0}+P_{0}+B^{2}_{0}\right)\delta\tilde{u}^{y}-ik_{z}B^{2}_{0}\delta\tilde{b}^{y}-ik_{x}B^{2}_{0}\delta\tilde{g}^{xy}=0,$
(38) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{z}+ik_{z}c_{s}^{2}\delta\tilde{\epsilon}+ik_{z}\delta\tilde{\Pi}+\frac{3ik_{z}k_{x}}{\omega}B^{2}_{0}\delta\tilde{u}^{x}+\frac{ik_{x}^{2}}{k_{z}}B^{2}_{0}\delta\tilde{g}^{xx}-ik_{x}B^{2}_{0}\delta\tilde{b}^{x}=0.$
(39)

We note that bulk stress in the present formulation is an independent
dynamical variable; linearising the bulk evolution equation we get :

$\displaystyle\delta\tilde{\Pi}\left(\frac{1}{\tau_{c}}+i\omega\right)+\beta_{\Pi}ik_{x}\delta\tilde{u}^{x}+\beta_{\Pi}ik_{z}\delta\tilde{u}^{z}=0.$
(40)

Now we have all the ingredients to derive the dispersion relations and study
the causality and stability of the system under consideration. For clarity, we
consider three different cases one by one, as discussed below.

###  V.1 Non-viscous fluid with Gravity

The first case we study is for non-viscous fluid with metric perturbation.
Here the independent variables are
$X\equiv\left(\delta\tilde{\epsilon},\delta\tilde{u}^{x},\delta\tilde{u}^{y},\delta\tilde{u}^{z},\delta\tilde{g}^{xx},\delta\tilde{g}^{xy}\right)^{T}$,
the superscript $T$ denotes the transpose. Using Eqs.((29),(30)) and
Eqs.(36,37,38,39) we obtain $\mathcal{A}X=0$ where the matrix $\mathcal{A}$ is
given below. We set $\Pi=0$, and $B_{0}=0$ while constructing $\mathcal{A}$.

$\displaystyle\mathcal{A}=\left[{\begin{array}[]{ccccccc}\begin{smallmatrix}i\omega&(\epsilon_{0}+P_{0})ik_{x}&0&(\epsilon_{0}+P_{0})ik_{z}&0&0\\\
ik_{x}c_{s}^{2}&(\epsilon_{0}+P_{0})i\omega&0&0&0&0\\\
0&0&(\epsilon_{0}+P_{0})i\omega&0&0&0\\\
ik_{z}c_{s}^{2}&0&0&(\epsilon_{0}+P_{0})i\omega&0&0\\\ 8\pi
Gc_{s}^{2}&0&0&0&\mathcal{M}&0\\\ 0&0&0&0&0&\mathcal{M}\\\
\end{smallmatrix}\end{array}}\right].$

Here we defined $\mathcal{M}=-\left(\Lambda+\frac{1}{2}(\omega^{2}-k^{2})-8\pi
GP_{0}\right)$ for brevity. Setting det($\mathcal{A}$)=0 we obtain the
dispersion relations. A careful analysis shows that in this case we have two
different modes of wave propagation namely (i) sound modes with the dispersion
relation $\omega=\pm c_{s}\sqrt{k_{x}^{2}+k_{z}^{2}}$, and (ii) the
gravitational modes with the dispersion relation
$\omega=\pm\sqrt{k^{2}-2\Lambda+16G\pi P_{0}}$ , where
$k=\sqrt{k_{x}^{2}+k_{z}^{2}}$. As expected the group velocity for the sound
mode is $v_{g}=\frac{\partial\omega}{\partial k}=\pm c_{s}$, whereas for the
gravitation mode $v_{g}=\frac{\partial\omega}{\partial
k}=\pm\frac{k}{\sqrt{k^{2}-2\Lambda+16G\pi P_{0}}}$ . Here we observe a few
points from the group velocity of the gravitational mode: since the term
inside the square root in the denominator may become negative for
$k<\pm\sqrt{(2\Lambda-16\pi GP_{0})}$; we do not have any propagating modes
for $k$ smaller than this cutoff and one of the roots will give rise to a
unstable mode. We also note that for $k$ just above this cutoff value may
results in $v_{g}>1.$ However, note that we should actually consider
$\lim_{k\rightarrow\infty}\frac{d\omega}{dk}=1$ for checking the causality
[23]; which for this case shows that the gravitational wave travels with speed
of light in vacuum. These are expected results as it is known that in non-
viscous fluid the disturbance propagates with the speed of sound $c_{s}$ and
the gravitational waves propagates with the speed of light. Now let us turn to
the second case where we include bulk viscosity in the fluid.

###  V.2  Bulk viscosity with Gravity

When we add bulk viscosity to the fluid we have the following set of
independent variables
$X\equiv\left(\delta\tilde{\epsilon},\delta\tilde{u}^{x},\delta\tilde{u}^{y},\delta\tilde{u}^{z},\delta\tilde{\Pi},\delta\tilde{g}^{xx},\delta\tilde{g}^{xy}\right)^{T}$.
In this case, the matrix $\mathcal{A}$ takes the following form:

$\mathcal{A}=\left[{\begin{array}[]{cccccccc}\begin{smallmatrix}i\omega&(\epsilon_{0}+P_{0})ik_{x}&0&(\epsilon_{0}+P_{0})ik_{z}&-i\omega&0&0\\\
ik_{x}c_{s}^{2}&(\epsilon_{0}+P_{0})i\omega&0&0&ik_{x}&0&0\\\
0&0&(\epsilon_{0}+P_{0})i\omega&0&0&0&0\\\
ik_{z}c_{s}^{2}&0&0&(\epsilon_{0}+P_{0})i\omega&ik_{z}&0&0\\\
0&ik_{x}\beta_{\Pi}&0&ik_{z}\beta_{\Pi}&(\frac{1}{\tau_{c}}+i\omega)&0&0\\\
8\pi Gc_{s}^{2}&0&0&0&8\pi G&\mathcal{M}&0\\\ 0&0&0&0&0&0&\mathcal{M}\\\
\end{smallmatrix}\end{array}}\right].$ (42)

Here $\mathcal{M}=-\left(\Lambda+\frac{1}{2}(\omega^{2}-k^{2})-8\pi
GP_{0}\right)$. Like the previous case here also we notice that the
gravitational modes are decoupled from the hydrodynamic modes. The dispersion
relation for the gravitational waves are
$\omega=\pm\sqrt{k^{2}-2\Lambda+16G\pi P_{0}}$. The full dispersion relations
for the hydro modes are quite long and hence we show here only the small and
large $k$ limits.
For small $k$ limit we have :
$\omega=\begin{cases}\frac{i}{\tau_{c}}+\sqrt{\frac{\beta_{\Pi}}{2\epsilon_{0}}}k+\mathcal{O}(k^{2}),\\\
\pm c_{s}k+\mathcal{O}(k^{2}).\end{cases}$

Here the first expression corresponds to a non-hydro mode and the second
expression corresponds to the usual sound mode with speed of sound $c_{s}$. On
the other hand, for the large k limit we have $\omega=\pm v_{g}k$ where
$v_{g}=\pm\sqrt{\frac{\beta_{\Pi}+P_{0}}{\epsilon_{0}}}$. Thus for causality
to hold we must have $\beta_{\Pi}\leq(1-c_{s}^{2})\epsilon_{0}$. Since for
this case $\Im{(\omega)}>0$ all modes are stable. For the gravitational
sector, we have exactly the same conclusion as was derived in the previous
sub-section(V.1).

###  V.3  Bulk viscosity , magnetic field, and Gravity

Finally, in this section, we consider a fluid with finite bulk viscosity in a
magnetic field and we allow small metric perturbation. A straightforward
extension of the previous cases give the following characteristic matrix:

$\mathcal{A}=\left[{\begin{array}[]{cccccccccc}\begin{smallmatrix}i\omega&(\epsilon_{0}+P_{0})ik_{x}&0&(\epsilon_{0}+P_{0})ik_{z}&-i\omega&0&0&0&0\\\
ik_{x}c_{s}^{2}&(\epsilon_{0}+P_{0}+B_{0}^{2})i\omega-i\frac{k_{x}^{2}}{\omega}B_{0}^{2}&0&0&ik_{x}&-ik_{z}B_{0}^{2}&0&-ik_{x}B_{0}^{2}\frac{k_{x}^{2}+k_{z}^{2}}{k_{z}^{2}}&0\\\
0&0&i\omega(\epsilon_{0}+P_{0}+B_{0}^{2})&0&0&0&-ik_{z}B^{2}_{0}&0&-ik_{x}B^{2}_{0}\\\
ik_{z}c_{s}^{2}&i\frac{3k_{x}k_{z}}{\omega}B_{0}^{2}&0&i\omega(\epsilon_{0}+P_{0})&ik_{z}&-ik_{x}B^{2}_{0}&0&-ik_{z}B^{2}_{0}+ik_{z}B_{0}^{2}\frac{k^{2}}{k_{z}^{2}}&0\\\
0&ik_{x}\beta_{\Pi}&0&ik_{z}\beta_{\Pi}&(\frac{1}{\tau_{c}}+i\omega)&0&0&0&0\\\
0&ik_{z}&0&0&0&-i\omega&0&0&0\\\ 0&0&ik_{z}&0&0&0&-i\omega&0&0\\\ 8\pi
Gc_{s}^{2}&\frac{-8\pi Gk_{x}}{\omega}B_{0}^{2}&0&0&8\pi
G&0&0&\mathcal{N}&0\\\ 0&0&0&0&0&0&0&0&\mathcal{S}\\\
\end{smallmatrix}\end{array}}\right],$

where $\mathcal{N}=-\left(\Lambda+\frac{\omega^{2}-k^{2}}{2}-8\pi
G\left(P_{0}+B_{0}^{2}-\frac{B_{0}^{2}}{2}\frac{k^{2}}{k_{z}^{2}}\right)\right)$,
$\mathcal{S}=-\left(\Lambda+\frac{\omega^{2}-k^{2}}{2}-8\pi GP_{0}-8\pi
GB_{0}^{2}/2\right)$, for this case the vector $X$ is
$\left(\delta\tilde{\epsilon},\delta\tilde{u}^{x},\delta\tilde{u}^{y},\delta\tilde{u}^{z},\delta\tilde{\Pi},\delta\tilde{b}^{x},\delta\tilde{b}^{y},\delta\tilde{g}^{xx},\delta\tilde{g}^{xy}\right)^{T}$.
We checked the dispersion relation and found that the usual Alfven modes due
to finite magnetic field has the following form
$\omega=\pm\frac{B_{0}k_{z}}{\sqrt{B_{0}^{2}+\epsilon_{0}+P_{0}}}$ , where
$k_{z}=kcos\theta$, $\theta$ is the angle between $\vec{k}$ and $\vec{B}$. The
corresponding group velocity for the Alfven wave is
$v_{g}=\pm\frac{B_{0}cos\theta}{\sqrt{B_{0}^{2}+\epsilon_{0}+P_{0}}}=\pm
v_{A}cos\theta$, it is clear that the Alfven waves are unaffected from the
metric perturbation. Interestingly, the gravitational sector shows explicit
dependence on the magnetic field; the dispersion relation for gravitational
modes are :

$\omega=\pm\sqrt{k^{2}+8B_{0}^{2}G\pi-2\Lambda+16G\pi P_{0}}.$ (43)

Corresponding group velocity is
$v_{g}=\pm\frac{k}{\sqrt{k^{2}+8B_{0}^{2}G\pi-2\Lambda+16G\pi P_{0}}}$. As it
is clear form the expression for $v_{g}$ that it may become purely imaginary
if the denominator becomes negative, which corresponds to the cutoff
$k_{c}=\pm\sqrt{-8B_{0}^{2}G\pi+2\Lambda-16G\pi P_{0}}$. Here we observe that
instability arises for $k\leq k_{c}$ and the magnetic field plays an important
role. Also, It seems that zones with varying magnetic fields gives rise to
dispersive gravitational waves.

The different hydrodynamical modes in the small $k$ limits are the following:
For small k limit we have :
$\omega=\begin{cases}\frac{i}{\tau_{c}}+\mathcal{O}(k^{2}),\\\
\pm\sqrt{2}\sqrt{4B_{0}^{2}\pi G-\Lambda+8\pi GP_{0}-4B_{0}^{2}\pi
G\tan^{2}{\theta}}+\mathcal{O}(k^{2}),\\\ \pm c_{s}k,\text{for }\theta=0,\\\
\pm v_{A}k,\text{for }\theta=0,\\\ \pm v_{M}k,\text{for
}\theta=\frac{\pi}{4}.\end{cases}$ where the expression for the magneto-sonic
wave $v_{M}$ is cumbersome and hence is not given here (see Eq.(76) in
appendix (C)), we also note that $v_{M}$ diverges for $\theta=\pi/2$. Here we
have found additional non-hydro modes arising due to gravity; these
gravitational modes are also unstable when
$tan^{2}\theta\geq\pi\left(1+\frac{2P_{0}}{B_{0}^{2}}\right)-\frac{\Lambda}{4B_{0}^{2}G}$.

In the large $k$ limit we have the following expression for $v_{g}$:

$\displaystyle v_{g}^{2}$ $\displaystyle=$
$\displaystyle\left[\frac{2v_{A}^{2}\beta_{\Pi}}{4\epsilon_{0}}+\frac{\epsilon_{0}}{4\epsilon_{0}(\epsilon_{0}+P_{0}+B_{0}^{2})}(B^{2}_{0}(2+2c_{s}^{2}cos^{2}\theta)+2(1+c_{s}^{2})(\beta_{\Pi}+P_{0}))\right]$
$\displaystyle\pm\frac{1}{4\epsilon_{0}(\epsilon_{0}+P_{0}+B_{0}^{2})}\sqrt{\left[\left[2\beta_{\Pi}\left(B^{2}_{0}cos^{2}\theta+\epsilon_{0}+P_{0}\right)+\epsilon_{0}\left(B^{2}_{0}(2+2c_{s}^{2}cos^{2}\theta)+2(P_{0}+c_{s}^{2}P_{0})\right)\right]^{2}\right]}$
$\displaystyle\pm\frac{1}{4\epsilon_{0}(\epsilon_{0}+P_{0}+B_{0}^{2})}\sqrt{\left[16B^{2}_{0}(-2+cos2\theta)cos^{2}\theta\epsilon_{0}(\beta_{\Pi}+P_{0})(B_{0}^{2}+\epsilon_{0}+P_{0})\right]}.$

Here in the large $k$ limit we have found the fast and slow magneto-sonic
modes with causality bound as a function of $\theta$ and $B_{0}$. For
consistency check, we set $\theta=0$ which gives
$v_{g}^{2}=\begin{cases}v_{A}^{2},\text{ for slow magneto-sonic},\\\
c_{s}^{2}+\frac{\beta_{\Pi}}{\epsilon_{0}}.\text{ for fast magneto-
sonic}.\end{cases}$. These expressions are identical to the ones found
in.[27]. Clearly, for $\theta=0$, the causality holds if $v_{A}^{2}\leq 1$ and
$\beta_{\Pi}\leq(1-c_{s}^{2})\epsilon_{0}$. Similarly for
$\theta=\frac{\pi}{2}$; $v_{g}^{2}=\begin{cases}0,\text{ for slow magneto-
sonic},\\\
\frac{B_{0}^{2}+(1+c_{s}^{2})(\beta_{\Pi}+P_{0})}{B^{2}_{0}+\epsilon_{0}+P_{0}},\text{for
fast magneto-sonic}.\end{cases}$
In this case, the causality requires
$c_{s}^{2}\leq\frac{\epsilon_{0}-\beta_{\Pi}}{P_{0}+\beta_{\Pi}}$.

Figure 1: (Color online) Top panel: Gravitational modes $\omega^{\pm}$ vs $k$
for different pressure values (0,0.5,1.0) (in unit of $m_{\pi}^{4}$) for
$B_{0}=1$ (in unit of $m_{\pi}^{2}$) , G=1 , $\Lambda=0$, the one passing
through the origin corresponds to ($B_{0}=P_{0}=0$). Bottom panel: Polar plot
of the non-hydro modes due to gravity for different values of magnetic field
(in unit of $m_{\pi}^{2}$) for $G=1,\Lambda=0,P_{0}=1$.

Top panel of Fig.(1) shows the dispersion relation for the gravitational modes
Eq.(43) for different values of pressure and magnetic field. The gravitational
waves travel on the light cone for $B_{0}=P_{0}=0$. But for non-zero $B_{0}$
and/or $P_{0}$ non-hydrodynamical modes arise as can be seen from Fig.(1) (top
panel), the blue, yellow, and green lines correspond to $B_{0}=1$,
$P_{0}=0-1.0$ in steps of 0.5 respectively. The bottom panel of Fig.(1) is the
polar plot of $\sqrt{2}\sqrt{4B_{0}^{2}\pi G-\Lambda+8\pi GP_{0}-4B_{0}^{2}\pi
G\tan^{2}{\theta}}$ (the non-hydro modes due to gravity) as a function of
$\theta$. As can be seen from the figure that in the small $k$ limit, for
$B_{0}$=0, $\omega$ remains real $\forall$ $\theta$. Whereas, if we have
finite magnetic field, then $\omega$ becomes complex for some ranges of
$\theta$, and this angular dependence is a function of $B_{0}$; for example
the yellow and green lines correspond to $B_{0}=1,2$ (in unit of
$m_{\pi}^{2}$).

##  VI Summary and conclusion

In this work, we formulate the first and second-order viscous
magnetohydrodynamics for an arbitrary smooth (no-singularity) space-time
geometry in the ideal-MHD limit (large magnetic Reynolds number). We found no
explicit curvature dependence in the evolution of dissipative stresses as was
found for the holographic calculation at finite temperature [1] (without
electromagnetic fields). This might not be surprising as it was mentioned in
ref [1] that for ordinary (i.e., without electromagnetic fields) relativistic
hydrodynamics derived from kinetic theory using RTA approximation should not
contain any explicit curvature dependence because the collision integral only
allows symmetric tensors (for more details see ref [1]). We further
investigate wave propagation within the framework of linearised perturbation
of fluid, field, and space-time metric. In this approach, for the case of a
non-viscous fluid, we get the normal sound mode and the gravitational modes,
which in the large $k$ limit gives rise to causal propagation of
perturbations. We show that one of these gravitational modes leads to
instability. After adding the bulk viscosity, we get the usual gravitational
modes, which were found in the non-viscous case, with the causality and
stability conditions remaining the same.

In the small $k$ limit, we found the usual non-hydro and the sound modes
(arising due to bulk viscosity). All the hydrodynamic modes for this case are
stable. On the other hand, for the large k limit, we found the causality
condition to be $\beta_{\Pi}\leq(1-c_{s}^{2})\epsilon_{0}$. Lastly, the usual
Alfven waves appear along with the magnetic field-dependent gravitational
waves when we add the magnetic field. For this case, instability arises when
$k$ becomes less than a cutoff value
$k_{c}=\pm\sqrt{-8B_{0}^{2}G\pi+2\Lambda-16G\pi P_{0}}$. In the small $k$
limit, we found the usual non-hydro mode (due to the bulk viscosity), sound
modes, and the Alfven modes. Interestingly, the magneto-sonic modes turned out
to be dependent on the gravitational field and the angle between $\vec{k}$ and
$\vec{B}$. We found additional non-hydro modes due to gravity which blows up
at $\theta=\pi/2$. Whereas, in the large $k$ limit, only slow and fast
magneto-sonic modes appear for which we have shown the corresponding causality
conditions. 111While finishing our work we came across a similar study by Noh
et al [41] where they study the gravito-magnetic instabilities of a static
homogeneous medium with magnetic field in the contexts of MHD with post-
Newtonian (PN) corrections, and, special relativistic (SR) MHD with weak
gravity. In future, we plan to investigate the causality using real space
analysis which takes into account non-linear effects.

###### Acknowledgements.

We would like to acknowledge Rajesh Biswas for his valuable inputs. AP
acknowledges the CSIR-HRDG financial support. VR acknowledge support from the
DAE, Govt. of India and financial support from the DST Inspire faculty
research grant (IFA-16-PH-167), India.

##  Appendix A Useful definitions

In this appendix we give details of a few definitions used in the main text.
The relativistic thermodynamic integrals are defined as:

$\displaystyle J^{(m)\pm}_{nq}$ $\displaystyle=$
$\displaystyle\frac{1}{(2q+1)!!}\int\frac{\mathcal{G}_{s}d^{3}p}{(2\pi)^{3}p_{0}}\sqrt{-g}(g_{\mu\nu}u^{\mu}p^{\nu})^{n-2q-m}\left(\Delta_{\alpha\beta}p^{\alpha}p^{\beta}\right)^{q}\left(f_{0}\tilde{f}_{0}\pm\bar{f}_{0}\tilde{\bar{f}}_{0}\right),$

here $\mathcal{G}_{s}$ is the spin degeneracy, where $n,q$ are integers with
$n\geq 1$ and $q\geq 0$. The $n$-th moments integral for the distribution
function is defined as:

$\displaystyle
I_{\mu_{1}\mu_{2}\cdots\mu_{n}}^{(m)\pm}=\int\frac{\mathcal{G}_{s}d^{3}p}{(2\pi)^{3}p_{0}\left(u\cdot
p\right)^{m}}\sqrt{-g}p_{\mu_{1}}p_{\mu_{2}}\cdots
p_{\mu_{n}}\left(f_{0}\pm\bar{f_{0}}\right),$ (44)

which can be docomposed as:

$\displaystyle
I_{\mu_{1}\mu_{2}\dots\mu_{n}}^{(m)\pm}=I_{n0}^{(m)\pm}u_{\mu_{1}}\cdots
u_{\mu_{n}}+I_{n1}^{(m)\pm}\left(\Delta_{\mu_{1}\mu_{2}}u_{\mu_{3}}\cdots
u_{\mu_{n}}+\text{perm.}\right)+\cdots$
$\displaystyle\cdots+I_{nq}^{(m)\pm}\left(\Delta_{\mu_{1}\mu_{2}}\Delta_{\mu_{3}\mu_{4}}\cdots\Delta_{\mu_{n-1}\mu_{n}}+\text{perm.}\right).$
(45)

where $n\geq 2q$.

We also note that $J$ is connected to the I integral via the following
relation ,

$J^{(0)\pm}_{nq}=\frac{1}{\beta}\left[-I^{(0)\pm}_{n-1,q-1}+(n-2q)I^{(0)\pm}_{n-1,q}\right].$
(46)

Equilibrium number density , energy density and pressure are given as
$I_{10}^{(0)-}$, $I_{20}^{(0)+}$, $-I_{21}^{(0)+}$ respectively. Next we
discuss the explicit form of $\delta f^{(1)}$ given in Eq.(17). For that we
need the following derivatives:

$\displaystyle\partial_{\mu}f_{0}$ $\displaystyle=$ $\displaystyle-
f_{0}\tilde{f}_{0}\left(g_{\alpha\beta}u^{\alpha}p^{\beta}\partial_{\mu}\beta+\beta\partial_{\mu}\left(g_{\alpha\beta}u^{\alpha}p^{\beta}\right)-\partial_{\mu}\alpha\right),$
(47) $\displaystyle\frac{\partial f_{0}}{\partial p^{\mu}}$ $\displaystyle=$
$\displaystyle-f_{0}\tilde{f}_{0}\beta u_{\sigma}\delta^{\sigma}_{\mu}.$ (48)

Here $\delta^{s}_{\mu}$ is the Kronecker-Delta.We have also used the following
formulas in the main text as and whenever required:

$\displaystyle\partial_{\mu}(\sqrt{-g}p^{\mu}f)$ $\displaystyle=$
$\displaystyle\sqrt{-g}p^{\mu}\partial_{\mu}f+\sqrt{-g}p^{\mu}\Gamma^{\alpha}_{\alpha\mu}f,$
$\displaystyle\partial_{\mu}\left(2\theta(p^{0})\delta(p^{2}-m^{2})\right)$
$\displaystyle=$ $\displaystyle
2\theta(p^{0})\delta^{{}^{\prime}}(p^{2}-m^{2})2p^{\alpha}p_{\beta}\Gamma^{\beta}_{\alpha\mu},$
$\displaystyle 2p_{\beta}\delta^{{}^{\prime}}(p^{2}-m^{2})$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial p^{\beta}}\delta(p^{2}-m^{2}).$

##  Appendix B $\delta f$ calculation

In this appendix we give the detail derivation of first and second order
corrections to the single particle distribution functions used in the main
text. First, let us discuss $\delta f^{(1)}$; using Eq.(16)), Eq.(47), and
Eq.(48) we have the following expression for $\delta f^{(1)}$:

$\delta
f^{(1)}=f_{0}\tilde{f}_{0}\frac{\tau_{c}}{u.p}\left((u.p)p^{\mu}\partial_{\mu}\beta+\beta
p^{\beta}p^{\mu}D_{\mu}u_{\beta}-p^{\mu}\partial_{\mu}\alpha\right)+f_{0}\tilde{f}_{0}\frac{\tau_{c}}{u.p}{\beta
qF^{\mu\nu}p_{\nu}\delta^{\sigma}_{\mu}u_{\sigma}}.$ (49)

The covariant derivative of $D_{\mu}u_{\beta}$ in the above expression can be
decomposed in the following manner

$D_{\mu}u_{\beta}=u_{\mu}\dot{u}_{\beta}+\nabla_{\mu}u_{\beta},$ (50)

where $\dot{u}_{\beta}$ is the proper acceleration of fluid element. In a
similar manner we can calculate the second-order $\delta f$ from Eq.(16)),
that has the following form:

$\delta f^{(2)}=-\frac{\tau_{c}}{u.p}\left(p^{\alpha}\partial_{\alpha}\delta
f^{(1)}+qF^{\alpha\beta}p_{\beta}\frac{\partial\delta f^{(1)}}{\partial
p^{\alpha}}-\Gamma^{\alpha}_{\nu\rho}p^{\nu}p^{\rho}\frac{\partial\delta
f^{(1)}}{\partial p^{\alpha}}\right).$ (51)

Using Eq.(49), Eq.(50) in Eq.(51) we have:

$\displaystyle\delta f^{(2)}$ $\displaystyle=$
$\displaystyle-\frac{\tau_{c}}{u.p}p^{\alpha}\left(\tau_{c}(\partial_{\mu}\beta)p^{\mu}\left[-f_{0}\tilde{f}_{0}\left(\tilde{f}_{0}-rf_{0}\right)\left((u.p)\partial_{\alpha}\beta+\beta
p^{\beta}\left(D_{\alpha}u_{\beta}\right)-\partial_{\alpha}\alpha\right)\right]\right)$
$\displaystyle-\frac{\tau_{c}^{2}}{(u.p)^{2}}p^{\alpha}(\partial_{\mu}\alpha)p^{\mu}f_{0}\tilde{f}_{0}\left(\left[\left(\tilde{f}_{0}-rf_{0}\right)\left((u.p)\partial_{\alpha}\beta+\beta
p^{\beta}\left(D_{\alpha}u_{\beta}\right)-\partial_{\alpha}\alpha\right)\right]+\frac{1}{(u.p)}p^{\beta}\left(D_{\alpha}u_{\beta}\right)\right)$
$\displaystyle-\frac{\tau_{c}}{u.p}p^{\alpha}f_{0}\tilde{f}_{0}\left(p^{\mu}D_{\alpha}\left(\tau_{c}\partial_{\mu}\beta\right)+\tau_{c}\beta
qE^{\nu}p_{\nu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\left((u.p)\partial_{\alpha}\beta+\beta
p^{\beta}\left(D_{\alpha}u_{\beta}\right)-\partial_{\alpha}\alpha\right)}{u.p}\right]\right)$
$\displaystyle+\frac{\tau_{c}^{2}\beta}{(u.p)^{2}}p^{\alpha}f_{0}\tilde{f}_{0}\left(-qE^{\nu}p_{\nu}\frac{1}{(u.p)}p^{\beta}\left(D_{\alpha}u_{\beta}\right)+(D_{\mu}u_{\beta})p^{\beta}p^{\mu}\left[\left(\tilde{f}_{0}-rf_{0}\right)\left((u.p)\partial_{\alpha}\beta+\beta
p^{\beta}\left(D_{\alpha}u_{\beta}\right)-\partial_{\alpha}\alpha\right)\right]\right)$
$\displaystyle+\frac{\tau_{c}}{(u.p)^{2}}p^{\alpha}p^{\mu}f_{0}\tilde{f}_{0}\left(\beta\tau_{c}\left(D_{\mu}u_{\beta}\right)p^{\beta}\frac{1}{(u.p)}p^{\nu}\left(D_{\alpha}u_{\nu}\right)+D_{\alpha}\left(\tau_{c}\partial_{\mu}\alpha\right)+qD_{\alpha}\left[\tau_{c}\beta
E_{\mu}\right]\right)$
$\displaystyle-\frac{\tau_{c}f_{0}\tilde{f}_{0}p^{\alpha}p^{\beta}p^{\mu}}{(u.p)^{2}}\left[\left(D_{\mu}u_{\beta}\right)\partial_{\alpha}\left(\beta\tau_{c}\right)+\beta\tau_{c}D_{\alpha}\left(D_{\mu}u_{\beta}\right)\right]$
$\displaystyle-\frac{\tau_{c}^{2}}{(u.p)}qF^{\alpha\sigma}p_{\sigma}\left(\left[f_{0}\tilde{f}_{0}\partial_{\alpha}\beta-p^{\mu}(\partial_{\mu}\beta)f_{0}\tilde{f}_{0}\beta
u_{\alpha}\left(\tilde{f}_{0}-rf_{0}\right)\right]+\beta\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\left(p^{\mu}D_{\alpha}u_{\mu}+p^{\beta}D_{\alpha}u_{\beta}\right)\right]\right)$
$\displaystyle+\frac{\beta\tau_{c}^{2}}{(u.p)}qF^{\alpha\sigma}p_{\sigma}p^{\beta}p^{\mu}\left(\sigma_{\mu\beta}+\omega_{\mu\beta}+\frac{\Delta_{\mu\beta}\theta}{3}+u_{\mu}\dot{u}_{\beta}\right)\left[\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(\tilde{f}_{0}-rf_{0}\right)}{u.p}\right]$
$\displaystyle+\frac{\tau_{c}^{2}}{(u.p)}qF^{\alpha\sigma}p_{\sigma}\left(\partial_{\mu}\alpha\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\delta^{\mu}_{\alpha}-p^{\mu}\left[\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(\tilde{f}_{0}-rf_{0}\right)}{u.p}\right]\right]+\beta
qE_{\alpha}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\right]\right)$
$\displaystyle-\frac{\beta
q^{2}\tau_{c}^{2}}{(u.p)}F^{\alpha\sigma}p_{\sigma}E^{\nu}p_{\nu}\left[\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(\tilde{f}_{0}-rf_{0}\right)}{u.p}\right].$

###  B.1 Shear viscosity

We decompose the second-order expression for the shear stress for the sake of
convenient as:

$\displaystyle\pi^{\lambda\chi}_{(2)}$ $\displaystyle=$ $\displaystyle
I+II+III,$

where :

$\displaystyle I$ $\displaystyle=$
$\displaystyle-\frac{\Delta^{\lambda\chi}_{\alpha\beta}}{3}\tau_{c}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2})\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\partial_{\kappa}\delta
f^{(1)},$ (52) $\displaystyle II$ $\displaystyle=$
$\displaystyle-\frac{\Delta^{\lambda\chi}_{\alpha\beta}}{3}\tau_{c}qF^{\mu}_{\quad\nu}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2})\frac{p^{\alpha}p^{\beta}p^{\nu}}{(u.p)}\frac{\partial\delta
f^{(1)}}{\partial p^{\mu}},$ $\displaystyle III$ $\displaystyle=$
$\displaystyle\frac{\Delta^{\lambda\chi}_{\alpha\beta}}{3}\tau_{c}\Gamma^{\mu}_{\nu\rho}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2})\frac{p^{\alpha}p^{\beta}p^{\nu}p^{\rho}}{(u.p)}\frac{\partial\delta
f^{(1)}}{\partial p^{\mu}}.$

Let us evaluate each term one-by-one, for convenience let us define
$\mathcal{I}^{\lambda\chi}_{\alpha\beta}\equiv\frac{\Delta^{\lambda\chi}_{\alpha\beta}}{3}\tau_{c}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2})$,
then we have

$\displaystyle I$ $\displaystyle=$
$\displaystyle-\mathcal{I}^{\lambda\chi}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\left(\tau_{c}\partial_{\mu}\beta
p^{\mu}\left[-f_{0}\tilde{f}_{0}\left(\tilde{f}_{0}-rf_{0}\right)\left((u.p)\partial_{\kappa}\beta+\beta
p^{\beta}\left(D_{\kappa}u_{\beta}+\Gamma^{\nu}_{\kappa\beta}u_{\nu}\right)-\partial_{\kappa}\alpha\right)\right]\right)$
$\displaystyle+\mathcal{I}^{\lambda\chi}_{\alpha\beta}\left(\tau_{c}(\partial_{\mu}\alpha)p^{\mu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\kappa}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\kappa}u_{\beta}+\Gamma^{\nu}_{\kappa\beta}u_{\nu}\right)\right]\right)$
$\displaystyle-\mathcal{I}^{\lambda\chi}_{\alpha\beta}\left(-\tau_{c}\beta
qE^{\nu}p_{\nu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\kappa}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\kappa}u_{\beta}+\Gamma^{\nu}_{\kappa\beta}u_{\nu}\right)\right]+f_{0}\tilde{f}_{0}p^{\mu}\partial_{\kappa}\left(\tau_{c}\partial_{\mu}\beta\right)\right)$
$\displaystyle-\mathcal{I}^{\lambda\chi}_{\alpha\beta}\left(\beta\tau_{c}(D_{\mu}u_{\beta})p^{\beta}p^{\mu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\kappa}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\kappa}u_{\beta}+\Gamma^{\nu}_{\kappa\beta}u_{\nu}\right)\right]\right)$
$\displaystyle+\mathcal{I}^{\lambda\chi}_{\alpha\beta}\left(\left(\frac{f_{0}\tilde{f}_{0}p^{\mu}}{u.p}\right)\partial_{\kappa}\left(\tau_{c}\partial_{\mu}\alpha\right)+q\left(f_{0}\tilde{f}_{0}\frac{p_{\nu}}{u.p}\right)\left[(\tau_{c}\beta)\partial_{\kappa}E^{\nu}+E^{\nu}\partial_{\kappa}(\tau_{c}\beta)\right]\right)$
$\displaystyle-\mathcal{I}^{\lambda\chi}_{\alpha\beta}\left(\left(\frac{f_{0}\tilde{f}_{0}p^{\beta}p^{\mu}}{u.p}\right)\left[\left(D_{\mu}u_{\beta}\right)\partial_{\kappa}\left(\beta\tau_{c}\right)+\beta\tau_{c}\partial_{\kappa}\left(D_{\mu}u_{\beta}\right)\right]\right),$

$\displaystyle II$ $\displaystyle=$
$\displaystyle-\mathcal{I}^{\lambda\chi}_{\alpha\beta}qF^{\mu}_{\quad\nu}\frac{p^{\alpha}p^{\beta}p^{\nu}}{(u.p)}\tau_{c}f_{0}\tilde{f}_{0}\partial_{\mu}\beta$
$\displaystyle-\mathcal{I}^{\lambda\chi}_{\alpha\beta}qF^{\mu}_{\quad\nu}\left(\left[\tau_{c}p^{\kappa}(\partial_{\kappa}\beta)f_{0}\tilde{f}_{0}\beta
u_{\mu}\left(-\tilde{f}_{0}+rf_{0}\right)\right]+\beta\tau_{c}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}p^{\kappa}\left(p^{\mu}D_{\mu}u_{\alpha}+p^{\kappa}D_{\alpha}u_{\kappa}\right)\right]\right)$
$\displaystyle-\mathcal{I}^{\lambda\chi}_{\alpha\beta}qF^{\mu}_{\quad\nu}\left(\beta\tau_{c}\left(\sigma_{\kappa\beta}+\omega_{\kappa\beta}+\frac{\Delta_{\kappa\beta}\theta}{3}+u_{\kappa}\dot{u}_{\beta}\right)p^{\beta}p^{\kappa}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right)$
$\displaystyle+\mathcal{I}^{\lambda\chi}_{\alpha\beta}qF^{\mu}_{\quad\nu}\left(\tau_{c}\partial_{\kappa}\alpha\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\delta^{\kappa}_{\mu}+p^{\kappa}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right]+\tau_{c}\beta
qE_{\mu}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\right]\right)$
$\displaystyle+\mathcal{I}^{\lambda\chi}_{\alpha\beta}qF^{\mu}_{\quad\nu}\left(\tau_{c}\beta
qE^{\nu}p_{\nu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right),$

$\displaystyle III$ $\displaystyle=$
$\displaystyle\mathcal{I}^{\lambda\chi}_{\alpha\beta}\Gamma^{\mu}_{\nu\rho}\frac{p^{\alpha}p^{\beta}p^{\nu}p^{\rho}}{(u.p)}\left(\tau_{c}\left[f_{0}\tilde{f}_{0}\partial_{\mu}\beta\right]\right)$
$\displaystyle+\mathcal{I}^{\lambda\chi}_{\alpha\beta}\Gamma^{\mu}_{\nu\rho}\left(\left[\tau_{c}p^{\kappa}(\partial_{\kappa}\beta)f_{0}\tilde{f}_{0}\beta
u_{\mu}\left(-\tilde{f}_{0}+rf_{0}\right)\right]+\beta\tau_{c}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}p^{\kappa}\left(p^{\mu}D_{\mu}u_{\alpha}+p^{\kappa}D_{\alpha}u_{\kappa}\right)\right]\right)$
$\displaystyle+\mathcal{I}^{\lambda\chi}_{\alpha\beta}\Gamma^{\mu}_{\nu\rho}\left(\beta\tau_{c}\left(\sigma_{\kappa\beta}+\omega_{\kappa\beta}+\frac{\Delta_{\kappa\beta}\theta}{3}+u_{\kappa}\dot{u}_{\beta}\right)p^{\beta}p^{\kappa}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right)$
$\displaystyle+\mathcal{I}^{\lambda\chi}_{\alpha\beta}\Gamma^{\mu}_{\nu\rho}\left(-\tau_{c}\partial_{\kappa}\alpha\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\delta^{\kappa}_{\mu}+p^{\kappa}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right]-\tau_{c}\beta
qE_{\mu}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\right]\right)$
$\displaystyle+\mathcal{I}^{\lambda\chi}_{\alpha\beta}\Gamma^{\mu}_{\nu\rho}\left(-\tau_{c}\beta
qE^{\nu}p_{\nu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right).$

Combining the above expressions for I,II,III, and Eq.(21) we get Eq.(24).

###  B.2 Bulk viscosity

We follow exactly same procedure to calculate $\Pi_{(2)}$ i.e.,

$\displaystyle\Pi_{(2)}$ $\displaystyle=$ $\displaystyle I+II+III,$

where

$\displaystyle I$ $\displaystyle=$
$\displaystyle\frac{\Delta_{\alpha\beta}}{3}\tau_{c}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2})\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\partial_{\kappa}\delta
f^{(1)},$ $\displaystyle II$ $\displaystyle=$
$\displaystyle\frac{\Delta_{\alpha\beta}}{3}\tau_{c}qF^{\mu}_{\quad\nu}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2})\frac{p^{\alpha}p^{\beta}p^{\nu}}{(u.p)}\frac{\partial\delta
f^{(1)}}{\partial p^{\mu}},$ $\displaystyle III$ $\displaystyle=$
$\displaystyle-\frac{\Delta_{\alpha\beta}}{3}\tau_{c}\Gamma^{\mu}_{\nu\rho}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2})\frac{p^{\alpha}p^{\beta}p^{\nu}p^{\rho}}{(u.p)}\frac{\partial\delta
f^{(1)}}{\partial p^{\mu}}.$

For brevity let us define

$\displaystyle\mathcal{I}_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{\Delta_{\alpha\beta}}{3}\tau_{c}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2}).$

Now putting $\delta f^{(1)}$ and $\mathcal{I}_{\alpha\beta}$ in I,II, and III
we get :

$\displaystyle I$ $\displaystyle=$
$\displaystyle\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\left(\tau_{c}(\partial_{\mu}\beta)p^{\mu}\left[-f_{0}\tilde{f}_{0}\left(\tilde{f}_{0}-rf_{0}\right)\left((u.p)\partial_{\kappa}\beta+\beta
p^{\beta}\left(D_{\kappa}u_{\beta}+\Gamma^{\nu}_{\kappa\beta}u_{\nu}\right)-\partial_{\kappa}\alpha\right)\right]\right)$
$\displaystyle-\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\tau_{c}(\partial_{\mu}\alpha)p^{\mu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\kappa}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\kappa}u_{\beta}+\Gamma^{\nu}_{\kappa\beta}u_{\nu}\right)\right]$
$\displaystyle-\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\left(\tau_{c}\beta
qE^{\nu}p_{\nu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\kappa}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\kappa}u_{\beta}+\Gamma^{\nu}_{\kappa\beta}u_{\nu}\right)\right]-f_{0}\tilde{f}_{0}p^{\mu}\partial_{\kappa}\left(\tau_{c}\partial_{\mu}\beta\right)\right)$
$\displaystyle+\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\left(\beta\tau_{c}(D_{\mu}u_{\beta})p^{\beta}p^{\mu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\kappa}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\kappa}u_{\beta}+\Gamma^{\nu}_{\kappa\beta}u_{\nu}\right)\right]\right)$
$\displaystyle-\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\left(\left(\frac{f_{0}\tilde{f}_{0}p^{\mu}}{u.p}\right)\partial_{\kappa}\left(\tau_{c}\partial_{\mu}\alpha\right)+q\left(f_{0}\tilde{f}_{0}\frac{p_{\nu}}{u.p}\right)\left[(\tau_{c}\beta)\partial_{\kappa}E^{\nu}+E^{\nu}\partial_{\kappa}(\tau_{c}\beta)\right]\right)$
$\displaystyle+\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\kappa}}{(u.p)}\left(\left(\frac{f_{0}\tilde{f}_{0}p^{\beta}p^{\mu}}{u.p}\right)\left[\left(D_{\mu}u_{\beta}\right)\partial_{\kappa}\left(\beta\tau_{c}\right)+\beta\tau_{c}\partial_{\kappa}\left(D_{\mu}u_{\beta}\right)\right]\right),$

$\displaystyle II$ $\displaystyle=$ $\displaystyle
qF^{\mu}_{\quad\nu}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}}{(u.p)}\tau_{c}f_{0}\tilde{f}_{0}\partial_{\mu}\beta$
$\displaystyle+qF^{\mu}_{\quad\nu}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}}{(u.p)}\left(\left[\tau_{c}p^{\kappa}(\partial_{\kappa}\beta)f_{0}\tilde{f}_{0}\beta
u_{\mu}\left(-\tilde{f}_{0}+rf_{0}\right)\right]+\beta\tau_{c}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\left(p^{\mu}D_{\mu}u_{\alpha}+p^{\kappa}D_{\alpha}u_{\kappa}\right)\right]\right)$
$\displaystyle+qF^{\mu}_{\quad\nu}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}}{(u.p)}\left(\beta\tau_{c}\left(\sigma_{\kappa\beta}+\omega_{\kappa\beta}+\frac{\Delta_{\kappa\beta}\theta}{3}+u_{\kappa}\dot{u}_{\beta}\right)p^{\beta}p^{\kappa}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right)$
$\displaystyle-
qF^{\mu}_{\quad\nu}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}}{(u.p)}\left(\tau_{c}\partial_{\kappa}\alpha\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\delta^{\kappa}_{\mu}+p^{\kappa}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right]+\tau_{c}\beta
qE_{\mu}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\right]\right)$ $\displaystyle-
qF^{\mu}_{\quad\nu}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}}{(u.p)}\left(\tau_{c}\beta
qE^{\nu}p_{\nu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right),$

$\displaystyle III$ $\displaystyle=$
$\displaystyle-\Gamma^{\mu}_{\nu\rho}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}p^{\rho}}{(u.p)}\tau_{c}f_{0}\tilde{f}_{0}\partial_{\mu}\beta$
$\displaystyle-\Gamma^{\mu}_{\nu\rho}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}p^{\rho}}{(u.p)}\left(\left[\tau_{c}p^{\kappa}(\partial_{\kappa}\beta)f_{0}\tilde{f}_{0}\beta
u_{\mu}\left(-\tilde{f}_{0}+rf_{0}\right)\right]+\beta\tau_{c}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}p^{\kappa}\left(p^{\mu}D_{\mu}u_{\alpha}+p^{\kappa}D_{\alpha}u_{\kappa}\right)\right]\right)$
$\displaystyle-\Gamma^{\mu}_{\nu\rho}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}p^{\rho}}{(u.p)}\left(\beta\tau_{c}\left(\sigma_{\kappa\beta}+\omega_{\kappa\beta}+\frac{\Delta_{\kappa\beta}\theta}{3}+u_{\kappa}\dot{u}_{\beta}\right)p^{\beta}p^{\kappa}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right)$
$\displaystyle+\Gamma^{\mu}_{\nu\rho}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}p^{\rho}}{(u.p)}\left(\tau_{c}\partial_{\kappa}\alpha\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\delta^{\kappa}_{\mu}+p^{\kappa}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right]+\tau_{c}\beta
qE_{\mu}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\right]\right)$
$\displaystyle+\Gamma^{\mu}_{\nu\rho}\mathcal{I}_{\alpha\beta}\frac{p^{\alpha}p^{\beta}p^{\nu}p^{\rho}}{(u.p)}\left(\tau_{c}\beta
qE^{\nu}p_{\nu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\mu}+f_{0}\tilde{f}_{0}\beta
u_{\mu}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right).$

Using these above expressions for I, II, III, and Eq.(19) we have Eq.(22).

###  B.3  Particle Diffusion

The calculation of the diffusion four-current is exactly same:

$\displaystyle V^{\mu}_{(2)}$ $\displaystyle=$ $\displaystyle I+II+III,$

where

$\displaystyle I$ $\displaystyle=$
$\displaystyle-\tau_{c}\Delta^{\mu}_{\beta}\int
dp\frac{p^{\beta}p^{\alpha}}{u.p}\partial_{\alpha}\delta f^{(1)},$
$\displaystyle II$ $\displaystyle=$
$\displaystyle-\tau_{c}qF^{\alpha}_{\quad\rho}\Delta^{\mu}_{\beta}\int
dp\frac{p^{\beta}p^{\rho}}{u.p}\frac{\partial\delta f^{(1)}}{\partial
p^{\alpha}},$ $\displaystyle III$ $\displaystyle=$
$\displaystyle\tau_{c}\Delta^{\mu}_{\beta}\Gamma^{\alpha}_{\nu\rho}\int
dp\frac{p^{\beta}p^{\nu}p^{\rho}}{u.p}\frac{\partial\delta f^{(1)}}{\partial
p^{\alpha}}.$

For convenience we define

$\displaystyle\mathcal{I}^{\mu}_{\beta}$ $\displaystyle=$
$\displaystyle\tau_{c}\Delta^{\mu}_{\beta}\int\frac{d^{4}p}{(2\pi)^{3}}\sqrt{-g}2\Theta(p^{0})\delta(p^{2}-m^{2}).$
$\displaystyle I$ $\displaystyle=$
$\displaystyle-\mathcal{I}^{\mu}_{\beta}\frac{p^{\alpha}p^{\beta}}{u.p}\tau_{c}(\partial_{\mu}\beta)p^{\mu}\left[-f_{0}\tilde{f}_{0}\left(\tilde{f}_{0}-rf_{0}\right)\left((u.p)\partial_{\alpha}\beta+\beta
p^{\beta}\left(D_{\alpha}u_{\beta}+\Gamma^{\nu}_{\alpha\beta}u_{\nu}\right)-\partial_{\alpha}\alpha\right)\right]$
$\displaystyle+\mathcal{I}^{\mu}_{\beta}\frac{p^{\alpha}p^{\beta}}{u.p}\tau_{c}(\partial_{\mu}\alpha)p^{\mu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\alpha}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\alpha}u_{\beta}+\Gamma^{\nu}_{\alpha\beta}u_{\nu}\right)\right]$
$\displaystyle+\mathcal{I}^{\mu}_{\beta}\frac{p^{\alpha}p^{\beta}}{u.p}\left(\tau_{c}\beta
qE^{\nu}p_{\nu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\alpha}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\alpha}u_{\beta}+\Gamma^{\nu}_{\alpha\beta}u_{\nu}\right)\right]-f_{0}\tilde{f}_{0}p^{\mu}\partial_{\alpha}\left(\tau_{c}\partial_{\mu}\beta\right)\right)$
$\displaystyle-\mathcal{I}^{\mu}_{\beta}\frac{p^{\alpha}p^{\beta}}{u.p}\beta\tau_{c}(D_{\mu}u_{\beta})p^{\beta}p^{\mu}\left[\frac{\left(\tilde{f}_{0}-rf_{0}\right)\partial_{\alpha}f_{0}}{u.p}-\frac{\left(f_{0}\tilde{f}_{0}\right)}{(u.p)^{2}}p^{\beta}\left(D_{\alpha}u_{\beta}+\Gamma^{\nu}_{\alpha\beta}u_{\nu}\right)\right]$
$\displaystyle+\mathcal{I}^{\mu}_{\beta}\frac{p^{\alpha}p^{\beta}}{u.p}\left(\left(\frac{f_{0}\tilde{f}_{0}p^{\mu}}{u.p}\right)\partial_{\alpha}\left(\tau_{c}\partial_{\mu}\alpha\right)+q\left(f_{0}\tilde{f}_{0}\frac{p_{\nu}}{u.p}\right)\left[(\tau_{c}\beta)\partial_{\alpha}E^{\nu}+E^{\nu}\partial_{\alpha}(\tau_{c}\beta)\right]\right)$
$\displaystyle-\mathcal{I}^{\mu}_{\beta}\frac{p^{\alpha}p^{\beta}}{u.p}\left(\frac{f_{0}\tilde{f}_{0}p^{\beta}p^{\mu}}{u.p}\right)\left[\left(D_{\mu}u_{\beta}\right)\partial_{\alpha}\left(\beta\tau_{c}\right)+\beta\tau_{c}\partial_{\alpha}\left(D_{\mu}u_{\beta}\right)\right],$

$\displaystyle II$ $\displaystyle=$ $\displaystyle-
qF^{\alpha}_{\quad\rho}\mathcal{I}^{\mu}_{\beta}\frac{p^{\beta}p^{\rho}}{u.p}\left(\tau_{c}\left[f_{0}\tilde{f}_{0}\partial_{\alpha}\beta+p^{\mu}(\partial_{\mu}\beta)f_{0}\tilde{f}_{0}\beta
u_{\alpha}\left(-\tilde{f}_{0}+rf_{0}\right)\right]+\beta\tau_{c}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\left(p^{\mu}D_{\mu}u_{\alpha}+p^{\kappa}D_{\alpha}u_{\kappa}\right)\right]\right)$
$\displaystyle-
qF^{\alpha}_{\quad\rho}\mathcal{I}^{\mu}_{\beta}\frac{p^{\beta}p^{\rho}}{u.p}\beta\tau_{c}\left(\sigma_{\mu\beta}+\omega_{\mu\beta}+\frac{\Delta_{\mu\beta}\theta}{3}+u_{\mu}\dot{u}_{\beta}\right)p^{\beta}p^{\mu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]$
$\displaystyle-
qF^{\alpha}_{\quad\rho}\mathcal{I}^{\mu}_{\beta}\frac{p^{\beta}p^{\rho}}{u.p}\left(-\tau_{c}\partial_{\mu}\alpha\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\delta^{\mu}_{\alpha}+p^{\mu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right]-\tau_{c}\beta
qE_{\alpha}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\right]\right)$
$\displaystyle+qF^{\alpha}_{\quad\rho}\mathcal{I}^{\mu}_{\beta}\frac{p^{\beta}p^{\rho}}{u.p}\tau_{c}\beta
qE^{\nu}p_{\nu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right],$

$\displaystyle III$ $\displaystyle=$
$\displaystyle\Gamma^{\alpha}_{\nu\rho}\mathcal{I}^{\mu}_{\beta}\frac{p^{\beta}p^{\nu}p^{\rho}}{u.p}\left(\tau_{c}\left[f_{0}\tilde{f}_{0}\partial_{\alpha}\beta+p^{\mu}(\partial_{\mu}\beta)f_{0}\tilde{f}_{0}\beta
u_{\alpha}\left(-\tilde{f}_{0}+rf_{0}\right)\right]+\beta\tau_{c}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\left(p^{\mu}D_{\mu}u_{\alpha}+p^{\kappa}D_{\alpha}u_{\kappa}\right)\right]\right)$
$\displaystyle+\Gamma^{\alpha}_{\nu\rho}\mathcal{I}^{\mu}_{\beta}\frac{p^{\beta}p^{\nu}p^{\rho}}{u.p}\left(\beta\tau_{c}\left(\sigma_{\mu\beta}+\omega_{\mu\beta}+\frac{\Delta_{\mu\beta}\theta}{3}+u_{\mu}\dot{u}_{\beta}\right)p^{\beta}p^{\mu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right)$
$\displaystyle-\Gamma^{\alpha}_{\nu\rho}\mathcal{I}^{\mu}_{\beta}\frac{p^{\beta}p^{\nu}p^{\rho}}{u.p}\left(\tau_{c}\partial_{\mu}\alpha\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\delta^{\mu}_{\alpha}+p^{\mu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right]\right]+\tau_{c}\beta
qE_{\alpha}\left[\frac{f_{0}\tilde{f}_{0}}{u.p}\right]\right)$
$\displaystyle-\Gamma^{\alpha}_{\nu\rho}\mathcal{I}^{\mu}_{\beta}\frac{p^{\beta}p^{\nu}p^{\rho}}{u.p}\tau_{c}\beta
qE^{\nu}p_{\nu}\left[-\frac{f_{0}\tilde{f}_{0}}{(u.p)^{2}}u_{\alpha}+f_{0}\tilde{f}_{0}\beta
u_{\alpha}\frac{\left(-\tilde{f}_{0}+rf_{0}\right)}{u.p}\right].$

Adding I, II, III, and Eq.(20) we have the desired result Eq.(23).

##  Appendix C Causality and stability

In this appendix we discuss details of causality and stability given in
Sec.(V) in the main text. For later use, we require the following expression
for the $\delta\tilde{\Gamma}^{\mu}_{\nu\alpha}$:

$\displaystyle\delta\tilde{\Gamma}^{\mu}_{\nu\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[i\delta\tilde{g}^{\mu}_{\alpha}k_{\nu}+i\delta\tilde{g}^{\mu}_{\nu}k_{\alpha}-i\delta\tilde{g}_{\nu\alpha}k^{\mu}\right].$
(53)

We consider three different scenarios and discuss them one-by-one.

###  C.1 Non-viscous fluid with metric perturbation

The first case we consider deals with non-viscous fluid in presence of metric
perturbation. Various equations pertaining to this are :

$\displaystyle
i\omega\delta\tilde{\epsilon}+\left(\epsilon_{0}+P_{0}\right)ik_{x}\delta\tilde{u}^{x}+ik_{z}\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{z}$
$\displaystyle=$ $\displaystyle 0,$ (54) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{x}+ik_{x}c_{s}^{2}\delta\tilde{\epsilon}$
$\displaystyle=$ $\displaystyle 0,$ (55) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{y}$ $\displaystyle=$
$\displaystyle 0,$ (56) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{z}+ik_{z}c_{s}^{2}\delta\tilde{\epsilon}$
$\displaystyle=$ $\displaystyle 0,$ (57) $\displaystyle 8\pi
Gc_{s}^{2}\delta\tilde{\epsilon}-\left(\Lambda+\frac{1}{2}\left(\omega^{2}-k^{2}\right)-8\pi
GP_{0}\right)\delta\tilde{g}^{xx}$ $\displaystyle=$ $\displaystyle 0,$ (58)
$\displaystyle-\left(\Lambda+\frac{1}{2}\left(\omega^{2}-k^{2}\right)-8\pi
GP_{0}\right)\delta\tilde{g}^{xy}$ $\displaystyle=$ $\displaystyle 0.$ (59)

These above set of equations are used in writing Matrix $\mathcal{A}$ in
sec.(V.1)

###  C.2  Viscous fluid with metric perturbation

Let us now turn to the second case where we consider a fluid with non-zero
bulk viscosity in presence of metric perturbations. The various equations
pertaining to this is :

$\displaystyle
i\omega\left(\delta\tilde{\epsilon}-\delta\tilde{\Pi}\right)+i\left(\epsilon_{0}+P_{0}\right)k_{x}\delta\tilde{u}^{x}+ik_{z}\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{z}$
$\displaystyle=$ $\displaystyle 0,$ (60) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{x}+ik_{x}\left(c_{s}^{2}\delta\tilde{\epsilon}+\delta\tilde{\Pi}\right)$
$\displaystyle=$ $\displaystyle 0,$ (61) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{y}$ $\displaystyle=$
$\displaystyle 0,$ (62) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{z}+ik_{z}\left(c_{s}^{2}\delta\tilde{\epsilon}+\delta\tilde{\Pi}\right)$
$\displaystyle=$ $\displaystyle 0,$ (63)
$\displaystyle\delta\tilde{\Pi}\left(\frac{1}{\tau_{c}}+i\omega\right)+\beta_{\Pi}ik_{x}\delta\tilde{u}^{x}+\beta_{\Pi}ik_{z}\delta\tilde{u}^{z}$
$\displaystyle=$ $\displaystyle 0,$ (64) $\displaystyle 8\pi
G\left(c_{s}^{2}\delta\tilde{\epsilon}+\delta\tilde{\Pi}\right)-\left(\Lambda+\frac{1}{2}\left(\omega^{2}-k^{2}\right)-8\pi
GP_{0}\right)\delta\tilde{g}^{xx}$ $\displaystyle=$ $\displaystyle 0,$ (65)
$\displaystyle-\left(\Lambda+\frac{1}{2}\left(\omega^{2}-k^{2}\right)-8\pi
GP_{0}\right)\delta\tilde{g}^{xy}$ $\displaystyle=$ $\displaystyle 0.$ (66)

These above set of equations are used in writing Matrix $\mathcal{A}$ in
sec.(V.2).

###  C.3  Viscous fluid with metric perturbation and magnetic field

The last case we consider is applicable for a viscous fluid (only bulk
viscosity) within the ideal MHD framework with metric perturbation. Here in
addition to the usual fluid variables we also need to consider magnetic field
as an additional parameter. The linearised equations take the following form:

$\displaystyle
i\omega\left(\delta\tilde{\epsilon}-\delta\tilde{\Pi}\right)+\left(\epsilon_{0}+P_{0}\right)ik_{x}\delta\tilde{u}^{x}+ik_{z}\left(\epsilon_{0}+P_{0}\right)\delta\tilde{u}^{z}$
$\displaystyle=$ $\displaystyle 0,$ (67) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}+B^{2}_{0}\right)\delta\tilde{u}^{x}+ik_{x}\left[c_{s}^{2}\delta\tilde{\epsilon}+\delta\tilde{\Pi}-\frac{k_{x}}{\omega}\delta\tilde{u}^{x}B^{2}_{0}\right]$
$\displaystyle-
ik_{z}B^{2}_{0}\delta\tilde{b}^{x}+ik_{x}\left(\delta\tilde{g}^{xx}-\frac{k_{x}^{2}+k_{z}^{2}}{k_{z}^{2}}\delta\tilde{g}^{xx}\right)B_{0}^{2}-ik_{x}\delta\tilde{g}^{xx}B^{2}_{0}$
$\displaystyle=$ $\displaystyle 0,$ (68) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}+B^{2}_{0}\right)\delta\tilde{u}^{y}-ik_{z}B^{2}_{0}\delta\tilde{b}^{y}-i\frac{k_{x}}{k_{z}}\delta\tilde{g}^{xy}k_{z}B^{2}_{0}$
$\displaystyle=$ $\displaystyle 0,$ (69) $\displaystyle
i\omega\left(\epsilon_{0}+P_{0}+B^{2}_{0}\right)\delta\tilde{u}^{z}+ik_{z}\left(c_{s}^{2}\delta\tilde{\epsilon}+\delta\tilde{\Pi}\right)-i\omega
B^{2}_{0}\delta\tilde{u}^{z}$
$\displaystyle+ik_{z}\left(\frac{3k_{x}}{\omega}\delta\tilde{u}^{x}-\delta\tilde{g}^{xx}+\frac{k_{x}^{2}+k_{z}^{2}}{k_{z}^{2}}\delta\tilde{g}^{xx}\right)B^{2}_{0}-ik_{x}B^{2}_{0}\delta\tilde{b}^{x}$
$\displaystyle=$ $\displaystyle 0,$ (70)
$\displaystyle\delta\tilde{\Pi}\left(\frac{1}{\tau_{c}}+i\omega\right)+i\beta_{\Pi}k_{x}\delta\tilde{u}^{x}+i\beta_{\Pi}k_{z}\delta\tilde{u}^{z}$
$\displaystyle=$ $\displaystyle 0,$ (71) $\displaystyle
ik_{z}B_{0}\delta\tilde{u}^{x}-i\omega B_{0}\delta\tilde{b}^{x}$
$\displaystyle=$ $\displaystyle 0,$ (72) $\displaystyle
ik_{z}B_{0}\delta\tilde{u}^{y}-i\omega B_{0}\delta\tilde{b}^{y}$
$\displaystyle=$ $\displaystyle 0,$ (73)
$\displaystyle-\Lambda\delta\tilde{g}^{xx}-\frac{1}{2}\delta\tilde{g}^{xx}\left(\omega^{2}-k^{2}\right)+8\pi
G\left(c_{s}^{2}\delta\tilde{\epsilon}+\delta\tilde{\Pi}-\frac{k_{x}}{\omega}\delta\tilde{u}^{x}B^{2}_{0}-\frac{\delta\tilde{g}^{zz}}{2}B^{2}_{0}\right)+8\pi
G\left(P_{0}+\frac{B^{2}_{0}}{2}\right)\delta\tilde{g}^{xx}$ $\displaystyle=$
$\displaystyle 0,$ (74)
$\displaystyle-\Lambda\delta\tilde{g}^{xy}-\frac{1}{2}\delta\tilde{g}^{xy}\left(\omega^{2}-k^{2}\right)+8\pi
G\left(P_{0}+\frac{B^{2}_{0}}{2}\right)\delta\tilde{g}^{xy}$ $\displaystyle=$
$\displaystyle 0.$ (75)

These above set of equations are used in writing Matrix $\mathcal{A}$ in
sec.(V.3).
The full expression for the $v_{M}$ which was not included in the main text
for $\theta=\pi/4$ is as follows :

$\displaystyle v_{M}$ $\displaystyle=$
$\displaystyle\pm\frac{1}{4\left(B^{2}_{0}+\epsilon_{0}+P_{0}\right)\left(8\pi
GP-\Lambda\right)}\left[-B^{2}_{0}(2+c_{s}^{2})\Lambda-2(1+c_{s}^{2})\Lambda
P_{0}+8\pi
G\left(B^{2}_{0}(2-c_{s}^{2})+B^{2}_{0}c_{s}^{2}(3\epsilon_{0}+2P_{0})+2(P_{0}^{2}+c_{s}^{2}P_{0}^{2})\right)\right]$
(76)
$\displaystyle\pm\frac{1}{4\left(B^{2}_{0}+\epsilon_{0}+P_{0}\right)}\left[\pm\sqrt{-16B^{2}_{0}c_{s}^{2}(B^{2}_{0}+\epsilon_{0}+P_{0})(\Lambda-8\pi
GP_{0})\left(-10B^{2}_{0}\pi G+\Lambda-8\pi
GP_{0}\right)+\mathcal{T}}\right],$

where

$\mathcal{T}=\left[8B^{4}_{0}(c_{s}^{2}-2)\pi
G+B^{2}_{0}(2+c_{s}^{2})\Lambda+2P_{0}\left(-4B^{2}_{0}(3+2c_{s}^{2})\pi
G+\Lambda+\Lambda c_{s}^{2}-8c_{s}^{2}(\epsilon_{0}+P_{0})\pi
G\right)\right]^{2}.$

## References

  * [1] Relativistic viscous hydrodynamics, conformal invariance, and holography, Baier, Rudolf and Romatschke, Paul and Son, Dam Thanh and Starinets, Andrei O. and Stephanov, Mikhail A. 10.1088/1126-6708/2008/04/100
  * [2] Magnetohydrodynamic Simulations of Binary Neutron Star Mergers in General Relativity: Effects of Magnetic Field Orientation on Jet Launching,Ruiz, Milton and Tsokaros, Antonios and Shapiro, Stuart L. Phys.Rev.D 101, 064042 (2020)
  * [3] GRMHD simulations of accreting neutron stars I: nonrotating dipoles,Çıkıntoğlu, Sercan and Ekşi, K. Yavuz and Rezzolla, Luciano. 2204.12275
  * [4] General relativistic resistive magnetohydrodynamics with robust primitive variable recovery for accretion disk simulations,Ripperda, Bart and Bacchini, Fabio and Porth, Oliver and Most, Elias R. and Olivares, Hector and Nathanail, Antonios and Rezzolla, Luciano and Teunissen, Jannis and Keppens, Rony. 10.3847/1538-4365/ab3922
  * [5] General-relativistic hydrodynamics of non-perfect fluids: 3+1 conservative formulation and application to viscous black hole accretion,Chabanov, Michail and Rezzolla, Luciano and Rischke, Dirk H. 10.1093/mnras/stab1384
  * [6] Relativistic magnetohydrodynamics in dynamical spacetimes: Numerical methods and tests,Matthew D. Duez, Yuk Tung Liu, Stuart L. Shapiro, and Branson C. Stephens. Phys. Rev. D 72, 024028 – Published 28 July 2005
  * [7] Classical and Quantum Gravity WhiskyMHD: a new numerical code for general relativistic magnetohydrodynamics.Bruno Giacomazzo and Luciano Rezzolla.  Class. Quantum Grav. 24 S235
  * [8] Long-term general relativistic magnetohydrodynamics simulations of magnetic field in isolated neutron stars,Ankan Sur, William Cook, David Radice, Brynmor Haskell, Sebastiano Bernuzzi.  https://doi.org/10.1093/mnras/stac353
  * [9] General relativistic magnetohydrodynamics simulations of prompt-collapse neutron star mergers: The absence of jets,Milton Ruiz and Stuart L. Shapiro. Phys. Rev. D 96, 084063 – Published 30 October 2017
  * [10] Dissipative magnetohydrodynamics for nonresistive relativistic plasmas: An implicit second-order flux-conservative formulation with stiff relaxation, Most, Elias R. and Noronha, Jorge. 10.1103/PhysRevD.104.103028
  * [11] Modeling general-relativistic plasmas with collisionless moments and dissipative two-fluid magnetohydrodynamics,Most, Elias R. and Noronha, Jorge and Philippov, Alexander A. 2111.05752
  * [12] Event-by-event distribution of magnetic field energy over initial fluid energy density in $\sqrt{s_{\rm NN}}$= 200 GeV Au-Au collisions,Roy, Victor and Pu, Shi. 10.1103/PhysRevC.92.064902
  * [13] Effect of intense magnetic fields on reduced-MHD evolution in $\sqrt{s_{\rm NN}}$ = 200 GeV Au+Au collisions,Roy, Victor and Pu, Shi and Rezzolla, Luciano and Rischke, Dirk H. 10.1103/PhysRevC.96.054909
  * [14] Estimate of the magnetic field strength in heavy-ion collisions,Skokov, V. and Illarionov, A. Yu. and Toneev, V. 10.1142/S0217751X09047570
  * [15] Particle production in strong electromagnetic fields in relativistic heavy-ion collisions,Tuchin, Kirill. 10.1155/2013/490495
  * [16] A short course on Relativistic Heavy Ion Collisions,Chaudhuri, A. K. 10.1088/978-0-750-31060-4
  * [17] Hydrodynamic Modeling of Heavy-Ion Collisions,Gale, Charles and Jeon, Sangyong and Schenke, Bjoern. 10.1142/S0217751X13400113
  * [18] Relativistic Fluid Dynamics In and Out of Equilibrium, Romatschke, Paul and Romatschke, Ulrike. 10.1017/9781108651998
  * [19] Collective flow and viscosity in relativistic heavy-ion collisions, Heinz, Ulrich and Snellings, Raimond. 10.1146/annurev-nucl-102212-170540
  * [20] Stability and causality in dissipative relativistic fluids,Hiscock, W. A. and Lindblom, L. 10.1016/0003-4916(83)90288-9
  * [21] J.M. Stewart Proc.Roy.Soc.Lond.A 357 (1977) 59
  * [22] Israel, Werner, and John M. Stewart. ”Transient relativistic thermodynamics and kinetic theory.” Annals of Physics 118.2 (1979): 341-372. 10.1016/0003-4916(79)90130-1
  * [23] Does stability of relativistic dissipative fluid dynamics imply causality?,Pu, Shi and Koide, Tomoi and Rischke, Dirk H. 10.1103/PhysRevD.81.114039
  * [24] Stability and Causality in relativistic dissipative hydrodynamics,Denicol, G. S. and Kodama, T. and Koide, T. and Mota, Ph. 10.1088/0954-3899/35/11/115102
  * [25] Linearized dispersion relations in viscous relativistic hydrodynamics.Perna, Guillermo and Calzetta, Esteban. 10.1103/PhysRevD.104.096005
  * [26] Causal Relativistic Hydrodynamics of Conformal Fermi-Dirac Gases, Aguilar, Milton and Calzetta, Esteban. 10.1103/PhysRevD.95.076022
  * [27] Causality and stability in relativistic viscous non-resistive magneto-fluid dynamics,Biswas, Rajesh and Dash, Ashutosh and Haque, Najmul and Pu, Shi and Roy, Victor. 10.1007/JHEP10(2020)171
  * [28] Nonresistive dissipative magnetohydrodynamics from the Boltzmann equation in the 14-moment approximation, Denicol, Gabriel S. and Huang, Xu-Guang and Molnár, Etele and Monteiro, Gustavo M. and Niemi, Harri and Noronha, Jorge and Rischke, Dirk H. and Wang, Qun”. 10.1103/PhysRevD.98.076009
  * [29] Relativistic non-resistive viscous magnetohydrodynamics from the kinetic theory: a relaxation time approach, Panda, Ankit Kumar and Dash, Ashutosh and Biswas, Rajesh and Roy, Victor. 10.1007/JHEP03(2021)216
  * [30] Resistive dissipative magnetohydrodynamics from the Boltzmann-Vlasov equation, Denicol, Gabriel S. and Molnár, Etele and Niemi, Harri and Rischke, Dirk H. 10.1103/PhysRevD.99.056017
  * [31] Relativistic resistive dissipative magnetohydrodynamics from the relaxation time approximation, Panda, Ankit Kumar and Dash, Ashutosh and Biswas, Rajesh and Roy, Victor. Phys. Rev. D 104, 054004
  * [32] Novel Relaxation Time Approximation to the Relativistic Boltzmann Equation,Rocha, Gabriel S. and Denicol, Gabriel S. and Noronha, Jorge. 10.1103/PhysRevLett.127.042301
  * [33] Relativistic hydrodynamics with momentum dependent relaxation time,Mitra, Sukanya. 10.1103/PhysRevC.103.014905
  * [34] The Relativistic Boltzmann Equation: Theory and Applications ,Carlo Cercignani ,Gilberto Medeiros Kremer. 2002, Volume 22 ISBN : 978-3-0348-9463-0
  * [35] Relativistic second-order dissipative fluid dynamics at finite chemical potential,Jaiswal, Amaresh and Friman, Bengt and Redlich, Krzysztof. 10.1051/epjconf/201612003008
  * [36] One particle distribution function and shear viscosity in magnetic field: a relaxation time approach, Mohanty, Payal and Dash, Ashutosh and Roy, Victor. 10.1140/epja/i2019-12705-7
  * [37] Relativistic second order dissipative hydrodynamics from an effective fugacity quasiparticle model,Mitra, Sukanya. 10.1103/PhysRevD.100.014012
  * [38] Second-order Hydrodynamics in Next-to-Leading-Order QCD , Ghiglieri, Jacopo and Moore, Guy D. and Teaney, Derek. 10.1103/PhysRevLett.121.052302
  * [39] Transport coefficients in high temperature gauge theories. 1. Leading log results,Arnold, Peter Brockway and Moore, Guy D. and Yaffe, Laurence G. 10.1088/1126-6708/2000/11/001
  * [40] Transport coefficients in high temperature gauge theories. 2. Beyond leading log, 10.1088/1126-6708/2003/05/051
  * [41] Gravitomagnetic Instabilities of Relativistic Magnetohydrodynamics, Noh, Hyerim and Hwang, Jai-chan. 10.3847/1538-4357/abc422