# Polarized anisotropic synchrotron emission and absorption and its
application to Black Hole Imaging

Alisa Galishnikova  [  Department of Astrophysical Sciences, Princeton
University, 4 Ivy Lane, Princeton, NJ 08544, USA  Alexander Philippov
Department of Physics, University of Maryland, College Park, MD 20742, USA
Eliot Quataert  Department of Astrophysical Sciences, Princeton University, 4
Ivy Lane, Princeton, NJ 08544, USA

###### Abstract

Low-collisionality plasma in a magnetic field generically develops anisotropy
in its distribution function with respect to the magnetic field direction.
Motivated by the application to radiation from accretion flows and jets, we
explore the effect of temperature anisotropy on synchrotron emission. We
derive analytically and provide numerical fits for the polarized synchrotron
emission and absorption coefficients for a relativistic bi-Maxwellian plasma
(we do not consider Faraday conversion/rotation). Temperature anisotropy can
significantly change how the synchrotron emission and absorption coefficients
depend on observing angle with respect to the magnetic field. The emitted
linear polarization fraction does not depend strongly on anisotropy, while the
emitted circular polarization does. We apply our results to black hole imaging
of Sgr A* and M87* by ray-tracing a GRMHD simulation and assuming that the
plasma temperature anisotropy is set by the thresholds of kinetic-scale
anisotropy-driven instabilities. We find that the azimuthal asymmetry of the
230 GHz images can change by up to a factor of $3$, accentuating
($T_{\perp}>T_{\parallel}$) or counteracting ($T_{\perp}<T_{\parallel}$) the
image asymmetry produced by Doppler beaming. This can change the physical
inferences from observations relative to models with an isotropic distribution
function, e.g., by allowing for larger inclination between the line of sight
and spin direction in Sgr A*. The observed image diameter and the size of the
black hole shadow can also vary significantly due to plasma temperature
anisotropy. We describe how the anisotropy of the plasma can affect future
multi-frequency and photon ring observations. In Appendices we calculate
kinetic anisotropy-driven instabilities (mirror, whistler, and firehose) for
relativistically hot plasmas.

<EMAIL_ADDRESS>
##  1 Introduction

Synchrotron emission produced by relativistic electrons in the presence of a
magnetic field appears in many astrophysical systems. It is the source of
emission across much of the electromagnetic spectrum in pulsar wind nebulae
and jets from neutron stars and black holes (BHs). Synchrotron emission is
also the source of the mm-wavelength radio emission observed on event-horizon
scales in M87* and Sgr A* by the Event Horizon Telescope (EHT) (Event Horizon
Telescope Collaboration et al., 2019a, 2022a).

Models of synchrotron emission from astrophysical plasmas typically assume
that the plasma has a thermal or power-law distribution function or a hybrid
of the two, such as a kappa distribution function. The latter two are
motivated by the power-law (non-thermal) synchrotron spectra often observed
from astrophysical sources. Another explicit assumption typically made is that
the electron distribution function is isotropic relative to the local magnetic
field, i.e., that the electrons have the same temperature or energy density in
all directions.111An exception to this is in very strongly magnetized plasmas
such as neutron star magnetospheres where the synchrotron cooling time is so
short that the perpendicular energy is nearly instantaneously radiated away.
In this paper we are focused on applications with weaker magnetic fields, such
as black hole accretion flows and jets.

In the presence of dynamically strong magnetic fields, the assumption of an
isotropic electron distribution function is not theoretically or
observationally well-motivated. By dynamically strong here, we mean an energy
density in the magnetic field similar to or larger than that in the plasma.
Such magnetized collisionless (and weakly collisional) plasmas can readily
depart from thermal equilibrium and develop anisotropies with respect to the
local magnetic field direction (Quataert et al., 2002). Although the
distribution function will in general be gyrotropic (isotropic in the plane
perpendicular to the magnetic field), it can have significant anisotropies
parallel and perpendicular to the local magnetic field (Kulsrud, 1983).

There is extensive observational evidence for such anisotropy in the solar
corona and solar wind (Bale et al., 2009). In the most extreme cases, Oxygen
ions in the solar corona have perpendicular temperatures that are a factor of
$\sim 10-100$ times that of their parallel temperature (Cranmer et al., 1999).
This anisotropy is in fact critical to interpreting spectroscopy of the solar
corona. By analogy, one might expect that anisotropy in the electron
distribution function could be important for interpreting synchrotron
radiation from astrophysical plasmas. This is particularly true in high
spatial resolution observations where our viewing angle relative to the local
magnetic field likely changes significantly across the image (e.g., the EHT or
radio interferometry more generally).

The anisotropy in a plasma’s distribution function cannot, however, grow
without bound. It is limited by kinetic-scale instabilities such as the
mirror, whistler, firehose, and ion cyclotron instabilities (Rosenbluth, 1956;
Southwood & Kivelson, 1993; Chandrasekhar et al., 1958; Rudakov & Sagdeev,
1961; Sudan, 1963; Gary, 1992). When the anisotropy in the distribution
function becomes too large (relative to the threshold of the
instability222Some instabilities, e.g., the ion cyclotron instability,
formally do not have a threshold, but their growth rate becomes sufficiently
small at low anisotropies that in practice they do.), such instabilities
rapidly grow, driving the anisotropy towards the instability threshold. This
endows the plasma with an effective collisionality that acts to partially
isotropize the distribution function. A very rough rule of thumb is that
instabilities set in vigorously when the fractional temperature anisotropy
satisfies $|\Delta T/T|\gtrsim\mathcal{O}(\beta^{-1})$ (where $\Delta T$ is
the temperature anisotropy and $\beta$ is the ratio of thermal to magnetic
energy). Anisotropy can thus be much larger in strongly magnetized plasmas
with $\beta\lesssim 1$. Anisotropy in the distribution function is thus
expected to be particularly important in jets and in models of accretion flows
with dynamically strong magnetic fields, such as the Magnetically Arrested
Disc (MAD) models favored by EHT observations of M87* (Event Horizon Telescope
Collaboration et al., 2021).

Observations of protons and electrons in the solar wind show that they obey
the expected anisotropy-driven instability thresholds and that the anisotropy
is larger at lower $\beta$ (Bale et al., 2009) (however, the measured
anisotropy is smaller than the instability thresholds at $\beta\lesssim 0.1$).
We expect that in accretion flows and jets, inflow, outflow, and heating of
the plasma will likewise drive temperature anisotropies to the point that
instabilities set in (Foucart et al., 2017). Global axisymmetric GR kinetic
simulations of collisionless plasma accreting onto a BH indeed find the growth
of the mirror and firehose instability and that they regulate the plasmas’s
temperature anisotropy (Galishnikova et al., 2023).

Motivated by the potential importance of an anisotropic distribution function
in synchrotron emitting plasmas, in this paper we theoretically calculate
emission and absorption of polarized synchrotron radiation for a physically
motivated gyrotropic distribution function. The study of polarized synchrotron
radiation dates back to the work of Westfold (1959), who studied emission from
an ultra-relativistically gyrating electron. General formulae for Stokes
parameters for ultra-relativistic synchrotron emission from an assemble of
electrons can be found in the review of Ginzburg & Syrovatskii (1965), who
noted that a substantial amount of circular polarization is present only in
the case of a highly anisotropic pitch-angle distribution. Melrose (1971)
presented the general equations for Stokes parameters for an arbitrary
anisotropic distribution function separable in momentum and pitch-angle, while
Sazonov (1972) focused on the case of a power-law momentum distribution with a
separable pitch-angle anisotropy.

In the last few decades, the study of synchrotron radiation was extended to a
broader range of validity and a number of different distribution functions via
numerical integration methods (Mahadevan et al., 1996; Shcherbakov, 2008;
Leung et al., 2011; Pandya et al., 2016, 2018; Dexter, 2016). This is useful
for improving analytical results at arbitrary frequency, emission direction
with respect to the magnetic field, and distribution function. These works
provide fits for the Stokes emissivities, absorptivities, and rotativities
that have been widely used in modeling polarized synchrotron radiation from
accreting black holes, particularly in the context of the EHT sources M87* and
Sgr A* [e.g., Dexter (2016); Mościbrodzka & Gammie (2018); White (2022) and
others; see also Gold et al. (2020)]. However, no pitch-angle anisotropy was
considered in these studies.

In this paper, we extend previous work on synchrotron radiation by studying
the intrinsic emission from an ensemble of electrons with an anisotropic
relativistic distribution function. We focus on the case of a relativistic
generalization of a bi-Maxwellian that has different temperatures
perpendicular and parallel to the local magnetic field (§2) and provide fits
for the polarized emissivity and absorption coefficients in Section 2.2. We
defer the case of an anisotropic power-law distribution function to future
work. We also defer the calculation of Faraday rotation and conversion
coefficients for an anisotropic distribution function to future work. We then
implement these expressions in a GR radiative transfer code to ray trace GRMHD
MAD simulations and study the impact of pitch-angle anisotropy on the
observable quantities (Section 3). Finally, in Section 4 we summarise the
application of our results to current and future EHT observations.

##  2 Synchrotron Emission From Gyrotropic Distribution Functions

In this section, we describe radiation transfer and emission produced by
electrons with a gyrotropic distribution function $f(\gamma,\xi)$ in the
presence of a background magnetic field ${\bf B}$; $\gamma$ and $\xi$ denote
the Lorentz factor of electrons and pitch angle with the magnetic field
respectively; we will use $\mu=\cos\xi$ and $\xi$ interchangeably in what
follows. Throughout the paper, $m_{e}$, $e$, and $c$ are constants that stand
for electron mass, electron charge, and speed of light. Therefore, the
momentum of a particle with velocity ${\bf v}$ is ${\bf p}=m_{e}\gamma{\bf v}$
and $\beta={\bf v}/c$. In what follows, we normalize the frequency of emission
$\nu$ by a non-relativistic cyclotron frequency given by $\nu_{c}=eB/2\pi
m_{e}c$. The angle between the propagation direction along the wavevector
${\bf k}$ and the background magnetic field ${\bf B}$ is set by $\theta_{B}$.

Polarized emission is described in the Stokes basis as
$I_{a}=\\{I,Q,U,V\\}^{T}$, where $I$ stands for intensity, $Q$ and $U$
describe linear polarization, and $V$ describes circular polarization. Given
emissivities $j_{a}=\\{j_{I},j_{Q},j_{U},j_{V}\\}^{T}$, absorption
coefficients
$\alpha_{a}=\\{\alpha_{I},\alpha_{Q},\alpha_{U},\alpha_{V}\\}^{T}$, and
Faraday rotativities $\rho_{a}=\\{\rho_{Q},\rho_{U},\rho_{V}\\}^{T}$, the
polarized emission can then be found using [see, e.g, Leung et al. (2011)]

$\frac{dI_{a}}{ds}=j_{a}-M_{ab}I_{b},$ (1)

where $M_{ab}$ is the Mueller matrix,

$M_{ab}=\begin{pmatrix}\alpha_{I}&\alpha_{Q}&\alpha_{U}&\alpha_{V}\\\
\alpha_{Q}&\alpha_{I}&\rho_{V}&-\rho_{U}\\\
\alpha_{U}&-\rho_{V}&\alpha_{I}&\rho_{Q}\\\
\alpha_{V}&\rho_{U}&-\rho_{Q}&\alpha_{I}\end{pmatrix},$ (2)

where $U$ components vanish if ${\bf B}$ is aligned with $U$: $j_{U}=0$,
$\alpha_{U}=0$, and $\rho_{U}=0$. Then $I$ components of $j_{a}$ and
$\alpha_{a}$ describe total emission, $Q$ components describe linearly
polarised emission and $V$ describe circularly polarised emission, while
$\rho_{Q}$ and $\rho_{V}$ account for Faraday conversion and rotation
respectively. In this work, we focus on emissivities $j_{a}$ and absorption
coefficients $\alpha_{a}$, while Faraday rotativities $\rho_{a}$ will be
studied in future work.

We need to evaluate $j_{a}$, $\alpha_{a}$, and $\rho_{a}$ through
$f(\gamma,\xi)$ to describe the radiation emission and transfer. In the Stokes
basis at frequency $\nu$ [see, e.g., Leung et al. (2011)]:

$\begin{split}j_{a}&=\frac{2\pi e^{2}\nu^{2}}{c}\int
d^{3}pf(\gamma,\xi)\sum_{n=1}^{\infty}\delta(y_{n})K_{a}(z),\\\
\alpha_{a}&=\frac{2\pi\nu}{m_{e}c^{2}}\int
d^{3}pDf(\gamma,\xi)\sum_{n=1}^{\infty}\delta(y_{n})K_{a}(z),\end{split}$ (3)

where $\delta(y_{n})$ is a delta function of argument
$y_{n}=n\nu_{c}/\gamma-\nu(1-\beta\cos\xi\cos\theta_{B})$,
$z=\nu\gamma\beta\sin\theta_{B}\sin\xi/\nu_{c}$, $d^{3}p=2\pi
m_{e}^{3}c^{3}\gamma^{2}\beta d\gamma d\cos\xi$ for a gyrotropic
$f(\gamma,\xi)$, and $Df$ is an operator that includes a full derivative of
the distribution function:

$\begin{split}&Df\equiv\left(k_{\parallel}\frac{\partial}{\partial
p_{\parallel}}+\frac{\omega-
k_{\parallel}v_{\parallel}}{v_{\perp}}\frac{\partial}{\partial
p_{\perp}}\right)f(\gamma,\xi)\\\
&=\frac{2\pi\nu}{m_{e}c^{2}}\left(\frac{\partial}{\partial\gamma}+\frac{\beta\cos\theta_{B}-\cos\xi}{\beta^{2}\gamma}\frac{\partial}{\partial\cos\xi}\right)f(\gamma,\xi),\end{split}$
(4)

In Equation 3, $K_{a}$ is defined as

$K_{a}=\begin{cases}M^{2}J_{n}^{2}(z)+N^{2}J^{\prime 2}_{n}(z),a=I,\\\
M^{2}J_{n}^{2}(z)-N^{2}J^{\prime 2}_{n}(z),a=Q,\\\ 0,a=U,\\\
2MNJ_{n}(z)J^{\prime}_{n}(z),a=V,\end{cases}$ (5)

where $J_{n}$ is a Bessel function of the first kind,
$M=(\cos\theta_{B}-\beta\cos\xi)/\sin\theta_{B}$, and $N=\beta\sin\xi$. Given
$f(\gamma,\xi)$, one can find $j_{a}$ and $\alpha_{a}$ through Equations 3, 4,
and 5.

###  2.1 Anisotropic electron distribution function

We will use an anisotropic distribution function $f(\gamma,\xi)$ for emitting
electrons, written in cgs units:

$\begin{split}&f(p_{\perp},p_{\parallel})=\frac{n_{e}\eta^{1/2}}{4\pi
m_{e}^{3}c^{3}\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}\times\\\
&\exp{(-\sqrt{1+(p_{\perp}/m_{e}c)^{2}+\eta(p_{\parallel}/m_{e}c)^{2}}/\epsilon_{\perp})},\end{split}$
(6)

where $n_{e}$ is the electron number density,
$\epsilon_{\perp}=kT_{\perp,e}/m_{e}c^{2}$ is the dimensionless perpendicular
electron temperature, $K_{2}$ is the modified Bessel function of the second
kind, $p_{\perp}$ and $p_{\parallel}$ stand for the relativistic momentum
perpendicular and parallel to the magnetic field direction. Here $\eta$ is a
measure of anisotropy, with $\eta=1$ corresponding to an isotropic
relativistic Maxwellian distribution function. In the non-relativistic limit,
$T_{\perp,e}/T_{\parallel,e}=\eta$, while
$T_{\perp,e}/T_{\parallel,e}\approx\eta^{0.8}$ in the ultra-relativistic limit
(see Appendix C for a detailed fit). Transforming $f(p_{\perp},p_{\parallel})$
to $\gamma-\xi$ variables:

$\begin{split}&f(\gamma,\xi)=\frac{n_{e}\eta^{1/2}}{4\pi
m_{e}^{3}c^{3}\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}\times\\\
&\exp{(-\sqrt{1+(\gamma^{2}-1)(\sin^{2}\xi+\eta\cos^{2}\xi)}/\epsilon_{\perp})}.\end{split}$
(7)

In the limit of high $\gamma$:

$\begin{split}f(\gamma,\xi)&=\frac{n_{e}\eta^{1/2}}{4\pi
m_{e}^{3}c^{3}\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}\times\\\
&\exp{(-\gamma\sqrt{1+(\eta-1)\cos^{2}\xi}/\epsilon_{\perp})}\\\
&=\frac{n_{e}\eta^{1/2}}{4\pi
m^{3}c^{3}\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}\times\exp{(-\gamma/\epsilon^{\ast}_{\perp})},\end{split}$
(8)

where

$\epsilon^{\ast}_{\perp}=\epsilon^{\ast}_{\perp}(\xi)=\frac{\epsilon_{\perp}}{\sqrt{1+(\eta-1)\cos^{2}\xi}}$
(9)

is the new renormalized temperature. Thus, in the high-$\gamma$ limit, the
temperature in the distribution function depends on both the anisotropy $\eta$
and pitch angle or $\mu=\cos\xi$. In the isotropic case $\eta=1$, the
temperature is described by $\epsilon=\epsilon_{\perp}=\epsilon_{\parallel}$
in all directions. Note that in the analytical fitting functions in the next
subsection, $\epsilon^{\ast}_{\perp}$ will be evaluated at $\xi=\theta_{B}$
because the radiation is beamed along the local direction of motion of
relativistic electrons (as is standard in synchrotron radiation, see Appendix
A for details). In our numerical evaluations, however, we integrate and sum
over $\xi$ and $\theta_{B}$ separately using Equation 3.

The total derivative $Df$ that is used in calculating $\alpha_{a}$ in Equation
3 contains

$\begin{split}\partial_{\gamma}f(\gamma,\xi)&=-\frac{\gamma}{\epsilon_{\perp}}\frac{1+\mu^{2}(\eta-1)}{\sqrt{\gamma^{2}+(\gamma^{2}-1)(\eta-1)\mu^{2}}}f(\gamma,\xi),\\\
\partial_{\mu}f(\gamma,\xi)&=-\frac{(\gamma^{2}-1)(\eta-1)\mu}{\epsilon_{\perp}\sqrt{\gamma^{2}+(\gamma^{2}-1)(\eta-1)\mu^{2}}}f(\gamma,\xi).\\\
\end{split}$ (10)

While $\partial_{\mu}f(\gamma,\xi)$ is non-zero, we find that the absorption
coefficients change negligibly if we include this term. This is due to the
prefactor it goes with in equation 4 since $\gamma\gg 1$ and the absorption is
mainly concentrated around $\xi\approx\theta_{B}$ (see Appendix A for
details).

###  2.2 Emissivities and absorption coefficients

Figure 1: Emissivity $j_{I}$ (a), absorption $\alpha_{I}$ (b), and emitted
polarization fractions (c) as functions of the angle between propagation
direction and magnetic field $\theta_{B}$ at different anisotropy values:
$T_{\perp}<T_{\parallel}$ ($\eta<1$, blue), isotropic ($\eta=1$, black), and
$T_{\perp}>T_{\parallel}$ ($\eta>1$, red). The free parameters are
$\nu/\nu_{c}=10^{3}$ and $\epsilon_{\perp}=10$ (near the peak of the optically
thin synchrotron spectrum for an isotropic distribution function). In (c) the
emitted linear $|j_{Q}|/j_{I}$ and circular $|j_{V}|/j_{I}$ polarization
fractions are shown by solid and dashed lines respectively. In (a)
$\sin^{2}\theta_{B}$ dependence is shown by a black dotted line.

We obtain the following fits for emissivities and absorption coefficients for
a relativistic plasma with an anisotropic bi-Maxwellian distribution function
(see Appendix A for details on the derivation):

$j_{a}=\begin{cases}\eta^{1/2}\mathcal{K}_{a}(\epsilon^{\ast}_{\perp},\epsilon_{\perp})j_{a,{\rm
iso}}(\epsilon=\epsilon_{\perp}^{\ast},\nu/\nu_{c},\theta_{B}),a=\\{I,Q\\}\\\
\eta^{3/2}\mathcal{K}_{a}(\epsilon^{\ast}_{\perp},\epsilon_{\perp})j_{a,{\rm
iso}}(\epsilon=\epsilon_{\perp}^{\ast},\nu/\nu_{c},\theta_{B}),a=V,\end{cases}$
(11)

where

$\mathcal{K}_{a}(\epsilon^{\ast}_{\perp}/\epsilon_{\perp})=\begin{cases}(\epsilon^{\ast}_{\perp}/\epsilon_{\perp})[K_{2}(1/\epsilon^{\ast}_{\perp})/K_{2}(1/\epsilon_{\perp})],a=\\{I,Q\\}\\\
(\epsilon^{\ast}_{\perp}/\epsilon_{\perp})^{3}[K_{2}(1/\epsilon^{\ast}_{\perp})/K_{2}(1/\epsilon_{\perp})],a=V,\end{cases}$
(12)

and $K_{2}(1/\epsilon_{\perp})\approx 2\epsilon_{\perp}^{2}$ when
$\epsilon_{\perp}\gg 1$. Here $\epsilon^{\ast}_{\perp}$ is evaluated at
$\xi=\theta_{B}$ and $j_{a,{\rm
iso}}(\epsilon=\epsilon_{\perp}^{\ast},\nu/\nu_{c},\theta_{B})$ and
$\alpha_{a,{\rm
iso}}(\epsilon=\epsilon_{\perp}^{\ast},\nu/\nu_{c},\theta_{B})$ correspond to
emission and absorption in the case of an isotropic relativistic Maxwellian.
Absorption coefficients $\alpha_{a}$ can be obtained via Kirchoff’s law for a
thermal distribution function:

$\displaystyle j_{a,{\rm iso}}-\alpha_{a,{\rm iso}}B_{\nu}=0,$ (13)

where
$B_{\nu}(T_{\perp,e})=\frac{2h\nu^{3}}{c^{2}}(\exp{(h\nu/kT_{\perp,e})-1})^{-1}$
for an anisotropic distribution has the same functional form as in the
isotropic case (with $T_{\perp,e}=T_{e}$ for an isotropic distribution).
Equations 11 correspond to total intensity and linearly polarized emissivities
with the same functional form as in an isotropic plasma but with a temperature
$\epsilon_{\perp}^{\star}$ that depends on the observer-angle $\theta_{B}$ due
to the anisotropy in the distribution function relative to the magnetic field
(the factors of
$\eta^{1/2}\mathcal{K}_{a}(\epsilon^{\ast}_{\perp},\epsilon_{\perp})\approx\eta^{1/2}(\epsilon^{\star}_{\perp}/\epsilon_{\perp})^{3}$
in Equation 11 reflect the change in normalization of the distribution
function due to the different number of particles whose radiation is beamed in
the direction of the observer). By contrast, the Stokes V (circular
polarization) emissivity in Equation 11 differs by a larger factor because of
a change in the efficiency of producing circularly polarized radiation for an
anisotropic distribution function (see Eq. A1 and A9 in Appendix A).

Given Equation 11, it is straightforward to derive the location of the peak of
the optically thin emission $\nu j_{I}(\epsilon,\nu,\theta_{B},\eta)$ as

$\nu_{\rm peak}\approx 36.7\nu_{c}\epsilon_{\perp}^{\ast
2}\sin\theta_{B}\simeq\frac{36.7\nu_{c}\epsilon_{\perp}^{2}\sin\theta_{B}}{1+(\eta-1)\cos^{2}\theta_{B}},$
(14)

which shifts to lower (higher) frequencies with increasing (decreasing) $\eta$
at a fixed $\epsilon_{\perp}$ and $\theta_{B}$ (though we show below in Figure
1 that $\eta$ changes the efficiency of producing radiation as a function of
$\theta_{B}$).

We find good agreement between our analytic expressions and numerical
calculation over a wide parameter range using a publicly available code
symphony, which integrates Equations 3 for a given distribution function
$f(\gamma,\xi)$. For a detailed derivation of the fits given by Equations 11,
full expressions for $j_{a,{\rm iso}}$ and $\alpha_{a,{\rm iso}}$, and the
comparison with numerical solutions see Appendix A. The fits presented here
become inaccurate for $\epsilon^{\star}_{\perp}\lesssim 3$ and low frequencies
$\nu/\nu_{c}\lesssim 10$, where the isotropic fits that we scale to in
Equation 11 themselves become inaccurate.

We demonstrate the resulting emission properties in Figure 1, where $j_{I}$
(a), $\alpha_{I}$ (b), and the emitted linear and circular polarization
fractions (solid and dashed lines in (c) respectively) are shown as functions
of the angle between propagation direction and magnetic field $\theta_{B}$ at
different values of anisotropy parameter $\eta$, represented by different
colors (we intentionally choose relatively large anisotropy to highlight the
large differences in synchrotron radiation possible in this limit). The
parameters used in this Figure are a high frequency of $\nu/\nu_{c}=10^{3}$
and temperature of $\epsilon_{\perp}=10$. The isotropic case (black line)
closely follows a $\sin^{2}\theta_{B}$ dependence of $j_{I}$ (dotted) in panel
(a) due to the frequency being near the peak of the synchrotron emissivity.

Figure 1 shows that there are significant differences in the synchrotron
emission/absorption for an anisotropic plasma distribution, compared to the
isotropic case. This change can be understood as a renormalization of the
number of relativistic particles emitting toward the observer at $\theta_{B}$.
In particular, the plasma is less prone to emitting along the magnetic field
at $\eta>1$ ($T_{\perp}>T_{\parallel}$, red lines), hence the rapid fall off
of $j_{\nu,I}$ with decreasing $\theta_{B}$, compared to $\eta\equiv 1$. That
is, the emission is even more concentrated towards $\theta_{B}=90^{\circ}$
when $\eta>1$. For the opposite anisotropy, $\eta<1$ and
$T_{\parallel}>T_{\perp}$ (blue lines), the number of particles capable of
emitting along the magnetic field direction increases. Thus, more emission can
be produced at smaller $\theta_{B}$ (along the magnetic field), relative to
the isotropic case with $\eta\equiv 1$. The unpolarised absorption coefficient
$\alpha_{I}$ (b) shows a similar but smaller dependence on $\eta$ as $j_{I}$.

The polarization fractions have a weaker dependence on $\eta$. This is shown
in Figure 1 (c) with solid and dashed lines for the intrinsic linear
$|j_{Q}|/j_{I}$ and circular $|j_{V}|/j_{I}$ polarization fractions,
respectively.

Quantitatively, both $|j_{Q}|/j_{I}$ and $|j_{V}|/j_{I}$ are higher for higher
$\eta$ but the change is particularly modest for the intrinsic linear
polarization $|j_{Q}|/j_{I}$. Since most of the emission comes from small
(large) angles for $\eta<1$ ($\eta>1$), the emitted circular polarisation
degree can significantly vary with $\eta$ due to the change in which pitch
angles dominate the emission. In particular, $\eta>1$ is significantly more
circularly polarised, and $\eta<1$ is less circularly polarized, compared to
emission from an isotropic plasma. This is because $\eta>1$ decreases the
effective temperature $\epsilon_{\perp}^{\ast}$ by suppressing the parallel
temperature at a fixed $\epsilon_{\perp}$.

Figure 2: Plasma-$\beta_{\rm th}$ (column 1, a and g), plasma temperature
$P/2\rho$ (column 2, b and h), anisotropy $\eta$ (columns 3 and 4 for mirror
and firehose respectively), and normalized electron perpendicular temperature
$\epsilon_{\perp}$ (columns 5 and 6 for mirror and firehose respectively) for
$a=0.98$ (at time of $14300r_{g}/c$, top row) and $a=0.5$ (at time of
$16000r_{g}/c$, bottom row).

##  3 Black Hole Imaging

In this section we study the observational implications of synchrotron
emission by a plasma with anisotropic temperatures in the context of black
hole accretion flows. Specifically, we focus on the application to the EHT
targets Sgr A* and M87* (Event Horizon Telescope Collaboration et al., 2019a,
2022a). Our goal in this initial study is to determine the rough magnitude of
the effect and which observables are most sensitive to electron temperature
anisotropy. The exact electron temperature anisotropy in the near-horizon
plasma is uncertain so we will use general stability arguments to bound the
anisotropy and thus the effect of anisotropy on the synchrotron radiation.

###  3.1 Method

We use a publicly available radiative transfer code blacklight (White, 2022)
to ray trace synchrotron emission in GRMHD simulations and study the resulting
intensity and polarization images. We implement the formulas for the
emissivity and absorption coefficients of hot electrons with an anisotropic
distribution function discussed in §2 (we use the limit of high temperature
such that $K_{2}(1/x)\approx 2x^{2}$). Since EHT observational constraints
favor highly magnetized models (Event Horizon Telescope Collaboration et al.,
2021), we restrict our study to a MAD simulation of plasma accreting onto a
spinning BH with dimensionless spin parameters of $a=0.98$ and $0.5$. Our
results are averaged over $100$ snapshots which span a time of $1000r_{g}/c$
when the accretion rate and magnetic flux on the horizon are in approximate
steady-state (see Appendix D for details on the simulation setup and choice of
this time period). Since the MHD method cannot handle vacuum, our GRMHD
simulations have a ceiling plasma magnetization parameter of
$\sigma=B^{2}/[4\pi\rho c]=100$, and we ignore emission from $\sigma>10$
regions.333The magnetization in the jet region can, in reality, be
significantly larger than the ceiling value set in our GRMHD simulations.
These low-density regions with $\sigma\gtrsim 10$ are, however, not expected
to contribute significantly to the observed flux at $230$ GHz. We choose a BH
mass and distance to the BH to match M87*, $M_{BH}=6.5\times 10^{9}M_{\odot}$
and $d=1.67\times 10^{7}$pc, unless otherwise specified. In the GRMHD
simulations, the plasma number density normalization is a free parameter,
which we choose such that the total flux of the image $F_{\nu}$ at $230$ GHz
matches EHT observations of M87*, i.e., $0.66$ Jy. The raytraced images have a
resolution of $128\times 128$ cells, with a point camera located at $100r_{g}$
and inclination [observing angle] of $\theta$. We consider both
$\theta=163\deg$, appropriate for M87*, as well as less face-on viewing angles
to demonstrate the change with viewing angle.

Since the GRMHD equations evolve a single fluid, while in the plasmas of
interest the electrons and ions likely have different temperatures, we have
the freedom to set the electron temperature. The heating of collisionless
electrons should depend on local plasma parameters, in particular, the
magnetic field strength (Quataert & Gruzinov, 1999) via $\beta_{\rm th}=P_{\rm
th}/P_{B}$ – the ratio of thermal pressure to magnetic pressure. To
parameterize the electron temperate, we use the widely-employed $R_{\rm
high}-R_{\rm low}$ model (Mościbrodzka et al., 2016). In this model, the ion-
to-electron temperature is set by

$R=\frac{T_{i}}{T_{e}}=\frac{\beta_{\rm th}^{2}R_{\rm high}+R_{\rm
low}}{1+\beta_{\rm th}^{2}},$ (15)

where $\beta_{\rm th}=P_{{\rm th}}/P_{B}$ is the plasma $\beta_{\rm th}$ for
an MHD fluid, $R_{\rm high}$ and $R_{\rm low}$ are ion-to-electron temperature
ratios in the high and low-$\beta_{\rm th}$ regions respectively. The fluid
GRMHD temperature is $T=(T_{i}+T_{e})/2$, and thus $T_{e}=2T/(1+R)$. In this
work, we explore three cases: $R_{\rm high}=1$, $10$ and $100$, while $R_{\rm
low}$ is set to $1$ always.

To study the effect of the anisotropy of the plasma on images, we also have
the freedom to set the anisotropy parameter $\eta\sim
T_{\perp,e}/T_{\parallel,e}$ since the GRMHD simulations have no information
about plasma anisotropy. The anisotropy of the plasma is limited by kinetic-
scale instability thresholds, which allow for a large anisotropy in
low-$\beta_{\rm th}$ regions. Ion-scale mirror and firehose instabilities are
clearly present in global GR kinetic simulations of collisionless plasma
accreting onto a BH (Galishnikova et al., 2023) and in kinetic shearing box
simulations (Kunz et al., 2014; Riquelme et al., 2015). The electrons also
contribute to driving mirror, firehose, and whistler instabilities, which are
important for setting the electron temperature anisotropy. Since the magnitude
of the electron temperature anisotropy in the near-horizon environment is not
fully understood, we consider all three limiting cases – where the plasma sits
at the mirror ($\eta>1$), whistler ($\eta>1$), or firehose ($\eta<1$)
instability thresholds everywhere. We then compare these limiting cases to the
usually considered isotropic plasma distribution. This should bracket the
magnitude of the effect introduced by an anisotropic electron distribution
function. We note that in the single-fluid global “extended GRMHD” simulations
of Foucart et al. (2017) in which the pressure anisotropy is a dynamical
variable, most of the plasma was near the mirror threshold. If generically
true, and applicable to electrons, this would suggest that the mirror and
whistler instability thresholds are the most important.

The microinstability thresholds can be expressed as
$T_{\perp,e}/T_{\parallel,e}=g(\beta_{e})$, where $g(\beta_{e})$ is a function
of election plasma$-\beta$, different for each of the anisotropy-driven
instability [see Appendix B for a derivation of relativistic mirror, firehose,
and whistler instabilities and Appendix C for additional details on
$g(\beta_{e})$ for each instability]. Therefore, our procedure for obtaining
$T_{\perp,e}$ and $\eta$ for each of the four instability cases (here and
after mirror, whistler, isotropic, firehose) is as follows. We first compute
electron-$\beta$ as $\beta_{e}=2\beta_{\rm th}/(1+R)$, according to Eq. 15.
Knowing $\beta_{e}$ allows us to calculate
$T_{\perp,e}/T_{\parallel,e}=g(\beta_{e})$ for each case of interest and thus
$\eta$ (we consider the relativistic limit, where
$T_{\perp,e}/T_{\parallel,e}\approx\eta^{0.8}$). The perpendicular temperature
can then be separately determined from the definition
$T_{e}=(T_{\parallel,e}+2T_{\perp,e})/3=T_{\perp,e}(\eta^{1/0.8}+2)/3=T_{i}/R$.
Now that we have $\eta$ and $\epsilon_{\perp}$ for electrons, we can then
calculate $j_{a}$ and $\alpha_{a}$, given by Equations 11 (with
$\epsilon_{\perp}=kT_{\perp,e}/m_{e}c^{2}$).

In Figure 2 we show an example of the inferred physical conditions in MAD
accretion flows from GRMHD simulations: MHD-$\beta_{\rm th}$ (first column),
plasma temperature $P/2\rho$ (second column), and the resulting $\eta$ (third
and fourth columns for mirror and firehose respectively) and electron
$\epsilon_{\perp}$ (fifth and sixth columns for mirror and firehose
respectively). The top row is for $a=0.98$ while the bottom row is for
$a=0.5$. Grey regions indicate $\sigma\geq\sigma_{\rm cut}=10$. Since MAD
simulations are highly magnetized with low plasma-$\beta$ in much of the
volume, the mirror (c and i, $\eta>1$) and firehose (d and j, $\eta<1$)
instabilities allow for a large temperature anisotropy in much of the volume.
The mirror case also results in a higher electron temperature
$\epsilon_{\perp}$, while the firehose case results in a lower and more
uniform $\epsilon_{\perp}$.

Future observations aim to probe not only the direct emission from the BH but
also the lensed emission associated with the “photon ring” (Johnson et al.,
2023). The latter can be decomposed into a series of sub-rings labeled by the
ray order $n$ – the number of half-orbits a photon traveled to the observer,
defined as $(\Delta\phi_{\rm ray}\mod\pi)$, where $\Delta\phi_{\rm ray}$ is
the change in the angular coordinate $\phi_{\rm ray}$ along the ray in the
plane of its orbit. To distinguish $n=0$ (direct) and $n=1$ (“photon ring” of
order $1$) in the ray tracing, we track $dz/d\lambda$ along each ray, where
$z$ and $\lambda$ are Cartesian Kerr Schild coordinate along the spin axis and
coordinate along the ray respectively. The number of times that $dz/d\lambda$
crosses zero for a particular ray defines the order of this ray $n$, allowing
us to approximately distinguish $n=0$ and $n=1$.

###  3.2 230 GHz images

Figure 3: Synchrotron emission of accreting plasma raytraced from a MAD
simulation for a BH with $a=0.98$ at $R_{\rm high}=10$ and inclination of
$\theta=163^{\circ}$. Each column represents different plasma anisotropy:
mirror instability threshold (column 1), whistler instability threshold
(column 2), isotropic plasma distribution (column 3), and firehose instability
threshold (column 4). The first row represents the full image blurred with
$20\mu$as FWHM Gaussian kernel on a linear scale (a-d), the second row shows a
full unblurred image on a logarithmic scale, from which $I_{\nu,0}$ ($n=0$)
and $I_{\nu,1}$ ($n=1$) are decoupled on the third (i-l) and forth (m-p) rows
respectively. Figure 4: Synchrotron emission of accreting plasma raytraced
from a MAD simulation at an inclination of $\theta=135^{\circ}$ for a BH with
$a=0.98$. As in Fig. 3, each column represents different plasma anisotropy.
The first and second rows represent $R_{\rm high}=1$ and $100$ models
respectively. Figure 5: Similar to Fig. 3 but for a BH with spin parameter of
$a=0.5$. Inclination of $\theta=163^{\circ}$ and $R_{\rm high}=10$ are
identical to Figure 3. The first row represents the full image blurred with
$20\mu$as FWHM Gaussian kernel on a linear scale (a-d), and the second row
shows a full unblurred image on a logarithmic scale.

Total intensity images observed at $\theta=163^{\circ}$, expected for M87*
(Mertens et al., 2016), with $R_{\rm high}=10$ are shown in Figure 3 for
$a=0.98$. The top row (a-d) shows the brightness blurred with $20\mu as$ FWHM
Gaussian kernel on a linear scale to match current EHT observations. Each
column represents a different anisotropy model: mirror, whistler, isotropic,
and firehose (from left to right, from largest to smallest $\eta$). The
density normalization is roughly the same for each of these cases at fixed
$R_{\rm high}$ and observing angle $\theta$, with the density in the firehose
model being larger than in the isotropic case by a factor of a few. The three
bottom rows in Figure 3 show unblurred full emission (second row), which is
decomposed into the direct emission ($n=0$, third row) and the $n=1$ photon
ring (fourth row) on a logarithmic scale (as appropriate for future higher
dynamic range measurements).

The azimuthal anisotropy in the images in Figure 3 is due to a combination of
two effects: Doppler beaming and differences in the angle $\theta_{B}$
relative to the local magnetic field that photons are emitted at, in order to
arrive at a given location in the observed image. Anisotropy in the electron
distribution function can significantly change the synchrotron emission as a
function of $\theta_{B}$, thus changing this second source of azimuthal image
anisotropy. Figure 3 shows that, compared to the isotropic case (c), the
mirror and whistler images [$\eta>1$, (a) and (b)] are more azimuthally
asymmetric, while plasma at the firehose instability threshold [$\eta<1$, (d)]
results in a more symmetric image. This is also noticeable in the unblurred
case, as well as separately in $n=0$ and $n=1$ images.

The dependence of the azimuthal image symmetry on plasma anisotropy is also
more apparent with increasing viewing angle, i.e. as we look more “edge-on”
instead of “face-on”. Additionally, the effect of the anisotropy of the
distribution function is more prominent for larger $R_{\rm high}$. This is
because larger $R_{\rm high}$ suppresses the emission from high-$\beta$
regions (where distribution function anisotropies are constrained to be
smaller) relative to low-$\beta$ regions (where distribution function
anisotropies can be larger). The more azimuthally asymmetric image at higher
inclination and higher $R_{\rm high}$ are demonstrated in Figure 4, where we
show the full intensity images at $R_{\rm high}=1$ (top) and $R_{\rm
high}=100$ (bottom) at a higher inclination relative to the spin axis of
$\theta=135^{\circ}$ for $a=0.98$.

Table 1 shows the image-averaged, emission-weighted ratio of the two
components of anisotropic temperature, $\langle
j_{\nu}T_{\perp}/T_{\parallel}\rangle/\langle j_{\nu}\rangle$, for $a=0.98$,
$\theta=163^{\circ}$, $R_{\rm high}=1$, $10$, and $100$. As $R_{\rm high}$
increases, the anisotropic temperature ratio approaches our maximum allowed
values of $10$ ad $0.1$ for mirror and firehose models respectively. The
significant changes in image morphology found here thus require large
temperature anisotropy in the emitting plasma.

$R_{\rm high}$ | mirror | whistler | firehose
---|---|---|---
1 | 3.13 | 1.47 | 0.52
10 | 5.3 | 2.0 | 0.1100
100 | 8.9 | 2.9 | 0.102
Table 1: Values of image-averaged emission-weighted temperature anisotropy
$T_{\perp}/T_{\parallel}$ [pixel-averaged $\langle
j_{\nu}T_{\perp}/T_{\parallel}\rangle/\langle j_{\nu}\rangle$] for our three
temperature anisotropy cases at a = 0.98, an inclination of $\theta=163\deg$,
and three $R_{\rm high}$ values.

To better understand the interplay between Doppler-induced asymmetry and
magnetic field viewing angle-induced asymmetry, we also consider the case of a
moderately spinning BH, $a=0.5$, shown in Figure 5, where the Doppler effect
is smaller than for $a=0.98$ studied above. This figure is organized
identically to Fig. 3 and the viewing angle relative to the spin axis and the
choice of $R_{\rm high}=10$ are the same. We find that the asymmetry of the
image due to the plasma temperature anisotropy is still pronounced, similar to
the case of a highly spinning BH. As in the $a=0.98$ case, mirror and whistler
anisotropies make the image more asymmetric, while temperature anisotropy near
the firehose boundary results in a more symmetric image.

Our calculations show that the anisotropic temperature distribution of plasma
sitting at the firehose and mirror thresholds leads to a more azimuthally
symmetric or asymmetric synchrotron image, respectively. At first glance, it
is not entirely obvious why the firehose sense of anisotropy (rather than the
mirror sense of anisotropy) should be associated with a more symmetric image.
Our interpretation of this is that if the rotation rate of the magnetic field
lines is small relative to the rotation rate of the plasma, then in ideal MHD
models, the plasma velocity is approximately parallel to the magnetic field
direction (see, e.g., eq. E148b of Chael et al. 2023b for a relativistic
version of this expression). For nearly (but not exactly) face-on viewing
angles, the Doppler effect and the effect of changing viewing angle relative
to the magnetic field are “in phase”: the brightening and dimming produced by
the two effects peak in roughly the same places in the image plane (this
follows, e.g., from the analytic model in Narayan et al. 2021). The firehose
instability sense of anisotropy counteracts this by making the emission a
significantly weaker function of angle relative to the magnetic field (Fig. 1)
thus making the overall emission more isotropic.

Another key difference between images with different electron temperature
anisotropy is the image diameter; this is noticeable at both spin values in
Figures 3 and 5: the size of the bright region in the image increases as
$\eta$ increases. Additionally, $a=0.5$ shows variations in the size of the BH
shadow between different models in Fig. 5. Both of these effects, as well as
the asymmetry of the images, are quantified below.

Figure 6: Average angle with the magnetic field along the ray, measured by
emission-weighted sine of $\theta_{B}$ for different spins and viewing angles.
Lower spin decreases the average angle of the emitted photons relative to the
magnetic field. This in turn enhances the effects of plasma anisotropy on the
observed image (Figures 3-5).

Figure 6 shows the emissivity-weighted angle between the magnetic field and
photon direction along the ray $\langle j_{\nu}\theta_{B}\rangle/j_{\nu}$.
This angle is larger for $a=0.98$ (a) than for $a=0.5$ (b,c) in the inner
region of the image. Physically for roughly face-on viewing angles, the
magnetic field in the accretion flow onto a BH with a smaller spin has a more
vertical field than onto a highly spinning BH (where the field is wrapped up
to be more azimuthal). This leads to the average angle between the propagation
direction and the local magnetic field decreasing for $a=0.5$ relative to
$a=0.98$. A less face-on viewing angle produces a similar effect (c). Figure 6
shows results for the isotropic emission model but can be used to gain insight
into why the central “shadow” is noticeably different in the firehose and
mirror cases in Figures 3 and 5. In particular, the lower average angle
between the photon and magnetic field in Figure 6 at lower spin and observer
viewing angle implies (via Figure 1) that in the mirror (firehose) case the
emission in the shadow should be suppressed (enhanced). This is exactly what
is seen in the images. Plasma anisotropy could thus have an effect on
observational efforts to infer physical properties of the black hole using the
“inner shadow” (Chael et al., 2021).

We now quantify the effects of changing image size and asymmetry for different
plasma anisotropy models. Following (Event Horizon Telescope Collaboration et
al., 2019b), we measure the image diameter $d$ as twice the distance from the
center of the image to the peak $I_{\nu}$ averaged over all directions and $w$
is the Full Width Half Maximum (FWHM) of $I_{\nu}$ averaged over all
directions. We can then infer $r_{\rm in}=(d-w)/2$ and $r_{\rm out}=(d+w)/2$ –
inner and outer radius of the image. The asymmetry parameter $A$ of the image,
defined in image plane coordinates $r_{\rm im}-\phi_{\rm im}$, is

$A=\left\langle\frac{\int_{0}^{2\pi}I(\phi_{\rm im})e^{i\phi_{\rm
im}}d\phi_{\rm im}}{\int_{0}^{2\pi}I(\phi_{\rm im})d\phi_{\rm
im}}\right\rangle_{r_{\rm im}\in[r_{\rm in},r_{\rm out}]},$ (16)

where $I(\phi_{\rm im})$ is the brightness profile across image coordinate
$\phi_{\rm im}$ at a fixed radial coordinate $r_{\rm im}$. A fully symmetric
image has $A=0$, while an antisymmetric image has $A=1$.

Figure 7: Asymmetry $A$ (a) and diameter $d$ (b) as functions of $R_{\rm
high}$ for $a=0.98$ and M87* observing angle, $\theta=163^{\circ}$, for images
at $230$ GHz blurred with $20\mu$as FWHM Gaussian kernel. (c): Asymmetry of
unblurred images at $\theta=163^{\circ}$ and $R_{\rm high}=10$ as a function
of observing frequency, $a=0.98$. (d): Asymmetry at $R_{\rm high}=10$ as a
function of observing angle $\theta$ for $a=0.98$ (solid lines) and $a=0.5$
(thin dotted lines) for Sgr A*. The green regions highlight EHT constraints
for M87* (a,b) and Sgr A* (d). In each panel, the color of the lines
represents 4 limiting cases: mirror instability, whistler instability,
isotropic plasma distribution, and firehose instability.

The asymmetry $A$ and diameter $d$ measured from raytraced images are shown in
Figure 7 with different models represented by different colors, identical
across all panels; $230$ GHz results and the variation with frequency are
shown in panels (a,b,d) and (c) respectively. Top panels show $A$ (a) and $d$
(b) for an M87* observing angle $\theta=163^{\circ}$ as functions of $R_{\rm
high}$. Panel (d) shows $A$ for $R_{\rm high}=10$ as a function of observing
angle $\theta$ for $a=0.98$ (solid lines) and $a=0.5$ (thin dotted lines) for
Sgr A*. Shaded regions indicate the allowed range as inferred from
observations for M87* (a-b) (Event Horizon Telescope Collaboration et al.,
2019b, c) and Sgr A* (d) (Event Horizon Telescope Collaboration et al.,
2022b). Panel (c) shows $A$ measured from unblurred images as a function of
frequency for an M87* viewing angle.

As expected, the difference in anisotropy $A$ between the models becomes
larger with increasing $R_{\rm high}$ (a) since larger $R_{\rm high}$
suppresses the emission from high-$\beta$ regions relative to low-$\beta$
regions, where the plasma can develop significant anisotropy (see also Table
1). The firehose case ($\eta<1$) always shows smaller $A$, consistent with the
images in Figures 3-5, while models with $\eta>1$ show higher $A$ compared to
the isotropic case. The firehose models typically have anisotropy $A$ up to
$\lesssim 3$ smaller than the mirror case, with the exact value depending on
$R_{\rm high}$ and viewing angle. As explained above, this is because plasma
at the firehose limit ($\eta<1$) emits more isotropically (over a wide range
of angles) with respect to the magnetic field direction, relative to the
mirror case which emits mostly at $\theta_{B}=90^{\circ}$. This leads to less
anisotropy in the image overall. In M87* the viewing angle is constrained to
be $\theta\approx 163^{\circ}$, while $A\approx 0.16-0.32$; thus, from Figure
7a, a better fit to the observed $A$ is obtained for $\eta<1$ at larger
$R_{\rm high}$ or $\eta>1$ with smaller $R_{\rm high}$.

We also show the diameter of the image in Figure 7b, calculated for the same
images as shown in panel (a), i.e., $a=0.98$, $\theta=163^{\circ}$, and
$R_{\rm high}=10$; the diameter is generally larger for models with larger
$\eta$. The shaded region indicates M87* constraint of $d=(42\pm 3)\mu$as
(Event Horizon Telescope Collaboration et al., 2019c). As was shown in Figure
2, the temperature $\epsilon_{\perp}$ for plasma at the firehose limit (f,l)
is smaller and varies less with radius than at the mirror limit (e,k) for both
spin parameters. This lower temperature leads to emission more concentrated
near the black hole and thus a smaller image diameter.

Figure 7d shows that the image becomes more asymmetric (d) as we look more
“edge-on” instead of “face-on”. Both spin values of $0.98$ (solid lines) and
$0.5$ (thin dotted lines) show similar behavior. The images used for panel (d)
are produced for Sgr A* with $M=4.3\times 10^{6}M_{\odot}$, and the density is
normalized such that the total flux matches EHT observations, i.e.,
$F_{\nu}=2.4$Jy at $230$ GHz and distance of $d=8178$pc. A quantitatively
similar trend, however, is also present for our M87* models. We also show the
EHT constraints on $A$ for Sgr A* by the shaded region in (d) [$A\approx
0-0.5$]. As before, plasma at the firehose anisotropy limit leads to a more
symmetric image, compared to mirror and whistler limits, at any observing
angle. Note the quite isotropic image (small A) at the firehose limit even at
$\theta=135^{\circ}$, especially for the lower spin case $a=0.5$. This effect
can significantly change the constraint on our viewing angle relative to Sgr
A* suggested by the EHT data, allowing for larger observing angles than for an
isotropic plasma.

Figure 8: Angular profiles of $n=0$ and $n=1$ brightness ($I_{0}$ [solid] and
$I_{1}$ [dotted] as functions of $\phi_{\rm im}$ in the image plane, top row)
and their ratio (bottom row) at observing angle of $\theta=163^{\circ}$ at
spin of $0.98$ (thick lines) and $0.5$ (thin lines). The first and second
columns represent $R_{\rm high}=1$ and $100$ respectively. The color of the
lines, as in Fig. 7, represents 4 limiting cases: mirror, whistler, isotropic,
and firehose.

We will now quantify the imprint of the anisotropy of the plasma distribution
function on the direct emission $I_{0}$ and the $n=1$ photon ring $I_{1}$
separately. As seen in Fig. 3 (i-p), both $n=0$ and $n=1$ emission have their
azimuthal asymmetry modified with varying $\eta$ in a way that is similar to
the full blurred image. Both are more symmetric at the firehose limit with
$\eta<1$ and more asymmetric at the mirror and whistler limits with $\eta>1$,
compared to the isotropic plasma distribution case. To distinguish the imprint
of the plasma anisotropy on the two components, we show angular profiles of
$n=0$ and $n=1$ emission ($I_{0}$ and $I_{1}$ as functions of the polar angle
in the image plane $\phi_{\rm im}$, top row, a and b) and their ratio
($I_{1}/I_{0}$, bottom row, c and d). The polar angle is plotted such that the
dimmest region of the image, $\phi_{\rm im}\sim 0$, is in the center of the
profile. This is for an observing angle of $\theta=163^{\circ}$ for both of
our spin values of $0.98$ and $0.5$ (thick and thin lines, respectively);
$R_{\rm high}=1$ and $100$ are shown in the left and right columns
respectively.

As expected, the $R_{\rm high}=100$ case shows a stronger dependence of
$I_{0}$ and $I_{1}$ on plasma anisotropy than $R_{\rm high}=1$ due to the
higher anisotropy in the low-$\beta_{\rm th}$ regions. The quantitative
dependence of the $n=0$ and $n=1$ intensities on plasma anisotropy differ
because the $n=0$ and $n=1$ photons at the same place in the image plane are
emitted at different directions relative to the local magnetic field. The
largest difference between $I_{0}$ and $I_{1}$ is reached in the case of
smaller electron temperatures at the mirror limit. In principle, measurements
of the azimuthal intensity profiles at $n=0$ and $n=1$ could thus be used to
constrain plasma anisotropy though it is unclear if this is feasible in
practice given uncertainties in black hole spin, the electron temperature,
degree of Doppler beaming, etc.

In addition to calculating the synchrotron emission and absorption produced by
an anisotropic distribution function, we have also calculated how the emitted
linear and circular polarization depends on plasma anisotropy. Because we do
not consider the impact of plasma anisotropy on Faraday rotation and
conversion in this paper, we defer a detailed discussion of the polarization
due to plasma anisotropy to future work. We can, however, quantify the change
in intrinsic linear and circular polarization, i.e. neglecting the effects of
Faraday rotation and conversion. We find that the image-averaged linear
polarisation fraction can change by up to roughly $+10$% or $-10$% for the
mirror and firehose limits respectively, compared to the isotropic case.
Circular polarization exhibits the same trend, but the mirror case can be $5$
times more circularly polarized compared to the isotropic case, at $R_{\rm
high}=100$. We also note that because models with plasma at the firehose
anisotropy have smaller $\epsilon_{\perp}$, a higher density is required to
match the observed EHT flux. This leads to an increase in pixel-averaged
optical depth, e.g.: $1.1\times 10^{-3}$, $1.2\times 10^{-3}$, $1.3\times
10^{-3}$, and $3.7\times 10^{-3}$ for the mirror, whistler, isotropic, and
firehose cases respectively at an inclination of $163^{\circ}$ and $R_{\rm
high}=10$. Thus, $\tau$ is by a factor of $3-4$ larger in the firehose case,
compared to other cases, which might also lead to a higher Faraday
depolarization.

###  3.3 Multi-wavelength observations

Figure 9: Synchrotron emission of accreting plasma raytraced from a MAD
simulation with $a=0.98$ at an inclination of $\theta=163^{\circ}$ and $R_{\rm
high}=10$ at frequencies of $345$ GHz (a-b) and $86$ GHz (c-d). The first and
second columns represent two limiting cases (mirror and firehose
respectively). Figure 10: Synchrotron emission spectra for BHs with $a=0.98$
(solid lines) and $a=0.5$ (thin dotted lines) viewed at inclinations of
$\theta=135^{\circ}$ (a,b) and $163^{\circ}$ (c,d) at $R_{\rm high}=10$ (a,c)
and $100$ (b,d). As in Fig. 7, different models are represented by different
colors.

Future mm interferometric observations will include 2 more frequencies, $345$
GHz and $86$ GHz (Johnson et al., 2023), with the latter (former) expected to
be more (less) optically thick (Chael et al., 2023a). We thus explore the
impact of an anisotropic plasma distribution function on observable images and
spectra at these frequencies. In Figure 9 we show intensity images for a BH
with $a=0.98$ at $345$ GHz on top (a-b) and $86$ GHz on the bottom (c-d), with
the parameters being identical to Figure 3 – $\theta=163^{\circ}$ and $R_{\rm
high}=10$. The mirror and firehose models are shown in the first (a,c) and
second (b,d) columns respectively. The respective images at $230$ GHz are
shown in Figure 3 for mirror (e) and firehose (h) cases.

The differences between the mirror and firehose in Figure 9 at $345$ GHz are
similar to the differences at $230$ GHz in Figure 3: the mirror case is more
azimuthally asymmetric than the firehose case. Images at $345$ GHz (a-b) are
particularly similar to their $230$ GHz counterparts because the emission is
predominantly optically thin in both cases. At lower frequency (c-d), however,
the higher synchrotron optical depth somewhat suppresses the differences
between the mirror and firehose limits and overall makes the emission more
azimuthally symmetric. Figure 7c quantifies the asymmetry $A$ as a function of
frequency for the 4 different distribution function models – the difference in
the asymmetry between the different distribution function models persists at
all frequencies though the overall asymmetry is largest at high frequencies.
In the firehose model at $86$ GHz, Figure 9d also shows that the photon ring
emission is much less evident. This is because the firehose model has a lower
temperature and higher density (at fixed $230$ GHz flux) than the other plasma
anisotropy models, and so the emission is optically thick at $86$ GHz. The
same trend, i.e. optically thin emission at high frequencies ($345$ GHz and
$230$ GHz) and optically thick emission at $86$ GHz in the firehose case,
persists at a lower spin parameter of $a=0.5$ (not shown here).

We also calculate the synchrotron emission spectra from $10^{10}$ to $10^{15}$
Hz, shown in Figure 10 at $135^{\circ}$ (left, a-b) and $163^{\circ}$ (right,
c-d) for a spin of $0.98$ (solid lines) and $0.5$ (dotted lines) at two
$R_{\rm high}$ values of $10$ (a,c) and $100$ (b,d) (shown on the left and
right side panels for each angle respectively). The different spectra for
different black hole spins are due to the higher temperatures found in more
rapidly spinning GRMHD simulations (Mościbrodzka et al., 2009). The color of
the lines is organized as in previous plots with different colors representing
different plasma anisotropy. The firehose case shows a significantly different
spectrum for both $a=0.98$ and $a=0.5$. The change is minor at low
frequencies, with firehose being slightly fainter than the other models. The
peak of the spectrum, however, can significantly shift to lower frequencies,
steepening the spectral slope just below the peak. At higher frequencies, the
emission in the firehose model is substantially fainter and the spectral slope
is steeper, compared to other cases. The qualitative results do not depend on
the value of $R_{\rm high}$. Our physical interpretation of this is that at
fixed GMRHD temperature, the firehose model (with
$T_{\parallel,e}>T_{\perp,e}$) has a lower value of $T_{\perp,e}$. This
suppresses the peak frequency of the synchrotron emission as given by Equation
14 leading to a more rapid decline in emission at high frequency.

##  4 Summary and Conclusions

Magnetized collisionless plasmas are prone to developing anisotropies in their
distribution function with respect to the magnetic field direction: the
distribution function is isotropic in the plane perpendicular to the magnetic
field because of rapid cyclotron motion (“gyrotropic”), but can be very
different along and perpendicular to the local magnetic field. In this work we
have calculated the synchrotron radiation from distribution functions with
anisotropy of this form. We are motivated by the application to low accretion
rate black holes such as those found in Sgr A* and M87* but we anticipate that
the synchrotron radiation calculations presented here will have broader
applicability.

First, we have derived and provided fits for synchrotron emissivities and
absorption coefficients for relativistic thermal electrons with an anisotropic
distribution function in (Eq. 11). The distribution function we choose (Eq. 7)
is a natural relativistic generalization of a non-relativistic bi-Maxwellian
and allows for arbitrary temperature anisotropies relative to the local
magnetic field $T_{\perp}/T_{\parallel}$ via a parameter $\eta$. The derived
fits we present are accurate to $\sim 10\%$ or better compared to numerical
solutions using the publicly available synchrotron code symphony (Pandya et
al., 2016) in the parameter range of interest (high frequency and high
temperature); the main source of error is the inaccuracy of the fits for
synchrotron emission and absorption for an isotropic thermal plasma, which our
fits are scaled to.

The change in synchrotron emission as the plasma transitions from an isotropic
to anisotropic distribution function at a fixed perpendicular temperature
$T_{\perp}$ can be understood as a renormalization of the number of particles
that emit toward the observer. The reason is that synchrotron emission emitted
at an angle $\theta_{B}$ relative to the local magnetic field is produced
primarily by particles whose vector momenta are in the same direction as
$\theta_{B}$, or, equivalently, the pitch angle of the emitting particles is
$\xi\approx\theta_{B}$. The emission thus depends on the distribution function
at pitch angle $\xi\approx\theta_{B}$. For a plasma with an isotropic
distribution function the temperature is independent of $\xi$ but temperature
anisotropy in the distribution function implies that the temperature is now
effectively a function of pitch angle $\xi$ and thus viewing angle
$\theta_{B}$ (Eq. 9). For an isotropic plasma, synchrotron emission is peaked
near $\theta_{B}\sim 90\deg$, i.e., orthogonal to the magnetic field. This
trend is enhanced for $T_{\perp}>T_{\parallel}$ (anisotropy parameter
$\eta>1$) while for $T_{\parallel}>T_{\perp}$ ($\eta<1$) the emission can peak
at significantly smaller observing angles, depending on the exact value of
$\eta$ (Fig. 1). The case of $\eta<1$ also shows more uniform emission across
observing angles than does $\eta>1$.

In addition to calculating the total emitted synchrotron radiation as a given
frequency, we have also calculated the emitted linear and circular
polarization fractions as a function of plasma anisotropy. We find that the
intrinsic linear polarization degree depends only weakly on the plasma
anisotropy $\eta$. On the other hand, circular polarization, which is very
weak in synchrotron emission from relativistic isotropic plasmas, increases
significantly for $T_{\parallel}<T_{\perp}$ at a fixed $T_{\perp}$ ($\eta>1$).
In addition, since most of the emission comes from large (small) angles
relative to the magnetic field for $\eta>1$ ($\eta<1$), the respective angle-
averaged circular polarization degree is higher (smaller).

We have employed the newly developed fits for synchrotron emission and
absorption by anisotropic electrons in a GR radiative transfer code
blacklight, capable of propagating synchrotron radiation in curvilinear space-
time. To assess how anisotropy of the accreting plasma affects mm-wavelength
observations of Sgr A* and M87*, we ray-trace GRMHD MAD simulations – the
accretion model most favored observationally (Event Horizon Telescope
Collaboration et al., 2021). Other accretion models, such as Standard and
Normal Evolution (SANE) models, are also possible. In such models the
plasma-$\beta$ is considerably higher. This suggests that the effect of plasma
anisotropies is relatively smaller in SANE models compared to MAD models, but
more detailed work in the future is required to assess this quantitatively.
Since the ideal MHD approach describes a collisional isotropic fluid, the main
source of uncertainty in this work is the temperature and temperature
anisotropy of the synchrotron-emitting electrons. In particular, the ion-to-
electron temperature ratio, which we approximate by the widely-used $R_{\rm
high}-R_{\rm low}$ model and the electron’s anisotropy $\eta$, are the main
free parameters in our study. Since $\eta$ is a prescribed quantity, absent in
our ideal GRMHD simulations, the conclusions of this work should be thought of
as qualitative rather than quantitative.

The temperature anisotropy in a collisionless plasma cannot grow without bound
because small-scale instabilities set in and limit the magnitude of the
temperature anisotropy. We thus examine the effect of an anisotropic
synchrotron emitting plasma on observed emission by considering 3 limiting
cases, defined by the anisotropy thresholds of three anisotropy-driven
instabilities: the mirror and whistler instabilities ($\eta>1$) and the
firehose instability ($\eta<1$). We present relativistic derivations of these
thresholds in Appendix B. In particular, we derive a fully kinetic mirror
instability threshold in the case of anisotropic relativistic electrons with
anisotropy parameter $\eta$ (and anisotropic non-relativistic ions with a
different anisotropy parameter $\eta_{i}$). The temperature anisotropy allowed
by kinetic-scale instabilities is larger for stronger magnetic fields, i.e.,
smaller $\beta$ (the ratio of thermal to magnetic energy density). The effects
of temperature anisotropy on observed synchrotron emission are thus likely to
be the largest when the emission is dominated by regions with $\beta\lesssim
1$, as is often the case in magnetically-arrested disk models favored on
theoretical and observational grounds.

We find that anisotropy in the accreting plasma can significantly modify the
observed synchrotron emission in horizon-scale images, including the azimuthal
asymmetry in the image plane and size of the image. This is primarily due to
the following two effects. The first effect is that the emission and
absorption for different distribution anisotropies are concentrated at
different observing angles with the magnetic field, with $\eta<1$ emitting
more uniformly across all angles as $\eta$ decreases, and $\eta>1$
emission/absorption being more concentrated near $\theta_{B}\sim 90\deg$ (Fig.
1). This can significantly modify the azimuthal asymmetry in the image plane
because different parts of the image contain radiation that was initially
emitted at different angles relative to the local magnetic field. The second
key effect is that the local perpendicular temperature $T_{\perp}$ of the
electrons changes with an assumed anisotropy $\eta$ at a given total fluid
temperature $T$ given by the GRMHD solution (Fig. 2). Models with $\eta>1$
($\eta<1$) have a larger (smaller) $T_{\perp}$, compared to the isotropic
case. Higher (lower) temperatures produce larger (smaller) 230 GHz images
because the emission at 230 GHz occurs over a larger (smaller) range of radii
(Fig. 3). Higher temperatures also lead to a smaller (higher) density of the
accreting plasma at a fixed $230$ GHz flux and thus more optically thin
(thick) emission; this is especially pronounced for $\eta>1$, i.e., the
firehose regime, in which the image-averaged optical depth can increase by a
factor of $3-4$.

More specifically we find that emission from plasma with $\eta<1$ ($\eta>1$)
produces a more azimuthally symmetric (asymmetric) image, up to a factor of
$3$ difference in the asymmetry parameter $A$. This result is of particular
interest in application to Sgr A*, where the observed EHT azimuthal asymmetry
is surprisingly modest given expectations for a random viewing angle. This
appears to suggest we are observing Sgr A* closer to face-on than not, which
is a priori surprising. Models with $\eta<1$ have significantly less variation
in the synchrotron emissivity with photon direction relative to the magnetic
field. This produces a more azimuthally symmetric image, alleviating the
restrictive constraints on viewing angle (Fig. 7d).

Anisotropy in the plasma distribution function also changes the image diameter
and the size of the central flux depression (or the observed “BH shadow”). The
smaller perpendicular temperature $T_{\perp}$ in $\eta<1$ firehose model
results in a reduced image diameter (Fig. 7b). At lower BH spins, the viewing
angle relative to the magnetic field is also smaller in the near-horizon
region. This suppresses (enhances) the emission in the image center interior
to the true photon ring (i.e., the critical curve). The BH “shadow” therefore
appears to be larger in low spin models with $\eta>1$ (Fig. 5). Chael et al.
(2021) showed that the size and shape of the “inner shadow” depend on BH spin
and our viewing angle relative to the BH spin, potentially providing a route
to measuring these quantities. Our results show that anisotropy in the
distribution function in this region close to the event horizon may be
important to consider as well.

In this paper we have not calculated the Faraday conversion coefficients for
an anisotropic plasma. We defer this to future work. We have, however,
calculated the emitted linear and circular polarization fractions and how they
depend on plasma anisotropy. We find that the imaged-averaged emitted linear
polarization fraction can increase (decrease) by up to $10\%$ in the mirror
and whistler (firehose) cases. The emitted circular polarization fraction
shows a similar trend, although the magnitude of the effect is much larger,
with the $T_{\perp}>T_{\parallel}$ regime showing an emitted circular
polarization in the mm that is up to $5$ times larger than in an isotropic
plasma.

The high frequency synchrotron emission is particularly sensitive to plasma
anisotropy. As a result, interpreting and modeling GRAVITY observations of Sgr
A* may require incorporating the effects of plasma anisotropy; this emission
is also likely non-thermal, however, so an extension of our results to non-
thermal distribution functions would be valuable.

We have also assessed how the anisotropy of the plasma affects future multi-
frequency and $n=1$ photon ring observations. We find that the effect of the
plasma distribution function on the azimuthal image asymmetry persists
throughout the frequencies of interest to future ngEHT observations, i.e. $86$
GHz and $345$ GHz, though the effect is more pronounced at higher frequencies
(Fig. 7c). We also find that the $n=1$ photon ring emission is even more
azimuthally asymmetric (symmetric) for $\eta>1$ ($\eta<1$) than the direct
$n=0$ emission, leading to an increased (decreased) ratio of photon ring to
direct emission brightness – up to a factor of $6$ in intensity ratio relative
to the isotropic distribution function case for the parameter range we
considered. Anisotropy in the distribution function has a particularly large
effect on the ratio of the $n=1$ to $n=0$ emission because plasma anisotropy
directly changes the emissivity as a function of viewing angle relative to the
magnetic field, and the $n=0$ and $n=1$ images contain emission emitted at
different angles relative to the local magnetic field.

The largest limitation of the present study as applied to modeling Sgr A*,
M87* and related sources is that the true electron temperature anisotropy in
the near-horizon environment is poorly constrained. In this work we have
attempted to bracket the magnitude of the effect that temperature anisotropy
can produce on near-horizon synchrotron radiation by considering the extreme
limit in which all of the plasma is at the temperature-anisotropy associated
with the instability thresholds for the mirror, whistler, or firehose
instabilities. The image-averaged emission-weighted electron temperature
anisotropies in these models are given in Table 1 and range from $\sim 0.1-9$.
Real systems likely do not follow just one of the limiting anisotropy models
considered here since different temperature anisotropy can co-exist in
different parts of the accretion flow. In magnetically dominated jet regions,
the plasma is in principle capable of developing large anisotropy in its
distribution function. This could occur due to differential parallel and
perpendicular heating and/or as a result of outflow-driven expansion of the
jet (as in the solar wind). Consequently, it would be interesting to apply the
methods developed here – likely extended to non-thermal distribution functions
– to model and interpret the emission from spatially extended jets (e.g., Lu
et al. 2023).

Fortunately, there is a clear path forward for improving our understanding of
the role of temperature anisotropy in the radiation from accretion flows and
jets. Global “extended” MHD models that evolve the pressure anisotropy as a
dynamical variable can predict $T_{\perp}/T_{\parallel}$ as a function of time
and space (Foucart et al., 2017), removing the need to specify the temperature
anisotropy in post-processing as we have done here; such models will, however,
needed to be extended to consider both electron and proton temperature
anisotropies. Global GRPIC simulations can go one step further and predict the
full distribution function in the accretion flow and outflow, including
temperature anisotropy, and deviations from a Maxwellian, etc. (Galishnikova
et al., 2023). One aspect that is important to account for in future modeling
is that in plasmas with $T_{p}>T_{e}$, the mirror and fluid firehose
instabilities are most sensitive to the proton temperature anisotropy (see
Appendix B). As a result it is plausible that the electron anisotropy is set
primarily by resonant instabilities such as the whistler and resonant firehose
instabilities.

We thank Alex Lupsasca for useful conversations and for sharing his analytical
calculations of source-plane emission angle for different photon sub-rings. We
also thank Chris White for useful conversations and for sharing the details of
his GRMHD MAD simulations, and Charles Gammie and Michael Johnson for useful
conversations. This research was facilitated by the Multimessenger Plasma
Physics Center (MPPC), NSF grant PHY-2206610. A.P. acknowledges support by
NASA ATP grant 80NSSC22K1054. This work was also supported in part by a Simons
Investigator Grant from the Simons Foundation (EQ), and was completed during
EQs stay at the Aspen Center for Physics, which is supported by National
Science Foundation grant PHY-2210452.

##  Appendix A Comparison of analytical synchrotron expressions with numerical
results

In this Appendix we analytically calculate the synchrotron emission and
absorption coefficients for our assumed gyrotropic distribution function and
compare the resulting analytic expressions to full numerical evaluations of
Equations 3. The analytic calculations are carried out in the limit of high
Lorentz factors for the emitting electrons, the same regime in which analytic
progress can be made for an isotropic distribution function [see, e.g.,
(Ginzburg & Syrovatskii, 1965; Melrose, 1971)].

###  A.1 Derivation of Analytical Fits for Total Intensity, Linear
Polarization, and Circular Polarization

Under the assumption of high Lorentz factor $\gamma$ (or energy $E$) for the
emitting electrons, the emission is predominantly concentrated in a narrow
cone around the pitch angle $\mu\simeq\cos(\theta_{B})$ where $\theta_{B}$ is
the viewing angle with respect to the magnetic field. Following Melrose
(1971), the electron emissivity in the Stokes basis from Equations 3 and 5 in
the main text can be expressed in tensor form as

$\begin{split}j^{\alpha\beta}&=\int
d^{3}pf(E,\xi)\eta^{\alpha\beta}=\int_{0}^{\infty}dEf(E,\theta_{B})\frac{\sqrt{3}e^{2}\nu_{c}\sin\theta_{B}}{8\pi
c}H^{\alpha\beta}(X),\\\
&H^{11}(X)=X\left[\int_{X}^{\infty}dtK_{5/3}(t)+K_{2/3}(X)\right],\\\
&H^{22}(X)=X\left[\int_{X}^{\infty}dtK_{5/3}(t)-K_{2/3}(X)\right],\\\
&H^{12}(X)=-H^{21}=-\frac{2i\cot\theta_{B}}{3\gamma}\left[(2+g(\theta_{B}))\int_{X}^{\infty}dtK_{1/3}(t)+2XK_{1/3}(X)\right],\\\
\end{split}$ (A1)

where $\nu_{c}=eB/2\pi m_{e}c$ is a non-relativistic cyclotron frequency,
$X=\nu/\nu_{cr}$ and $\nu_{cr}=(3/2)\nu_{c}\gamma^{2}\sin\theta_{B}$. The
first expression in Equation A1 is general while in the second expression we
have integrated over pitch angle $\xi$ by assuming $\xi\simeq\theta_{B}$. The
Stokes emissivities are related to Equation A1 as $j_{I}=j^{22}+j^{11}$,
$j_{Q}=j^{22}-j^{11}$, $j_{U}=j^{12}+j^{21}\equiv 0$, and
$j_{V}=i(j^{12}-j^{21})$. Here, unlike in Melrose (1971), we define
$g(\theta_{B})$ for a general non-separable gyrotropic distribution function,
which for our choice of the distribution (Equation 7 in the main text) is

$g(\theta_{B})=\tan\theta_{B}\left.\frac{df(E,\xi)}{d\xi}\right|_{\xi=\theta_{B}}\frac{1}{f(E,\theta_{B})}=\frac{\gamma}{\epsilon_{\perp}^{\ast}}\frac{(\eta-1)\sin^{2}\theta_{B}}{1+(\eta-1)\cos^{2}\theta_{B}}=\frac{\gamma}{\epsilon_{\perp}^{\ast}}\left(\frac{\epsilon_{\perp}^{\ast}}{\epsilon_{\perp}}\right)^{2}(\eta-1)\sin^{2}\theta_{B}\equiv
A\gamma.$ (A2)

In the last equality in Equation A2 we have defined the anisotropy parameter A
(a function of $\eta)$ that will appear below.

We now proceed analytically evaluating the emissivities $j_{I}$, $j_{Q}$ and
$j_{V}$, beginning with $j_{I}$. Equation A1 for $j_{I}$ can be rewritten as:

$j_{I}(\epsilon_{\perp},\nu/\nu_{c},\theta_{B},\eta)=\frac{\sqrt{3}Bm_{e}^{2}ce^{3}\sin(\theta_{B})}{8\pi}\int
d\gamma\gamma^{2}\beta
f(\gamma,\theta_{B})F(X)=\eta^{1/2}\frac{\sqrt{3}n_{e}Be^{3}\sin(\theta_{B})}{32\pi^{2}m_{e}c^{2}\epsilon_{\perp}K_{2}(\epsilon_{\perp})}\int
d\gamma\gamma^{2}\beta e^{-\gamma/\epsilon_{\perp}^{\ast}}F(X),$ (A3)

where

$F(X)=X\int_{X}^{\infty}dtK_{5/3}(t)=\begin{cases}2^{2/3}\Gamma(1/3)X^{1/3}+\mathcal{O}(X),\
X\ll 1,\\\ \sqrt{\frac{\pi}{2}X}e^{-X}(1+\mathcal{O}(1/X)),\ X\gg
1\end{cases}$ (A4)

is the asymptotic behavior of the synchrotron power at low and high
frequencies and $\Gamma(a)$ is the gamma-function. To express the emissivity
in terms of the new temperature
$\epsilon_{\perp}^{\ast}=\epsilon_{\perp}^{\ast}(\xi=\theta_{B})=\epsilon_{\perp}/\sqrt{1+(\eta-1)\cos^{2}\theta_{B}}$,
as given by distribution in Equation 8, we consider separately the low and
high frequency limits in Equation A4 applied to Equation A3. In the low-
frequency limit,

$j_{I}(\epsilon_{\perp},\nu/\nu_{c},\theta_{B},\eta)\propto\eta^{1/2}\int_{1}^{\infty}d\gamma\gamma^{2}\beta\frac{e^{-\gamma/\epsilon_{\perp}^{\ast}}}{\epsilon_{\perp}K_{2}(\epsilon_{\perp})}\gamma^{-2/3}\approx\frac{\eta^{1/2}}{\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}\int_{1}^{\infty}d\gamma\gamma^{4/3}e^{-\gamma/\epsilon_{\perp}^{\ast}}\approx\eta^{1/2}\frac{{\epsilon_{\perp}^{\ast}}^{7/3}}{\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}.$
(A5)

Therefore, the final expression for the emissivity is

$j_{I}(\epsilon_{\perp},\nu/\nu_{c},\theta_{B},\eta)\approx\frac{2^{4/3}\pi\eta^{1/2}}{3}\frac{n_{e}e^{2}\nu_{s}^{\ast}}{c{K_{2}(1/\epsilon_{\perp})}}\left(\frac{\nu}{\nu_{s}^{\ast}}\right)^{1/3}\left(\frac{\epsilon_{\perp}^{\ast}}{\epsilon_{\perp}}\right)=\eta^{1/2}\frac{\epsilon_{\perp}^{\ast}}{\epsilon_{\perp}}\frac{K_{2}(1/\epsilon_{\perp}^{\ast})}{K_{2}(1/\epsilon_{\perp})}j_{I,{\rm
iso}}(\epsilon=\epsilon_{\perp}^{\ast},\nu/\nu_{c},\theta_{B}),$ (A6)

where $\nu_{s}^{\ast}=\frac{2}{9}\nu_{c}\epsilon_{\perp}^{\ast
2}\sin\theta_{B}$ and
$\epsilon_{\perp}^{\ast}=\epsilon_{\perp}^{\ast}(\xi=\theta_{B})$. This
calculation was done in the limit of low $\nu$, but the same expression can be
obtained in the limit of high $\nu$ as well. The integral over Lorentz factor
in Equation A3 now becomes

$j_{I}(\epsilon_{\perp},\nu/\nu_{c},\theta_{B},\eta)\propto\eta^{1/2}\int_{1}^{\infty}d\gamma\gamma\beta\frac{e^{-\gamma/\epsilon_{\perp}^{\ast}-X}}{\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}.$
(A7)

The maximum of the exponent in Equation A7 occurs at
$\gamma_{0}=(2B\epsilon_{\perp}^{\ast})^{1/3}$, where
$B=(\nu/\nu_{cr})\gamma^{2}=2/3(\nu/\nu_{c})\sin^{-1}\theta_{B}\gg 1$. The
integral over $\gamma$ can then be carried out using the method of steepest
descent (as in the case of an isotropic distribution function), leading again
to

$j_{I}(\epsilon_{\perp},\nu/\nu_{c},\theta_{B},\eta)=\eta^{1/2}\frac{\epsilon_{\perp}^{\ast}}{\epsilon_{\perp}}\frac{K_{2}(1/\epsilon_{\perp}^{\ast})}{K_{2}(1/\epsilon_{\perp})}j_{I,{\rm
iso}}(\epsilon=\epsilon_{\perp}^{\ast},\nu/\nu_{c},\theta_{B}).$ (A8)

The fact that $j_{I}$ for the anisotropic relativistic Maxwellian can be
expressed as Equation A6 in both the low and high frequency limits motivates
our using this expression as the proposed fit in Equation 11 of the main text.
Physically, this corresponds to the total intensity emissivity just changing
due to a different effective distribution function in the angle $\theta_{B}$
towards the observer. Note as well that although we derived Equation A6 for
total intensity the same expression scaled to the isotropic distribution
function emissivity holds for the intrinsic linear polarization emissivity,
i.e., $j_{Q}$. This is because $K_{2/3}(X)$ has the same functional form as
$\int_{X}^{\infty}K_{5/3}(X)$ at both high and low frequencies.

Circular polarization, however, has a different functional form:

$j_{V}(\epsilon_{\perp},\nu/\nu_{c},\theta_{B},\eta)\propto\eta^{0.5}\cot\theta_{B}\sin\theta_{B}\int_{1}^{\infty}d\gamma\gamma\frac{e^{-\gamma/\epsilon_{\perp}^{\ast}}}{\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}\left[(g(\gamma,\theta_{B})+2)\int_{X}^{\infty}K_{1/3}(t)dt+2XK_{1/3}\left(X\right)\right].$
(A9)

Unlike in the case of total intensity and linear polarization, the circular
polarization emissivity requires expanding the distribution function in a
narrow cone around $\theta_{B}$; the resulting $j_{V}$ depends on the
derivative of the distribution function, included in $g(\gamma,\theta_{B})$.
To understand the origin of our fit for $j_{V}$ in Equations 11, we first
consider the high-frequency limit when both $\int_{X}^{\infty}K_{1/3}(t)dt$
and $K_{1/3}\left(X\right)$ scale as $e^{-X}/\sqrt{X}$ for $X\gg 1$. The
integrand in equation A9 can be written as $h(\gamma,A)e^{S(\gamma,A)}$, where

$S(\gamma,A)=-\gamma/\epsilon_{\perp}^{\star}-B/\gamma^{2}\ \ \ \ {\rm and}\ \
\ \ h(\gamma,A)=(A\gamma^{3}+2\gamma^{2}+2B)\approx(A\gamma^{3}+2B).$ (A10)

The exponential term $e^{S(\gamma)}$ is again maximum at
$\gamma_{0}\approx(2B\epsilon_{\perp}^{\ast})^{1/3}$, where
$B=(\nu/\nu_{cr})\gamma^{2}=2/3(\nu/\nu_{c})\sin^{-1}\theta_{B}\gg 1$.
Equation A9 can then be integrated via the method of steepest descent as:

$\int d\gamma
h(\gamma)e^{S(\gamma)}\approx\sqrt{\frac{2\pi}{-S^{\prime\prime}(\gamma_{0})}}e^{S(\gamma_{0})}h(\gamma_{0})=\sqrt{\frac{2\pi}{-S^{\prime\prime}(\gamma_{0})}}e^{S(\gamma_{0})}\times
2B(A\epsilon_{\perp}^{\ast}+1)=\sqrt{\frac{2\pi}{-S^{\prime\prime}(\gamma_{0})}}e^{S(\gamma_{0})}\times
2B\eta\left(\frac{\epsilon_{\perp}^{\ast}}{\epsilon_{\perp}}\right)^{2}.$
(A11)

As with $j_{I}$ and $j_{Q}$ we choose to express $j_{V}$ relative to the
result for an isotropic Maxwellian with temperature
$\epsilon^{\star}_{\perp}$. The latter can be derived in an identical manner
to Equation A11. We find that ratio of $j_{V}$ in the anisotropic case to
$j_{V,{\rm iso}}$ at a temperature of $\epsilon_{\perp}^{\ast}$ and $A=0$ has
two terms. One is the ratio of distribution function normalizations
$\eta^{1/2}(\epsilon^{\star}_{\perp}/\epsilon_{\perp})(K_{2}(1/\epsilon^{\star}_{\perp})/K_{2}(1/\epsilon_{\perp}))$
that appears in $j_{I}$ and $j_{Q}$. The other is the factor
$\eta(\epsilon^{\star}_{\perp}/\epsilon_{\perp})^{2}$ in Equation A11 –
present only in $j_{V}$ and not $j_{Q}$ and $j_{I}$ – that is due to the
presence of the distribution function derivative $g(\theta_{B})$ in the
circular polarization emissivity. The net result is

$\frac{j_{V}(\epsilon_{\perp},\nu/\nu_{c},\theta_{B},\eta)}{j_{V,{\rm
iso}}(\epsilon=\epsilon_{\perp}^{\ast},\nu/\nu_{c},\theta_{B})}=\eta\left(\frac{\epsilon_{\perp}^{\ast}}{\epsilon_{\perp}}\right)^{2}\times\eta^{1/2}\left(\frac{\epsilon_{\perp}^{\ast}}{\epsilon_{\perp}}\right)\left(\frac{K_{2}(1/\epsilon_{\perp}^{\ast})}{K_{2}(1/\epsilon_{\perp})}\right)=\eta^{3/2}\left(\frac{\epsilon_{\perp}^{\ast}}{\epsilon_{\perp}}\right)^{3}\left(\frac{K_{2}(1/\epsilon_{\perp}^{\ast})}{K_{2}(1/\epsilon_{\perp})}\right),$
(A12)

which gives the analytical fit given by Equation 11 in the main text. The same
result can be derived in the low frequency limit via suitable expansion of
Equation A9.

###  A.2 Comparison of Analytics and Numerics

Figure 11: Integrands $|K_{I}f(\gamma,\xi)|$ for $j_{I}$ (a) and
$|K_{I}Df(\gamma,\xi)|$ for $\alpha_{I}$ (b). Two parts of $Df$ that include
$\partial_{\gamma}f(\gamma,\xi)$ (c) and $\partial_{\mu}f(\gamma,\xi)$ (d).
The approximate location of the peak that corresponds to
$\cos\theta_{B}=\beta\mu$ is shown by dotted lines. The free parameters are
$\nu/\nu_{c}=10^{3}$, $\eta=10$, $\epsilon_{\perp}=10$. While the location of
the peak is still the same as in the $\eta=1$ case, a non-zero term with
$\partial_{\mu}f(\gamma,\mu)$ appears which, however, goes through zero at the
peak of the integrand. Note the saturated colorbar in (c) and (d).

We now solve Equations 3-5 in the main text numerically and check the validity
of the approximations used in the previous section for obtaining analytical
fits for the polarized synchrotron emissivity and absorption coefficients. To
do so, we use the publicly available code symphony to compare our theoretical
fits (Eq.11) with a numerical solution. We implemented an anisotropic
distribution function to calculate $j_{S}$ and $\alpha_{S}$. In particular, we
added the possibility for the distribution to depend on harmonic number $n$ as
well as a non-zero $\partial_{\mu}f$ term in the absorption coefficient
calculation (Eq. 4 that shows up in Eq. 3 includes $\partial_{\mu}f$) – both
were absent in symphony. The distribution function and analytical derivatives
$\partial_{\gamma}f$ and $\partial_{\mu}f$ can now depend on $\mu=\cos\xi$.
However, as described in Section 2 and below, the term with $\partial_{\mu}f$
in the absorption coefficient is negligible because it shows up proportional
to a term that vanishes when the pitch angle is approximately $\theta_{B}$.

The integrands in Equations 3 are $K_{a}f(\gamma,\xi)$ for $j_{a}$ and
$K_{a}Df(\gamma,\xi)$ for $\alpha_{a}$, where $\xi$ can be substituted for $n$
since at $y_{n}=0$ (as required by $\delta(y_{n})$):

$\cos\xi=\frac{1-\frac{n}{\gamma}\frac{\nu_{c}}{\nu}}{\beta\cos\theta_{B}}.$
(A13)

Thus, the integrands can be expressed as functions of $\gamma$ and $n$ and
integrated in $\gamma-n$ space.

In Figure 11 we show the integrands for $j_{I}$ (a) and $\alpha_{I}$ (b) and
the two terms from $Df$ that include $\partial_{\gamma}f(\gamma,\xi)$ and
$\partial_{\mu}f(\gamma,\xi)$ for $\nu/\nu_{c}=10^{3}$, $\epsilon_{\perp}=10$,
and $\eta=10$. The location of the sharp peak in the $\gamma-n$ plane is where
$\xi\simeq\theta_{B}$ (as in the isotropic distribution function case).
However, the exact harmonic at which the emission peaks moves along this line
depending on $\eta$. This is equivalent to the result in Figure 1 that
different $\theta_{B}$ dominate the emission as we vary $\eta$. Panels (c) and
(d) in Figure 11 show that, at the location in $\gamma-n$ space where the
absorption coefficient peaks (panel b), the first term in $Df$ due to
gradients in $\gamma$ is much larger than the second term due to gradients in
$\mu$. This is because most of the emission and absorption is coming from
pitch angles of $\xi\approx\theta_{B}$. Thus, the propagation is almost
parallel, and the term with $\beta\cos\theta_{B}-\mu$ shown in panel (d) does
not contribute significantly to $Df$. This implies that in practice the total
intensity emission and absorption coefficients for the anistropic distribution
function are equivalent to calculations for a thermal isotropic distribution
function at a new temperature $\epsilon^{\star}_{\perp}$. This allows us to
calculate $\alpha_{a}$ from $j_{a}$ via Kirchoff’s law even for our
anisotropic distribution function (at least in the limit of high $\gamma$
where $\xi\approx\theta_{B}$ is justified).

A number of fitting functions for $j_{a}$ and $\alpha_{a}$ are used in the
literature [see, e.g., Pandya et al. (2016); Dexter (2016)]. Here we compare
our results for the fits used by blacklight:

$j_{S,{\rm
iso}}(\epsilon,X,\theta_{B})=\frac{n_{e}e^{2}\nu_{c}}{c}e^{-X^{1/3}}\times\begin{cases}\frac{\sqrt{2}\pi}{27}\sin(\theta_{B})(X^{1/2}+2^{11/12}X^{1/6})^{2}\
[S=I],\\\
-\frac{\sqrt{2}\pi}{27}\sin(\theta_{B})\Big{(}X^{1/2}+\Big{(}\frac{7\epsilon^{24/25}+35}{10\epsilon^{24/25}+75}\Big{)}2^{11/12}X^{1/6}\Big{)}^{2}\
[S=Q],\\\ 0\ [S=U],\\\
\frac{\cos\theta_{B}}{\epsilon}\Big{(}\frac{\pi}{3}+\frac{\pi}{3}X^{1/3}+\frac{2}{300}X^{1/2}+\frac{2}{19}\pi
X^{2/3}\Big{)}\ [S=V],\end{cases}$ (A14)

where $X=\nu/\nu_{s}$. Absorption coefficients $\alpha_{S}$ for a thermal
distribution can be obtained via Kirchoff’s law.

Figure 12: Comparison of numerical results for $j_{I}$ (a-b), $\alpha_{I}$
(c-d), $j_{Q}$ (e-f), and $j_{V}$ (g,h) with the theoretical fits given by
Equations 11 and A14. Numerical results and theoretical fits are shown on the
left by solid and dotted lines respectively, and on the right – the relative
error is shown. The dashed gray line on the right shows a relative error of
$30\%$. The free parameters are $\epsilon_{\perp}=10$ and
$\nu/\nu_{c}=10^{3}$.

.

Figure 13: Same as Figure 12 but for $\epsilon_{\perp}=3$ and
$\nu/\nu_{c}=10$.

.

Figure 12 shows numerical integration results from symphony (solid lines)
along with their respective theoretical fits (dotted lines) for $j_{I}$ (a),
$\alpha_{I}$ (b), $j_{Q}$ (c), and $j_{V}$ (d) on the left. On the right,
their respective absolute errors are shown in panels (b,d,f,h). All results
are shown as a function of observing angle $\theta_{B}$ at different
anisotropy parameters $\eta$ represented by different colors, at
$\nu/\nu_{c}=10^{3}$ and for $\epsilon_{\perp}=10$. These are typical
parameters for application to Sgr A* and M87*. The agreement is excellent for
all $\eta$, with maximal errors $\lesssim 10\%$

We show a more challenging case of low temperature $\epsilon_{\perp}=3$ and
low frequency $\nu/\nu_{c}=10$ in Figure 13, which is organized identically to
Fig. 12. This case is more challenging for our analytic fits than Figure 12
because the emission for $\epsilon_{\perp}=3$ and low frequency
$\nu/\nu_{c}=10$ is dominated by much lower energy electrons. The errors in
our fits in Figure 13 are, not surprisingly, larger. Generally, $\eta<1$ has
smaller relative errors than $\eta>1$. This is because at a fixed observing
angle $\theta_{B}$ and $\epsilon_{\perp}$, the effective temperature
$\epsilon_{\perp}^{\ast}$ is larger than $\epsilon_{\perp}$ for $\eta<1$. By
contrast, $\epsilon_{\perp}^{\ast}<\epsilon_{\perp}$ for $\eta>1$, which can
start to approach the non-relativistic cyclotron limit for which our fits do
not apply. Figure 13 shows that the case of $\eta<1$ has a relative error of
$<30\%$ across all considered angles $\theta_{B}\in[5,85]^{\circ}$, and
$<10\%$ for most angles. The fits have a relative error of larger than $30\%$
for $\eta\geq 1$, however, at the largest and smallest angles. This is true
for the isotropic case as well at small angles. We note, though, that the
actual value of the emissivity and absorption coefficient are very small at
small angles for $\eta\gtrsim 1$ (Fig. 6) so that most of the emission and
absorption will arise at larger angles where the fits are better. In addition
at this low frequency, the emission will in most practical cases of interest
be self-absorbed and approximately a blackbody. Finally, we note that a
significant cause of error here is that our fits in Equation 11 for the
anisotropic emission and absorption coefficients are scaled to the isotropic
emissivity and absorption fits in Equation A14, which become inaccurate at low
temperatures, low frequencies, and small angles, as indicated by the large
fractional errors for the isotropic case in Figure 13. In practice we advise
caution in using the fits here if $\epsilon^{\star}_{\perp}\lesssim 3$ and the
frequency is low $\lesssim 10\nu_{c}$. The regime of most interest for our
applications is much higher frequencies where the analytic fits in Equation 11
are accurate.

##  Appendix B Anisotropy-driven instabilities in relativistic plasmas

###  B.1 Mirror instability

To calculate the kinetic threshold for the relativistic mirror instability, we
consider Vlasov and Maxwell’s equations:

$\displaystyle\frac{\partial f_{s}}{\partial t}+{\bf v}_{s}\cdot\nabla
f_{s}+q_{s}({\bf E}+\frac{{\bf v}_{s}}{c}\times{\bf B})\cdot\frac{\partial
f_{s}}{\partial{\bf p}}=0,$ (B1) $\displaystyle\frac{1}{c}\frac{\partial{\bf
E}}{\partial t}=\nabla\times{\bf B}-\frac{4\pi}{c}{\bf j},$ (B2)
$\displaystyle\frac{1}{c}\frac{\partial{\bf B}}{\partial t}=-\nabla\times{\bf
E},$ (B3)

where $s$ is the particle species (ions $i$ or electrons $e$) with mass
$m_{s}$ and charge $q_{s}$, ${\bf v}_{s}={\bf p}_{s}/m_{s}\gamma$, ${\bf E}$
is electric field (with ${\bf E}_{0}=0$ initially), ${\bf B}$ is the magnetic
field, and the axes are chosen such that ${\bf B}=B_{0}\hat{z}$. We now
consider a small perturbation in the form of displacement $\propto e^{i{\bf
kr}-i\omega t}$, where we consider ${\bf
k}=k_{\perp}\hat{x}+k_{\parallel}\hat{z}$. We will initially consider
electrons with an anisotropic distribution and ions with an isotropic
distribution, $\delta{\bf E}=\delta E_{y}\hat{y}$ and thus $\delta{\bf
B}=\delta B_{x}\hat{x}+\delta B_{z}\hat{z}$; the distribution function is
perturbed as $f_{s}+\delta f_{s}$. The linearized equations are then

$\displaystyle(-i\omega+i{\bf k}\cdot{\bf v}_{s})\delta f_{s}+q_{s}\frac{{\bf
v}_{s}}{c}\times{\bf B}_{0}\cdot\partial_{\bf p}\delta
f_{s}+q_{s}\left(\delta{\bf E}+\frac{{\bf v}_{s}}{c}\times\delta{\bf
B}\right)\cdot\partial_{\bf p}f_{s}=0,$ (B4)
$\displaystyle\frac{4\pi}{c^{2}}\omega\delta
j_{y}=i\left(\frac{\omega^{2}}{c^{2}}-k^{2}\right)\delta E_{y},\
\frac{4\pi}{c^{2}}\omega\delta
j_{x}=i\left(\frac{\omega^{2}}{c^{2}}-k_{\parallel}^{2}\right)\delta E_{x},$
(B5)

where we used $\delta{\bf B}=\frac{c}{\omega}{\bf k}\times\delta{\bf E}$. We
seek the solution of the linearized Vlasov equation for $\delta f_{s}$ and the
corresponding current via the method of characteristics (e.g., Mikhailovsky,
1976). The current response $\delta j_{y}$ due to $\delta E_{y}$ is:

$\delta j_{y,s}=-2\pi iq_{s}^{2}\int dpd\mu
p^{2}\sum_{n=-\infty}^{+\infty}\frac{v^{2}\sin^{2}\xi}{\omega-
k_{\parallel}v\cos\xi-n\Omega_{s}}\left[\frac{1}{v}\frac{\partial f}{\partial
p}-\frac{\cos\xi}{vp}\frac{\partial
f}{\partial\mu}+\frac{k_{\parallel}}{\omega}\frac{1}{p}\frac{\partial
f}{\partial\mu}\right]J^{{}^{\prime}2}_{n}\left(\frac{k_{\perp}v_{\perp}}{\Omega_{s}}\right)\delta
E_{y},$ (B6)

where $\Omega_{s}$ is the relativistic cyclotron frequency of species $s$. For
$\Omega_{s}\gg\omega$ and $\Omega_{s}\gg k_{\parallel}v_{\parallel}$, keeping
the leading terms $n=0,\pm 1$ of order $\Omega_{s}^{-2}$ and using
$J^{{}^{\prime}}_{0}(z)\approx-\frac{z}{2}$ and $J^{{}^{\prime}}_{\pm
1}(z)\approx\pm\frac{1}{2}$:

$\begin{split}n=0:\ &-\frac{\pi iq_{s}^{2}}{2}\int dpd\mu
p^{2}\frac{v^{4}\sin^{4}\xi}{\omega-
k_{\parallel}v\cos\xi}\left[\frac{1}{v}\frac{\partial f_{s}}{\partial
p}-\frac{\cos\xi}{vp}\frac{\partial
f_{s}}{\partial\mu}+\frac{k_{\parallel}}{\omega}\frac{1}{p}\frac{\partial
f_{s}}{\partial\mu}\right]\frac{k_{\perp}^{2}}{\Omega_{s}^{2}}\delta E_{y},\\\
n=\pm 1:\ &\frac{\pi iq_{s}^{2}}{2}\int dpd\mu
p^{2}v^{2}\sin^{2}\xi\frac{\omega-
k_{\parallel}v\cos\xi}{\Omega_{s}^{2}}\left[\frac{1}{v}\frac{\partial
f_{s}}{\partial p}-\frac{\cos\xi}{vp}\frac{\partial
f_{s}}{\partial\mu}+\frac{k_{\parallel}}{\omega}\frac{1}{p}\frac{\partial
f_{s}}{\partial\mu}\right]\delta E_{y}.\end{split}$ (B7)

For isotropic ions the terms with $\partial_{\mu}f$ can be dropped, resulting
in the following current response

$\delta j_{y,i}=\frac{\pi ic}{2B^{2}}\delta E_{y}\int dpd\mu
p^{3}\sqrt{m_{i}^{2}c^{2}+p^{2}}S(\mu)\frac{\partial f_{i}}{\partial
p}\left[\omega-k_{\parallel}v\mu-\frac{v^{2}S(\mu)}{\omega-
k_{\parallel}v\mu}k_{\perp}^{2}\right],$ (B8)

where the second term in brackets equals zero due to the odd function
$S(\mu)=1-\mu^{2}=\sin^{2}\xi$. For the mirror mode, we are interested in the
$k_{\parallel}v_{\parallel}/\omega\gg 1$ limit, which leaves only the resonant
term. Using

$\frac{v^{2}}{\omega}\int_{-1}^{1}d\mu\frac{S^{2}(\mu)}{1-\frac{k_{\parallel}v\mu}{\omega}}=-i\pi\frac{v}{k_{\parallel}}+\frac{16}{3}\frac{\omega}{k_{\parallel}^{2}}+\mathcal{O}((k_{\parallel}v_{\parallel}/\omega)^{-3}),$
(B9)

and considering non-relativistic ions, $p/m_{i}c\ll 1$, with a Maxwellian
distribution function $f_{i}(p)$ and number density $n$

$f_{i}(p)=\frac{n}{(2\pi
m_{i}^{2}c^{2}\epsilon_{i})^{3/2}}e^{-p^{2}/2m_{i}^{2}c^{2}\epsilon_{i}},$
(B10)

integration by parts of the third resonant term in Eq. B8 results in

$\delta j_{y,i}=\frac{2\pi^{2}c^{2}}{B_{0}^{2}}\delta
E_{y}\frac{k_{\perp}^{2}}{k_{\parallel}}\int_{0}^{\infty}dpp^{3}f_{s}(p)=\frac{\delta
E_{y}}{\sigma_{i}}\frac{c\epsilon_{i}^{1/2}}{4(\pi/2)^{1/2}}\frac{k_{\perp}^{2}}{k_{\parallel}}=\frac{\delta
E_{y}}{4\pi\sigma_{i}}\pi^{1/2}v_{th,i}\frac{k_{\perp}^{2}}{k_{\parallel}},$
(B11)

where $v_{th,i}=\sqrt{2\epsilon_{i}}c$ and $\sigma_{i}=B_{0}^{2}/4\pi
nm_{i}c^{2}=v_{a}^{2}/c^{2}$, $v_{a}^{2}=B_{0}^{2}/4\pi nm_{i}$ is the Alfven
speed.

We will now analyze the electron’s current in Equation B6, splitting it by the
three terms in the brackets $j_{y,e,1}$, $j_{y,e,2}$, and $j_{y,e,3}$. The
second term in Eq. B6 is $\mu\omega/k_{\parallel}v\ll 1$ times smaller than
the third term and thus $j_{y,e,2}$ is negligible. As with the ions,
considering the resonant term’s residue of $-i\pi\omega/k_{\parallel}v$ at
$\mu_{0}=\omega/k_{\parallel}v$

$\delta j_{y,e,1}\approx-\frac{\pi^{2}c^{2}}{2B^{2}}\delta
E_{y}\frac{k_{\perp}^{2}}{k_{\parallel}}\int dpp^{4}\frac{\partial
f_{e}(p,\mu_{0})}{\partial p}=\frac{\delta
E_{y}}{4\pi\sigma_{i}}\frac{2\pi^{2}}{m_{i}n}\frac{k_{\perp}^{2}}{k_{\parallel}}\int
dpp^{3}f_{e}(p,\mu_{0}).$ (B12)

The dispersion relation in the limit of $\omega\ll k_{\parallel}v\mu$ is
therefore

$-\omega\frac{k_{\perp}^{2}}{k_{\parallel}}\mathcal{J}=ik^{2}v_{a}^{2}+4\pi
k_{\parallel}c\sigma_{i}\omega\frac{\delta j_{y,e,3}}{\delta B_{x}},$ (B13)

where

$\mathcal{J}=\frac{2\pi^{2}}{m_{i}n}\int
dpp^{3}f_{e}(p,\mu_{0})+\pi^{1/2}v_{th,i}>0.$ (B14)

As in Osipov et al. (2017), the current response from anisotropic electrons,
which drives the mirror instability, can be calculated in the same form:

$\delta j_{y,e,3}=-i\frac{\pi c\delta
B_{x}}{2B_{0}^{2}}\frac{k_{\perp}^{2}}{k_{\parallel}}\int
dpp^{3}v\int_{-1}^{1}d\mu\frac{(1-\mu^{2})^{2}}{\mu}\frac{\partial
f_{e}(p,\mu)}{\partial\mu},$ (B15)

which we will calculate in two parts

$\int_{-1}^{1}d\mu\frac{(1-\mu^{2})^{2}}{\mu}\frac{\partial
f_{e}}{\partial\mu}=\underbrace{\int_{-1}^{1}d\mu(\mu^{3}-2\mu)\frac{\partial
f_{e}}{\partial\mu}}_{I_{1}}+\underbrace{\int_{-1}^{1}d\mu\frac{1}{\mu}\frac{\partial
f_{e}}{\partial\mu}}_{I_{2}}.$ (B16)

Integral $I_{1}$ can be calculated by parts and expressed through parallel and
perpendicular pressure $P_{\parallel,e}$ and $P_{\perp,e}$ since

$I_{1}=-2f_{e}(1)-\int_{-1}^{1}f_{e}(p,\mu)d(\mu^{3}-2\mu)=-2f_{e}(1)+{2\int_{-1}^{1}d\mu(1-\mu^{2})f_{e}(p,\mu)-\int_{-1}^{1}d\mu\mu^{2}f_{e}(p,\mu)},$
(B17)

and

$P_{\parallel,e}=2\pi\int_{0}^{\infty}dpp^{3}v\int_{-1}^{1}d\mu\mu^{2}f_{e}(p,\mu),\
P_{\perp,e}=\pi\int_{0}^{\infty}dpp^{3}v\int_{-1}^{1}d\mu(1-\mu^{2})f_{e}(p,\mu).$
(B18)

Therefore, the two integral terms in $I_{1}$ in Eq. B17 lead to

$2\pi\times
2\int_{0}^{\infty}dpp^{3}v\int_{-1}^{1}d\mu(1-\mu^{2})f_{e}(p,\mu)-2\pi\int_{0}^{\infty}dpp^{3}v\int_{-1}^{1}d\mu\mu^{2}f_{e}(p,\mu)=4P_{\perp,e}-P_{\parallel,e}=P_{\parallel,e}\left(4\eta^{\lambda}-1\right).$
(B19)

For the distribution function given by Equation 7 in the main text,
$f_{e}(p,\mu)\propto\exp{(-a\sqrt{1+b\mu^{2}})}$, where
$a=\gamma/\epsilon_{\perp}$ and $b=(1-1/\gamma^{2})(\eta-1)$. Then, the
boundary term in Eq. B17 can be expressed as

$-2\times
2\pi\int_{0}^{\infty}dpp^{3}vf_{e}(p,\mu=1)=-\frac{nm_{e}c^{2}\eta^{1/2}}{\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}e^{-a\sqrt{1+b}}=-P_{\parallel,e}\frac{\eta^{1/2+\lambda}}{\epsilon_{\perp}^{2}K_{2}(1/\epsilon_{\perp})}\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}e^{-a\sqrt{1+b}},$
(B20)

