2.1. Vlasov–Maxwell equations

The Vlasov–Maxwell equations on phase space $(x,k)$ are given by (Sarbach & Zannias Reference Sarbach and Zannias2013)

(2.3a) \begin{align} &k^{\mu } \frac {\partial f_{s}}{\partial x^{\mu }} + \left[{-} \varGamma ^{\alpha }_{\mu \nu } k^{\mu } k^{\nu } + \frac {q_{s}}{m_{s}} F^{\alpha }{}_{\mu } k^{\mu } \right] \frac {\partial f_{s}}{\partial k^{\alpha }} = \mathcal{C}_s\left [f_s\right ] ,\\[-9pt]\nonumber \end{align}

(2.3b) \begin{align} &\boldsymbol{\nabla} _{\mu } (\star F^{\mu \nu }) = 0, \quad \star F^{ \mu \nu } = \frac {1}{2} \epsilon ^{\mu \nu \kappa \lambda } F_{\kappa \lambda },\\[-10pt]\nonumber \end{align}

(2.3c) \begin{align} & \boldsymbol{\nabla} ^{\mu } F_{\mu \nu } = -4 \pi \sum _{s=1}^{S} q_{s} N_{\nu }^{s} , \end{align}

where $\varGamma ^{\alpha }_{\mu \nu }$ are the Christoffel symbols, $\mathcal{C}_s$ is the collision operator and $N^{\mu }_{s}$ is the number density current,

(2.4) \begin{align} N^{\mu }_{s} \equiv \left \langle k^{\mu }\right \rangle _{s} = \int d^4 k \sqrt {-g} f_{s} k^{\mu } = \int \frac {d^3k}{\left |k_0\right |} \sqrt {-g} f_{s} k^{\mu } . \end{align}

The first and second moments of the Vlasov equations yield the conservation laws for the number density and stress–tensor of species $s$ ,

(2.5a) \begin{align} &\boldsymbol{\nabla} _{\alpha } \left ( N_{s}^{\alpha } \right ) = 0 ,\\[-10pt]\nonumber \end{align}

(2.5b) \begin{align} &\boldsymbol{\nabla} _{\mu } T^{\mu \nu }_{s} = q_{s} F^{\nu}_{\mu } N_{s}^{\mu } , \end{align}

where $T^{\mu \nu }_{s} \equiv m_{s} \left \langle k^{\mu } k^{\nu } \right \rangle _{s}$ . Given a timelike vector $U^{\mu }$ , we decompose $N^{\mu }_{s}$ and $T^{\mu \nu }_{s}$ as

(2.6a) \begin{align} &N_{s}^{\mu } = n_{s} U^{\mu } + V^{\mu }_{s} ,\\[-10pt]\nonumber \end{align}

(2.6b) \begin{align} &T^{\mu \nu }_{s} = \mathcal{E}_{s} U^{\mu } U^{\nu } + \mathcal{P}_{s} \Delta ^{\mu \nu } + 2 Q^{(\mu }_{s} U^{\nu )} + \pi ^{\mu \nu }_{s} . \end{align}

In terms of the moments of the distribution function, the decomposition is given by

(2.7a) \begin{align} &n_{s} \equiv -N^{\mu }_{s} U_{\mu } = -\left \langle k_{\mu } U^{\mu }\right \rangle _{s} ,\\[-8pt]\nonumber \end{align}

(2.7b) \begin{align} &V^{\mu }_{s} \equiv \Delta ^{\mu }_{\nu } N^{\nu }_{s} = \Delta ^{\mu \nu }\left \langle k_{\nu } \right \rangle _{s} ,\\[-8pt]\nonumber \end{align}

(2.7c) \begin{align} &\mathcal{E}_{s} \equiv U_{\mu } U_{\nu } T^{\mu \nu }_{s} = m_s \left \langle (k_{\mu } U^{\mu })^2 \right \rangle _{s} ,\\[-8pt]\nonumber \end{align}

(2.7d) \begin{align} &\mathcal{P}_{s} \equiv \frac {1}{3} \Delta _{\mu \nu } T^{\mu \nu }_{s} = \frac {1}{3} m_{s} \left \langle \Delta _{\mu \nu } k^{\mu } k^{\nu }\right \rangle _{s},\\[-8pt]\nonumber \end{align}

(2.7e) \begin{align} &Q^{\mu }_{s} \equiv - \Delta ^{\mu \beta } U^{\alpha } T_{\alpha \beta }^{s} = - m_{s} \left \langle \Delta ^{\mu \beta } k_{\beta } U_{\alpha } k^{\alpha } \right \rangle _{s} ,\\[-8pt]\nonumber \end{align}

(2.7f) \begin{align} & \pi ^{\mu \nu }_{s} \equiv m_{s} \Delta ^{\mu \nu }_{\alpha \beta }\left \langle k^{\mu } k^{\nu }\right \rangle _{s} , \end{align}

where $n_{s}$ is the number density, $V_{s}^{\mu }$ is the particle diffusion current, $\mathcal{E}_{s}$ is the energy density, $\mathcal{P}_{s}$ is the pressure, $\mathcal{Q}^{\mu }_{s}$ is the heat vector and $\pi _{s}^{\mu \nu }$ is the shear tensor. The projection tensors $\Delta _{\mu \nu }$ and $\Delta ^{\alpha \beta }_{\mu \nu }$ are defined as

(2.7g) \begin{align} &\Delta _{\mu \nu } \equiv g_{\mu \nu } + U_{\mu } U_{\nu } ,\\[-10pt]\nonumber \end{align}

(2.7h) \begin{align} &\Delta ^{\alpha \beta }_{\mu \nu } \equiv \Delta ^{(\alpha }_{\mu }\Delta ^{\beta )}_{\nu } - \frac {1}{3}\Delta ^{ \alpha \beta } \Delta _{\mu \nu } . \end{align}

The fluid velocity vector $U^{\mu }$ , which was used in (2.7), has been arbitrary so far. In this paper, we adopt the Eckhart frame to define $U^{\mu }$ ,

(2.8) \begin{align} \sum _{s} m_{s} V^{\mu }_{s} = 0 . \end{align}

For a review of different hydrodynamic frame choices, please see Cercignani & Kremer (Reference Cercignani and Kremer2002), Bemfica et al. (Reference Bemfica, Disconzi and Noronha2022) and Kovtun (Reference Kovtun2019).

With the choice of the fluid frame specified, we can define the mass-averaged current and the stress-energy tensor of the system,

(2.9a) \begin{align} &m N^{\mu } = \sum _{s} m_s N_{s}^{\mu } = U^{\mu } \sum _{s} m_s n_s \equiv m n U^{\mu } ,\\[-10pt]\nonumber \end{align}

(2.9b) \begin{align} &T^{\mu \nu } = \mathcal{E} U^{\mu } U^{\nu } + \mathcal{P} \Delta ^{\mu \nu } + 2 Q^{(\mu } U^{\nu )} + \pi ^{\mu \nu } ,\\[-10pt]\nonumber \end{align}

(2.9c) \begin{align} &\mathcal{E} = \sum _{s} \mathcal{E}_{s} = \sum _{s} m_s \left \langle (k_{\mu } U^{\mu })^2 \right \rangle _{s} ,\\[-10pt]\nonumber \end{align}

(2.9d) \begin{align} &\mathcal{P} = \sum _{s} \mathcal{P}_{s} = \frac {1}{3} \sum _{s} m_{s} \left \langle \Delta _{\mu \nu } k^{\mu } k^{\nu }\right \rangle _{s},\\[-10pt]\nonumber \end{align}

(2.9e) \begin{align} &Q^{\mu } = \sum _{s} \mathcal{Q}^{\mu }_{s} = - \sum _{s} m_{s} \left \langle \Delta ^{\mu \beta } k_{\beta } U_{\alpha } k^{\alpha } \right \rangle _{s} ,\\[-10pt]\nonumber \end{align}

(2.9f) \begin{align} & \pi ^{\mu \nu } = \sum _{s} \pi ^{\mu \nu }_{s} = \sum _{s} m_{s} \Delta ^{\mu \nu }_{\alpha \beta }\left \langle k^{\mu } k^{\nu }\right \rangle _{s} , \end{align}

where $m \equiv \sum _{s} m_{s}$ is a reference mass scale, equal to the sum of the particle masses of the species included in the plasma. The evolution of the mean density current and stress energy tensor can be obtained by summing (2.5a ) and (2.5b ) over all the particles,

(2.10a) \begin{align} &\boldsymbol{\nabla} _{\alpha }\left (N^{\alpha } \right ) = \boldsymbol{\nabla} _{\alpha } \left ( n U^{\alpha } \right ) = 0,\\[-10pt]\nonumber \end{align}

(2.10b) \begin{align} &\boldsymbol{\nabla} _{\alpha } \left ( T^{\mu \nu } + T^{\mu \nu }_{\mathrm{EM}} \right ) = 0, \end{align}

where the electromagnetic stress energy tensor is

(2.11) \begin{align} T^{\mu \nu }_{\mathrm{EM}} \equiv \frac {1}{4\pi } \bigg (F^{\mu \lambda } F^{\nu }{}_{\lambda } - \frac {1}{4} g^{\mu \nu } F^{\alpha \beta } F_{\alpha \beta } \bigg ) . \end{align}

Observe that (2.9a ) and (2.10a ) are only valid in the Eckhart frame where particle diffusion over all the species in the plasma vanishes (2.8).

2.2. Kinetic MHD

Let us consider the characteristic scaling of each of the terms in (2.3a ) in a locally Minkowskian frame. Let $L_{c}$ denote the characteristic scale of spatiotemporal variation at the macroscopic level, $\tau _{c}$ the collision time scale, $R_{c}$ the characteristic radius of particle motion and $\omega _{c}$ the characteristic frequency. With these definitions, we can obtain scaling of the different terms in (2.3a ),

(2.12a) \begin{align} &\left |k^{\mu } \frac {\partial f_{s}}{\partial x^{\mu }} \right | \sim \frac {|k| |f_s|}{L_{c}} \sim \frac {R_c \omega _c |f_s|}{L_{c}} ,\\[-10pt]\nonumber \end{align}

(2.12b) \begin{align} &\left [\frac {q_{s}}{m_{s}} F^{\alpha }{}_{\mu } k^{\mu }\right ]\frac {\partial f_{s}}{\partial k^{\alpha }} \sim \frac {\left |a_{\mathrm{particle}}\right |}{|k|} |f_s| \sim \frac {R_c \omega _c^2 |f_s| }{|k|} \sim \omega _c |f_s| ,\\[-10pt]\nonumber \end{align}

(2.12c) \begin{align} &\mathcal{C}_s[f_s] \sim \tau _{c}^{-1} |f_s| . \end{align}

In this paper, we are interested in capturing the dynamics of diffuse accretion flow around an SMBH. Consider an accretion flow around an SMBH of mass $M_{\mathrm{BH}}$ , with a characteristic magnetic field $B_c$ . For simplicity, we assume the charged particles in the plasma are electrons and protons. In typical hot and diffusive black hole accretion flows, $L_c \sim M_{\mathrm{SMBH}}$ ; the characteristic radius for particle motion is the Larmor radius, the characteristic frequency is the cyclotron frequency and the collisional time scale is set by Coulomb collisions. For an M87-like system (Akiyama et al. Reference Akiyama2021), we obtain