resulting in the following contribution to $\delta j_{y,e,3}$ from $I_{1}$:

$2\pi\int_{0}^{\infty}dpp^{3}vI_{1}=P_{\parallel,e}\left[4\eta^{\lambda}-1-\frac{\eta^{1/2+\lambda}}{\epsilon_{\perp}^{2}K_{2}(1/\epsilon_{\perp})}\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}e^{-a\sqrt{1+b}}\right].$
(B21)

For calculating $I_{2}$, we find the derivative of the distribution function
as $\partial
f_{e}(p,\mu)/\partial\mu=-ab\frac{\mu}{\sqrt{1+b\mu^{2}}}f_{e}(p,\mu)$. Since
the integrand is an even function of $\mu$, the contribution of $I_{2}$ to the
current is

$\begin{split}2\pi\int_{0}^{\infty}dpp^{3}vI_{2}=-\frac{nm_{e}c^{2}\eta^{1/2}}{\epsilon_{\perp}K_{2}(1/\epsilon_{\perp})}\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}ab\int_{0}^{1}d\mu\frac{e^{-a\sqrt{1+b\mu^{2}}}}{\sqrt{1+b\mu^{2}}}\\\
=-P_{\parallel,e}\frac{\eta^{1/2+\lambda}}{\epsilon_{\perp}^{2}K_{2}(1/\epsilon_{\perp})}\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}ab\int_{0}^{1}d\mu\frac{e^{-a\sqrt{1+b\mu^{2}}}}{\sqrt{1+b\mu^{2}}}.\end{split}$
(B22)

Therefore, the relevant current can be expressed as

$\begin{split}\delta j_{y,e,3}&=-i\frac{c\delta
B_{x}}{4B_{0}^{2}}\frac{k_{\perp}^{2}}{k_{\parallel}}P_{\parallel,e}\left[4\eta^{\lambda}-1-\frac{\eta^{1/2+\lambda}}{\epsilon_{\perp}^{2}K_{2}(1/\epsilon_{\perp})}\left(\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}e^{-a\sqrt{1+b}}+\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}ab\int_{0}^{1}d\mu\frac{e^{-a\sqrt{1+b\mu^{2}}}}{\sqrt{1+b\mu^{2}}}\right)\right]\\\
&\equiv-i\frac{c\delta
B_{x}}{4B_{0}^{2}}\frac{k_{\perp}^{2}}{k_{\parallel}}P_{\parallel,e}\left[4\eta^{\lambda}-1-\mathcal{I}\right]\end{split}$
(B23)

which defines the integral $\mathcal{I}$.

Therefore, the final dispersion relation is

$-i\omega\frac{k_{\perp}^{2}}{k_{\parallel}}\mathcal{J}=k^{2}v_{a,i}^{2}-k_{\perp}^{2}\frac{P_{\parallel,e}}{4\rho_{i}}(4\eta^{\lambda}-1-\mathcal{I}),$
(B24)