(2.13a) \begin{align} &L_c \sim M_{\mathrm{SMBH}} \sim 1.5 \times 10^{11}\,\mathrm{m} \left (\frac {M}{10^8\, M_{\odot }} \right ) ,\\[-10pt]\nonumber \end{align}

(2.13b) \begin{align} &R_c \sim \frac {m_{\kern-1pt p} v_{\mathrm{th}}}{{e\kern-1.1pt B}_c} \sim 521\,\mathrm{m} \,\left (\frac {v_{\mathrm{th}}}{c}\right ) \left (\frac {30 \, \mathrm{G}}{B_c} \right ),\\[-10pt]\nonumber \end{align}

(2.13c) \begin{align} &\omega _c \sim \frac {{e\kern-1.1pt B}_c}{m_{\kern-1.1pt p}} \sim 2.8 \times 10^5\,\mathrm{Hz} \left (\frac {B_c}{30 \, \mathrm{G}}\right ) ,\\[-10pt]\nonumber \end{align}

(2.13d) \begin{align} &\tau _c^{-1} \sim \frac {v_{\mathrm{th}}}{n_{\mathrm{e,i}} \lambda _{D}^4} \sim 1.32 \times 10^{-6}\,\mathrm{Hz} \, \left (\frac {v_{\mathrm{th}}}{c}\right ) \left (\frac {n_{e,i}}{10^{7} \,\mathrm{cm}^{-3}}\right ) \left (\frac {T_{e,i}}{10^{10}\, \mathrm{K}}\right )^{-2} , \end{align}

where $\lambda _D$ is the Debye length. Hence, for diffusive accretion flows around an SMBH, we can safely assume that

(2.14) \begin{align} \omega _c \gg \frac {R_c \omega _c}{L_c} \gt \tau _c^{-1} , \end{align}

which implies that

(2.15) \begin{align} \left |\left [\frac {q_{s}}{m_{s}} F^{\alpha }{}_{\mu } k^{\mu }\right ]\frac {\partial f_{s}}{\partial k^{\alpha }} \right | \gg \left |k^{\mu } \frac {\partial f_{s}}{\partial x^{\mu }} \right | \gt \mathcal{C}_s[f_s] . \end{align}

We now implement the ordering obtained in (2.15) formally. Consider an expansion of the form (Kulsrud Reference Kulsrud1983)

(2.16) \begin{align} f_{s} = \sum _{n=0}^{\infty } f_{n,s} \epsilon ^{-n} , \end{align}

where $\epsilon \sim \mathcal{O}(|L_c/R_c|)$ . At leading order, the Vlasov equation simplifies to

(2.17) \begin{align} F^{\alpha }{}_{\mu } k^{\mu } \frac {\partial f_{s,0}}{\partial k^{\alpha }} = 0. \end{align}

We show later that this condition implies the distribution function is gyrotropic, as in the non-relativistic limit (Kulsrud Reference Kulsrud1983). Next, we analyse the stress-energy conservation equation (2.5b ) for particle $s$ . To leading order, the Lorentz-force term on the right-hand side of (2.5b ) dominates, so that each species satisfies

(2.18) \begin{align} F^{\nu}{}_{\mu }N^{\mu }_{s} = 0 . \end{align}

Multiplying this equation by $m_s$ , summing over species and using (2.9a ), we obtain

(2.19) \begin{align} F^{\nu}{}_{\mu } \sum _s m_s N^{\mu }_{s} = mnF^{\nu}{}_{\mu }U^{\mu } = 0 . \end{align}

Thus,

(2.20) \begin{align} F^{\nu}{}_{\mu }U^{\mu } = 0 \end{align}

is a leading-order degeneracy condition of the KMHD expansion. It does not imply that the comoving electric field vanishes at all orders. In particular, the subleading parallel electric field $E_{\parallel }$ enters the gyroaveraged drift-kinetic equation and the corresponding moment equations; see (2.25) and (2.35).

Finally, the electromagnetic field equations (2.3c ) imply charge neutrality,

(2.21) \begin{align} \sum _s q_s N^{\mu }_{s} = 0 \implies \sum _s q_s n_s = 0 , \qquad \sum _s q_s V^{\mu }_{s} = 0 . \end{align}

To show that the plasma is gyrotropic (2.24), we need to perform a coordinate transformation for the velocity variables. In the non-relativistic theory (Kulsrud Reference Kulsrud1983), one decomposes the velocity $k^{\mu }$ parallel and perpendicular to the magnetic field. In relativity, we use the fact that electrogmagnetic field tensor is magnetically dominated (2.18) to pick an orthonormal basis $(U^{\mu }, \hat {X}^{\mu }, \hat {Y}^{\mu }, \hat {b}^{\mu })$ and decompose the velocity vector $k^{\mu }$ as [see Appendix B.6 of Chael et al. (Reference Chael, Lupsasca, Wong and Quataert2023)]

(2.22) \begin{align} k^{\mu } = W_{U} U^{\mu } + v_{\perp } \cos (\vartheta ) \hat {X}^{\mu } + v_{\perp } \sin (\vartheta ) \hat {Y}^{\mu } + v_{\parallel } \hat {b}^{\mu }, \end{align}

where $b^{\mu }$ is the magnetic field in the fluid rest frame, $\vartheta$ is the gyroangle, $W_U$ is the Lorentz factor, $v_{\parallel }$ is the component of the velocity parallel to the magnetic field and $v_{\perp }$ is the component of the velocity perpendicular to the magnetic field. The covariant decomposition of these quantities is

(2.23a) \begin{align} &b^{\mu } \equiv U_{\nu } \left (\star F\right )^{\nu \mu } = \frac {1}{2} U_{\delta } \epsilon ^{\delta \mu \alpha \beta } F_{\alpha \beta } , \quad \hat {b}^{\mu } \equiv \frac {b^{\mu }}{|b|} ,\quad W_{U} \equiv - k_{\mu } U^{\mu } = \sqrt {1 + v_{\perp }^2 + v_{\parallel }^2} , \end{align}

(2.23b) \begin{align} &v_{\parallel } \equiv k_{\mu } \hat {b}^{\mu } ,\quad v_{\perp } = \sqrt {\big(k_{\mu } \hat {Y}^{\mu }\big)^2 + \big(k_{\mu } \hat {X}^{\mu }\big)^2} ,\quad \tan \vartheta = \frac {\big(\hat {Y}^{\mu } k_{\mu } \big)}{\big(\hat {X}^{\mu } k_{\mu } \big)}. \end{align}

As we show in Appendix A, one can change variables from $(x^{\mu },k) \to (x,v_{\parallel }, v_{\perp }, \vartheta )$ and show that (2.17) reduces to

(2.24) \begin{align} \frac {\partial f_{0,s}}{\partial \vartheta } = 0 , \end{align}

signifying that the plasma is gyrotropic.Footnote ¹

To understand the evolution of the distribution function in $(x,v_{\parallel }, v_{\perp })$ space, we continue the formal expansion of the distribution function (2.16) to the next order and average the Vlasov equation over $\vartheta$ (see Appendix A for details). The result is the relativistic generalisation of the gyroaveraged kinetic equation

(2.25) \begin{align} &W_{U} D_{\parallel } f_{0,s} + \frac {W_{U} v_{\perp }D_{\parallel } b }{2 b} \frac {\partial f_{0,s}}{\partial v_{\perp }} + \frac {\partial f_{0,s}}{\partial v_{\parallel }} \bigg [ \frac {q_s}{m_s} E_{\parallel } W_{U} - W_{U}^2 \hat {b}^{\beta } D_{\parallel } U_{\beta } - \frac {\hat {b}^{\alpha } v_{\perp }^2 \boldsymbol{\nabla} _{\alpha }b}{2 b} \bigg ] \nonumber \\[4pt] &= \mathcal{C}_{s}[f_{s}] , \end{align}

where $E_{\parallel } \equiv F^{\mu \nu } U_{\nu } \hat {b}_{\mu }$ is the parallel electric field and $D_{\parallel }$ is the derivative along the gyroaveraged particle trajectory,

(2.26) \begin{align} D_{\parallel } \equiv D + \frac {v_{\parallel } \hat {b}^{\alpha }}{W_U}\boldsymbol{\nabla} _{\alpha } = U^{\alpha }\boldsymbol{\nabla} _{\alpha } + \frac {v_{\parallel } \hat {b}^{\alpha }}{W_U}\boldsymbol{\nabla} _{\alpha } . \end{align}

To obtain the non-relativistic limit, we can set $W_{U} = 1$ and see that (2.25) reduces to (20.90) of Schekochihin (Reference Schekochihin2015).

The gyrotropic nature of the distribution function simplifies the number density current and the stress-energy tensor. Using (2.24) in (2.6), one can show that for a gyrotropic plasma,

(2.27a) \begin{align} &n_{s} = \left \langle W_{U} \right \rangle _{s} ,\\[-10pt]\nonumber \end{align}

(2.27b) \begin{align} &V^{\mu }_{s} = \langle v_{\parallel } \rangle _{s}\ \hat {b}^{\mu } \equiv V_{s} \hat {b}^{\mu } ,\\[-10pt]\nonumber \end{align}

(2.27c) \begin{align} &\mathcal{E}_{s} = m_{s} \left \langle W^2_{U} \right \rangle _{s} ,\\[-10pt]\nonumber \end{align}

(2.27d) \begin{align} &\mathcal{P}_{s} = \frac {m_{s}}{3} \left \langle v_{\perp }^2 + v_{\parallel }^2\right \rangle _{s} = \frac {2 p_{\perp ,s} + p_{\parallel ,s}}{3} ,\\[-10pt]\nonumber \end{align}

(2.27e) \begin{align} &\mathcal{Q}^{\mu }_{s} = m_{s} \left \langle W_{U} v_{\parallel } \right \rangle _{s} \hat {b}^{\mu } \equiv \mathcal{Q}_{s} \hat {b}^{\mu } ,\\[-10pt]\nonumber \end{align}

(2.27f) \begin{align} &\pi ^{\mu \nu }_{s} = m_{s} \left \langle v_{\parallel }^2 - \frac {v_{\perp }^2}{2} \right \rangle _{s} \hat {b}^{\lt \mu } \hat {b}^{\nu \gt } = \big(p_{\parallel ,s} - p_{\perp ,s}\big)\, \hat {b}^{\lt \mu } \hat {b}^{\nu \gt } , \end{align}

where the pressures parallel and perpendicular to the magnetic field are

(2.28) \begin{align} p_{\parallel ,s} \equiv m_s \langle v_{\parallel }^2\rangle _s , \quad p_{\perp ,s} \equiv \frac {m_s \langle v_{\perp }^2\rangle _s}{2}. \end{align}

Physically, (2.27) shows that particle and heat transport in the plasma ( $V^{\mu }_{s}$ and $\mathcal{Q}^{\mu }_{s}$ ) are always along the magnetic field in the fluid rest frame. The difference between $p_{\parallel ,s}$ and $p_{\perp ,s}$ drives a shear stress in the plasma through $\pi ^{\mu \nu }_{s}$ . In the ideal MHD limit, $p_{\parallel ,s} = p_{\perp ,s}$ and there is no anisotropic stress in the plasma.

2.3. Moments of the gyroaveraged equations

Before giving the evolution equations for the moments of the gyrotropic distribution function, we introduce some notation. We denote the most general moment of the gyrotropic distribution function as