or, for $k_{\perp}\gg k_{\parallel}$

$-i\frac{\omega}{k_{\parallel}}\rho_{i}\mathcal{J}=\frac{B_{0}^{2}}{4\pi}-P_{\parallel,e}(4\eta^{\lambda}-1-\mathcal{I})/4.$
(B25)

The growth rate of the instability is positive when

$2\beta_{\parallel,e}^{-1}-(4\eta^{\lambda}-1-\mathcal{I})/4>0.$ (B26)

Therefore, the threshold can be expressed as

$\beta_{\perp,e}<\frac{8\eta^{\lambda}}{4\eta^{\lambda}-1-\mathcal{I}},$ (B27)

where

$\mathcal{I}=\frac{\eta^{1/2+\lambda}}{\epsilon_{\perp}^{2}K_{2}(1/\epsilon_{\perp})}\left(\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}e^{-a\sqrt{1+b}}+\int_{1}^{\infty}d\gamma(\gamma^{2}-1)^{3/2}ab\int_{0}^{1}d\mu\frac{e^{-a\sqrt{1+b\mu^{2}}}}{\sqrt{1+b\mu^{2}}}\right).$
(B28)

In the ultra-relativistic limit, $\sqrt{\gamma^{2}-1}\approx\gamma$ and
$b\approx\eta-1$, $\mathcal{I}$ reduces to