(2.29) \begin{align} &\mathcal{I}^{(P,Q,R)}_{s} \equiv \left \langle v_{\perp }^{P} v_{\parallel }^{Q} W_{U}^{R} \right \rangle _{s} \nonumber \\[4pt] & \quad = \int _{v_{\perp }=0}^{\infty }\int _{v_{\parallel }=-\infty }^{\infty }\int _{\vartheta =0}^{2\pi } \frac {v_{\perp } {\rm d}v_{\perp } {\rm d}v_{\parallel } {\rm d}\vartheta }{\sqrt {1+v_{\perp }^2+v_{\parallel }^2}} f_{0,s} v_{\perp }^{P} v_{\parallel }^{Q} \left ( 1+v_{\perp }^2+v_{\parallel }^2\right )^{R/2} . \end{align}

For example, using the above-mentioned notation, the number density, energy density and pressure (2.7) are

(2.30) \begin{align} n_{s} = - \mathcal{I}^{(0,0,1)}_{s} ,\quad \mathcal{E}_{s} = m_s \mathcal{I}^{(0,0,2)}_{s} ,\quad \mathcal{P}_{s} = \frac {m_s}{3} \left (\mathcal{I}^{(2,0,0)}_{s} + \mathcal{I}^{(0,2,0)}_{s}\right ) . \end{align}

We also introduce short hands for certain common moments of the distribution function that appear in our equations [for definitions in the non-relativistic limit, see below (13) of Snyder et al. (Reference Snyder, Hammett and Dorland1997)],

(2.31a) \begin{align} &\mathcal{Q}_{\parallel ,s}^{({-}k)} \equiv m_s \mathcal{I}_{s}^{(0,3,-1-k)} , \quad \mathcal{Q}_{\perp ,s}^{({-}k)} \equiv \frac {1}{2} m_s \mathcal{I}_{s}^{(2,1,-1-k)},\\[-10pt]\nonumber \end{align}

(2.31b) \begin{align} &p_{\perp ,s}^{({-}k)} \equiv \frac {m_s \mathcal{I}_{s}^{(2,0,-k)}}{2} , \quad p_{\parallel ,s}^{({-}k)} \equiv m_s \mathcal{I}_{s}^{(0,2,-k)} ,\\[-10pt]\nonumber \end{align}

(2.31c) \begin{align} &r_{\parallel ,\perp ,s}^{({-}k)} \equiv \frac {1}{2}m_s \mathcal{I}_{s}^{(2,2,-k)} , \quad r_{\parallel ,\parallel ,s}^{({-}k)} \equiv m_s \mathcal{I}_{s}^{(0,4,-k)} , \quad r_{\perp ,\perp ,s}^{({-}k)} \equiv \frac {1}{4} m_s \mathcal{I}_{s}^{(4,0,-k)}. \end{align}

In the above-mentioned equations, when $k=0$ , we drop the superscript to simplify notation. For example, we denote $\mathcal{Q}_{\parallel ,s}^{({-}0)}$ as $\mathcal{Q}_{\parallel ,s}$ . Physically, $\mathcal{Q}_{\parallel ,s}^{({-}k)}$ is the (Lorentz factor weighted) heat flux parallel to the magnetic field, $\mathcal{Q}_{\perp ,s}^{({-}k)}$ is the (Lorentz factor weighted) heat flux perpendicular to the magnetic field and $(r_{\parallel ,\perp ,s}^{({-}k)}, r_{\parallel ,\parallel ,s}^{({-}k)}, r_{\perp ,\perp ,s}^{({-}k)})$ are Lorentz factor weighted fourth-order moments of the distribution function.

Let us now introduce operators that will be used to write down the evolution equations. Define

(2.32a) \begin{align} &O_{\parallel } \left [\mathcal{I}^{(P,Q,R)}_s\right ] \equiv \left \langle v_{\perp }^{P} v_{\parallel }^{Q} W_{U}^{R} \frac {\partial \log f_{0,s} }{\partial {v_{\parallel }}} \right \rangle _{s} = - Q \mathcal{I}^{(P,Q-1,R)}_s - (R-1) \mathcal{I}^{(P,Q+1,R-2)}_s ,\nonumber\\[4pt] \end{align}

(2.32b) \begin{align} &O_{\perp } \left [\mathcal{I}^{(P,Q,R)}_s\right ] \equiv \left \langle v_{\perp }^{P} v_{\parallel }^{Q} W_{U}^{R} \frac {\partial \log f_{0,s} }{\partial {v_{\perp }}} \right \rangle _{s} = - (P+1) \mathcal{I}^{(P-1,Q,R)}_s - (R-1) \mathcal{I}^{(P+1,Q,R-2)}_s ,\\[-10pt]\nonumber \end{align}

(2.32c) \begin{align} &O_{c}\left [\mathcal{I}_s^{(P,Q,R)}\right ] \equiv \int \frac {v_{\perp } {\rm d}v_{\perp } {\rm d}v_{\parallel } {\rm d}\vartheta }{\sqrt {1+v_{\perp }^2+v_{\parallel }^2}} \mathcal{C}_{s}[f_{0,s}] v_{\perp }^{P} v_{\parallel }^{Q} \left ( 1+v_{\perp }^2+v_{\parallel }^2\right )^{R/2} . \end{align}

The second equality in (2.32a ) and (2.32b ) is obtained by integrating by parts in $v_{\parallel }$ and $v_{\perp }$ , respectively. To obtain the evolution equation for the moment $\mathcal{I}_{s}^{(P,Q,R)}$ , we multiply the drift-kinetic equation (2.25) by $ v_{\perp }^{P} v_{\parallel }^{Q} W_{U}^{R-1}$ and integrate over velocity space. Before we present the results of this calculation, let us define

(2.33a) \begin{align} &I_0 \equiv \left \langle D_{\parallel } f_{0,s} v_{\perp }^{P} v_{\parallel }^{Q} W_{U}^{R} \right \rangle = D \mathcal{I}^{(P,Q,R)}_{s} + \hat {b}^{\alpha } \boldsymbol{\nabla} _{\alpha } \mathcal{I}_{s}^{(P,Q+1,R-1)} ,\\[-7pt]\nonumber \end{align}

(2.33b) \begin{align} &I_1 \equiv \left \langle \frac {v_{\perp }^{P +1} v_{\parallel }^{Q} W_{U}^{R}D_{\parallel } b }{2 b} \frac {\partial f_{0,s}}{\partial v_{\perp }} \right \rangle = \frac {D b}{2 b} O_{\perp } \left [\mathcal{I}^{(P+1,Q,R)}_{s}\right ] + \frac {\hat {b}^{\alpha } \boldsymbol{\nabla} _{\alpha } b}{2 b} O_{\perp } \left [\mathcal{I}^{(P+1,Q+1,R-1)}_{s}\right ] ,\\[-7pt]\nonumber \end{align}

(2.33c) \begin{align} &I_2 \equiv \left \langle v_{\perp }^{P} v_{\parallel }^{Q} W_{U}^{R} \frac {\partial f_{0,s}}{\partial v_{\parallel }} \bigg [ \frac {q_s}{m_s} E_{\parallel } - W_{U} \hat {b}^{\beta } D_{\parallel } U_{\beta } - \frac {\hat {b}^{\alpha } v_{\perp }{}^2 \boldsymbol{\nabla} _{\alpha }b}{2 b W_{U}} \bigg ] \right \rangle \nonumber \\[6pt] &= \frac {q_s}{m_s} E_{\parallel } O_{\parallel } \left [\mathcal{I}^{(P,Q,R)}_s\right ] - \hat {b}^{\beta } D U_{\beta } O_{\parallel } \left [\mathcal{I}^{(P,Q,R+1)}_s\right ] - \hat {b}^{\alpha } \hat {b}^{\beta } \boldsymbol{\nabla} _{(\alpha } U_{\beta )} O_{\parallel } \left [\mathcal{I}^{(P,Q+1,R)}_s\right ] \nonumber \\[6pt] & \quad - \frac {{\hat {b}}^{\alpha } \boldsymbol{\nabla} _{\alpha } b}{2 b} O_{\parallel } \left [\mathcal{I}^{(P+2,Q,R-1)}_s\right ] , \end{align}

where $I_0$ represents the advection of the moments parallel to the velocity and magnetic fields, and $I_1$ , $I_2$ are the contributions from the velocity derivative of the distribution function perpendicular and parallel to the magnetic field. The moment equations are the combined sum of $I_0$ , $I_1$ , $I_s$ and the contributions obtained from collisions,

(2.34) \begin{align} &I_0 + I_1 + I_2 = O_{c}\left [\mathcal{I}^{(P,Q,R-1)}\right ]. \end{align}

The previous equations describe the evolution equation for the most general moment of the gyrotropic distribution function. In this paper, we focus only on the evolution equation for the first few moments: the number density, the stress-energy tensor, the parallel pressure and the perpendicular pressure. These equations are

(2.35a) \begin{align} &\boldsymbol{\nabla} _{\alpha } \left (n_{s} U^{\alpha } + V^{\alpha }_{s} \right ) = 0,\\[-10pt]\nonumber \end{align}

(2.35b) \begin{align} &\boldsymbol{\nabla} _{\mu } T^{\mu \nu }_{s} = q_s F^{\nu}_{\mu } N^{\mu }_{s} ,\\[-10pt]\nonumber \end{align}

(2.35c) \begin{align} & D p_{\parallel ,s} + \hat {b}^{\alpha } \boldsymbol{\nabla} _{\alpha } \mathcal{Q}_{\parallel ,s} -\frac {\hat {b}^{\alpha }\boldsymbol{\nabla} _{\alpha }b \left (\mathcal{Q}_{\|, s}-2 \mathcal{Q}_{\perp ,s}\right )}{b} + \frac {Db r_{\|,\perp ,s}^{({-}2)}}{b} -\hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta } r_{\|,\|,s}^{({-}2)} \nonumber \\[3pt] &+E_{\|} q_s \left (\frac {\mathcal{Q}_{\|,s}^{({-}1)}}{m_s}-2 V_s\right ) +p_{\|,s} \left (3 \hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta }-\frac {Db}{b}\right ) + 2 Q_s \hat {b}^{\alpha }DU_{\alpha } = m_s\kern-1.3pt O_{c}\left [ \mathcal{I}^{(0,2,-1)}_s \right ]\! ,\\[-10pt]\nonumber \end{align}

(2.35d) \begin{align} &D p_{\perp ,s} + \hat {b}^{\alpha } \boldsymbol{\nabla} _{\alpha }\mathcal{Q}_{\perp ,s} -\hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta } r_{\|,\perp ,s}^{({-}2)} + \frac {Db r_{\perp ,\perp ,s}^{({-}2)}}{b}- \frac {2 \hat {b}^{\alpha }\boldsymbol{\nabla} _{\alpha }b \mathcal{Q}_{\perp ,s}}{b} +\frac {E_{\|} q_s \mathcal{Q}_{\perp ,s}^{({-}1)}}{m_s} \nonumber \\ &+p_{\perp ,s} \left (\hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta }-\frac {2 Db}{b}\right ) = \frac {m_s}{2}O_{c}\left [ \mathcal{I}^{(2,0,-1)}_s \right ] . \end{align}

The moments required for closure of the above-mentioned equations are listed in table 2. We emphasise that (2.35b ) retains the Lorentz force on each species. Although the leading-order KMHD ordering implies $F^{\nu}{}_{\mu }U^{\mu }=0$ and a vanishing leading-order perpendicular current, it does not imply that the full current vanishes. Indeed, summing (2.35b ) over species gives