$\begin{split}\mathcal{I}&=\frac{\eta^{1/2+\lambda}}{\epsilon_{\perp}^{2}K_{2}(1/\epsilon_{\perp})}\left(\int_{1}^{\infty}d\gamma\gamma^{3}e^{-\frac{\gamma}{\epsilon_{\perp}}\sqrt{\eta}}+\frac{\eta-1}{\epsilon_{\perp}}\int_{0}^{1}\frac{d\mu}{\sqrt{1+(\eta-1)\mu^{2}}}\int_{1}^{\infty}d\gamma\gamma^{4}e^{-\frac{\gamma}{\epsilon_{\perp}}\sqrt{1+(\eta-1)\mu^{2}}}\right)\\\
&=\frac{3\eta^{1/2+\lambda}\epsilon_{\perp}^{2}}{K_{2}(1/\epsilon_{\perp})}\left(\frac{2}{\eta^{2}}+\left[3-\frac{1}{\eta}-\frac{2}{\eta^{2}}+3\sqrt{\eta-1}\tan^{-1}\sqrt{\eta-1}\right]\right)\rightarrow\frac{9\pi\eta^{1+\lambda}\epsilon_{\perp}^{2}}{2K_{2}(1/\epsilon_{\perp})},\
\eta\rightarrow\infty,\end{split}$ (B29)