(2.36) \begin{align} \boldsymbol{\nabla} _{\mu }T^{\mu \nu } = F^{\nu}{}_{\mu }J^{\mu } , \qquad J^{\mu } = \sum _s q_s N_s^{\mu } . \end{align}

The first non-vanishing perpendicular current arises at subleading order in the KMHD expansion and provides the usual MHD Lorentz force, see (2.40).

Table 2.

List of the moments needed to close (2.35).

Table 2. Long description

A table with two columns and four rows. The first column is labeled ‘Description’ and the second column is labeled ‘Moment’. The rows list different descriptions and their corresponding moments. Row 1: Description: Heat flux and particle diffusion, Moment: Q s and V s. Row 2: Description: Parallel heat flux, Moment: {Q parallel,s, Q parallel,-1,s}. Row 3: Description: Perpendicular heat flux, Moment: {Q perpendicular,s, Q perpendicular,-1,s}. Row 4: Description: Lorentz factor weighted pressure moments, Moment: {r perpendicular,-2,s, r parallel,-2,s, r parallel,-2,s}.

In the non-relativistic limit, the characteristic particle velocity $v_{\textrm {ch}}$ is much less than the speed of light $c$ . The higher-order moments, such as $r^{({-}2)}_{\parallel ,\parallel ,s}$ , that appear in (2.35c ) and (2.35d ) are suppressed relative to the pressure moments by additional powers of $v_{\textrm {ch}}/c$ , where $v_{\textrm {ch}}$ is the characteristic flow velocity. For example, restoring factors of $c$ ,

(2.37) \begin{align} r^{({-}2)}_{\parallel ,\parallel ,s} = m_s \left \langle \frac {v_{\parallel }^4}{c^2} W_U^{-2} \right \rangle _s \sim m_s \left \langle \frac {v_{\parallel }^4}{c^2} \right \rangle _s , \end{align}

where we have used $W_U^{-2}=1+\mathcal{O}(v_{\textrm {ch}}^2/c^2)$ . Since

(2.38) \begin{align} p_{\parallel ,s} \sim m_s \big\langle v_{\parallel }^2 \big\rangle _s , \end{align}

we have

(2.39) \begin{align} \frac {r^{({-}2)}_{\parallel ,\parallel ,s}}{p_{\parallel ,s}} \sim \mathcal{O} \left ( \frac {v_{\textrm {ch}}^2}{c^2} \right ) . \end{align}

Thus, the appearance of higher powers of the particle velocity in the negative-Lorentz-factor-weighted moments makes these moments higher order in the non-relativistic expansion. Neglecting such higher-order moments and setting $\mathcal{Q}=V=0$ in (2.35c ) results in equations identical to (16) and (17) of Snyder et al. (Reference Snyder, Hammett and Dorland1997).

We obtain the evolution equation for the species-averaged flow from (2.35) by summing over $s$ ,

(2.40a) \begin{align} \kern-25pt\boldsymbol{\nabla} _{\alpha } \left (n U^{\alpha } \right ) = 0, \kern260pt \end{align}

(2.40b) \begin{align} \kern-180pt\boldsymbol{\nabla} _{\mu } \left (T^{\mu \nu } + T^{\mu \nu }_{\mathrm{EM}}\right ) = 0 ,\kern110pt\end{align}

(2.40c) \begin{align} &D p_{\parallel } + \hat {b}^{\alpha } \boldsymbol{\nabla} _{\alpha } \mathcal{Q}_{\parallel } -\frac {\hat {b}^{\alpha }\boldsymbol{\nabla} _{\alpha }b \left (\mathcal{Q}_{\|}-2 \mathcal{Q}_{\perp }\right )}{b} + \frac {Db r_{\|,\perp }^{({-}2)}}{b} -\hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta } r_{\|,\|}^{({-}2)} \nonumber \\ &+E_{\|} \sum _{s} q_s \left (\frac {\mathcal{Q}_{\|,s}^{({-}1)}}{m_s}\right ) +p_{\|} \left (3 \hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta }-\frac {Db}{b}\right ) + 2 \mathcal{Q} \hat {b}^{\alpha }DU_{\alpha } = \sum _{s} m_s O_{c}\left [ \mathcal{I}^{(0,2,-1)}_s \right ] , \end{align}

(2.40d) \begin{align} &D p_{\perp } + \hat {b}^{\alpha } \boldsymbol{\nabla} _{\alpha }\mathcal{Q}_{\perp } -\hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta } r_{\|,\perp }^{({-}2)} + \frac {Db r_{\perp ,\perp }^{({-}2)}}{b}- \frac {2 \hat {b}^{\alpha }\boldsymbol{\nabla} _{\alpha }b \mathcal{Q}_{\perp }}{b} +E_{\|} \sum _{s}\frac { q_s \mathcal{Q}_{\perp ,s}^{({-}1)}}{m_s} \qquad\qquad \nonumber \\ &+p_{\perp } \left (\hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta }-\frac {2 Db}{b}\right ) = \sum _{s}\frac {m_s}{2}O_{c}\left [ \mathcal{I}^{(2,0,-1)}_s \right ] . \end{align}

We now analyse the previous equations to understand their implications. Equations (2.10a ) and (2.40b ) provide the conservation laws for the number density and the stress-energy tensor averaged over all the species in the plasma. The stress-energy tensor (2.27) depends on the parallel pressure, perpendicular pressure and heat flux. The evolution of the parallel and perpendicular pressures is given by (2.40c ) and (2.40d ), which contain unknown moments listed in table 2. In § 4, we analyse possible closures to consistently evolve the above-mentioned equations in the presence of heat fluxes. For now, we ignore the heat fluxes and show that our equations reduce to the familiar CGL equations.

2.4. Double adiabatic limit

Consider a simple closure for (2.40) (Chew et al. Reference Chew, Goldberger and Low1956), and ignore the contributions of collisions, heat fluxes and Lorentz factor weighted pressure moments in table 2. With these assumptions, the stress-energy tensor does not contain contributions from the heat flux, and (2.40c) and (2.40d) reduce to

(2.41a) \begin{align} &D p_{\parallel } +p_{\|} \left (3 \hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta }-\frac {Db}{b}\right ) = 0 ,\\[-10pt]\nonumber \end{align}

(2.41b) \begin{align} &D p_{\perp } +p_{\perp } \left (\hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta }-\frac {2 Db}{b}\right ) = 0. \end{align}

Using $b_{\beta } \boldsymbol{\nabla} _{\alpha } (\star F^{\alpha \beta } ) = 0$ and $\boldsymbol{\nabla} _{\mu }(n U^{\mu }) = 0$ , one can show that

(2.42) \begin{align} \hat {b}^{\alpha }\hat {b}^{\beta }\boldsymbol{\nabla} _{\alpha }U_{\beta } = \frac {Db}{b} - \frac {Dn}{n} . \end{align}

Hence, (2.41) reduces to

(2.43a) \begin{align} &\frac {n^3}{b^2} D \left (\frac {p_{\parallel } b^2}{n^3}\right ) = 0 ,\\[-10pt]\nonumber \end{align}

(2.43b) \begin{align} &n b D \left (\frac {p_{\perp }}{n b} \right ) = 0. \end{align}

The above-mentioned equations generalise the CGL equations of non-relativistic MHD to relativity (Gedalin & Oiberman Reference Gedalin and Oiberman1995). However, there is an important assumption underlying their derivation. While one can neglect the heat fluxes in table 2, one cannot always ignore the Lorentz factor weighted pressure moments for all equations of state. For non-relativistic flows and equations of state, these contributions can be safely ignored. For ultra-relativistic flows, however, these moments are of the same order as the parallel pressure moments and cannot be ignored (see § 4).

2.5. Summary of the important equations

The evolution of gyrotropic plasma in general relativity is described by the following equations.

1. (i) Eckhart frame choice (2.8),

   (2.44) \begin{align} \sum _{s} m_{s} V_{s} = 0. \end{align}

2. (ii) Conservation of laws for the averaged number density, stress-energy tensor across all species (2.40a ) and (2.40b ), and Gauss–Farday law (2.3b ),

   (2.45) \begin{align} \boldsymbol{\nabla} _{\alpha } \left (n U^{\alpha } \right ) = 0,\quad \boldsymbol{\nabla} _{\mu } \left (T^{\mu \nu } + T^{\mu \nu }_{\mathrm{EM}}\right ) = 0 , \quad \boldsymbol{\nabla} _{\mu } (\star F^{\mu \nu }) = 0 . \end{align}

3. (iii) A kinetic equation describing the evolution of the gyrotropic distribution function (2.25),Footnote ²

   (2.46) \begin{align} &W_{U} D_{\parallel } f_{s} + \frac {W_{U} v_{\perp }D_{\parallel } b }{2 b} \frac {\partial f_{s}}{\partial v_{\perp }} + \frac {\partial f_{s}}{\partial v_{\parallel }} \bigg [ \frac {q_s}{m_s} E_{\parallel } W_{U} - W_{U}^2 \hat {b}^{\beta } D_{\parallel } U_{\beta } - \frac {\hat {b}^{\alpha } v_{\perp }^2 \boldsymbol{\nabla} _{\alpha }b}{2 b} \bigg ] = \mathcal{C}_{s}[f_{s}]. \end{align}

4. (iv) Charge neutrality (2.21),

   (2.47) \begin{align} \sum _{s} q_{s} n_s = 0 , \quad \sum _{s} q_{s} V_{s} = 0. \end{align}

The kinetic equation derived directly couples to the macroscopic fluid equations. The general evolution equations for the moments of the gyroaveraged distribution function are listed in (2.34). Simplified evolution equations for the parallel and perpendicular pressures for species $s$ are given in (2.35c ) and (2.35d ). The evolution equations for the parallel and perpendicular moments, averaged across species and given in (2.40c ) and (2.40d ), provide a practical approach to evolving multi-species relativistic plasmas. The higher-order moments needed to close these equations are provided in table 2.

3.1. Equilibrium distribution function

Let us suppose that we are studying a gyrotropic plasma at rest, threaded by a constant magnetic field and a vanishing electric field in Minkowski space. We set up our coordinate system so the fluid velocity aligns with the time direction and the magnetic field points along the $z$ -direction,

(3.3) \begin{align} U^{\mu } = (1,0,0,0), \quad \hat {b}^{\mu } = (0,0,0,1). \end{align}

A global equilibrium state of the KMHD equations (2.46) without a parallel electric field is any function $f_{s}(v_{\parallel }, v_{\perp })$ . In non-relativistic theory, a common choice for the equilibrium distribution $f_{0,s}(v_{\parallel }, v_{\perp })$ is a bi-Maxwellian distribution function [see (23) of Snyder et al. (Reference Snyder, Hammett and Dorland1997)]

(3.4) \begin{align} f_{s}^{\mathrm{non-rel}} \propto \exp \bigg [ -\frac {m_s v_{\parallel }^2}{2 T_{\parallel ,s}} -\frac {m_s v_{\perp }^2}{2 T_{\perp ,s}} \bigg ], \end{align}

where $T_{\parallel ,s}$ and $T_{\perp ,s}$ are the temperatures parallel and perpendicular to the magnetic field. A bi-Maxwellian distribution function captures pressure anisotropies and one can analyse the linearised equations to obtain a closure for the heat fluxes (Snyder et al. Reference Snyder, Hammett and Dorland1997).

However, it is important to note that a bi-Maxwellian distribution function in relativity does not reduce to the equilibrium Maxwell-Juttner distribution function in the absence of pressure anisotropies. Therefore, we choose a distribution function inspired by anisotropic relativistic hydrodynamics (Alqahtani, Nopoush & Strickland Reference Alqahtani, Nopoush and Strickland2018),

(3.5) \begin{align} f_{s}(v_{\parallel }, v_{\perp }) = \alpha _s \exp \Big({-}z_s \sqrt {1 + v_{\parallel }^2 \delta _{1,s} + v_{\perp }^2 \delta _{2,s}} \Big) , \end{align}

where $\alpha _{s}$ is a constant, $z_{s} \sim m_{s}/T_{s}$ quantifies the ratio of the rest mass energy to the thermal energy $(T_s)$ of a particle and $\delta _{1,s}\sim T_{s}/T_{\parallel ,s}$ ( $\delta _{2,s}\sim T_{s}/T_{\perp ,s}$ ) quantifies the ratio of the total temperature to the temperature parallel (perpendicular) to the magnetic field. In the non-relativistic limit, (3.5) reduces to the bi-Maxwellian distribution function (3.4).

Apart from the moments $\mathcal{I}_{s}^{(P,Q,R)}$ (2.34), we will also need to evaluate the following moments of $f_s$ :

(3.6) \begin{align} \mathcal{I}^{({-}1)}_{(P,Q,R),s} \equiv \left \langle \frac {1}{\tilde {W}} v_{\perp }^{P} v_{\parallel }^{Q} W^{R} \right \rangle _s, \end{align}

where $\tilde {W} \equiv \sqrt {1 + v_{\parallel }^2 \delta _{1,s} + v_{\perp }^2 \delta _{2,s}}$ . Strategies for numerically obtaining the moments $\mathcal{I}^{(P,Q,R)}_s$ and $\mathcal{I}^{({-}1)}_{(P,Q,R),s}$ for general plasma temperature are in Appendix B. In the high-temperature limit, the moments can be evaluated analytically; see Appendix B.1 for explicit expressions for the first few moments of the distribution function.

3.2. Linearised equations

The linearised gyroaveraged kinetic equation (2.46) in the absence of collisions is

(3.7) \begin{align} D_{\parallel } \delta f_s + \frac {v_{\perp } D_{\parallel } \delta b }{2 b} \frac {\partial f_s}{\partial v_{\perp }} + \frac {\partial f_s}{\partial v_{\parallel }} \bigg [ \frac {q_s}{m_s} \delta E_{\parallel } - W_{U} \hat {b}^{\beta } D_{\parallel } \delta U_{\beta } - \frac {v_{\perp }^2 \hat {b}^{\alpha } \boldsymbol{\nabla} _{\alpha } \delta b}{2 b W_{U}} \bigg ] = 0 . \end{align}

Here, $\delta U^{\mu }$ denotes the dynamical fluid four-velocity perturbation defined by our choice of Eckhart fluid frame. This arises due to the changes in the bulk-motion of the fluid 4-velocity due to inertia and magnetic pressure. Its equilibrium value is provided in (3.3). To prevent clutter, we shall use $U^{\mu }$ to refer to its equilibrium value and $\delta U^{\mu }$ to denote the perturbation. Note that the velocity-space variables $v_{\parallel }$ and $v_{\perp }$ are defined with respect to the corresponding local fluid tetrad and remain well defined under this perturbation.

Now, we Fourier transformFootnote ³ in space and time, and replace

(3.8) \begin{align} \partial _t \to - i \omega , \partial _x \to i k_x , \partial _y \to i k_y , \partial _z \to i k_z . \end{align}

In this section, $k_{x,y,z}$ represent the wave numbers for the perturbations and should not be confused with the particle $4$ -velocity in the Vlasov–Maxwell equations. In Fourier space, we simplify (3.7) as

(3.9) \begin{align} \delta f_{s} = - \frac {v_{\perp } \delta b}{2 b} \frac {\partial f_s}{\partial v_{\perp }} + W_U \hat {b}^{\beta } \delta U_{\beta } \frac {\partial f_s}{\partial v_{\parallel }} - \frac {1}{D_{\parallel }} \frac {\partial f_s}{\partial v_{\parallel }} \bigg [ \frac {q_s}{m_s} \delta E_{\parallel } \bigg ] + \frac { v_{\perp }^2 }{ W_{U} D_{\parallel }} \frac {i k_z \delta b }{2 b} \frac {\partial f_s}{\partial v_{\parallel }} , \end{align}

where in Fourier space,

(3.10) \begin{align} D_{\parallel } = i \left ({-} \omega + \frac {v_{\parallel }}{W_U} k_z\right ) . \end{align}

Integrating (3.9) over velocity space gives the linearised moments of the system

(3.11) \begin{align} \delta \mathcal{I}^{(M,N,Q)}_s &= O_{\parallel }\left [ \mathcal{I}^{(M,N,Q+1)}_s \right ] \hat {b}^{\beta } \delta U_{\beta } -\frac {\delta b}{2 b} O_{\perp }\left [ \mathcal{I}^{(M+1,N,Q)}_s \right ] - \frac {q_s \mathrm{sgn}^{N}(k_z)}{m_s (i k_z)} \delta E_{\parallel } \mathcal{H}^{(M,N,Q)}_{s} \nonumber \\ & \quad + \frac {\mathrm{sgn}^N(k_z) \delta b}{2 b} \mathcal{H}^{(M+2,N,Q-1)}_{s} , \end{align}

where $\mathrm{sgn}(\boldsymbol{\cdot })$ is the sign function

(3.12) \begin{align} \zeta \equiv \frac {\omega }{|k_z|} , \end{align}

and

(3.13) \begin{align} \mathcal{H}^{(P,Q,R)}_{s}(\zeta ) \equiv \left \langle \frac {v_{\perp }^{P} v_{\parallel }^{Q} W_U^{R} }{v_{\parallel } W_{U}^{-1} - \zeta }\frac {\partial \log (f_s)}{\partial v_{\parallel }} \right \rangle _s = -z_s \delta _{1,s} \left \langle \frac {v_{\perp }^{P} v_{\parallel }^{Q+1} W_U^{R} }{ \left (v_{\parallel } W_{U}^{-1} - \zeta \right ) \tilde {W}} \right \rangle _s . \end{align}

We have used (3.5) to obtain the second equality. We note that the linearised plasma must have zero particle diffusion (3.2),

(3.14) \begin{align} \delta \mathcal{I}_{s}^{(0,1,0)} & = O_{\parallel }\left [ \mathcal{I}^{(0,1,1)}_s \right ] \hat {b}^{\beta } \delta U_{\beta } -\frac {\delta b}{2 b} O_{\perp }\left [ \mathcal{I}^{(1,1,0)}_s \right ] - \frac {q_s \mathrm{sgn}(k_z)}{m_s (i k_z)} \delta E_{\parallel } \mathcal{H}^{(0,1,0)}_{s} \nonumber \\ & \quad + \frac {\mathrm{sgn}(k_z) \delta b}{2 b} \mathcal{H}^{(2,1,-1)}_{s} = 0 . \end{align}

We eliminate $\hat {b}^{\beta } \delta U_{\beta }$ using (3.14) to get

(3.15) \begin{align} \delta \mathcal{I}^{(M,N,Q)}_s &= -\frac {q_s \delta E_{\parallel } \mathrm{sgn}(k_z)}{m_s (ik_z)} \bigg \{ \frac {\mathcal{H}_s^{(0,1,0)} }{n_s} O_{\parallel }\left [ \mathcal{I}^{(M,N,Q+1)}_s \right ] + \mathrm{sgn}^{N-1}(k_z) \mathcal{H}^{(M,N,Q)}_{s} \bigg \} \nonumber \\ & \quad + \frac {\delta b}{2b} \bigg \{ \mathrm{sgn}^N(k_z) \mathcal{H}^{(M+2,N,Q-1)}_{s} + \frac {\mathrm{sgn}(k_z) }{n_s} \mathcal{H}^{(2,1,-1)}_{s}O_{\parallel }\left [ \mathcal{I}^{(M,N,Q+1)}_s \right ] \nonumber \\ & \qquad\qquad\qquad\qquad\qquad\quad\qquad - O_{\perp }\left [ \mathcal{I}^{(M+1,N,Q)}_s \right ] \bigg \} . \end{align}

This equation will be used as a starting point in § 4 to obtain the closures for the moments in table 2. To evaluate the coefficients proportional to $\delta E_{\parallel }$ and $\delta b$ , we analyse (3.13). In the non-relativistic limit, (3.13) can be rewritten in terms of the plasma response function (Snyder et al. Reference Snyder, Hammett and Dorland1997). The same simplification applies in relativity. For integer $Q\gt 0$ , one can show that (see Appendix C)

(3.16) \begin{align} \mathcal{H}^{(P,Q,R)}_s = -z_s \delta _{1,s} \sum _{k=0}^{Q} \zeta ^k \, \mathcal{I}^{({-}1)}_{(P,Q-k,R+1+k),s} - z_s \delta _{1,s} \mathcal{Z}^{(P,Q+R+1)}_s \zeta ^{Q+1} , \end{align}

where $\mathcal{Z}^{(P,R)}_s(\zeta )$ is a relativistic generalisation of the plasma response function

(3.17) \begin{align} \mathcal{Z}^{(P,R)}_s(\zeta ) \equiv \left \langle \frac {v_{\perp }^P W^{R}}{\left (v_{\parallel } W_{U}^{-1} - \zeta \right ) \tilde {W}} \right \rangle _s . \end{align}

Given a value for $\zeta$ , we evaluate all terms proportional to $\delta E_{\parallel }$ and $\delta b/b$ in (3.15) using (3.16), (2.32a ) and (2.32b ). Useful identities for numerically evaluating $\mathcal{Z}^{(P,R)}_s(\zeta )$ , including analytical results in the high-temperature limit, are discussed in Appendix C.

3.3. High-temperature and low-frequency limit

In this section, we provide analytical expressions for the first few moments of the linearised distribution function (3.15) in the high-temperature $z_{s}\ll 1$ and low-frequency limit $\zeta \ll 1$ . Analytical expressions for the moments of the distribution function in this limit are given in (B7). The linearised moments from the KMHD equations in this limit can be obtained by using (3.15),

(3.18a) \begin{align} & \frac {\delta n_s}{n_s} = \frac {\delta b}{b } \left [-\tilde {\mu }-\frac {3}{4} i \pi \zeta (\tilde {\mu }+1)^{3/2}\right ] \nonumber \\[5pt] & + \frac {q_s z_s \delta E_{\|}}{m_s k_z} \left [\frac {1}{2} \pi \zeta (\tilde {\mu }+1)^{3/2} \sqrt {\delta _2}-\frac {i \sqrt {\tilde {\mu }+1} \sqrt {\delta _2} \big(\sqrt {\tilde {\mu }}+(\tilde {\mu }+1) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)\big)}{2 \sqrt {\tilde {\mu }}}\right ]\! ,\nonumber\\[5pt]\end{align}

(3.18b) \begin{align} &\frac {\delta p_{\parallel ,s}}{p_{\parallel ,s}} = \frac {\delta b}{b } \left [ \frac {\left ({-}\tilde {\mu }^2+2 \tilde {\mu }+3\right ) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)-\sqrt {\tilde {\mu }} (\tilde {\mu }+3)}{2 (\tilde {\mu }+1) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)-2 \sqrt {\tilde {\mu }}} \right ] \nonumber \\[5pt] & + \frac {q_s z_s \delta E_{\|}}{m_s k_z} \left [ -\frac {2 i \tilde {\mu }^{3/2} \sqrt {\tilde {\mu }+1} \sqrt {\delta _2}}{3 (\tilde {\mu }+1) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)-3 \sqrt {\tilde {\mu }}}\right ] ,\\[-5pt]\nonumber \end{align}