where the boundary term (first term in Eq. B28) cancels. Therefore, the mirror
instability threshold in the ultra-relativistic limit is

$\beta_{\parallel,e}<\frac{8}{4\eta^{\lambda}-1-1.5\eta^{1/2+\lambda}\left[3+3\sqrt{\eta-1}\tan^{-1}\sqrt{\eta-1}-1/\eta\right]}.$
(B30)

In the non-relativistic limit, the instability threshold defined by Equations
B27 and B28 reduces to $\beta_{\perp,e}<1/(\eta-1)$, consistent with previous
work. This is because for a non-relativistic distribution
$P_{\perp,e}/P_{\parallel,e}=\eta^{1}$, i.e., $\lambda=1$. The contribution to
$\delta j_{y,e,3}$ from the boundary term in $I_{1}$ (first term in Eq. B28)
is $3P_{\perp,e}/\eta^{2}$ and the contribution from $I_{2}$ (second term in
Eq. B28) is $P_{\perp,e}(\eta-1)(8\eta^{2}+4\eta+3)/\eta^{2}$. Thus,
$4\eta-1-\mathcal{I}=8\eta(\eta-1)$ in a non-relativistic limit.

In the case of non-relativistic anisotropic ions with anisotropy parameter
$\eta_{i}$, the same derivation leads to a threshold

$\beta_{\perp,e}<\frac{8\eta^{\lambda}}{4\eta^{\lambda}-1-\mathcal{I}+8\frac{P_{\parallel,i}}{P_{\parallel,e}}\eta_{i}(\eta_{i}-1)}.$
(B31)

Therefore, the threshold condition is defined by the ions as
$\beta_{\perp,i}<1/(\eta_{i}-1)$ when $P_{\parallel,i}\gg P_{\parallel,e}$.

The assumption of a zero parallel current $j_{z}$ holds only when
plasma-$\beta$ of at least one species is $\ll 1$. Consequently, when this
assumption is invalid, a non-zero $\delta E_{z}$ also leads to a more complex
$j_{y}$. In this case, both $j_{y}$ and $j_{z}$ contain terms proportional to
$\delta E_{y}$ and $\delta E_{z}$ through Vlasov’s equation. In a non-
relativistic limit, the threshold is modified by an additional stabilizing
term, which depends on plasma-$\beta$ of all species [see, e.g., Hall (1979);
Hellinger (2007)]. A similar stabilizing relation can be obtained in a
relativistic limit. The dispersion relation in the long-wavelength limit with
$\Omega_{e}\gg\omega$ and $\Omega_{e}\gg kv$ is then defined by the following
two inseparable equations:

$\frac{4\pi}{c^{2}}\omega\delta
j_{y}=i\left(\frac{\omega^{2}}{c^{2}}-k^{2}\right)\delta E_{y},\
\frac{4\pi}{c^{2}}\omega\delta
j_{z}=i\left(\frac{\omega^{2}}{c^{2}}-k_{\perp}^{2}\right)\delta E_{z},$ (B32)

which is equivalent to writing it in terms of plasma dielectric tensor
$\mathcal{E}_{\alpha\beta}$:

$\mathcal{E}_{22}-\frac{\mathcal{E}_{23}\mathcal{E}_{32}}{\mathcal{E}_{33}-\frac{k_{\perp}^{2}c^{2}}{\omega^{2}}}=\frac{k^{2}c^{2}}{\omega^{2}},$
(B33)

where $\mathcal{E}_{22}$ has already been calculated as the current response
along $\hat{y}$ due to $\delta E_{y}$:

$\mathcal{E}_{22}=1-\frac{ik_{\perp}^{2}}{\omega
k_{\parallel}^{2}\sigma_{i}}\mathcal{J}+\frac{\pi
c^{2}k_{\perp}^{2}}{B_{0}^{2}\omega^{2}}P_{\parallel,e}(4\eta-1-\mathcal{I}).$
(B34)

The general relations for relevant dielectric tensor components are:

$\begin{split}\mathcal{E}^{s}_{23}&=-\mathcal{E}^{s}_{32}=\frac{4\pi
iq_{s}^{2}}{\omega}\int
p_{\perp}v_{\perp}v_{\parallel}dp_{\perp}dp_{\parallel}\sum_{n=-\infty}^{+\infty}\frac{1}{\omega-
k_{\parallel}v_{\parallel}-n\Omega_{s}}\left[\frac{1}{v_{\perp}}\frac{\partial
f_{s}}{\partial p_{\perp}}+\frac{k_{\parallel}}{\omega}\left(\frac{\partial
f_{s}}{\partial p_{\parallel}}-\frac{v_{\parallel}}{v_{\perp}}\frac{\partial
f_{s}}{\partial
p_{\perp}}\right)\right]J_{n}\left(\frac{k_{\perp}v_{\perp}}{\Omega_{s}}\right)J^{\prime}_{n}\left(\frac{k_{\perp}v_{\perp}}{\Omega_{s}}\right)\\\
\mathcal{E}^{s}_{33}&=\frac{4\pi q_{s}^{2}}{\omega}\int
p_{\perp}v_{\parallel}dp_{\perp}dp_{\parallel}\sum_{n=-\infty}^{+\infty}\frac{1}{\omega-
k_{\parallel}v_{\parallel}-n\Omega_{s}}\left[\frac{\partial f_{s}}{\partial
p_{\parallel}}-\frac{n\Omega_{s}}{\omega}\left(\frac{\partial f_{s}}{\partial
p_{\parallel}}-\frac{v_{\parallel}}{v_{\perp}}\frac{\partial f_{s}}{\partial
p_{\perp}}\right)\right]J_{n}^{2}\left(\frac{k_{\perp}v_{\perp}}{\Omega_{s}}\right),\end{split}$
(B35)

where
$\mathcal{E}_{\alpha\beta}=\delta_{\alpha\beta}+\sum_{s}\mathcal{E}^{s}_{\alpha\beta}$.
Keeping the leading terms of order $\Omega_{s}^{-1}$ for $\delta j_{z,\delta
E_{y}}$ and $j_{y,\delta E_{z}}$ and terms of order $\Omega_{s}^{0}$ for
$\delta j_{z,\delta E_{z}}$

$\begin{split}\mathcal{E}^{s}_{23}=-\mathcal{E}^{s}_{32}&=2i\pi\frac{k_{\perp}}{\omega^{2}}\frac{q_{s}cm_{s}}{B_{0}}\int
dp_{\perp}dp_{\parallel}\frac{v_{\parallel}p^{2}_{\perp}}{\omega-
k_{\parallel}v_{\parallel}}\left[k_{\parallel}v_{\perp}\frac{\partial
f_{s}}{\partial p_{\parallel}}+(\omega-
k_{\parallel}v_{\parallel})\frac{\partial f_{s}}{\partial p_{\perp}}\right]\\\
&=-2i\pi\frac{k_{\perp}k_{\parallel}}{\omega^{2}}\frac{q_{s}cm_{s}}{B_{0}}\int
dp_{\perp}dp_{\parallel}p_{\perp}^{2}f_{s}\frac{\partial}{\partial
p_{\parallel}}\left(\frac{v_{\parallel}v_{\perp}}{\omega-
k_{\parallel}v_{\parallel}}\right)\\\ &=-2i\pi\frac{k_{\perp}}{\omega
m_{s}k_{\parallel}}\frac{q_{s}c}{B_{0}}\int
dp_{\perp}dp_{\parallel}\frac{f_{s}p_{\perp}^{3}}{\gamma^{4}(\omega/k_{\parallel}-v_{\parallel})^{2}}\left(\gamma^{2}-2\frac{p_{\parallel}^{2}}{m_{s}^{2}c^{2}}+\frac{k_{\parallel}v_{\parallel}}{\omega}\frac{p_{\parallel}^{2}}{m_{s}^{2}c^{2}}\right),\\\
\mathcal{E}^{s}_{33}&=\frac{4\pi q_{s}^{2}}{\omega}\int
dp_{\perp}dp_{\parallel}\frac{p_{\perp}v_{\parallel}}{\omega-
k_{\parallel}v_{\parallel}}\frac{\partial f_{s}}{\partial
p_{\parallel}}=-\frac{4\pi e^{2}}{m_{s}}\int
dp_{\perp}dp_{\parallel}\frac{p_{\perp}f_{s}}{\gamma^{2}(\omega-
k_{\parallel}v_{\parallel})^{2}}\left(\gamma-
p_{\parallel}\frac{\partial\gamma}{\partial p_{\parallel}}\right)\\\
&=-\frac{4\pi q_{s}^{2}}{k_{\parallel}^{2}m_{s}}\int
dp_{\perp}dp_{\parallel}\frac{p_{\perp}f_{s}}{\gamma^{3}(\omega/k_{\parallel}-v_{\parallel})^{2}}\left(1+\frac{p_{\perp}^{2}}{m^{2}c^{2}}\right).\end{split}$
(B36)

In a non-relativistic limit for an isotropic distribution, this reduces to:

$\begin{split}\mathcal{E}^{s}_{23}&=-\frac{2\pi k_{\perp}}{\omega
k_{\parallel}}\frac{q_{s}c}{B_{0}}\int\frac{dp_{\perp}dp_{\parallel}p_{\perp}^{2}v_{\perp}f_{s}}{(\omega/k_{\parallel}-v_{\parallel})^{2}},\\\
\mathcal{E}^{s}_{33}&=-\frac{4\pi
e^{2}}{m_{s}k_{\parallel}^{2}}\int\frac{dp_{\perp}dp_{\parallel}p_{\perp}f_{s}}{(\omega/k_{\parallel}-v_{\parallel})^{2}}\approx\frac{\omega_{p,s}^{2}}{v_{th,s}^{2}k_{\parallel}^{2}}\left(1+i\sqrt{\frac{\pi}{2}}\frac{\omega}{k_{\parallel}v_{th,s}}\right).\end{split}$
(B37)

In the dispersion relation, the imaginary term in $\mathcal{E}_{33}^{s}$ above
will group with other imaginary terms in $\mathcal{E}_{22}$ as a coefficient
in front of the growth rate $-i\omega$. Thus, the threshold condition,
considering $k_{\parallel}v_{\parallel}/\omega\gg 1$ and $k_{\perp}\gg
k_{\parallel}$, to leading order in $\Omega_{s}$, is

$\beta_{e,\parallel}\frac{4\eta-1-\mathcal{I}}{8}+\beta_{\parallel,i}\eta_{i}(\eta_{i}-1)>1-\frac{\pi}{B_{0}^{2}}\frac{\left(\sum_{s}\frac{q_{s}}{m_{s}}\int
p_{\perp}dp_{\perp}dp_{\parallel}\frac{f_{s}p_{\perp}^{2}(\gamma^{2}-2p_{\parallel}^{2}/m_{s}^{2}c^{2})}{\gamma^{4}(\omega/k_{\parallel}-v_{\parallel})^{2}}\right)^{2}}{\sum_{s}\frac{q_{s}^{2}}{m_{s}}\int
p_{\perp}dp_{\perp}dp_{\parallel}\frac{f_{s}(1+p_{\perp}^{2}/m_{s}^{2}c^{2})}{\gamma^{3}(\omega/k_{\parallel}-v_{\parallel})^{2}}}.$
(B38)

Figure 14: (a): Relativistic mirror instability thresholds (anisotropy $\eta$
as a function of $\beta_{\parallel,e}$) defined by Equations B27-B28
calculated for different $\epsilon_{\perp}$ represented by different colors
from darkest ($\epsilon_{\perp}=0.01$) to brightest ($\epsilon_{\perp}=100$).
Dotted black line represents a non-relativistic limit. (b): Ratio of
relativistic mirror instability threshold value $\beta_{\parallel,e}$ and its
non-relativistic limit $\beta_{\parallel,e,{\rm nr}}$ as a function of $\eta$.
(c): Ratio of relativistic mirror instability threshold value
$\beta_{\parallel,e}$ and its ultra-relativistic limit
$\beta_{\parallel,e,{\rm ur}}$, defined by Equation B30 as a function of
$\eta$.

In Figure 14 we show a comparison of the non-relativistic mirror threshold
$\beta_{\parallel,e,{\rm nr}}=1/(\eta-1)$ (thick black line) and the
numerically calculated relativistic electron mirror threshold from Equations
B27 and B28. In panel (a), $\eta$ is shown as a function of
$\beta_{\parallel,e}$ at different temperatures. Panel (b) shows the ratio of
the relativistic threshold value $\beta_{\parallel,e}$ and the non-
relativistic threshold $\beta_{\parallel,e,{\rm nr}}$ as a function of $\eta$.
Deviations are small for $\epsilon_{\perp}\lesssim 0.1$ and small values of
anisotropy parameter $\eta\lesssim 10$. At high temperatures, the numerical
solution is well approximated by the ultra-relativistic limit
$\beta_{\parallel,e,{\rm ur}}$ given by Equation B30. We show their ratio,
$\beta_{\parallel,e,{\rm ur}}/\beta_{\parallel,e}$, as a function of $\eta$ in
panel (c). Due to the large uncertainties in the electron anisotropy in
accretion flows and since we limit the anisotropy
$T_{\parallel,e}/T_{\perp,e}$ to be $\leq 10$ where the non-relativistic and
relativistic mirror thresholds are similar, we chose to use the analytically
simple non-relativistic mirror threshold for our application to BH images in §
3.

###  B.2 Parallel firehose instability

To calculate the relativistic firehose threshold, we linearize Equations B3
for ${\bf k}=k\hat{z}$, $\delta{\bf B}=\delta B_{y}\hat{y}$, and $\delta{\bf
E}=\delta E_{x}\hat{x}$ [see, e.g., Barnes & Scargle (1973)]. The resulting
equations for the ion and relativistic electron’s currents are:

$\begin{split}\delta j_{x,e}=-\frac{i\pi\delta E_{x}}{B_{0}^{2}c^{2}}\int
dpp^{3}vd\mu S(\mu)(kc\mu-\omega)\left(p\frac{\partial f_{e}}{\partial
p}-\mu\frac{\partial f_{e}}{\partial\mu}+\frac{kc}{\omega}\frac{\partial
f_{e}}{\partial\mu}\right),\\\ \delta j_{x,i}=-\frac{i\pi\delta
E_{x}}{B_{0}^{2}}\int dpd\mu
p^{3}\sqrt{m_{i}^{2}+(p/c)^{2}}S(\mu)(kv\mu-\omega)\frac{\partial
f_{i}}{\partial p}=-4i\omega\pi m_{i}c^{2}\frac{\delta E_{x}}{B_{0}^{2}}\int
dpp^{2}\frac{4(p/m_{i}c)^{2}/3+1}{\sqrt{(p/m_{i}c)^{2}+1}}f_{i}(p).\end{split}$
(B39)

Solving for sub-relativistic isotropic ions with distribution function B10 and
dropping the terms with odd $\mu$-integrands in $\delta j_{x,e}$ we find:

$\delta j_{x}=-\frac{i\omega\delta
E_{x}}{4\pi\sigma_{i}}\left[1+\frac{5}{2}\epsilon_{i}+\frac{\pi}{nm_{i}c}\int
dpd\mu
p^{3}(3-\mu^{2})f_{e}(p,\mu)+\frac{1}{nm_{i}c^{2}}(P_{\parallel,e}-P_{\perp,e})\frac{c^{2}k^{2}}{\omega^{2}}\right].$
(B40)

This results in the following dispersion relation

$\begin{split}\omega^{2}=k^{2}c^{2}\frac{\sigma_{i}-(P_{\parallel,e}-P_{\perp,e})/(nm_{i}c^{2})}{\mathcal{F}},\\\
\mathcal{F}=\sigma_{i}+1+\frac{5}{2}\epsilon_{i}+\frac{2\pi}{nm_{i}c}\int
dpd\mu p^{3}f_{e}(p,\mu)+\frac{P_{\perp,e}}{nm_{i}c^{2}}>0,\end{split}$ (B41)

which gives the usual non-relativistic firehose threshold

$\frac{P_{\perp,e}}{P_{\parallel,e}}<1-\frac{2}{\beta_{\parallel,e}},$ (B42)

where $\beta_{\parallel,e}=8\pi P_{\parallel,e}/B_{0}^{2}$. Note that if the
ions are anisotropic as well the relevant firehose threshold becomes
$P_{\perp,e}+P_{\perp,i}<P_{\parallel,e}+P_{\parallel,i}-B^{2}/4\pi$. Thus if
$P_{\perp,i}>P_{\perp,e}$ the ion anisotropy will in general be more important
than the electron anisotropy in setting stability to the fluid firehose
instability.