(3.18c) \begin{align} & \frac {\delta p_{\perp ,s}}{p_{\perp ,s}} = \frac {\delta b}{b } \bigg [ \frac {\left ({-}3 \tilde {\mu }^2+2 \tilde {\mu }-3\right ) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)-3 (\tilde {\mu }-1) \sqrt {\tilde {\mu }}}{2 \big(\sqrt {\tilde {\mu }}+(\tilde {\mu }-1) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)\big)} \nonumber \\[5pt] & \qquad\qquad\qquad -\frac {2 i \pi \zeta \tilde {\mu }^{3/2} (\tilde {\mu }+1)}{\sqrt {\tilde {\mu }}+(\tilde {\mu }-1) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)} \bigg ] \nonumber \\[5pt] &+ \frac {q_s z_s \delta E_{\|}}{m_s k_z} \left [ \frac {\pi \zeta \tilde {\mu }^{3/2} (\tilde {\mu }+1) \sqrt {\delta _2}}{\sqrt {\tilde {\mu }}+(\tilde {\mu }-1) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)}-\frac {4 i \tilde {\mu }^{3/2} \sqrt {\tilde {\mu }+1} \sqrt {\delta _2}}{3 \big(\sqrt {\tilde {\mu }}+(\tilde {\mu }-1) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)\big)} \right ] ,\nonumber\\ \\[-10pt]\nonumber \end{align}

(3.18d) \begin{align} & \frac {\delta \mathcal{Q}_{\parallel ,s}}{\mathrm{sgn}(k_z)} = \frac {\delta b}{b }\frac {\pi \alpha \zeta m_s}{z_s^4} \left [ \frac {6 \left (5 \tilde {\mu }^2+14 \tilde {\mu }+9\right ) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)}{\tilde {\mu }^{5/2} \delta _2{}^2}-\frac {6 (11 \tilde {\mu }+9)}{\tilde {\mu }^2 \delta _2{}^2} \right ] \nonumber \\[5pt] & + \frac {\pi \alpha \zeta q_s z_s E_{\|}}{k_z} \bigg [ -\frac {12 i \sqrt {\tilde {\mu }+1} \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)}{\tilde {\mu }^{3/2} \delta _2{}^{3/2} z_s^4}+\frac {18 i (\tilde {\mu }+1)^{5/2} \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)^2}{\tilde {\mu }^3 \delta _2{}^{3/2} z_s^4} \nonumber \\[5pt] & \qquad\qquad\qquad\qquad -\frac {2 i \big(4 \tilde {\mu }^2+15 \tilde {\mu }+9\big)}{\tilde {\mu }^2 \sqrt {\tilde {\mu }+1} \delta _2{}^{3/2} z_s^4} \bigg ] ,\\[-10pt]\nonumber \end{align}

(3.18e) \begin{align} & \frac {\delta \mathcal{Q}_{\perp ,s}}{\mathrm{sgn}(k_z)} = \frac {\delta b}{b }\frac {\pi \alpha \zeta m_s}{z_s^4} \left [ -\frac {3 (\tilde {\mu }-9) (\tilde {\mu }+1)}{\tilde {\mu }^2 \delta _2{}^2}-\frac {3 \left (\tilde {\mu }^2+2 \tilde {\mu }+9\right ) (\tilde {\mu }+1) \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)}{\tilde {\mu }^{5/2} \delta _2{}^2} \right ] \nonumber \\[5pt] & + \frac {\pi \alpha \zeta q_s z_s E_{\|}}{k_z} \bigg [ \frac {6 i (\tilde {\mu }+1)^{3/2} \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)}{\tilde {\mu }^{3/2} \delta _2{}^{3/2} z_s^4}+\frac {3 i (\tilde {\mu }-3) (\tilde {\mu }+1)^{5/2} \tan ^{-1}\big(\sqrt {\tilde {\mu }}\big)^2}{\tilde {\mu }^3 \delta _2{}^{3/2} z_s^4} \nonumber \\[5pt] & +\frac {i \left ({-}8 \tilde {\mu }^2+3 \tilde {\mu }+9\right ) \sqrt {\tilde {\mu }+1}}{\tilde {\mu }^2 \delta _2{}^{3/2} z_s^4} \bigg ] ,\\[-10pt]\nonumber \end{align}

(3.18f) \begin{align} &\delta \mathcal{Q}_{s} = \delta \mathcal{Q}_{\parallel ,s} + 2\delta \mathcal{Q}_{\perp ,s} . \end{align}

The quantities $\tilde {\mu }$ and $\tilde {\Delta }$ are defined in (B6). The above-mentioned expressions will be analysed later to provide a Landau-fluid closure for the heat fluxes; observe that we have already obtained a closure for $\mathcal{Q}_s$ in (3.18f ).

4.1. Closure strategy

If collisions dominate a kinetic system, then the distribution function evolves towards a local equilibrium characterised by macroscopic properties such as fluid velocity, energy density and pressure. For small departures from local equilibrium in a collisional system, there are well-developed strategies such as the Chapman–Enskog approximation (Chapman & Cowling Reference Chapman and Cowling1991) or Grad’s moment approximation (Grad Reference Grad1949) that allow one to systematically understand the evolution of the moments of the distribution functions and calculate non-equilibrium contributions such as viscous stresses in the system. For a general discussion of these techniques in the relativistic context, see Denicol et al. (Reference Denicol, Niemi, Molnár and Rischke2012, Reference Denicol, Huang, Molnár, Monteiro, Niemi, Noronha, Rischke and Wang2018) and Denicol & Rischke (Reference Denicol and Rischke2021); for numerical simulations of black hole accretion flow in the collisional Braginskii limit, see Chandra et al. (Reference Chandra, Gammie, Foucart and Quataert2015), Foucart et al. (Reference Foucart, Chandra, Gammie and Quataert2015, Reference Foucart, Chandra, Gammie, Quataert and Tchekhovskoy2017) and Dhruv et al. (Reference Dhruv, Prather, Chandra, Joshi and Gammie2025a ); for discussions of collisional closure in multi-species magnetised plasma, see Most & Noronha (Reference Most and Noronha2021) and Most et al. (Reference Most, Noronha and Philippov2022).

The plasmas studied in this paper are weakly collisional (2.14). To obtain a closure for the heat fluxes for these systems, we adopt the Landau fluid closure introduced by Hammett & Perkins (Reference Hammett and Perkins1990), Snyder et al. (Reference Snyder, Hammett and Dorland1997) and Goswami et al. (Reference Goswami, Passot and Sulem2005).

1. (i) First, define a set of master variables. In this paper, these variables will be $(\mathcal{E}_s/n_s, p_{\perp ,s} - p_{\parallel ,s}, b, n_s )$ .

2. (ii) Choose a moment of the distribution function that requires a closure. For illustration, suppose that this is $\mathcal{Q}_{\parallel ,s}$ .

3. (iii) Obtain the linearised expression for $\delta \mathcal{Q}_{\parallel ,s}$ using (3.15) in the low-frequency limit. Schematically, this is of the form

   (4.1) \begin{align} \delta \mathcal{Q}_{\parallel ,s} = \left [s_1 + s_2 \zeta \right ] \delta E_{\parallel } + \left [s_1 + s_2 \zeta \right ] \delta b , \end{align}
   where the coefficients $s_{i}$ depend on the equilibrium distribution function (3.5).

4. (iv) Consider a closure of the form

   (4.2) \begin{align} \delta \mathcal{Q}_{\parallel ,s} = i \mathrm{sgn}(k_z) \bigg [r_1 n_s \delta \left (\frac {\mathcal{E}_s}{n_s}\right ) + r_2 \delta \left (p_{\perp ,s} - p_{\parallel ,s}\right ) + r_3 \frac {\delta b}{b} \mathcal{E}_s + r_4 \mathcal{E}_s \frac {\delta n_s}{n_s} \bigg ] . \end{align}
   Equate (4.1) and (4.2), and determine $r_{i}$ .

In the non-relativistic limit, one can interpret this strategy as a procedure to build Padé approximations to the plasma response function, see Snyder et al. (Reference Snyder, Hammett and Dorland1997), Goswami et al. (Reference Goswami, Passot and Sulem2005) and Hunana et al. (Reference Hunana, Tenerani, Zank, Goldstein, Webb, Khomenko, Collados, Cally, Adhikari and Velli2019). This interpretation works at the linear level, mainly due to the choice of a bi-Maxwellian distribution function which allows one to rewrite all the linearised moments in terms of the plasma response function, see (36)–(38) of Snyder et al. (Reference Snyder, Hammett and Dorland1997). In relativity, we interpret this approach as a way to construct an effective field theory for heat fluxes that closely reproduces the kinetic dispersion relation and captures Landau damping.

4.2. Analytical closure in the ultra-relativistic limit

If we implement the Landau fluid closure strategy outlined previously to (3.18), we see that the heat fluxes satisfy

(4.3a) \begin{align} \delta \mathcal{Q}_{s} &= \delta \mathcal{Q}_{\parallel ,s} + 2\delta \mathcal{Q}_{\perp ,s} ,\\[-10pt]\nonumber \end{align}

(4.3b) \begin{align} \delta \mathcal{Q}_{\parallel ,s} &= i \mathrm{sgn}(k_z) \bigg [ n_s \left ( -\frac {8}{5\pi } + f_1 \right ) \delta \left (\frac {\mathcal{E}_s}{n_s}\right ) + \left ( \frac {28}{15 \pi } + f_2 \right ) \delta (p_{\perp ,s}-p_{\parallel ,s}) \nonumber \\[7pt] & \quad + f_3 \frac {\delta b}{b} \mathcal{E}_s + f_4 \mathcal{E}_s \frac {\delta n_s}{n_s} \bigg ] ,\\[-10pt]\nonumber \end{align}

(4.3c) \begin{align} \delta \mathcal{Q}_{\perp ,s} &= i \mathrm{sgn}(k_z) \bigg [ n_s \left ( -\frac {8}{15\pi } + g_1 \right ) \delta \left (\frac {\mathcal{E}_s}{n_s}\right ) + \left ( -\frac {4}{15 \pi } + g_2 \right ) \delta (p_{\perp ,s}-p_{\parallel ,s}) \nonumber \\[7pt] & \quad + g_3 \frac {\delta b}{b} \mathcal{E}_s + g_4 \mathcal{E}_s \frac {\delta n_s}{n_s} \bigg ] .\\[-10pt]\nonumber \end{align}

The coefficients $f_{1,2,3,4}$ and $g_{1,2,3,4}$ are complicated functions of $\tilde {\mu }$ and vanish when the initial anisotropy vanishes. We list these coefficients in the supplementary Mathematica file. One could repeat the Landau fluid procedure to obtain an accurate closure for $r_{\perp ,\perp }^{({-}2)} , r_{\parallel ,\parallel }^{({-}2)}$ and $r_{\parallel ,\perp }^{({-}2)}$ . However, we find that it is easier to use the equilibrium relationship for these functions. In the high-temperature limit $(z\to 0 \implies m \ll T)$ , these are given by

(4.4) \begin{align} r_{\perp ,\perp ,s}^{({-}2)} &= \mathcal{E}_s \left [ \frac {2}{15} + h_1 \right ] , \quad r_{\parallel ,\parallel ,s}^{({-}2)} = \mathcal{E}_s \left [ \frac {1}{5} + h_2 \right ] , \quad r_{\parallel ,\perp ,s}^{({-}2)} = \mathcal{E}_s \left [ \frac {1}{15} + h_3 \right ] , \end{align}

where $h_{1,2,3}$ are complicated functions of $\tilde {\mu }$ that vanish when the pressure anisotropy is zero. These coefficients are provided in the supplementary Mathematica file.

Equations (4.3), (4.4) and (3.2) provide the closure relationship for all the moments listed in table 2 in the ultra-relativistic limit for small pressure anisotropy. These Landau fluid closures can be used to evolve the fluid equations (2.40) instead of the kinetic equations to capture weakly collisional effects; see Sharma et al. (Reference Sharma, Hammett, Quataert and Stone2006), Sharma et al. (Reference Sharma, Quataert, Hammett and Stone2007), Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) for numerical simulations in the non-relativistic limit.

Before we proceed, let us briefly compare the heat flux closure in the ultra-relativistic limit (4.3) with the non-relativistic limit from (48) and (49) of Snyder et al. (Reference Snyder, Hammett and Dorland1997),

(4.5a) \begin{align} &\delta \mathcal{Q}_{\parallel }^{\mathrm{non-rel}} = - \sqrt {\frac {8}{\pi }} \frac {n_s i \mathrm{sgn}(k_z) }{\sqrt {z_{\parallel ,s}} } \delta \left (\frac { p_{\parallel ,s}}{n_s}\right ) , \quad z_{\parallel } \equiv n_s/p_{\parallel ,s} ,\\[-10pt]\nonumber \end{align}

(4.5b) \begin{align} &\delta \mathcal{Q}_{\perp }^{\mathrm{non-rel}} = \ - \sqrt {\frac {2}{\pi }} \frac {n_s i \mathrm{sgn}(k_z) }{ \sqrt {z_{\parallel ,s}} } \bigg [ \delta \left (\frac { p_{\perp }}{n_s}\right ) + \frac {m_s (p_{\perp ,s} - p_{\parallel ,s})}{p_{\parallel ,s}} \frac {\delta b}{b z_{\perp ,s}} \bigg ] , z_{\perp } \equiv n_s/p_{\perp ,s} . \end{align}

Observe that the non-relativistic parallel (perpendicular) heat flux is driven by the gradients in the parallel (perpendicular) pressure. However, in the relativistic heat fluxes, the anisotropic distribution is spheroidal (3.5), and both parallel and perpendicular pressure gradients contribute.

In relativistic flows applicable to SMBHs, the electrons are often ultra-relativistic and one can use the ultra-relativistic closure relation in (4.4) to evolve their heat fluxes. For the ions, which are typically non-relativistic or only mildly relativistic, one can instead use the covariant form of the non-relativistic closure in (4.5). More generally, the same matching procedure used in this section can, in principle, be repeated at finite temperature (Denicol, Jeon & Gale Reference Denicol, Jeon and Gale2014). In that case, the relevant kinetic response functions and moment integrals depend on the dimensionless temperature $z\sim m/T$ , and the closure coefficients need not take the simple analytic form found in the ultra-relativistic limit. They may instead be obtained by constructing asymptotic expansions about the non-relativistic and ultra-relativistic limits (Denicol et al. Reference Denicol, Jeon and Gale2014), or by evaluating the finite-temperature response functions numerically. Hence, the ultra-relativistic and non-relativistic closures should be viewed as limiting cases of a more general finite-temperature relativistic Landau-fluid closure.

Figure 1.

Comparison between the linear number density (top panels) and energy response (bottom panels) between the KMHD (solid blue line), CGL (dashed purple line), MHD (dash dotted green line) and the Landau fluid closure (bold orange line). as a function of the driving speed $\zeta = \omega /|k_z|$ . The left and right columns compare the real and imaginary parts of the response. Observe from the right column that the Landau fluid closure qualitatively models the Landau damping predicted by KMHD equations. This damping is not present in the CGL or the MHD models.

Figure 1. Long description

Panel A: The top left graph is a line graph showing the real part of the linear number density response. The x-axis is labeled with the Greek letter zeta, and the y-axis is labeled Re[delta log(ns)/delta log(b)]. Four lines represent different models: KMHD (solid blue line), CGL (dashed purple line), MHD (dash-dotted green line), and Landau fluid closure (bold orange line). Panel B: The top right graph is a line graph showing the imaginary part of the linear number density response. The x-axis is labeled with the Greek letter zeta, and the y-axis is labeled -Im[delta log(ns)/delta log(b)]. The same four lines represent the different models. Panel C: The bottom left graph is a line graph showing the real part of the energy response. The x-axis is labeled with the Greek letter zeta, and the y-axis is labeled Re[delta log(es)/delta log(b)]. The same four lines represent the different models. Panel D: The bottom right graph is a line graph showing the imaginary part of the energy response. The x-axis is labeled with the Greek letter zeta, and the y-axis is labeled -Im[delta log(es)/delta log(b)]. The same four lines represent the different models.

4.3. Comparison between kinetic and fluid response

In this section, we compare the density and energy density responses obtained from the Landau fluid closure and with the full KMHD calculation. For simplicity, as in the non-relativistic calculation (Snyder et al. Reference Snyder, Hammett and Dorland1997), we choose to compare the responses when the initial distribution function is isotropic. In this case, all the $\Delta p$ dependent terms drop out from the Landau fluid closure relations. The linearised Landau fluid system can be obtained from (2.35) to obtain a response equation similar to the KMHD result, see (3.18). The exact results for the fluid response from the MHD, CGL and Landau fluid closures are listed in Appendix D. We also include a comparison between the fluid response and the kinetic response in case of an initially anisotropic distribution in Appendix E.

The top and bottom panels of figure 1 compare the KMHD response results with the Landau fluid closure for the number density response and the energy density, respectively, as a function of the driving speed $\zeta = \omega /|k_z|$ . The left and right panels of each row plot the real and imaginary part of the response functions. Observe from the left column that the real part of the response functions for the number density and the energy density match the KMHD response function in the low frequency regime. The MHD (green dash dotted line) and CGL (purple dashed line) models have resonances at sound speeds $\zeta = \sqrt {1/3}$ and $\zeta = \sqrt {3/5}$ , respectively. This resonant behaviour is absent in the KMHD model (blue solid line) and the Landau fluid closure (bold orange line). From the right column of the figure, we see that the Landau fluid closure approximately captures the Landau damping predicted by the KMHD equations, a feature that this not present in either the CGL or the MHD equations. However, the Landau fluid closure over damps the system at $\zeta \sim 0.65$ . This behaviour is also observed in the non-relativistic theory; compare figure 1 with figure 3 of Snyder et al. (Reference Snyder, Hammett and Dorland1997). For an initially anisotropic distribution function, we observe results qualitatively similar to figure 1, see figure 2 in Appendix E.

The results from figure 1 show that the Landau fluid closure obtained in the previous section reproduces the dissipative KMHD response in the low-frequency limit and qualitatively captures the Landau damping predicted by the kinetic equations. At intermediate and high phase speeds, however, quantitative differences between the Landau fluid closure and the exact KMHD response remain. The agreement is expected to improve as more moments are retained in the closure relationship. While this has been demonstrated in the non-relativistic limit [see figure 3 of Snyder et al. (Reference Snyder, Hammett and Dorland1997)], here, we do not explore such higher-order closures, as it is unclear at present how to incorporate them consistently into a nonlinear relativistic simulation. Moreover, in the absence of a precise order-counting scheme, improved convergence in linear theory does not necessarily imply improved behaviour in a nonlinear simulation. We therefore choose a simple closure that captures some of the essential kinetic physics while keeping the resulting fluid model tractable.

4.4. Collisional operator and phenomenological closure

Let us now include the effects of collisions using a simplified model. We ignore inter-species collisions and use a simple linear operator operator introduced by Rocha, Denicol & Noronha (Reference Rocha, Denicol and Noronha2021) to model the collisions,

(4.6) \begin{align} \mathcal{C}_s^{\mathrm{RTA}}[f_s] &= -\frac {W_U}{\tau _{R,s}} \bigg [ \boldsymbol{\delta } f_s - f_{\mathrm{eq}}^{s} \frac {\left \langle W_U \delta f\right \rangle }{\left \langle W_U f_{\mathrm{eq}}^{s} \right \rangle } - P_{1} \frac {f_{\mathrm{eq}}^{s} \left \langle W_U P_{1} \boldsymbol{\delta } f_s \right \rangle }{\left \langle W_U P_{1} P_{1} f_{\mathrm{eq}}^{s} \right \rangle } \nonumber \\[4pt] & \quad - k^{\lt \mu \gt } \frac { f_{\mathrm{eq}}^{s} \left \langle W_U k_{\lt \mu \gt } \boldsymbol{\delta } f_s \right \rangle }{(1/3) \left \langle W_U k_{\lt \nu \gt } k^{\lt \nu \gt } f_{\mathrm{eq}}^{s} \right \rangle } \bigg ] , \end{align}

where $f_{\mathrm{eq}}^{s}$ is the Maxwell–Juttner distribution function,

(4.7) \begin{align} f_{\mathrm{eq}}^{s} = \alpha _{s,\mathrm{eq}} \exp \left({-} z_{s,\mathrm{eq}} \sqrt {1 + v_{\parallel }^2 + v_{\perp }^2}\right), \end{align}

$\boldsymbol{\delta } f_s = f_s - f_{\mathrm{eq}}^{s}$ , $\tau _{R,s}$ is a relaxation time and

(4.8) \begin{align} P_1 &\equiv 1-\frac {\left \langle W_U f_{\mathrm{eq}}^{s} \right \rangle }{\left \langle W_U^2 f_{\mathrm{eq}}^{s} \right \rangle } W_U . \end{align}

Suppose that the plasma is ultra relativistic and choose the equilibrium function so that it satisfies

(4.9) \begin{align} \left \langle W_U f_{\mathrm{eq}}^{s} \right \rangle = \left \langle W_U f_s\right \rangle , \left \langle W_U^2 f_{\mathrm{eq}}^{s} \right \rangle = \left \langle W_U^2 f_s\right \rangle \implies \alpha _{s,\mathrm{eq}} = \frac {\alpha _s}{\sqrt {\delta _1} \delta _2} , z_{s,\mathrm{eq}} = z . \end{align}

These assumptions simplify the collisional operator to

(4.10) \begin{align} \mathcal{C}_s^{\mathrm{RTA}}[f_s] &= -\frac {W_U}{\tau _{R,s}} \bigg [ \boldsymbol{\delta } f_s - k^{\lt \mu \gt } \frac { f_{\mathrm{eq}}^{s} \left \langle W_U k_{\lt \mu \gt } \boldsymbol{\delta } f_s \right \rangle }{(1/3) \left \langle W_U k_{\lt \nu \gt } k^{\lt \nu \gt } f_{\mathrm{eq}}^{s} \right \rangle } \bigg ] . \end{align}

For the gyro-averaged kinetic equation, we need the gyro-averaged collisional operator which takes the form

(4.11) \begin{align} \mathcal{C}_s[f_s] = \frac {1}{2\pi } \int {\rm d}\vartheta \, \mathcal{C}_s^{\mathrm{RTA}}[f_s] = -\frac {W_U}{\tau _{R,s}} \bigg [ \boldsymbol{\delta } f_s - v_{\parallel } \frac { f_{\mathrm{eq}}^{s} \mathcal{Q}_s}{p_s } \bigg ] . \end{align}