The calculation presented here focuses on the fluid parallel firehose
instability. There are also resonant parallel and oblique firehose
instabilities: Larmor-scale resonant instabilities destabilized by cyclotron
interaction. The resonant instabilities typically have faster growth rates and
somewhat lower anisotropy thresholds than the fluid firehose instability (Gary
et al., 1998; Hellinger & Matsumoto, 2000). Calculations of electron-scale
resonant firehose instabilities for relativistically hot electrons with
$T_{p}\gtrsim T_{e}$ would be valuable but we leave this to future work.

###  B.3 Whistler instability

The electron whistler instability, first noted in Sudan (1963) and followed by
a relativistic derivation (Sudan, 1965), is an instability of circularly
polarized electron waves propagating along the magnetic field direction
$B_{0}\hat{z}$. Considering a wavevector $k$ and fluctuating electric field
$\delta E_{x}$ and $\delta E_{y}$, the dispersion relation can be written as
(Gladd, 1983)

$\frac{\epsilon_{\parallel}k^{2}c^{2}}{\beta_{\parallel,e}}-\frac{\epsilon_{\perp}\omega^{2}}{\beta_{\perp,e}}+\pi
m_{e}^{2}c^{4}\Omega_{e,0}^{2}\int\frac{p_{\perp}^{2}v_{\perp}dp_{\perp}dp_{\parallel}}{kv_{\parallel}-(\omega-\Omega_{e})}\left[\frac{\partial
f_{e}}{\partial p_{\perp}^{2}}(\omega-kv_{\parallel})+\frac{\partial
f}{\partial p_{\parallel}^{2}}kv_{\parallel}\right]=0,$ (B43)

where $\epsilon_{\parallel}=T_{\parallel,e}/m_{e}c^{2}$, $f_{e}$ is defined by
Equation 7, $p_{\parallel}$ and $p_{\perp}$ are relativistic parallel and
perpendicular momentum, respectively, $\Omega_{e}$ and $\Omega_{e,0}$ are the
relativistic and non-relativistic electron cyclotron frequencies. The whistler
instability, like the ion cyclotron instability and unlike the mirror and
firehose instabilities considered in the previous section, typically does not
have a formal threshold but the growth rate becomes negligible for decreasing
anisotropy. This dispersion relation can thus be solved numerically to find
the target growth rate for a fixed $\beta_{\perp,e}$ and varying $\eta$.

The threshold for the relativistic whistler instability can be parameterized
as (Lynn, 2014)
$P_{\perp,e}/P_{\parallel,e}=1+S(\epsilon_{\perp})/\beta_{\perp,e}^{\alpha}$,
where $S(\epsilon_{\perp})=0.265-0.165(1+\epsilon_{\perp}^{-1})$ and
$\alpha=0.58-0.043\log\Gamma$ , where $\Gamma\sim 10^{-6}|eB/m_{e}c|$ is the
assumed growth rate. Since $S(\epsilon_{\perp})$ is a slowly varying and
monotonic function of temperature, $S(\epsilon_{\perp})\approx[0.1-0.25]$ for
$\epsilon_{\perp}=[10^{-2},10^{2}]$, we choose to use
$S(\epsilon_{\perp}=1)=0.183$.

##  Appendix C Anisotropy model for GR radiative transfer

Figure 15: Numerically calculated $\lambda$ as a function of
$\epsilon_{\perp}$ for $T_{\perp,e}/T_{\parallel,e}=\eta^{\lambda}$ for an
anisotropic relativistic bi-Maxwellian distribution function. A simple fit for
$\lambda$ is given in Equation C1.

For our choice of the distribution function (Eq. 7), the ratio of
perpendicular and parallel temperatures
$T_{\perp,e}/T_{\parallel,e}=\eta^{\lambda}$. The value of $\lambda$ is in
turn a function of temperature $\epsilon_{\perp}$, which we show in Figure 15.
The function $\lambda(\epsilon_{\perp}$) can be well-approximated by

$\lambda=-0.08\tanh{\left(1.5(\log_{10}\epsilon_{\perp}+0.5)\right)}+0.92,$
(C1)

In the non-relativistic limit, when $\epsilon_{\perp}\ll 1$, this gives
$T_{\perp,e}/T_{\parallel,e}\approx\eta$, while in the ultra-relativistic
limit, $T_{\perp,e}/T_{\parallel,e}\approx\eta^{0.8}$.

In our modeling of black hole accretion images we consider three limiting
cases for the anisotropy of the distribution function
$T_{\perp,e}/T_{\parallel,e}$, intended to bracket the magnitude of the effect
that an anisotropic distribution function can introduce:

$\begin{split}&\left(T_{\perp,e}/T_{\parallel,e}\right)_{\rm
mirror}=1+1/\beta_{\perp,e},\\\ &\left(T_{\perp,e}/T_{\parallel,e}\right)_{\rm
whistler}=1+S/\beta_{\perp,e}^{\alpha},\\\
&\left(T_{\perp,e}/T_{\parallel,e}\right)_{\rm isotropic}\equiv 1,\\\
&\left(T_{\perp,e}/T_{\parallel,e}\right)_{\rm
firehose}=1-2/\beta_{\parallel,e},\end{split}$ (C2)

where we take $S=0.183$ and $\alpha=0.838$ as in Appendix §B.3. Since the
firehose threshold is undefined at small electron-$\beta_{\parallel,e}$, we
choose to set the threshold to a constant value of
$T_{\perp,e}/T_{\parallel,e}=0.1$ at low $\beta_{\parallel,e}$. This is
motivated by local simulations (Riquelme et al., 2015). Likewise, for the
mirror instability, we limit $T_{\perp,e}/T_{\parallel,e}<10$. In reality, the
temperature anisotropy at low $\beta$ will depend on the heating, expansion
and contraction of the plasma, which is what drives the temperature anisotropy
in the first place.

It is useful to re-express Equations C2 in terms of the total electron
temperature

$T_{e}=\frac{1}{3}(2T_{\perp,e}+T_{\parallel,e}).$ (C3)

Using Equation C3, Equations C2 for the instability thresholds can be
rewritten as

$\begin{split}\beta_{\perp,e,{\rm
mirror}}&=\frac{\beta_{e}}{2}-\frac{1}{3}+\frac{1}{2}\sqrt{\frac{4}{9}+\frac{8}{3}\beta_{e}+\beta_{e}^{2}},\\\
\beta_{\parallel,e,{\rm firehose}}&=\beta_{e}+\frac{4}{3},\\\
\beta_{\perp,e,{\rm
whistler}}&=\begin{cases}0.141\beta_{e}^{3}-0.26\beta_{e}^{2}+1.171\beta_{e}+0.005,\beta_{e}<1\\\
1.054\beta_{e}+0.012,\beta_{e}>1\\\ \end{cases},\end{split}$ (C4)

where $\beta_{\parallel,e,{\rm firehose}}$ and $\beta_{\parallel,e,{\rm
mirror}}$ are exact solutions and $\beta_{\perp,e,{\rm whistler}}$ is a
polynomial fit to a numerical solution with growth rate $\Gamma\sim
10^{-6}|eB/m_{e}c|$ and $\epsilon_{\perp}=1$. The thresholds in Equation C4
can then be used in Equation C2, thus providing expressions for the threshold
temperature anisotropy in terms of $\beta_{e}$. This is a variable accessible
to a simulation that does not evolve temperature anisotropy, such as those
that we used in § 3. The threshold conditions are shown in Fig. 16 as a
function of electron $\beta_{e}$, which is extracted from the MHD
plasma-$\beta_{\rm th}$ via $\beta_{e}=2\beta_{\rm th}/(R+1)$.

Figure 16: Models used in the radiative transfer of GRMHD simulation,
represented by different colors: mirror and whistler
($T_{\perp,e}>T_{\parallel,e}$) and firehose ($T_{\perp,e}<T_{\parallel,e}$).

j

##  Appendix D GRMHD simulations

The GRMHD simulations used in §3 were performed using the publicly available
code Athena++ in spherical Kerr-Schild coordinates with a logarithmically
stretched grid in the radial direction $r$. The setup is identical to White et
al. (2019) with the outer radius of $1000r_{g}$ and the inner radius being
inside the horizon. The grid is refined with the level $0$ grid being
$N_{r}\times N_{\xi}\times N_{\phi}=64\times 32\times 64$. A total of 3
refinement levels are concentrated around the midplane, $\theta=\pi/2$,
resulting in an effective resolution of $512\times 256\times 512$ in $r$,
$\theta$, and $\phi$. We initialize a Fishbone torus Fishbone & Moncrief
(1976) with a purely poloidal magnetic field with mean plasma-$\beta$ of
$100$. We study two different spin values of the BH: $a=0.98$ and $0.5$. Each
of the two simulations is run up to a steady state and several eruption events
for a total simulation time of more than $15000r_{g}/c$.

Figure 17 shows the time evolution of the accretion rate $\dot{M}$ in code
units (a), magnetic flux $\Phi=0.5\int d\theta d\phi\sqrt{-4\pi g}|B^{r}|$
though a hemisphere (b), and dimensionless magnetic flux $\phi_{\rm
BH}=\Phi/\sqrt{\dot{M}r_{g}^{2}c}$ (c) measured at $2r_{g}$ as functions of
time, starting from $8000r_{g}/c$. Here $g$ is the determinant of spherical
Kerr-Schild metric. Spins of $0.5$ and $0.98$ are shown by blue and black
lines respectively. The time periods chosen for the GR radiative transfer in §
3 (shown by shaded blue and grey regions for $a=0.5$ and $a=0.98$
respectively) are such that the accretion rate is almost constant and no
magnetic flux eruptions occur. We have also performed the same analysis for
different quiescent time periods and found no qualitative difference in the
obtained results. The time interval we use to calculate average images is
relatively short but we do not analyze the time variability properties of our
results so this modest time interval is sufficient for our purposes.

Figure 17: Time evolution of accretion rate $\dot{M}$ (a), magnetic flux
$\Phi$ (b), and normalized magnetic flux $\Phi_{\rm
BH}/\sqrt{\dot{M}r_{g}^{2}c}$ after $8000r_{g}/c$ measured at $2r_{g}$ for two
spin parameters: $a=0.98$ (black) and $a=0.5$ (blue). Shaded regions
correspond to quiescent periods of $1000r_{g}/c$ chosen for GR radiative
transfer analysis: $14300-15300r_{g}/c$ for $a=0.98$ and $16000-17000r_{g}/c$
for $a=0.5$.

## References

  * Bale et al. (2009) Bale, S. D., Kasper, J. C., Howes, G. G., et al. 2009, Phys. Rev. Lett., 103, 211101, doi: 10.1103/PhysRevLett.103.211101
  * Barnes & Scargle (1973) Barnes, A., & Scargle, J. D. 1973, ApJ, 184, 251, doi: 10.1086/152324
  * Chael et al. (2023a) Chael, A., Issaoun, S., Pesce, D. W., et al. 2023a, ApJ, 945, 40, doi: 10.3847/1538-4357/acb7e4
  * Chael et al. (2021) Chael, A., Johnson, M. D., & Lupsasca, A. 2021, ApJ, 918, 6, doi: 10.3847/1538-4357/ac09ee
  * Chael et al. (2023b) Chael, A., Lupsasca, A., Wong, G. N., & Quataert, E. 2023b, arXiv e-prints, arXiv:2307.06372, doi: 10.48550/arXiv.2307.06372
  * Chandrasekhar et al. (1958) Chandrasekhar, S., Kaufman, A. N., & Watson, K. M. 1958, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 245, 435.  http://www.jstor.org/stable/100290
  * Cranmer et al. (1999) Cranmer, S. R., Kohl, J. L., Noci, G., et al. 1999, ApJ, 511, 481, doi: 10.1086/306675
  * Dexter (2016) Dexter, J. 2016, MNRAS, 462, 115, doi: 10.1093/mnras/stw1526
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a, ApJ, 875, L1, doi: 10.3847/2041-8213/ab0ec7
  * Event Horizon Telescope Collaboration et al. (2019b) —. 2019b, ApJ, 875, L4, doi: 10.3847/2041-8213/ab0e85
  * Event Horizon Telescope Collaboration et al. (2019c) —. 2019c, ApJ, 875, L6, doi: 10.3847/2041-8213/ab1141
  * Event Horizon Telescope Collaboration et al. (2021) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021, ApJ, 910, L13, doi: 10.3847/2041-8213/abe4de
  * Event Horizon Telescope Collaboration et al. (2022a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022a, ApJ, 930, L12, doi: 10.3847/2041-8213/ac6674
  * Event Horizon Telescope Collaboration et al. (2022b) —. 2022b, ApJ, 930, L15, doi: 10.3847/2041-8213/ac6736
  * Fishbone & Moncrief (1976) Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962, doi: 10.1086/154565
  * Foucart et al. (2017) Foucart, F., Chandra, M., Gammie, C. F., Quataert, E., & Tchekhovskoy, A. 2017, MNRAS, 470, 2240, doi: 10.1093/mnras/stx1368
  * Galishnikova et al. (2023) Galishnikova, A., Philippov, A., Quataert, E., et al. 2023, Phys. Rev. Lett., 130, 115201, doi: 10.1103/PhysRevLett.130.115201
  * Gary (1992) Gary, S. P. 1992, Journal of Geophysical Research: Space Physics, 97, 8519, doi: https://doi.org/10.1029/92JA00299
  * Gary et al. (1998) Gary, S. P., Li, H., O’Rourke, S., & Winske, D. 1998, Journal of Geophysical Research: Space Physics, 103, 14567, doi: https://doi.org/10.1029/98JA01174
  * Ginzburg & Syrovatskii (1965) Ginzburg, V. L., & Syrovatskii, S. I. 1965, ARA&A, 3, 297, doi: 10.1146/annurev.aa.03.090165.001501
  * Gladd (1983) Gladd, N. T. 1983, The Physics of Fluids, 26, 974, doi: 10.1063/1.864249
  * Gold et al. (2020) Gold, R., Broderick, A. E., Younsi, Z., et al. 2020, ApJ, 897, 148, doi: 10.3847/1538-4357/ab96c6
  * Hall (1979) Hall, A. N. 1979, Journal of Plasma Physics, 21, 431–443, doi: 10.1017/S0022377800022005
  * Hellinger (2007) Hellinger, P. 2007, Physics of Plasmas, 14, 082105, doi: 10.1063/1.2768318
  * Hellinger & Matsumoto (2000) Hellinger, P., & Matsumoto, H. 2000, J. Geophys. Res., 105, 10519, doi: 10.1029/1999JA000297
  * Johnson et al. (2023) Johnson, M. D., Akiyama, K., Blackburn, L., et al. 2023, Galaxies, 11, 61, doi: 10.3390/galaxies11030061
  * Kulsrud (1983) Kulsrud, R. M. 1983, Handbook of plasma physics. Vol. 1: Basic plasma physics I. edited by Galeev, A. A. and Sudan, R. N.
  * Kunz et al. (2014) Kunz, M. W., Schekochihin, A. A., & Stone, J. M. 2014, Phys. Rev. Lett., 112, 205003, doi: 10.1103/PhysRevLett.112.205003
  * Leung et al. (2011) Leung, P. K., Gammie, C. F., & Noble, S. C. 2011, ApJ, 737, 21, doi: 10.1088/0004-637X/737/1/21
  * Lu et al. (2023) Lu, R.-S., Asada, K., Krichbaum, T. P., et al. 2023, Nature, 616, 686, doi: 10.1038/s41586-023-05843-w
  * Lynn (2014) Lynn, J. W. 2014, PhD thesis, University of California, Berkeley
  * Mahadevan et al. (1996) Mahadevan, R., Narayan, R., & Yi, I. 1996, ApJ, 465, 327, doi: 10.1086/177422
  * Melrose (1971) Melrose, D. B. 1971, Ap&SS, 12, 172, doi: 10.1007/BF00656148
  * Mertens et al. (2016) Mertens, F., Lobanov, A. P., Walker, R. C., & Hardee, P. E. 2016, A&A, 595, A54, doi: 10.1051/0004-6361/201628829
  * Mikhailovsky (1976) Mikhailovsky, A. B. 1976, Electromagetic instabilities of non-Maxwellian plasma, ed. M. A. Leontovich (English translation by Herbert Lashinsey. Plenum Publishing Corporation (Consultants Bureau), N.Y. Pp. xiv 331.)
  * Mościbrodzka et al. (2016) Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A, 586, A38, doi: 10.1051/0004-6361/201526630
  * Mościbrodzka & Gammie (2018) Mościbrodzka, M., & Gammie, C. F. 2018, MNRAS, 475, 43, doi: 10.1093/mnras/stx3162
  * Mościbrodzka et al. (2009) Mościbrodzka, M., Gammie, C. F., Dolence, J. C., Shiokawa, H., & Leung, P. K. 2009, ApJ, 706, 497, doi: 10.1088/0004-637X/706/1/497
  * Narayan et al. (2021) Narayan, R., Palumbo, D. C. M., Johnson, M. D., et al. 2021, ApJ, 912, 35, doi: 10.3847/1538-4357/abf117
  * Osipov et al. (2017) Osipov, S. M., Bykov, A. M., Ellison, D. C., & Gladilin, P. E. 2017, Journal of Physics: Conference Series, 929, 012006, doi: 10.1088/1742-6596/929/1/012006
  * Pandya et al. (2018) Pandya, A., Chandra, M., Joshi, A., & Gammie, C. F. 2018, ApJ, 868, 13, doi: 10.3847/1538-4357/aae7d5
  * Pandya et al. (2016) Pandya, A., Zhang, Z., Chandra, M., & Gammie, C. F. 2016, ApJ, 822, 34, doi: 10.3847/0004-637X/822/1/34
  * Quataert et al. (2002) Quataert, E., Dorland, W., & Hammett, G. W. 2002, The Astrophysical Journal, 577, 524, doi: 10.1086/342174
  * Quataert & Gruzinov (1999) Quataert, E., & Gruzinov, A. 1999, ApJ, 520, 248, doi: 10.1086/307423
  * Riquelme et al. (2015) Riquelme, M. A., Quataert, E., & Verscharen, D. 2015, ApJ, 800, 27, doi: 10.1088/0004-637X/800/1/27
  * Rosenbluth (1956) Rosenbluth, M. 1956.  https://www.osti.gov/biblio/4329910
  * Rudakov & Sagdeev (1961) Rudakov, L. I., & Sagdeev, R. Z. 1961, Soviet Physics Doklady, 6, 415
  * Sazonov (1972) Sazonov, V. N. 1972, Ap&SS, 19, 25, doi: 10.1007/BF00643165
  * Shcherbakov (2008) Shcherbakov, R. V. 2008, ApJ, 688, 695, doi: 10.1086/592326