We can use (4.11) to simply the $O_{c}[ \mathcal{I}^{(0,2,-1)}_s ]$ and $O_{c}[ \mathcal{I}^{(2,0,-1)}_s ]$ appearing in (2.40c) and (2.40d ),

(4.12a) \begin{align} O_{c}\left [ \mathcal{I}^{(0,2,-1)}_s \right ] &= -\frac {2}{3 \tau _{R,s} m_s} \left (p_{\parallel ,s} - p_{\perp ,s}\right ) ,\\[-10pt]\nonumber \end{align}

(4.12b) \begin{align} O_{c}\left [ \mathcal{I}^{(2,0,-1)}_s \right ] &= \frac {2}{3 \tau _{R,s} m_s} \left (p_{\parallel ,s} - p_{\perp ,s}\right ) . \end{align}

These equations are exactly equivalent to their non-relativistic counter parts, see (16) and (17) of Snyder et al. (Reference Snyder, Hammett and Dorland1997).

It is possible to include the collisional operator in the kinetic equation to obtain a collisional closure. However, performing this exactly is not feasible because the collisional operator introduces poles in the Green’s function that are hard to treat. In the non-relativistic limit, Snyder et al. (Reference Snyder, Hammett and Dorland1997) followed an alternate approach: first deriving a collisionless closure for higher order moments of the distribution function, then taking the low-frequency limit with collisions included to obtain the closure for the heat fluxes. This procedure can also be carried out in relativity, but one typically needs to close more moments than in the non-relativistic limit. The result of this calculation will be a closure of the form

(4.13a) \begin{align} \delta \mathcal{Q}_{\parallel ,s} &= i \frac {|k_z|}{k_z + \tau _{R,s}^{-1} \boldsymbol{a}_{\parallel ,s}} \bigg [ n_s \left ( -\frac {8}{5\pi } + f_1 \right ) \delta \left (\frac {\mathcal{E}_s}{n_s}\right ) + \left ( \frac {28}{15 \pi } + f_2 \right ) \delta (p_{\perp ,s}-p_{\parallel ,s}) \nonumber \\[5pt]& \quad + f_3 \frac {\delta b}{b} \mathcal{E}_s + f_4 \mathcal{E}_s \frac {\delta n_s}{n_s} \bigg ] ,\\[-5pt]\nonumber \end{align}

(4.13b) \begin{align} \delta \mathcal{Q}_{\perp ,s} &= i \frac {|k_z|}{k_z + \tau _{R,s}^{-1} \boldsymbol{a}_{\perp ,s}} \bigg [ n_s \left ( -\frac {8}{15\pi } + g_1 \right ) \delta \left (\frac {\mathcal{E}_s}{n_s}\right ) + \left ( -\frac {4}{15 \pi } + g_2 \right ) \delta (p_{\perp ,s}-p_{\parallel ,s}) \nonumber \\[5pt] & \quad + g_3 \frac {\delta b}{b} \mathcal{E}_s + g_4 \mathcal{E}_s \frac {\delta n_s}{n_s} \bigg ] , \end{align}

where $\boldsymbol{a}_{\parallel ,s}$ and $\boldsymbol{a}_{\perp ,s}$ are $\mathcal{O}(1)$ coefficients. Notice that when $\tau _{R,s} \to \infty$ , the above-mentioned equation reduces to the collisionless closure (4.3). In a numerical implementation, one can simply treat $\boldsymbol{a}_{\parallel ,s}$ and $\boldsymbol{a}_{\perp ,s}$ as phenomenological parameters that decrease the pressure anisotropy (Sharma et al. Reference Sharma, Hammett, Quataert and Stone2006; Squire et al. Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023).

Before proceeding, it is useful to comment on the locality, stability and causality properties of the closure in (4.13). The closure was derived from the linearised equations in Fourier space. In the collisionless limit, the appearance of factors such as $\mathrm{sgn}(k_z)$ implies that, when transformed back to physical space, the closure is non-local along the magnetic field. This is the familiar Hilbert-transform-like non-locality of Landau-fluid closures (Hammett & Perkins Reference Hammett and Perkins1990; Snyder et al. Reference Snyder, Hammett and Dorland1997).

This non-locality should not be interpreted as a fundamental inconsistency of the underlying relativistic kinetic theory. Rather, it reflects the fact that the closure is a low-frequency representation of collisionless phase mixing. In practical fluid simulations, one often replaces the non-local operator by a local phenomenological approximation involving a characteristic parallel wavenumber $k_z$ , see Sharma et al. (Reference Sharma, Hammett, Quataert and Stone2006).

Even after such a replacement, the resulting equations should not be viewed as a manifestly causal, strictly hyperbolic relativistic fluid theory. This is expected, since the closure was constructed in a low-frequency approximation rather than from a relaxation-type causal theory. A practical way to handle this limitation is to use an order-reduction scheme: time derivatives appearing in the heat-flux sector are replaced using the lower-order CGL equations, thereby expressing the heat-flux corrections in terms of spatial derivatives of the leading-order variables. Such an order-reduced system can be evolved stably provided the heat-flux corrections remain perturbative, i.e. small compared with the leading CGL terms.

4.5. Comparisons with models in the literature

Before we compare our model to other relativistic flow models, let us briefly place it in the context of non-relativistic models. Historically, the description of magnetised, weakly collisional plasmas relied on the Braginskii fluid equations (Braginskii Reference Braginskii1958), which requires a viscous closure for the heat-fluxes. However, in many astrophysical environments – such as the solar wind, the intracluster medium and accretion flows – the mean free path of particles exceeds the scale lengths of interest, requiring collisionless models (Hammett & Perkins Reference Hammett and Perkins1990; Snyder et al. Reference Snyder, Hammett and Dorland1997; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009).

In the context of accretion flows around SMBHs, several studies have used collisionless Landau-fluid models to study the impact of pressure anisotropy in a local shearing-box simulation (Sharma et al. Reference Sharma, Hammett and Quataert2003, Reference Sharma, Hammett, Quataert and Stone2006). The model developed in this paper will help generalise these studies to black hole accretion flows to understand the impact of pressure anisotropy on relativistic jets and near-horizon flow dynamics. Beyond basic transport, the model studied here should allow the investigation of plasma dynamics absent in standard ideal GRMHD, such as the firehose and mirror instabilities (Kunz et al., Reference Kunz, Schekochihin and Stone2014) and magnetic immutability (Squire et al. Reference Squire, Schekochihin, Quataert and Kunz2019, Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023; Majeski, Kunz & Squire Reference Majeski, Kunz and Squire2024).

4.5.1. Comparison to EMHD

In this section, we briefly compare the Landau fluid model presented here with the EMHD models developed by Chandra et al. (Reference Chandra, Gammie, Foucart and Quataert2015), Most & Noronha (Reference Most and Noronha2021) and Most et al. (Reference Most, Noronha and Philippov2022). Both approaches yield an identical structural form for the stress-energy tensor, where non-equilibrium contributions are driven by pressure anisotropy and magnetic-field-aligned heat fluxes (2.27). However, the underlying closure physics is fundamentally different. The original model of Chandra et al. (Reference Chandra, Gammie, Foucart and Quataert2015) was a phenomenological model for a plasma in which micro-instabilities have already relaxed plasma to a local equilibrium state. Our model, however, does not make this assumption and instead tries to capture general pressure anisotropies.

More specifically, the form of the heat flux (4.3) and the pressure evolution equations (2.40c ) and (2.40d ) differ significantly between the two models. In the EMHD approach, the system is collision-dominated and the distribution function is assumed to remain close to a local Maxwell-Juttner equilibrium. This assumption restricts the magnitude of the pressure anisotropy and heat flux, forcing these out-of-equilibrium corrections to relax towards a Braginskii-like form [see (45) and (50) of Chandra et al. (Reference Chandra, Gammie, Foucart and Quataert2015)]:

(4.14a) \begin{align} &\mathcal{Q}_{\mathrm{EMHD}} \sim - \rho \chi \hat {b}^{\mu } \left [ \boldsymbol{\nabla} _{\mu } \varTheta + \varTheta a_{\mu } \right ] ,\\[-10pt]\nonumber \end{align}

(4.14b) \begin{align} &(p_{\perp } - p_{\parallel })_{\mathrm{EMHD}} \sim 3 \rho \nu \left (\hat {b}^{\mu } \hat {b}^{\nu } \boldsymbol{\nabla} _{\mu } U_{\nu } - \frac {1}{3} \boldsymbol{\nabla} _{\mu } U^{\mu } \right ) , \end{align}

where $\chi$ , $\nu$ and $\varTheta$ are the thermal transport coefficient, the viscous transport coefficients and the temperature. Unlike the EMHD model, the initial distribution function in the Landau fluid is assumed to be anisotropic and non-Maxwellian. This allows both parallel and perpendicular pressure variables to evolve dynamically. Moreover, the sources of dissipation in the two models differ. The EMHD models dissipation through local collisions, whereas in the Landau fluid model, dissipation arises from kinetic Landau damping and phase-mixing.

4.5.2. Comparison to recent weakly collisional models

During the final stages of preparing this manuscript, two recent studies by Wierzchucka et al. (Reference Wierzchucka, Bilbao, Thomas, Uzdensky and Schekochihin2026) and Ley et al. (Reference Ley, Tran and Zweibel2026) were submitted. These studies generalise the CGL equations to relativistic distribution functions with the assumptions of spatial homogeneity and non-relativistic bulk flow. Here, we briefly describe how to obtain the results from these papers using our more general approach. Consider a collisionless, homogeneous plasma with non-zero parallel electric, non-relativistic bulk flow ( $W_U \sim 1$ ) in Minkowski space and assume that $\partial _t f \gg \partial _{x,y,z} f$ . Using these assumptions, (2.46) reduces to (2.2) of Ley et al. (Reference Ley, Tran and Zweibel2026),

(4.15) \begin{align} \partial _t f_s + \frac {v_{\perp } \partial _t b }{2 b} \frac {\partial f_{s}}{\partial v_{\perp }} + v_{\parallel } \frac {\partial f_{s}}{\partial v_{\parallel }} \left (\frac {\partial _t n_s}{n_s} - \frac {\partial _t b}{b} \right ) = 0. \end{align}

Equation (4.15) can be solved using the method of lines to obtain (2.5), (2.11) and (2.12) of Ley et al. (Reference Ley, Tran and Zweibel2026). Once the evolution of the distribution function is obtained, Ley et al. (Reference Ley, Tran and Zweibel2026) derive the evolution of the parallel and perpendicular pressures for different initially isotropic distribution functions such as the Maxwell–Boltzmann distribution and the Maxwell–Juttner distribution. The results of Wierzchucka et al. (Reference Wierzchucka, Bilbao, Thomas, Uzdensky and Schekochihin2026) are similar but they use general symmetry principles in phase space to obtain these results.

In contrast, our model preserves full generality by deriving the evolution equations directly from the covariant Vlasov–Maxwell system without assuming homogeneity or non-relativistic flow. Our closure relationship relies on an anisotropic distribution function and we employ relativistic linear response theory to self-consistently incorporate heat fluxes driven by Landau damping. This allows our framework to describe the dissipative kinetic effects that are not analysed by Ley et al. (Reference Ley, Tran and Zweibel2026) and Wierzchucka et al. (Reference Wierzchucka, Bilbao, Thomas, Uzdensky and Schekochihin2026).