# On the Comparison of AGN with GRMHD Simulations: II. M87

Richard Anantua${}^{1,2,3,4,5\kern 0.56001pt\href
https://orcid.org/0000-0003-3457-7660}$,  Angelo Ricarte${}^{2,3\kern
0.56001pt\href https://orcid.org/0000-0001-5287-0452}$, George
Wong${}^{6,7\kern 0.56001pt\href https://orcid.org/0000-0001-6952-2147}$,
Razieh Emami,${}^{3\kern 0.56001pt\href
https://orcid.org/0000-0002-2791-5011}$ Roger Blandford${}^{4\kern
0.56001pt\href https://orcid.org/0000-0002-1854-5506}$, Lani Oramas${}^{1\kern
0.56001pt\href https://orcid.org/0009-0006-4575-5126}$,  Hayley
West${}^{1\kern 0.56001pt\href https://orcid.org/0009-0003-1798-8406}$,
Joaquin Duran${}^{1\kern 0.56001pt\href
https://orcid.org/0009-0003-6244-0271}$ and Brandon Curd${}^{1,2,3\kern
0.56001pt\href https://orcid.org/0000-0002-8650-0879}$

1 Department of Physics $\&$ Astronomy, The University of Texas at San
Antonio, One UTSA Circle, San Antonio, TX 78249, USA
2 Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA
02138, USA
3 Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
4Kavli Institute for Particle Astrophysics and Cosmology, Stanford University
P.O. Box 20450, MS 29, Stanford, CA 94309, USA
5Astronomy Department, University of California, Berkeley, 601 Campbell Hall,
Berkeley, CA 94720, USA
6 School of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive,
Princeton, NJ 08540
7 Princeton Gravity Initiative, Princeton University, Princeton, New Jersey
08544, USA E-mail<EMAIL_ADDRESS>
(Accepted ????. Received ; in original form )

###### Abstract

Horizon-scale observations of the jetted active galactic nucleus M87 are
compared with simulations spanning a broad range of dissipation mechanisms and
plasma content in three-dimensional general relativistic flows around spinning
black holes. Observations of synchrotron radiation from radio to X-ray
frequencies can be compared with simulations by adding prescriptions
specifying the relativistic electron-plus-positron distribution function and
associated radiative transfer coefficients. A suite of time-varying
simulations with various spins and plasma magnetizations is chosen to
represent distinct possibilities for the M87 jet/accretion flow/black hole
(JAB) system. We then input turbulent heating and equipartition-based emission
prescriptions (and piecewise combinations thereof) in the time-dependent 3D
simulations, in which jet morphology, polarization and variation are
“observed” and compared with real observations so as to try to infer the rules
that govern the polarized emissivity. The models in this paper support a
magnetically arrested disk (MAD) with several possible spin/emission model
combinations supplying the jet in M87, whose inner jet and black hole shadow
have been observed down to the photon ring at 230 GHz by the Event Horizon
Telescope (EHT). We also show that some MAD cases that are dominated by
intrinsic circular polarization have near-linear $V/I$ dependence on unpaired
electron or positron content while SANE polarization exhibits markedly greater
positron-dependent Faraday effects – future probes of the SANE/MAD dichotomy
and plasma content with the EHT. This is the second work in a series also
applying the observing simulations methodology to near-horizon regions of
supermassive black holes in Sgr A* and 3C 279.

###### keywords:

active galactic nucleus (AGN), general relativistic magnetohydrodynamic
(GRMHD) simulation, relativistic jet, very long baseline interferometry (VLBI)

††pagerange: On the Comparison of AGN with GRMHD Simulations: II. M87–B.2
Positron Distribution from Geometry††pubyear: 2023

##  1 Introduction

Over a century ago, M87 was described as a “curious straight ray” by Heber
Curtis (Curtis, 1918) due to its relativistic jet, and nearly three quarters
of a century ago, it was identified as a discrete radio source (Bolton et al.,
1949). M87 is now the best studied jet/accretion flow/black hole (JAB) system,
and the first to be imaged down to the horizon scale by the Event Horizon
Telescope (Event Horizon Telescope Collaboration et al., 2019a, b, c, d, e,
f). Throughout the years in which the giant elliptical galaxy M87 has been
observed from its lobes to its core, we have learned that it is one of the
closest examples of a common physical phenomenon, the production of twin,
relativistic jets by accreting, spinning, massive black holes. In recent
years, our understanding of how jets form, propagate and radiate has advanced
considerably. Much of this progress can be credited to advances in
observational capability, throughout the electromagnetic spectrum. In
particular, the technique of VLBI (including polarimetry) has been extended to
higher frequency, where the angular resolution is finer and depolarization is
smaller. Gamma-ray observations have also contributed much. Equivalent
progress can be seen in the numerical simulation of non-axisymmetric, general
relativistic, hydromagnetic flows where specific models can be evolved
dynamically with numerical confidence. The challenge today is to reconcile
these two approaches.

This reconciliation needs to take place at several levels. Radio/mm/submm jets
can be imaged down to merely tens of gravitational radii (defined henceforth
to be $M\equiv GM_{\mathrm{H}}/c^{2}$, where $M_{\mathrm{H}}$ is the mass of
the hole). Structure has been identified down to $\sim 10M$ (Fish et al.,
2016)) and beyond by the Event Horizon Telescope (EHT) project, which has made
linearly polarized images with resolution limit $\sim 5M$ (Event Horizon
Telescope Collaboration et al., 2021a, b). Jet structure has now been
connected with the emitting ring in 3.5mm intensity maps (Lu et al., 2023).

The general relativistic regime has often been simulated under a variety of
dynamical assumptions and initial conditions. There has also been progress in
adding radiative transfer to these codes to take account of absorption,
scattering and Faraday rotation (Dexter et al., 2012). However, in order to
achieve the ultimate goal of elucidating jet launch and collimation (Lu et
al., 2014) and to measure the spin of the black hole, it is necessary to have
a better understanding of how high energy electrons are accelerated.

Most jets are observed on the larger scale, where general relativity is
unimportant. In addition to radio observations, optical and X-ray observations
extend down to $\sim~{}0.1$”, $\sim~{}1$” respectively. Long term monitoring,
using VLBI, has taught us much about the apparent motion of emission sites
within jets and accounting for this is a goal too. It is also necessary to
uncover the character of the medium through which the jet is propagating–
inflow, outflow or thick, orbiting torus; fluid or hydromagnetic– and the
interaction with it. In addition, it is necessary to understand the remarkably
rapid variability, notably at $\gamma$-ray energies, which cannot originate
very close to the black hole and where the jets are totally unresolved.
Answering these questions presents an even greater challenge to our current
understanding of particle acceleration.

Relativistic jets are associated with a large and heterogeneous sample of
sources and the distribution of their properties is largely determined by the
statistics of relativistic beaming. (We are primarily concerned with AGN here,
but the problems we are addressing are also features of gamma-ray bursts and
galactic superluminal sources, from which much can also be learned.) The
orientation of a specific source is a parameter which can be adjusted to match
a simulation to a particular source. However, we also know that black hole
spin axes are isotropically distributed in space. Furthermore, the fluxes and
images should scale with the jet powers and hole masses. All of this implies
specific contributions to the overall distributions of total fluxes, apparent
expansion speeds, polarization and so on in a complete sample of sources
selected according to well-defined criteria. The overall nonthermal radiative
efficiency can be determined observationally (Soltan, 1982). This, too,
relates strongly to the particle acceleration.

There have been many proposals as to how particles are accelerated under these
conditions. Strong shocks are commonly invoked, but these are not efficient
accelerators in magnetically-dominated flows and may be too slow to account
for many observations. Supersonic jets are surely very noisy, and the
associated hydromagnetic turbulence can promote second order acceleration. The
very existence of fast variability, especially at $\gamma$-ray energy,
suggests that unstable electromagnetic configurations lead to very rapid
particle acceleration where electromagnetic energy is efficiently converted
locally to high energy particles by a process we have called
“magnetoluminescence” (Blandford et al., 2017).

Clearly the program that has just been sketched is a massive undertaking and
it is premature to try to execute it in full. In this paper, we limit
ourselves to a smaller exercise designed to link plasma microphysics to
discrete observable AGN features. Observing JAB simulations has already been
carried out for Sagittarius A* in Anantua et al. (2020)b and Anantua et al.
(2023), where a turbulent heating model exponentially suppressing emission
from high-gas-to-magnetic pressure regions outperformed other distinct
phenomenological model classes with respect to image size and spectral
constraints anticipating EHT results (Event Horizon Telescope Collaboration et
al., 2022)111 The Critical Beta two-parameter model was not among the $6\%$ of
models in Sgr A* Paper I passing non-EHT constraints on 86 GHz flux– with the
important caveat that it was explored at a single point in model parameter
space for limited inclinations relative to the fiducial models.  . Here, we
consider a few simulations and apply it to another well-observed source, M87 –
with the added flexibility of comparing two turbulent heating prescriptions,
using a separate presription for jet regions, and including positrons. We will
pay most attention to varying the particle acceleration/emissivity
prescription and explore how to see which choices come closest to reproducing
the observations and calculating the underlying physical properties. There is
no reason to expect the match to be especially good right now. However, better
simulations and observations, both imminent, should allow this approach to be
followed more productively, now armed with an arsenal of distinguishable
models ranging from M87-like ordered electromagnetic jet and magnetosphere to
dense, Faraday-thick turbulent plasma more suitable for other AGN.

In Section 2, we review observations of the M87 hole, disk and jet,
emphasizing those that are most directly connected to the synchrotron
emissivity. Section 3 presents GRMHD simulations with varying spin and
accretion mode and describes commonalities and differences in their plasma
flow structures. In Section 4, we introduce self-consistent prescriptions for
emission (and absorption), particle acceleration and dissipation, including
the essential effects of positron physics. In Section 5, we describe the
global properties of our GRMHD simulations. In Section 6, we apply our
emission prescriptions with positrons to the time-dependent simulations to
“observe” M87. Our general conclusions and plans for further investigations
are collected in Section 7. Synchrotron radiation theory calculations for
alternative emission model prescriptions including positrons, are expounded in
the Appendix.

##  2 Observations of M87

Located at the heart of the Virgo Cluster $d=16.7\pm 0.6{\rm\ Mpc}$ away
(Blakeslee et al., 2009) (cf. Table 1), the bright active galaxy M87 (3C 274)
serves as an exemplary laboratory for the investigation of black hole jets.
Observations of the jet on all scales suggest that M87 is a FRI misaligned BL-
Lac blazar. M87 has a remarkably prominent jet, with an equally remarkably
faint disk centered on a large black hole.

###  2.1 Black Hole

We adopt a black hole mass of $(6.5\pm 0.7)\times 10^{9}{\rm M}_{\odot}$
(Event Horizon Telescope Collaboration et al., 2019a), corresponding to
length, time, angular and energy scales $10^{13}{\rm\ m},9{\rm\ hr},4\ \mu{\rm
as}\mathrm{\ and\ }1.2\times 10^{57}{\rm\ J}$, respectively. The associated
Eddington luminosity is $\sim 8\times 10^{40}{\rm\ W}$. A lower bound of
$a>0.2M_{\mathrm{M87}}$ has been derived for the spin of the hole (Li et al.,
2009), assuming M87’s SMBH is surrounded by a prograde, radiatively
inefficient accretion flow (RIAF). Hybrid jet/advection-dominated accretion
flow (ADAF) models Feng & Wu (2017) have provided an estimate as high as
$a=0.98M_{\mathrm{H}}$. We first adopt an intermediate angular frequency of
$\Omega_{\mathrm{H}}=0.35/(GM_{\mathrm{H}}/c^{3})$ of $10^{-5}{\rm s}^{-1}$
for M87’s black hole. Using the relations $J_{H}=GM_{\mathrm{H}}^{2}a/c$ and
$\Omega_{\mathrm{H}}=\frac{a}{r_{+}^{2}+a^{2}}$, where
$r_{+}=\frac{1}{2}(r_{S}+\sqrt{r_{S}^{2}-4(J_{H}/(M_{\mathrm{H}}c))^{2}})$ is
the radius of the outer Kerr horizon, the corresponding dimensionless black
hole spin is $a/M_{\mathrm{H}}=0.94$. We also consider lower spin prograde
$a/M=0.5$ ($\Omega_{\mathrm{H}}=0.13/(GM_{\mathrm{M87}}/c^{3})$) and
retrograde $a/M=-0.5$ cases. Observations of the jets on all scales suggest
that they, and by hypothesis, the spin of the hole, are inclined to the line
of sight at an angle $\theta=20^{\circ}$ (Wang & Zhou, 2009; Prieto et al.,
2016). The extractable rotational energy and angular momentum are $\sim
2\times 10^{56}{\rm J}$ and $\sim 4\times 10^{61}{\rm kg\ m}^{2}{\rm s}^{-1}$,
respectively – ample to power the jets observed today for a Hubble time
without any accretion. Henceforth, we measure all lengths, angles and times in
units $M$ set by M87’s black hole, i.e., $GM_{\mathrm{M87}}/c^{2}$,
$(GM_{\mathrm{M87}}/c^{2})/d$, $GM_{\mathrm{M87}}/c^{3}$, respectively (confer
Table 2). When the spin direction is along the general direction of the
angular velocity of the orbiting gas (Walsh et al., 2013), it is aligned with
the receding jet.

###  2.2 Radio-mm-submm Observations

Following the pioneering observations of Junor & Biretta (1995), there have
been many impressive high resolution observations of the inner jet of M87.

  * •

A 15 GHz map with beam $600\ \mu{\rm as}\times 1300\ \mu{\rm as}\equiv\
150M\times\ 325M$ measured out to a projected length $Y\sim 2\times 10^{4}M$
(Kovalev et al., 2007) and therefore a length along the deprojected jet $z\sim
6\times 10^{4}M$.

  * •

Pilot monitoring at 22 GHz with resolution $\sim 250M\times 250M$ extending
out to projected length $\sim 6000M$ and showing superluminal motion with
speed up to $\sim 1.6c$ (Hada et al., 2017)

  * •

A 43 GHz VLBI time sequence– 11 maps made over 210 d with a beam FWHM $\sim
55M\times 115M$, which extends out to $z\sim 6000M$ in projected radius
(Walker et al., 2008, 2016; Mertens et al., 2016) (Fig. 2). Apparent speeds up
to $\sim 2c$ are observed.

  * •

An 86 GHz image with resolution $\sim 20M\times 60M$ extending out to $z\sim
3000M$ and exhibiting $\sim 20$ percent linear polarization and strong Faraday
rotation (Hada et al., 2016), and a Global mm-VLBI Array (GMVA) 86 GHz M87
observation exhibiting a limb-brightened jet base (see Kim et al. (2018b) Fig.
4).

  * •

Event Horizon Telescope (EHT) observational data made with effective beam of
size $\sim 10M$ (Akiyama et al., 2015). Increased coverage with next
generation facilities may further bridge the gap between this jet structure
and the ringlike EHT 2019 observation around the central supermassive black
hole with greater precision, building off of Lu et al. (2023).

These observations confirm that the approaching radio jet is modestly
relativistic and collimated within $z\sim 30M$. They are also broadly
consistent with there being a self-absorbed innermost jet called the core with
constant brightness temperature $\sim 2\times 10^{10}{\rm\ K}$ and flux
density $S_{\nu}\sim 1{\rm\ Jy}$ for $3{\rm\ GHz}\lesssim\nu\lesssim 300{\rm\
GHz}$. At higher frequency, the entire jet appears to be optically thin. The
resolved jet structure accounts for a minority of the flux at all frequencies
but is strongly edge-brightened and the mean intensity decays as an inverse
square with distance from the axis. There are indications that the southern
edge is brighter than the northern edge. The jet is quite variable. The shape
of the jet is roughly parabolic for $30M\lesssim z\lesssim 7000M$ with the
separation of the brightened edges roughly given by $\sim 6z^{1/2}$. At larger
radii, the jet expansion is closer to linear.

####  2.2.1 Observational Model

We have made a simple analytical model which captures many of the features of
the time-averaged observation over the range of frequencies where the inner
jet has been resolved. While this does not do full justice to the
observations, it is sufficient for our purpose.

We introduce a Cartesian coordinate system on the sky with $Y$ measuring
distance along the jet from the hole and $X$ distance across it (in units of
$M$). The intensity satisfies

$I=\frac{I_{0}e^{\xi(1-\xi)}}{1+\nu_{300}^{8/3}Y^{2}},$ (1)

where $\xi=X^{2}/Y$.

Figure 1: VLBA 2 cm image of M87. The swirling jet substructure suggests
magnetic Kelvin-Helmholtz instability– a feature also seen in the simulation
movie of the corkscrewing jet http://richardanantua.com/sample-
page/jetaccretion-diskblack-hole-movies/. Image adapted with permission of Dan
Homan, Yuri Kovalev, Matt Lister and Ken Kellermann.

(2)

(3)

Figure 2: M87 VLA observation on linear intensity scale (Top) and
observational movie snapshots at $t_{\mathrm{Obs0}}$ (Bottom Left) side-by-
side with a snapshot at $t_{\mathrm{Obs0}}+10t_{\mathrm{Step}}$ (Bottom Right)
for $t_{\mathrm{Step}}=21\mathrm{days}\approx 56M_{\mathrm{M87}}$
monotonically transformed by $(\cdot)^{1/4}$ for visual clarity on the bottom
row. The beam dimensions are $4.3\cdot 10^{-4}\mathrm{arcs}\cdot 2.1\times
10^{-4}\mathrm{arcs}$. These images can be viewed as a movie sans
transformation, courtesy of Craig Walker and his collaborators, here:
http://www.aoc.nrao.edu/~cwalker/M87/.

The observed radio-mm-submm isotropic power is dominated by mm observations
and is $\sim 10^{34}$ W (Prieto et al., 2016). A beaming corrected guess might
be as high as $\sim 10^{35}$ W.

The spectral energy distribution of the central 0.4 arcs (32 pc) of M87 from
Prieto et al. (2016) is $\nu F_{\nu}\approx 10^{11.6}$ Jy Hz for $10^{10.5}$
Hz $<\nu<$ $10^{14.5}$ Hz. The radio-to-UV bolometric luminosity is $3.6\cdot
10^{-6}L_{\mathrm{Edd}}=2.7\cdot 10^{35}$ W (Prieto et al., 2016).

####  2.2.2 Event Horizon Telescope Observations

In April 2019, the Event Horizon Telescope released the first images resolving
the boundary of a black hole (Event Horizon Telescope Collaboration et al.,
2019a), ushering in the age of direct observation of horizons. The results
have already resolved a wide discrepancy in the black hole mass for M87*– from
stellar dynamical measurements of $6.6\times 10^{9}M_{\odot}$ from Gebhardt et
al. (2011), compared to gas dynamical measurements from Walsh et al. (2013)
who measured half this mass. Simulations concordant with EHT M87 observations
require that the central black hole have nonzero spin in order to explain the
presence of the jet powered by the Blandford-Znajek mechanism (Blandford &
Znajek, 1977), and polarized observations (Event Horizon Telescope
Collaboration et al., 2021a, b) indicating the hole is supplied vertical
magnetic flux further support this interpretation.

(4) Figure 3: EHT M87 2017 observational campaign: intensity map on (Left) and
linear polarization map (Right) (Event Horizon Telescope Collaboration et al.,
2021a).

###  2.3 Optical-Infrared Observations

The Hubble Space Telescope has provided us with stunning optical band
observations of the M87 jets, including knots with superluminal components and
flatter spectrum than the rest of the jet (Perlman et al., 2001). The most
famous feature, HST-1 $\sim 80\ \mathrm{pc}$ from the nucleus, produces blobs
that appear to move up to $6c$ on the observer plane (Biretta et al., 1999).
$\mathit{HST}$-1 has exhibited 40$\%$ variability between 1993 and 1997
(Perlman et al., 2001).

The isotropic radiant power of M87 has a bolumetric luminosity given by
$L_{\mathrm{bol}}\sim 3\times 10^{35}\mathrm{W}$ Prieto et al. (2016). Using
the quiescent spectral energy distribution from Prieto et al. (2016), it is
inferred that the upper limit to disk power is $L_{\mathrm{disk}}\leq 3.4\cdot
10^{41}$erg/s. At 10$\%$ efficiency, $L=\eta\dot{m}c^{2}$ implies an upper
bound to the mass accretion rate of $3.8\cdot 10^{21}$ g/s = $6\times
10^{-5}M_{\odot}/\mathrm{yr}$. The twofold change in M87’s observed bolometric
luminosity from its quiescent state value to ($L=5.4\times 10^{42}$ergs/s)
during its 2005 outburst Prieto et al. (2016) suggests a Doppler boosting
factor of 8-16 given the mass accretion rate upper limit.

###  2.4 X-ray Observations

M87 has been observed at X-ray wavelengths by $\mathit{Chandra}$ (Wilson &
Yang, 2002). The equipartition magnetic field value for the knot HST-1 was
found to be $\sim 3\times 10^{-4}$ G (Owen et al., 1989).

The steady jet isotropic luminosity in 2-10 keV X-rays is $\sim 3\times
10^{34}$ W (Prieto et al., 2016). However, the variable source HST-1 can be
roughly ten times brighter including observations at optical wavelengths. As
HST-1 also displays features moving towards us exhibiting apparent speeds
$\sim 6c$, we suppose that this is a small part of the flow that, unlike the
main body of the jet, is directed along our line of sight. It may not
contribute significantly to the true integrated jet bolometric power.

M87 resides at the center of the $\sim 10^{15}{\rm M}_{\odot}$ Virgo cluster
of galaxies, surrounding M87 in a cooling flow region with X-ray luminosity
$10^{36}$ W (Churazov et al., 2001). The cooling time is short compared with
the flow time, suggesting that the hot gas is maintained in rough dynamical
equilibrium by mechanical heating associated with the jets (though other
possibilities have been widely discussed). If this is the case, then each jet
must carry a total power that is significantly larger than this and which is
mostly carried off by buoyant bubbles. On this rather uncertain basis, we
estimate a total jet power of $L_{\rm jet}=5\times 10^{36}$ W, noting that M87
could be in a relatively dormant state right now.

###  2.5 Gamma-ray Observations

The $\mathit{Fermi}$ Large Area Telescope (LAT) has seen M87 as a gamma ray
point source with variable TeV emission on yearly timescales (Aharonian, 2006)
and no observed yearly or decadal variability at MeV-GeV ranges (Abdo et al.,
2009).

The observed $>100$ MeV flux is $2.54\cdot 10^{-8}$ photons cm-2 $s^{-1}$ and
the corresponding luminosity is $4.9\cdot 10^{41}$ erg $s^{-1}$ (Abdo et al.,
2009).

###  2.6 Galaxy and Cluster

The largest scale observations (Owen et al., 2000; de Gasperin et al., 2012;
Forman et al., 2017) show that the M87 jets interact strongly with the
surrounding medium. The jet orientation changes at a radius $\sim 2{\rm\
kpc}$. The jets have inflated large, buoyant bubbles that probably produce
enough dissipation to balance radiative cooling loss. M87 has been active but
quite variable for several Gyr, although the underlying mass supply rate may
have varied significantly over this time. Analyzing the most recent activity
leads to a conservative estimate of the current power per jet $\sim
10^{37}{\rm\ W}$.

###  2.7 Summary

M87 is fairly extreme in many respects. It has the most massive black hole we
can study in detail. Its disk luminosity is $\lesssim 3\times 10^{-7}L_{\rm
Edd}$ while it creates jets with total power $\gtrsim 3\times 10^{-4}L_{\rm
Edd}$ which have been shown to be collimated on scales $\lesssim 10M$. Most
importantly, we have recently begun to learn much more about the region within
$\sim 100M$ from EHT observations. This makes it an excellent source to model.

##  3 3D GRMHD Simulations

###  3.1 Historical Overview

We now give a brief synopsis of the development of GRMHD simulations to
contextualize the principal simulation used in this work.

Koide and Meier pioneered GRMHD simulations (Koide et al., 2000) of jet
formation in the magnetosphere of rapidly rotating Kerr black holes. Their
code solves the GRMHD equations for conservation of particle number, momentum
and energy in a Kerr metric for $0.75r_{S}<x^{1}<20r_{S}$ (where $x^{1}$ is
the radius $r$ in Boyer-Lindquist coordinates).

The next major advance was the high-accuracy relativistic magnetohydrodynamics
(HARM) code of Gammie McKinney and Toth (Gammie et al., 2003). This
conservative numerical scheme for integrating the GRMHD equations is
guaranteed to obey shock jump conditions at discontinuities in the fluid
variables.

HARM led to a number of applications, such as Dexter and Fragile’s disk
simulations of Sgr A* Dexter et al. (2012) and Mishra et al. (2016) and
McKinney & Blandford (2009) simulations of stable relativistic jets.

Later simulations by Farris and Gold account for strong gravitational
curvature near black hole binaries (Farris et al., 2012) or describe models
with magnetorotational (MRI) disk instability and magnetically arrested disks
(MAD) (Gold et al., 2016). The simulations in Gold et al. (2016) that model
the disk – which is governed by distinct emission mechanisms from the jet –
require an evolution equation for proton temperature $T_{p}$ and electron
temperature $T_{e}$. Later simulations by Aloy et al. (Aloy et al., 2015) have
merged the Multi-Scale Fluid-Kinetic Simulation Suite with the high resolution
3D RMHD code MR-GENESIS. HARM continues to be the benchmark against which the
accuracy of GRMHD simulations is measured (Porth et al., 2016).

Table 1: M87 Geometry

Distance from Earth1 | Schwarzschild radius2 | Apparent angular width3 | Jet opening angle4 | Jet viewing angle5
---|---|---|---|---
$(16.7\pm 0.6)$MPc |  $(6.35\times 10^{-4}\pm 3.69\times 10^{-5}$)pc | 3.9$\mu$as |  $5^{\circ}$ (at 100pc) | $10^{\circ}-19^{\circ}$
($5.15\times 10^{25}\pm 2.78\times 10^{24}$)cm | ($1.94\times 10^{15}\pm 2.09\times 10^{14}$)cm |  |  $60^{\circ}$ (at the core) |

1Blakeslee et al. (2009)
2Event Horizon Telescope Collaboration et al. (2019a)
3Gebhardt et al. (2011)
4Doeleman et al. (2012)
5Wang & Zhou (2009)

Table 2: Code scale $M$ to physical units for M87 Unit type | $M$
---|---
Mass | $(6.6\pm 0.4)\times 10^{9}M_{\odot}$
Length |  $(3.2\times 10^{-4}\pm 1.8\times 10^{-5})$pc
Angular width at $d_{\mathrm{M87}}$ |  $(3.7\pm 0.3)\mu$as
Time |  $(9.1\pm 0.8)$h

###  3.2 Fiducial Simulations

####  3.2.1 Overview

In this work, we use a set of three numerical GRMHD simulations of black hole
accretion. The fluid simulations were produced with the KHARMA code, a GPU-
based descendant of iharm, a conservative second-order explicit shock-
capturing finite-volume code for arbitrary stationary spacetimes (Gammie et
al., 2003; Prather et al., 2021). The governing equations of ideal GRMHD can
be written as a set of conservation laws; in a coordinate basis, they are

$\displaystyle\partial_{t}\big{(}\sqrt{-g}\rho u^{t}\big{)}$
$\displaystyle=-\partial_{i}\big{(}\sqrt{-g}\rho u^{i}\big{)},$ (5)
$\displaystyle\partial_{t}\big{(}\sqrt{-g}T^{t}_{~{}\nu}\big{)}$
$\displaystyle=-\partial_{i}\big{(}\sqrt{-g}T^{i}_{~{}\nu}\big{)}+\sqrt{-g}T^{k}_{~{}\lambda}\Gamma^{\lambda}_{~{}\nu
k},$ (6) $\displaystyle\partial_{t}\big{(}\sqrt{-g}B^{i}\big{)}$
$\displaystyle=-\partial_{j}\big{(}\sqrt{-g}(b^{j}u^{i}-b^{i}u^{j})\big{)},$
(7)

along with a no-monopoles constraint
$\partial_{i}\left(\sqrt{-g}B^{i}\right)=0$. Here, the rest-mass density of
the fluid is $\rho$, $u^{\mu}$ is the fluid four-velocity, $b^{\mu}$ is the
magnetic induction four-vector, and the magnetohydrodynamic stress–energy
tensor is

$\displaystyle
T^{\mu\nu}=\left(\rho+u+P+b^{2}\right)u^{\mu}u^{\nu}+\left(P+b^{2}/2\right)g^{\mu\nu}-b^{\mu}b^{\nu},$
(8)

where $u$ and $P$ are the internal energy of the fluid and its pressure, which
is related to the internal energy via an ideal gas law equation of state with
constant adiabatic index $\hat{\gamma}$ via $P=\left(\hat{\gamma}-1\right)u$.
The effects of spcaetime are accounted for in the usual way, with
$g=\sqrt{-{\rm det}g_{\mu\nu}}$, the determinant of the covariant metric and
$\Gamma$ a Christoffel symbol encapsulating derivatives of the metric. More
detail about the end-to-end simulation procedure can be found in Wong et al.
(2022)

The simulations all used outflow boundary conditions at both the inner and
outer radial edges, located within the event horizon and at
$1,000GM_{\mathrm{H}}/c^{2}$ respectively. Each simulation was run from
$t=0\,GM_{\mathrm{H}}/c^{3}$ until $30,000\,GM_{\mathrm{H}}/c^{3}$ in order to
provide a converged characterization of the source, although we use snapshots
from the latter $25,000\,GM_{\mathrm{H}}/c^{3}$ of each simulation. In
particular, we consider two MAD simulations with dimensionless black hole
spins $a/M_{\mathrm{H}}=-0.5$ and $+0.94$. We also consider one SANE
simulation with spin $a/M_{\mathrm{H}}=-0.5$. In MAD accretion, magnetic flux
carried by the accretion flow builds up near the event horizon until the
magnetic pressure near the hole is large enough to counterbalance the inward
ram pressure of the accreting material. MAD accretion thus proceeds in chaotic
bursts of isolated, thin plasma streams beginning far from the hole, with the
overall flow characterized by occasional violent magnetic eruption events. In
contrast, SANE accretion proceeds in a turbulent but consistent, disk-like
flow. Further information about the KHARMA GRMHD library can be found in Event
Horizon Telescope Collaboration et al. (2022).

Figure 4 shows xz- (poloidal) slices of the fiducial simulations snapshots at
$T=25,000M$ for MAD $a/M_{H}=0.94$, $a/M_{H}=-0.5$ and SANE $a/M_{H}=-0.5$ for
the key quantity magnetization $\sigma=\frac{b^{2}}{\rho}$. Expectedly, the
MAD simulations are more highly magnetized– particularly along the polar
regions where there is a relativistic outflow. The columns in Figure 5
representing electron number density, internal energy and magnetic field
strength for the same simulations show that turbulence in the equatorial
inflow– such as that driven by the magnetorotational instability– is
particularly prominent for the SANE case. The MAD/SANE magnetic substructure
and field strength dichotomy is also apparent in slices of plasma $\beta$ in
Fig. 6.

Figure 4: Magnetization $\sigma\equiv b^{2}/\rho$ in azimuthal slices for
fiducial simulations MAD $a=0.94$ (Left) MAD $a=-0.5$ (Middle) and SANE
$a=-0.5$ (Right). Figure 5: Vertical slices for electron number density
$N_{e}$, internal energy $U$, and magnetic field strength $|B|$ in cgs units
for fiducial simulations MAD $a=0.94$ (Top) MAD $a=-0.5$ (Middle) and SANE
$a=-0.5$ (Bottom).

####  3.2.2 Mass Accretion Rate

The mass accretion rate in the code is an adjustable parameter with which the
flux scales. In this work, our target flux for synthetic images is 0.5 Jy. Now
that we have described the simulations we are using to test the radiative
properties of M87’s JAB system, we add the key physics governing energy
transfer from the GRMHD plasma to the high energy particles responsible for
the observed emission.

##  4 Model Prescriptions for Emission, Particle Acceleration and Dissipation

###  4.1 General Considerations

We suppose that the radio and mm emission is synchrotron radiation due to
particle acceleration arising from a number of different mechanisms discussed
here and in the Appendix. We expand upon the synchrotron prescriptions
implemented for jet models in Anantua et al. (2018, 2020)a, where the
relativistic electron distribution function is a power law with slope 2, which
implies that the emissivity in the comoving (primed) frame
$j_{\nu^{\prime},\Omega^{\prime}}^{\prime}(\nu^{\prime},{\bf
n}^{\prime})\propto{\tilde{P}}_{e}|{\bf B}^{\prime}\times{\bf
n}^{\prime}|^{3/2}\nu^{\prime-1/2}$, where ${\tilde{P}}_{e}$ is the (presumed
isotropic) partial pressure of the electrons emitting at the frequency
$\nu^{\prime}$. The choice to scale near-horizon jet emissivity in terms of
magnetic pressure is motivated by the observation that the jet becomes
increasingly simple and electromagnetic — high $\sigma$ — as the horizon is
approached and as exhibited by the RMHD simulations. However, we note the
assumptions underlying some of the previous simulations make no provision for
the particle acceleration and transport resulting in the spatial and temporal
variation of ${\tilde{P}}_{e}$. In this work, we use a more generic electron
pressure formalism distinguishing the highly magnetized outflow from the less
magnetized inflow to prescribe emission models. It is the primary objective of
this investigation to explore how much this matters and if indeed any
prescription for the particle acceleration is compatible with existing and
anticipated VLBI imaging.

Perfect MHD defines a set of reference frames in which the plasma, treated as
a fluid, is at rest. The motion perpendicular to the magnetic field has a
velocity ${\bf v}_{\perp}={\bf E}\times{\bf B}/B^{2}$ in the simulation frame.
In general, the component of the fluid velocity resolved along the field ${\bf
v}_{\parallel}$ is problematic when the inertia of the plasma is ignorable as
we are implicitly assuming here. It should be emphasized that the minimum
charged particle density needed to support the space charge and current
(Goldreich & Julian, 1969) is orders of magnitude smaller than what is needed
to account for the radio and $\gamma$-ray emission. It should also be stressed
that efficient and progressive pair production and particle acceleration is to
be expected in AGN jets as modeled here. The motional potential difference
across an electromagnetic jet near the black hole should be $\sim 1-300\,{\rm
EV}$, many orders of magnitude greater than the $\sim 1\,{\rm MV}$ minimum
needed to create positrons or the $\sim 1\,{\rm GV}$ needed to accelerate
electrons to the $\lesssim 1\,{\rm GeV}$ energies associated with the mm
emission. The numerical simulations express none of this physics and, in any
case, introduce a “floor” to the electron density for purely numerical
reasons. The simulation particle density should not be trusted within the
inner jet.

The composition of the plasma is also uncertain. Close enough to the event
horizon, the plasma must fall inward and be connected to a source as jets are
outflowing at larger radius. The simplest assumption, which we shall adopt, is
that pairs are continuously produced in the inner magnetosphere. Plasma can
also be entrained from the surrounding medium. This is expected to play a
large role in the dynamical evolution and emission of the jet at large radii.
The simulations suggest that it is unimportant within, say, $\sim 1000M$ and
we shall suppose that the associated radio emission is a locally accelerated
pair plasma where $\gamma$-ray production is balanced by annihilation and
other quantum electrodynamical processes.

###  4.2 Phenomenological Jet Models

Given all this uncertainty, we adopt a phenomenological approach where we
adopt a set of simple prescriptions for electron and positron pressure and
temperature that can lead to quite different simulated observations. We now
develop a compendium of formulas linking plasma variables to radiation
phenomenology, starting with jet models for the electromagnetically regions of
a force-free, relativistic plasma.

####  4.2.1 Constant Electron $\beta$ Model

The simplest prescription is the Constant $\beta_{e}$ Model, that the pressure
of synchrotron radiating electrons is a constant fraction of the magnetic
pressure

$P_{e}=\beta_{e0}P_{B}$ (11)

Close to the hole, the total pressure $P$ is dominated by the magnetic
pressure $P_{B}=b^{2}/2\mu_{0}$ but at large radius, there is gas pressure
contributed by the entrained gas. If we adopt $\beta\lesssim 0.01$, close to
the hole, then this is the sub-“equipartition” prescription that is often
invoked when interpreting jet observations e.g.,  Anantua et al. (2020)b.

Figure 6: Simulation plasma $\beta$ parameter $\equiv P_{\rm gas}/P_{\rm mag}$
on azimuthal slices for fiducial models MAD $a=0.94$ (Left) MAD $a=-0.5$
(Middle) and SANE $a=-0.5$ (Right).

####  4.2.2 Current Density ($j^{2}$) Model

Another type of model by which electromagnetic jets can dissipate power
employs currents. This type of model can be implemented using the field
gradient — the current density $j^{\prime}$ — and introducing a resistivity
$\eta=\mu_{0}L_{j}$ where $L_{j}$ is a length scale which we choose to be a
fixed fraction of the jet width. The dissipation rate is then $W^{\prime}=\eta
j^{\prime 2}$. This approach is partly motivated by particle-in-cell (PIC)
simulations of relativistic reconnection. This model has been compared to the
43 GHz M87 jet in Anantua et al. (2018) as a mechanism for generating limb
brightening. Though we do not make images of the Current Density Model here,
we note the physical significance of jet currents as a spatially inhomogeneous
source of dissipation.

###  4.3 JAB System Models

We now consider the entire inflow-outflow structure governed by the
supermassive black hole. We have previously described the relativistic polar
outflow as a Blandford-Znajek jet. Beyond the outflow or jet funnel,
astrophysical plasmas experience discontinuities in pressure and density. When
there is a sufficient velocity gradient, Kelvin-Helmholtz instabilities
produce gaseous swirling features– notable in M87 up to 1 kpc from the black
hole (Pasetto et al., 2021). The enveloping corona is loosely bound to the JAB
system. The inflowing disk is supported against its own inertia by magnetic
and thermal pressure and momentum transport– the latter of which may lead to
the magneto-rotational instability.

The property of turbulent heating to preferentially energize electrons in
magnetically dominated regions and protons in gas pressure dominated regions
has been explored in Howes (2010). There are several ways of parameterizing
this behavior (Mościbrodzka et al., 2011; Anantua et al., 2020)b.

####  4.3.1 $R-\beta$ Model

We start our JAB emission modeling by noting the tendency of plasma turbulence
to preferentially heat electrons at low $\beta$ and ions at high $\beta$, as
was originally conceptualized in the context of the solar corona (Quataert &
Gruzinov, 1999; Howes, 2010). Applied to JAB systems, the $R-\beta$ turbulent
heating model takes the form

$R=\frac{T_{i}}{T_{e}}=\frac{\beta^{2}}{1+\beta^{2}}R_{\mathrm{high}}+\frac{1}{1+\beta^{2}}R_{\mathrm{low}}$
(12)

It is the primary model used by the Event Horizon Telescope (Event Horizon
Telescope Collaboration et al., 2019e, 2021b) and developed by Moscibrodzka et
al. (2016).

####  4.3.2 Critical $\beta$ Model

The Critical $\beta$ Model is an alternative turbulent heating model with an
exponential parameter $\beta_{c}$ controlling the transition between electron-
and ion-dominated heating

$\frac{T_{e}}{T_{e}+T_{i}}=fe^{-\beta/\beta_{c}}$ (13)

This model was developed in Anantua et al. (2020)b. The models are compared
for reasonable parameter values in Fig. 7. We see at low $\beta$, the
electron-to-ion temperature ratios in the $R-\beta$ and the Critical Beta
Models have similar asymptotic behavior. However, at high $\beta$, the
Critical $\beta$ Model $T_{e}/T_{i}$ exponentially falls to 0, which by
preliminary indications may reduce the bremsstrahlung contribution (though we
reserve a more extensive investigation modeling emission processes beyond
synchrotron for future work). The Critical $\beta$ Model also has transition
behavior at intermediate $\beta$’s controlled by an exponential parameter
$\beta_{c}$, leading to greater variety of intermediate $\beta$ emission
regions probed between the same range of electron-to-ion temperature ratios
compared to the $R-\beta$ Model.

Figure 7: Comparison of R-$\beta$ (solid lines) and Critical Beta (dashed
lines) models for reasonable parameter values.

####  4.3.3 Multi-Zone Emission Models

We have seen how plasma inertial and electromagnetic properties within and
across our JAB simulations differ by orders of magnitude through $\beta$ and
$\sigma$. Moreover, differences in the plasma velocity towards and away from
the black hole leads to plasma mixing regions and whose behavior may not be
amenable to the existing smooth emission parametric models. We thus refine our
JAB system emission models by combining the $\beta$-dependent turbulent
heating models with jet regions radiating by conversion of magnetic to
particle energy. We define the jet region using a transitional value of
magnetization $\sigma$ of 1/2. Armed with these models of particle
thermodynamics, we turn to the emitted radiation.

###  4.4 Radiation

An electromagnetically-dominated jet outflow must be continually converting
electromagnetic energy into high energy pairs through a $-{\bf E}\cdot{\bf j}$
interaction. These particles will radiate through the synchrotron and Compton
processes. We can use the observations to draw some inferences about the
variation of the distribution function along the jet.

At sufficiently small radius, the jet will become optically thick to
synchrotron self-absorption at a given frequency given by

$\nu_{\rm SA}=\\\ 22\left(\frac{L_{\rm EM}}{10^{44}{\rm erg\
s}^{-1}}\right)^{0.1}\left(\frac{S}{\rm
1Jy}\right)^{0.4}\left(\frac{r}{50M}\right)^{-0.8}\left(\frac{\Gamma}{3}\right)^{-0.4}{\rm
GHz}$ (15)

Also, at a given radius there is a characteristic frequency where the
synchrotron cooling time of the emitting electrons is equal to the expansion
time scale

$\nu_{\rm cool}=\left(\frac{L_{\rm EM}}{10^{44}{\rm erg\
s}^{-1}}\right)^{-1.5}\left(\frac{r}{50M}\right)^{-0.2}\left(\frac{\Gamma}{3}\right)^{6}{\rm
THz}$ (16)

If a power law is accelerated, the local spectrum should break by
$\Delta\alpha=0.5$ at this frequency. However, note the extreme sensitivity to
the bulk Lorentz factor $\Gamma$. This probably controls what is actually
observed.

The entire nuclear spectrum within $r\sim 3\times 10^{5}M$ has been carefully
determined by Prieto et al. (2016). They find a sharp break at
$\nu=\nu_{b}\sim 150$ GHz. Presumably the flatter $\alpha\sim 0.2$ spectrum at
$\nu\lesssim\nu_{b}$ is attributable to the superposition of a radial sequence
of spectra with a frequency to radius mapping as defined above. The lowest
frequency considered, $\sim 3\ {\rm GHz}$, should originate at $r\sim 1000M$.

###  4.5 Emission Modeling with Positrons

####  4.5.1 Particle Production and Acceleration

Though the vast majority of GRMHD simulations consider ionic plasma without an
explicit fluid for electron-positron pairs, it is unknown whether jet matter
is typically dominated by an ion-electron plasma or a pair one. In the latter
case, black hole-powered jets generally get their initial mass-loading through
$\gamma-\gamma$ pair production. The means of particle acceleration within AGN
jets, however, is not understood. There have been several mechanisms invoked,
including those involving strong shock fronts, both non-relativistic and
relativistic, magnetic reconnection, stochastic acceleration though wave-
particle interaction and electrostatic acceleration – either along a magnetic
field or perpendicular to the field as a consequence of drift motion. Recent
observations, especially when seen in the context of observations of extreme
acceleration in pulsar wind nebulae and Gamma Ray Bursts, point to the need
for new and more rapid approaches. The need is best exemplified by
$\gamma$-ray observations which can show that electric fields as strong as
magnetic fields (setting $c=1$) may be needed in order for electrons to attain
the required energies in the face of strong radiative loss through synchrotron
emission and Compton scattering.

####  4.5.2 Positron Production Modeling and Radiative Transfer

Much of the plasma in the accreting component of the RIAF system is likely a
mixture of ionized hydrogen and helium from stellar winds and the interstellar
medium. Since the conductivity is high near the event horizon, particles in
the plasma are forced to follow magnetic field lines, so the jet, which is
canonically magnetically disconnected from the accretion disk, cannot be
directly supplied with plasma from the disk. If electron–positron pairs are
produced in these regions, then they may be the dominant matter source.
Electron–positron pairs may be produced by pair cascades in charge-starved
magnetospheres (like in evacuated jet funnels) or in the disk coronae.

In the systems we study, electron–positron pairs are produced mainly via the
Briet–Wheeler process, i.e., as a result of photon–photon collisions. In order
to create a pair, the center-of-momentum energy of the photons must exceed the
rest-mass energy of a pair $\approx 1\,$MeV $\approx 2\times\left(1.2\times
10^{20}\right)\,$Hz. The cross-section peaks near this threshold value, and
the participating photon couples lie over a spectrum of energy ratios: some
pairs of photons having approximately the same frequency while others are
matched with low/high frequencies. Pair-producing processes are often
differentiated based on the photon source and whether the newly created pairs
radiate and contribute (non-negligibly) as a new photon source.

_Pair drizzle_ occurs when the pairs are produced by photons from the
background radiation field (due to synchrotron and bremsstrahlung emission and
Compton upscattering) and typically exhibits variation on timescales
associated with the plasma fluid. Drizzle pair production has been studied in
a variety of scenarios ranging between stellar-mass to supermassive black hole
accretion Mościbrodzka et al. (2011); Laurent & Titarchuk (2018); Kimura &
Toma (2020); Wong et al. (2021); Yao et al. (2021). In the alternative
scenario, high energy photons with frequencies $\gg 10^{20}\,$Hz can interact
with background (low energy) photons from the disk to undergo pair production.
Here, the high-energy photons are produced when unscreened electric fields
accelerate stray charges (Beskin et al., 1992); when the acceleration is
large, the leptons radiate the requisite high energy photons. Often, the newly
created pairs are born in the same region with unscreened fields and are thus
themselves accelerated, restarting the process in a cascade of pair creation.
The short timescales associated with pair cascades means they may explain the
ultra-rapid high-frequency radio emission flares from AGN jets. Pair cascades
have been studied with a variety of analytic, semi-analytic, and numerical
methods, e.g., Fragile & Meier (2009); Broderick & Tchekhovskoy (2015);
Parfrey et al. (2019). Positrons not only effect images at the level of
emission, but also through radiative cooling, e.g., Fragile & Meier (2009);
Yoon et al. (2020).

Given the uncertainty in jet positron fraction, we focus on the special cases
of a sparse ionic jet with electron number density $n_{e0}$ and a jet where
all sources of pair production result in a plasma with an equal number density
of ionic and pair plasma ($f_{\mathrm{pos}}\equiv
n_{\mathrm{pairs}}/n_{e0}=1$).

There are benchmarks for positron content in the Literature, such as the
estimate by Ghisellini (2012) of the fraction of a jet (opening angle $\psi$,
distance from black hole $R_{0}$) converted to positrons, as
$f=0.1\min\\{1,\frac{\ell_{0}}{60}\\}$ where the compactness is
$\ell=\frac{\sigma_{T}L_{0}}{\psi R_{0}\mathcal{D}^{4}m_{e}c^{3}}$. More
possibilities for painting positrons on jets can be found in Appendix.

##  5 Observing a Time-dependent Simulation

###  5.1 Anatomy of a Time-Dependent KHARMA Jet Simulation

We have outlined a methodology for combining jet emission prescriptions with
detailed, 3D time-dependent simulations. To emphasize the 3D nature, we take
transverse slices in the equatorial plane in Figure 8. There, we see even
among two MADs (spins $a/M=0.94,-0.5$) there are large, azimuthal symmetry-
breaking patches of high magnetization emanating in different directions from
the black hole. Electron number density exhibits a similar pattern of
asymmetry, however, internal energy, magnetic field strength and plasma
$\beta$ are relatively azimuthally symmetric.

We now describe the process of ray tracing the resulting emission.

###  5.2 Radiative transfer with IPOLE: Azimuthal and Polar Variation

GRMHD simulations can be ray-traced using general relativistic radiative
transfer (GRRT) codes to simulate surveys of sources throughout the sky. In
this work, we use the GRRT code ipole (Mościbrodzka & Gammie, 2018, see also
Wong et al. 2022) to produce polarimetric images of the GRMHD simulations
described in Section 3.2.1. The ipole code solves for the evolution of the
polarized intensities at each point along a geodesic with a two-stage operator
splitting method. In the first stage, the covariant coherency tensor (which
can be written in terms of the invariant Stokes parameters and Pauli matrices)
is parallel transported along the direction of the geodesic. In the second
stage, the Stokes parameters are updated using the analytic solution to the
explicit, general polarized transport with constant emission, absorption, and
rotation coefficients, which are computed in the local orthonormal tetrad
defined by the fluid and magnetic field.

Each image comprises a square grid of $N$x$N$ pixels over a $160\,\mu$as field
of view, with each pixel reporting the Stokes intensities for $I,Q,U$, and
$V$. Since producing images requires evaluating transfer coefficients in
physical units, it is necessary to specify scales for the mass-density of the
accreting plasma and the size of the black hole as well as the orientation of
the observer (i.e., the software camera) with respect to the black hole. We
list the physical M87 black hole parameters (and references) in Table 1, and
in Table 2 we report the “code scale” values corresponding to the black hole
mass identified above.

Note that GRMHD codes are generally unable to accurately evolve the fluid
state in regions with high magnetization $\sigma\equiv b^{2}/\rho$ and
artificially inject mass and energy in these regions. Since the plasma density
(and temperature) are therefore typically unphysically high in regions of
large $\sigma$, ray-tracing codes like ipole normally introduce a so-called
$\sigma$ cutoff, where the plasma density in regions with $\sigma>\sigma_{\rm
cutoff}$ is explicitly set to zero before computing the transfer coefficients.
In this work, we set  $\sigma_{\rm cutoff}=2$, consistent with the typical
values used for such flows (see, e.g., Event Horizon Telescope Collaboration
et al. (2019e).

(17)

(18)

(19)

Figure 8: Equatorial slices for fiducal MAD $a/M=0.94,-0.5$ and SANE
$a/M=-0.5$ of $\sigma$, $N_{e}$, internal energy, $|B|$ and $\beta$.

##  6 COMPARISON OF SIMULATIONS WITH OBSERVATIONS

We now present a suite of 3D-GRMHD simulation images spanning various plasma
compositions and prescriptions for electron-positron thermodynamics. JAB
systems are often modeled as electron-proton plasmas with a single function
linking electron temperature to plasma variables such as $\beta$ throughout
the inflow/outflow system Event Horizon Telescope Collaboration et al. (2019e,
2021b). We start with models using this approach with the $R-\beta$ and
Critical $\beta$ turbulent heating prescriptions, and then we refine the
models by prescribing jet funnel emission in a region between
$\sigma_{\mathrm{min}}=1/2$ and $\sigma_{\mathrm{max}}=2$ and adding pairs.
Unless otherwise stated, the images are raytraced at 230 GHz and the
inclination angle is 17∘. Note that in the tables referenced in this section,
comparisons are made with snapshots in the GRMHD space. Due to this, model
performance with respect to observations may not hold when comparisons to
windows of simulations are made such as those in Event Horizon Telescope
Collaboration et al. (2019a).

The image library presented here greatly expands the $a=-0.5M$ MAD and SANE
snapshots from Anantua et al. (2023) to include the highly prograde spin
$a=0.94$, temporal evolution, extreme positron fractions
$n_{\mathrm{pairs}}/n_{e0}=50,100$, varying frequency from 230 GHz to 86 GHz
and varying inclination up to $40^{\circ}$ viewing angle. Here, the SANE-MAD
dichotomy manifest in the image library is also made quantitative by the
tabulation of Faraday conversion and rotation depths and comparison to M87
linear polarization data.

###  6.1 SANE Positron Effects

####  6.1.1 SANE R-$\beta$

Fig. 9 shows intensity with electric vector polarization angle (EVPA) and
circular polarization maps for the SANE $a=-0.5$ simulation in the $R-\beta$
model. This model has asymptotic ion-to-electron temperature ratios
$R_{\mathrm{low}}=1$ and $R_{\mathrm{high}}=20$. Here, the total flux is
greatest near the photon ring immediately surrounding the central depression
and slowly decreases radially becoming broadly distributed through the
equatorial annulus. Polarization oriented at the EVPA is spread throughout the
equatorial annulus in the $40\mu$as x $40\mu$as field of view. To the ionic
plasma in the Top Panels, an equal number density of pairs as the original
electron number density are added in the Bottom Panels (while renormalizing
$m_{\mathrm{unit}}$ to maintain a 0.5 Jy image flux). The added positrons
significantly rotate EVPAs, as the Faraday rotation measure depends
sensitively on the positron fraction for SANEs.

Figure 9: For the $a=-0.5$ SANE at $T=25,000M$: Top panel: R-Beta at 230 GHz
without positrons. Bottom panel: R-Beta at 230 GHz with for an even mix of
pair and ionic plasma with ions and positrons each accounting for 1/4 the
plasma number density and electrons accounting for the remaining 1/2.

####  6.1.2 SANE Critical $\beta$

In Fig. 10 we show our other turbulent heating model– the Critical Beta Model.
This model controls the transition from preferential electron heating to
preferential ion heating through the exponential parameter $\beta_{c}$, which
for higher values smooths the transition by allowing a larger range of betas
to include radiating electrons. For model parameters temperature ratio
prefactor $f=0.5$ and $\beta_{c}=0.5$, total flux is concentrated in a ring at
$\sim 20\mu$as and regions along lines of sight close to the polar axis,
though polarization morphology trends remain similar to the $R-\beta$ case.
Note our $R-\beta$ models are more linearly depolarized than Critical $\beta$
Models even with lower contributions to intrinsic emission at high $\beta$ in
the latter models. This is one of several examples of Faraday effects we will
see in this Section.

Figure 10: For the $a=-0.5$ SANE: Top panel: Critical Beta at 230 GHz without
positrons Bottom panel: Critical Beta at 230 GHz for an even mix of pair and
ionic plasma.

####  6.1.3 SANE R-$\beta$ with Constant $\beta_{E}$ Jet

In Fig. 11 we add a jet region where the energy of relativistic electrons is
directly derived from the magnetic pressure to the $R-\beta$ model. The
emission is extended more broadly and evenly throughout the field of view as
it is projected from a broader region of the outflow paraboloid governed by
the transitional value of $\sigma$ separating the constant $\beta_{e}$ jet
from the turbulently heated plasma.

Figure 11: For the $a=-0.5$ SANE at $T=25,000M$: Top panel: R-Beta with
$\beta_{e0}=0.01$ jet at 230 GHz without positrons Bottom panel: R-Beta at 230
GHz for an even mix of pair and ionic plasma.

####  6.1.4 SANE Critical $\beta$ with Constant $\beta_{E}$ Jet

In Fig. 12 we add a jet region of magnetic-to-particle energy conversion to
the Critical Beta model. In the SANE case, the jet does not appreciably change
the image morphology. Moreover, polarization does not vary monotonically with
the addition of positrons across different emission models.

Figure 12: For the $a=-0.5$ SANE at $T=25,000M$: Top panel: Critical Beta with
$\beta_{e0}=0.01$ jet at 230 GHz without positrons. Bottom panel: Critical
Beta at 230 GHz for an even mix of pair and ionic plasma.

###  6.2 MAD Positron Effects

The MAD images from our fiducial time are nearly indistinguishable when
Faraday effects are turned off. In these MAD images, we see a prominent flux
tube in a loop extending towards the lower left. In these particular MAD
images, whose circular polarization is dominated by the intrinsic, we see
another dramatic polarization effect: linear increase in the magnitude of
$V/I$ (confer Section 6.5 for polarimetric quantity definitions) as a function
of synchrotron emitters not in pairs (which is maximal for the ionic plasma
case).

####  6.2.1 MAD R-$\beta$

Starting in Fig. 13 with the R-$\beta$ model, the linear polarization ticks
oriented at the EVPA for the 0-positron case remain in their orientations with
minimal angular displacement when positrons are added to form the mixed plasma
in the ray tracing step. However, in radiative transfer using coefficients for
the mixed plasma in which 1/2 the particles are electrons, 1/4 of the
particles are positive ions and 1/4 of the particles are positrons (i.e, 2/3
of the synchrotron emitting leptons are paired), we have the degree of
circular polarization $V/I$ diminishing to 1/3 of the positron-free value. The
addition of positrons also reverses the polarity of the bottom left portion of
the flux eruption loop.

Figure 13: For the $a=-0.5$ MAD at $T=25,000M$: Top panel: R-Beta at 230 GHz
without positrons Bottom panel: R-Beta at 230 GHz for an even mix of pair and
ionic plasma.

####  6.2.2 MAD Critical $\beta$

In Fig. 14, the Critical $\beta$ image and $V/I$ map mirror the global
structure in the $R-\beta$ case in Fig. 13. They also share similar dependence
of the circular polarization dependence on the free electron-positron
fraction, and partial reversal of circular polariaztion sense in the flux
eruption loop.

Figure 14: For the $a=-0.5$ MAD at $T=25,000M$: Top panel: Critical Beta at
230 GHz without positrons Bottom panel: Critical Beta at 230 GHz for an even
mix of pair and ionic plasma.

####  6.2.3 MAD R-$\beta$ with Constant $\beta_{E}$ Jet

In Fig. 15, the R-Beta model with jet maintains the prominent flux eruption
loop as the above models. The presence of the Constant $\beta_{e}$ jets
slightly reduced the circular polarization degree both with and without
positrons.

Figure 15: For the $a=-0.5$ MAD at $T=25,000M$: Top panel: R-Beta with
$\beta_{e0}=0.01$ jet at 230 GHz without positrons Bottom panel: R-Beta at 230
GHz for an even mix of pair and ionic plasma.

####  6.2.4 MAD Critical $\beta$ with Constant $\beta_{E}$ Jet

In Fig. 16 the Critical Beta model with jet exhibits the same trends as its
$R-\beta$ counterpart above.

Figure 16: For the $a=-0.5$ MAD at $T=25,000M$: Top panel: Critical Beta with
$\beta_{e0}=0.01$ jet at 230 GHz without positrons Bottom panel: Critical Beta
with $\beta_{e0}=0.01$ jet at 230 GHz for an even mix of pair and ionic
plasma.

###  6.3 Comparison of R-Beta and Critical Beta Models

Figure 17: R-Beta vs. Critical Beta comparison for MAD $a=+0.94$. Top panel:
R-Beta at 230 GHz. Bottom panel: Critical Beta

In Fig. 17 we compare Critical Beta Model with $R-\beta$ Model for MAD
$a/M=0.94$. Like comparisons between Figs. 10 and 11 and between 13 and 14
show, these models are quite degenerate at the level of total intensity
morphology. However, the circular polarization of the Critical Beta $V/I$ map
does exhibit more scrambling near the photon ring– where a broad range of
$\beta$’s may contribute given our shallow exponential parameter
$\beta_{c}=0.5$.

###  6.4 Extreme Positron Fractions

As mentioned in Sec. 6 increasing the pair fraction causes intrinsically
emitted circular polarization to decrease, Faraday rotation to decrease, and
Faraday conversion to increase. Faraday rotation is essential for depolarizing
these accretion flows. Thus, dramatic effects occur when the pair fraction is
raised high enough to turn an model from Faraday thick to Faraday thin.

In Fig. 18 we show the effects of raising $n_{\mathrm{pairs}}/n_{0}=100$ for a
SANE $a=-0.5$ simulation and contrast with MAD $a=-0.5$ and $a=+0.94$ in Figs.
19 and 20, respectively (using the R-$\beta$ Model). The effect on the MAD
simulation is subtle, characterized by a decrease in the intrinsically emitted
circular polarization that dominates on large scales. Note that this is model-
specific, and in fact Faraday conversion is the only source of circular
polarization in pair plasma jet models in Anantua et al. (2020)a Meanwhile,
the effect is much stronger for the SANE simulation, which is intrinsically
more Faraday thick. After removing Faraday rotation, the simulation acquires a
much more ordered linear polarization pattern. In addition, very large
circular polarization fractions are produced in the absence of depolarization
by Faraday rotation.

Figure 18: Extreme positron variation comparison for SANE $a=-0.5$ R-Beta at
230 GHz: fPos = 0 (Top Panel) vs. fPos = 100 (Bottom Panel).

Figure 19: Extreme positron variation comparison for MAD $a=-0.5$ R-Beta at
230 GHz: fPos = 0 (Top Panel) vs. fPos = 100 (Bottom Panel).

Figure 20: MAD $a=+0.94$ R-Beta at 230 GHz.Top panel: fPos = 0 Bottom panel:
fPos = 100

###  6.5 Comparison of Models with Polarization Constraints

In Table 3, we compare observations from Event Horizon Telescope Collaboration
et al. (2021a) against fiducial model linear polarization: both summed from
net (unresolved) $Q$, $U$ and $I$ across the image plane,

$|m|_{\mathrm{net}}=\frac{\sqrt{\left(\sum_{i}{Q_{i}}\right)^{2}+\left(\sum_{i}{U_{i}}\right)^{2}}}{\sum_{i}{I_{i}}}$
(20)

and its average local (resolved) magnitude

$<|m|>=\frac{\sum_{i}I_{i}P_{i}}{\sum_{i}{I_{i}}}=\frac{\sum_{i}\sqrt{Q_{i}^{2}+U_{i}^{2}}}{\sum_{i}{I_{i}}}.$
(21)

SANE models tend to be less linearly polarized on net and MAD models more
linearly polarized on net than the constraint– though all models exceed the
beam/resolution-dependent averaged polarization magnitude constraint. The
fiducial models satisfying the net linear polarization constraint for MAD
$a=-0.5$ are the R-Beta, R-Beta with Jet and Critical Beta with Jet, all with
maximal positron fractions. For MAD $a=+0.94$, only the positron-free Critical
Beta moddels (with and without jet) satisfy the $|m|_{\mathrm{net}}$
constraint.

In Table 4, we project forward anticipating circular polarization measurements
will be performed in the near future to compare with our models. All of our
models satisfy the preliminary $V_{\mathrm{net}}$ constraint in Goddi et al.
(2021). For the structure parameter $\beta_{2}$, however, Table 5 shows only
MAD $a=0.94M$ pure turbulent heating models pass.

Table 3: Linear polarization $|m|_{\mathrm{net}}$ and $<|m|>$ for fiducial
models at $T=25,000M$. The observational constraints from EHT M87 Paper VII
take the form of the polarization ranges $0.01\leq|m|_{\mathrm{net}}\leq
0.037$ and $0.057<<|m|><0.107$.  Note that the bold values refer to fiducial
models which satisfy the net linear polarization constraints.

. $\mathrm{SANE}$ ($a=-0.5$) $\mathrm{MAD}$ ($a=-0.5$) $R\mathrm{-Beta}$
$R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\ Beta}$
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$\mathrm{w./\ Jet}$ $|m|_{\mathrm{net}}(f_{\mathrm{pos,min}})$ $3.93\cdot
10^{-3}$ $3.97\cdot 10^{-3}$ $8.51\cdot 10^{-3}$ $2.26\cdot 10^{-3}$
$5.51\cdot 10^{-2}$ $4.07\cdot 10^{-2}$ $6.93\cdot 10^{-2}$ $4.86\cdot
10^{-2}$ $|m|_{\mathrm{net}}(f_{\mathrm{pos,max}})$ $1.62\cdot 10^{-3}$
$2.88\cdot 10^{-3}$ $2.73\cdot 10^{-3}$ $2.50\cdot 10^{-3}$ $\mathbf{3.67\cdot
10^{-2}}$ $\mathbf{3.10\cdot 10^{-2}}$ $5.21\cdot 10^{-2}$ $\mathbf{3.55\cdot
10^{-2}}$ $<|m|>(f_{\mathrm{pos,min}})$ $1.29\cdot 10^{-1}$ $1.38\cdot
10^{-1}$ $1.35\cdot 10^{-1}$ $1.44\cdot 10^{-1}$ $3.49\cdot 10^{-1}$
$3.30\cdot 10^{-1}$ $2.71\cdot 10^{-1}$ $2.53\cdot 10^{-1}$
$<|m|>(f_{\mathrm{pos,max}})$ $1.31\cdot 10^{-1}$ $1.45\cdot 10^{-1}$
$1.40\cdot 10^{-1}$ $1.47\cdot 10^{-1}$ $4.20\cdot 10^{-1}$ $3.60\cdot
10^{-1}$ $3.83\cdot 10^{-1}$ $3.46\cdot 10^{-1}$ $\mathrm{MAD}$ ($a=+0.94$)
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$|m|_{\mathrm{net}}(f_{\mathrm{pos,min}})$ $4.91\cdot 10^{-2}$ $4.56\cdot
10^{-2}$ $\mathbf{3.35\cdot 10^{-2}}$ $\mathbf{3.67\cdot 10^{-2}}$
$|m|_{\mathrm{net}}(f_{\mathrm{pos,max}})$ $5.17\cdot 10^{-2}$ $5.06\cdot
10^{-2}$ $4.59\cdot 10^{-2}$ $4.82\cdot 10^{-2}$ $<|m|>(f_{\mathrm{pos,min}})$
$5.76\cdot 10^{-1}$ $5.18\cdot 10^{-1}$ $5.18\cdot 10^{-1}$ $4.81\cdot
10^{-1}$ $<|m|>(f_{\mathrm{pos,max}})$ $5.86\cdot 10^{-1}$ $5.25\cdot 10^{-1}$
$5.81\cdot 10^{-1}$ $5.26\cdot 10^{-1}$

Table 4: Circular polarization $|V|_{\mathrm{net}}$ and $<|V|>$ for fiducial
models at $T=25,000M$. Note that all models satisfy the
$V_{\mathrm{net}}$=0.008 EHT bound (Event Horizon Telescope Collaboration et
al., 2021a; Goddi et al., 2021)

. $\mathrm{SANE}$ ($a=-0.5$) $\mathrm{MAD}$ ($a=-0.5$) $R\mathrm{-Beta}$
$R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\ Beta}$
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$\mathrm{w./\ Jet}$ $|V|_{\mathrm{net}}(f_{\mathrm{pos,min}})$ $-5.60\cdot
10^{-3}$ $-1.65\cdot 10^{-3}$ $-7.77\cdot 10^{-3}$ $-1.72\cdot 10^{-3}$
$-4.50\cdot 10^{-3}$ $-3.54\cdot 10^{-3}$ $-4.24\cdot 10^{-3}$ $-3.59\cdot
10^{-3}$ $|V|_{\mathrm{net}}(f_{\mathrm{pos,max}})$ $-2.81\cdot 10^{-3}$
$2.67\cdot 10^{-3}$ $3.15\cdot 10^{-3}$ $7.44\cdot 10^{-4}$ $-1.39\cdot
10^{-3}$ $-8.99\cdot 10^{-4}$ $-1.19\cdot 10^{-3}$ $-8.45\cdot 10^{-4}$
$<|V|>(f_{\mathrm{pos,min}})$ $6.75\cdot 10^{-3}$ $5.14\cdot 10^{-3}$
$8.94\cdot 10^{-3}$ $6.62\cdot 10^{-3}$ $3.98\cdot 10^{-3}$ $2.96\cdot
10^{-3}$ $3.82\cdot 10^{-3}$ $3.03\cdot 10^{-3}$ $<|V|>(f_{\mathrm{pos,max}})$
$1.28\cdot 10^{-2}$ $7.65\cdot 10^{-3}$ $2.32\cdot 10^{-2}$ $9.62\cdot
10^{-3}$ $9.04\cdot 10^{-4}$ $5.75\cdot 10^{-4}$ $8.57\cdot 10^{-4}$
$5.66\cdot 10^{-4}$ $\mathrm{MAD}$ ($a=+0.94$) $R\mathrm{-Beta}$
$R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{w./\
Jet}$ $\mathrm{w./\ Jet}$ $|V|_{\mathrm{net}}(f_{\mathrm{pos,min}})$
$-3.36\cdot 10^{-3}$ $2.33\cdot 10^{-3}$ $-3.02\cdot 10^{-3}$ $-2.41\cdot
10^{-3}$ $|V|_{\mathrm{net}}(f_{\mathrm{pos,max}})$ $-1.13\cdot 10^{-2}$
$7.37\cdot 10^{-4}$ $-1.19\cdot 10^{-3}$ $-7.29\cdot 10^{-4}$
$<|V|>(f_{\mathrm{pos,min}})$ $3.18\cdot 10^{-3}$ $2.17\cdot 10^{-3}$
$2.85\cdot 10^{-3}$ $2.19\cdot 10^{-3}$ $<|V|>(f_{\mathrm{pos,max}})$
$7.62\cdot 10^{-4}$ $5.02\cdot 10^{-4}$ $7.74\cdot 10^{-4}$ $4.86\cdot
10^{-4}$

Table 5: Azimuthal structure mode $\beta_{2}$ for fiducial models at
$T=25,000M$. The observational constraints from EHT M87 Paper VII are in the
range $0.04\leq|\beta_{2}|\leq 0.07$. Note that the bold values refer to
fiducial models which satisfy the observational constraints.

. $\mathrm{SANE}$ ($a=-0.5$) $\mathrm{MAD}$ ($a=-0.5$) $R\mathrm{-Beta}$
$R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\ Beta}$
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$\mathrm{w./\ Jet}$ $\beta_{2}(f_{\mathrm{pos,min}})$ $3.51\cdot 10^{-3}$
$2.82\cdot 10^{-3}$ $3.84\cdot 10^{-3}$ $3.80\cdot 10^{-3}$ $6.22\cdot
10^{-3}$ $6.31\cdot 10^{-3}$ $1.17\cdot 10^{-2}$ $6.18\cdot 10^{-3}$
$\beta_{2}(f_{\mathrm{pos,max}})$ $2.96\cdot 10^{-3}$ $3.03\cdot 10^{-3}$
$4.24\cdot 10^{-3}$ $9.86\cdot 10^{-4}$ $1.47\cdot 10^{-2}$ $9.83\cdot
10^{-3}$ $2.15\cdot 10^{-2}$ $1.32\cdot 10^{-2}$ $\mathrm{MAD}$ ($a=+0.94$)
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$\beta_{2}(f_{\mathrm{pos,min}})$ $3.23\cdot 10^{-2}$ $2.42\cdot 10^{-2}$
$\mathbf{3.58\cdot 10^{-2}}$ $2.77\cdot 10^{-2}$
$\beta_{2}(f_{\mathrm{pos,max}})$ $\mathbf{3.93\cdot 10^{-2}}$ $2.72\cdot
10^{-2}$ $\mathbf{3.66\cdot 10^{-2}}$ $2.72\cdot 10^{-2}$

###  6.6 Faraday Effects

As linear polarization travels through a magnetized plasma, its EVPA is
rotated by Faraday rotation, interchanging Stokes $Q$ and $U$. Similarly,
Faraday conversion exchanges linear and circular polarization, interchanging
Stokes $U$ and $V$. Both of these effects can be significant in accreting
black hole systems. In particular, Faraday rotation is believed to be
extremely important for reducing the linear polarization fraction in models of
M87* to the observed values. Typically, SANEs have larger Faraday rotation and
conversion depths than MADs. This is largely because SANE models require
larger mass densities to match the observed flux of M87*. They also have lower
temperatures, which increases the efficiency of Faraday effects.

####  6.6.1 Faraday Rotation

Table 6 of fiducial model Faraday rotation depths shows a pronounced gap
between a marginal effect in MAD simulations relative to the corresponding
effect which is 3 orders of magnitude larger in SANE simulations. Varying
positron content even at the percent level leads to large EVPA rotational
swings for SANE plasmas due to the large absolute response of the Faraday
rotation measure to the increased fraction of positrons. This naturally leads
to a profoundly discriminating probe of plasma magnetic inflow properties in
regions of changing positron fraction. Even when the plasma composition is in
steady state, we may identify the rapid spatial variation of circular
polarization as a signature of high Faraday rotation depth characteristic of
SANE accretion flows.

Table 6: Faraday rotation depth at min and max fpos for fiducial models at
$T=25,000M$:

. $\mathrm{SANE}$ ($a=-0.5$) $\mathrm{MAD}$ ($a=-0.5$) $R\mathrm{-Beta}$
$R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\ Beta}$
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$\mathrm{w./\ Jet}$ $\tau_{V}(\mathrm{f_{pos,min}})$ $2.73\cdot 10^{3}$
$1.30\cdot 10^{3}$ $7.74\cdot 10^{3}$ $3.40\cdot 10^{3}$ 1.88 1.55 4.64 4.38
$\tau_{V}(\mathrm{f_{pos,max}})$ $5.92\cdot 10^{2}$ $2.01\cdot 10^{2}$
$1.52\cdot 10^{3}$ $6.19\cdot 10^{2}$ $3.31\cdot 10^{-2}$ $2.44\cdot 10^{-1}$
$9.93\cdot 10^{-1}$ $7.00\cdot 10^{-1}$ $\mathrm{MAD}$ ($a=+0.94$)
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$\tau_{V}(\mathrm{f_{pos,min}})$ $6.48\cdot 10^{-1}$ $4.78\cdot 10^{-1}$ 2.66
2.22 $\tau_{V}(\mathrm{f_{pos,max}})$ $1.16\cdot 10^{-1}$ $9.6\cdot 10^{-2}$
$5.25\cdot 10^{-1}$ $3.79\cdot 10^{-1}$

Table 7: Faraday conversion depth at min and max fpos for fiducial models at
$T=25,000M$:

. $\mathrm{SANE}$ ($a=-0.5$) $\mathrm{MAD}$ ($a=-0.5$) $R\mathrm{-Beta}$
$R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\ Beta}$
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$\mathrm{w./\ Jet}$ $\tau_{Q}(\mathrm{f_{pos,min}})$ $4.30\cdot 10^{0}$
$4.66\cdot 10^{0}$ $2.11\cdot 10^{0}$ $9.22\cdot 10^{-1}$ $3.88\cdot 10^{-2}$
$2.12\cdot 10^{-2}$ $2.84\cdot 10^{-2}$ $2.08\cdot 10^{-2}$
$\tau_{Q}(\mathrm{f_{pos,max}})$ $4.19\cdot 10^{0}$ $4.54\cdot 10^{0}$
$1.75\cdot 10^{0}$ $5.41\cdot 10^{-1}$ $3.54\cdot 10^{-2}$ $1.49\cdot 10^{-2}$
$3.22\cdot 10^{-2}$ $1.60\cdot 10^{-2}$ $\mathrm{MAD}$ ($a=+0.94$)
$R\mathrm{-Beta}$ $R\mathrm{-Beta}$ $\mathrm{Crit.\ Beta}$ $\mathrm{Crit.\
Beta}$ $\mathrm{w./\ Jet}$ $\mathrm{w./\ Jet}$
$\tau_{Q}(\mathrm{f_{pos,min}})$ $4.03\cdot 10^{-2}$ $2.06\cdot 10^{-2}$
$2.73\cdot 10^{-2}$ $1.84\cdot 10^{-2}$ $\tau_{Q}(\mathrm{f_{pos,max}})$
$3.77\cdot 10^{-2}$ $2.04\cdot 10^{-2}$ $2.90\cdot 10^{-2}$ $1.56\cdot
10^{-2}$

####  6.6.2 Faraday Conversion

Table 7 shows fiducial model Faraday conversion depths for SANEs are 1-2
orders greater than those for MADs. Faraday conversion depths tend to be lower
than Faraday rotation depths: by 3 orders of magnitude for SANEs and 1-2
orders for MADs. However, because Faraday conversion results in the direct
production of circular polarization (from linear), it may result in a
significant contribution of $V$. Faraday conversion can produce CP even in a
pure pair plasma as long as the magnetic field twists along the line of sight
(Wardle & Homan, 2003; Ricarte et al., 2021).

###  6.7 Frequency and Inclination Dependence

In Figs. 21 and 22, we search for extended structure at 86 GHz in the
R-$\beta$ with $\beta_{e0}=0.01$ Model in the MAD and SANE cases,
respectively. In our 86 GHz images, we use the same Munit that normalized the
230 GHz images to .5 Jy, though now we have more flux with a larger field of
view and shifting emitting regions at low frequency. The SANE Fig. 22 in
particular shows an upwardly extended feature reminiscent of the M87 jet base
in Lu et al. (2023).

Changing observer inclination induces further distinctive image morphological
features, as shown in the 86 GHz maps in Fig. 23 with inclination angle
$40^{\circ}$ (instead of the default M87 inclination $17^{\circ}$ used
throughout this work). In Fig. 23, the $R-\beta$ model with $\beta_{e0}=0.01$
jet has image plane projection horizontally elongated in the MAD case and
vertically elongated in the SANE case relative to the default orientation. The
jet collimation profile generally broadens as the jet inclination tilts away
from the line of sight; the broader jet is more like observed in Lu et al.
(2023). More edge-on morphologies are expected to break degeneracies in face-
on images due to general relativistic strong lensing effects.

Figure 21: For the MAD $a=-0.5$ at T=25,000M at 86 GHz for the R-Beta model
with $\beta_{e0}=0.01$ jet: Top panel: Ionic plasma (fPos = 0); Bottom panel:
fPos = 50.

Figure 22: For the SANE $a=-0.5$ at T=25,000M at 86 GHz for the R-Beta model
with $\beta_{e0}=0.01$ jet: Top panel: Ionic plasma (fPos = 0); Bottom panel:
fPos = 50.

Figure 23: For $a=-0.5$ at $T=25,000M$ and 86 GHz for the R-Beta model with
$\beta_{e0}=0.01$ jet: a MAD (Top) vs. SANE (Bottom) comparison now with
$40^{\circ}$ inclination.

###  6.8 Temporal Evolution

Figs. 24 and 25 show temporal variation of the GRMHD simulations for the MAD
and SANE cases, respectively. Both intensity and polarization morphology vary
noticeably over timescales of thousands of gravitational radii (years), as
expected from observations Wielgus & et al. (2020) in which bright spots
appear at various azimuths throughout the M87 emitting ring. The flux eruption
feature outside the photon ring at $T=25,000M$ is partially replaced by a
spiral arm at $T=20,000M$ and a smaller extrusion at $T=30,000M$ in the
circular polarization maps.

Figure 24: For the $a=-0.5$ MAD at $T=25,000M$: R Beta at 230 GHz without
positrons at $T=20,000M$ (Top), $T=25,000M$ (Middle) and $T=30,000M$ (Bottom).
Figure 25: For the $a=-0.5$ SANE at $T=25,000M$: R Beta at 230 GHz without
positrons at $T=20,000M$ (Top), $T=25,000M$ (Middle) and $T=30,000M$ (Bottom).

Figure 26 compares the evolution of circular polarization fraction with
positron fraction at fiducial time $T=25,000M$ with that of a later snapshot
of the simulation at $30,000M$. Fifty-one different positron fractions are
used in the series of frames producing these curves representing the variation
of $V/I$ with positron fraction. The fiducial time with a prominent flux
ejection loop has $V/I$ monotonically going from the most negatively polarized
to approaching 1/3 this value linearly in the fraction of unpaired emitters
(slightly slower than linearly around
$\frac{(n_{e})_{0}}{(n_{e})_{0}+2n_{\mathrm{pairs}}}=1/3$ where the plasma is
an even mix of electrons, positrons and protons).

Figure 26:  Degree of circular polarization as a function of
$\frac{n_{\mathrm{pairs}}}{(n_{e})_{0}}$ (Top Panel) and unpaired synchrotron
emitter fraction $\frac{(n_{e})_{0}}{(n_{e})_{0}+2n_{\mathrm{pairs}}}$ (Bottom
Panel). The fiducial time is $T=25,000M$ and has characteristically different
monotonicity behavior $5,000M$ after.

The flux loop occurrence at fiducial time $T=25,000M$ of the MAD $a=-0.5$
simulation may be representative of a broader episodic phenomenon occurring
throughout the evolution of the flow. The Fig 27 time series of mass accretion
rate $\dot{m}$ and horizon-threading flux $\phi$ from $10,000M<T<30,000M$
reveals at $T=25,000M$ a sharp rise in $\phi$ accompanied by a sharp decrease
in $\dot{m}$, which accords with the flux eruption scenario where a highly
polarized magnetic flux loop is added to a magnetically arrested disk. Similar
loop morphologies were observed for $T=17,730M$ and $27,110M$ where the
simulation time series have similar peaks in $\phi$ and troughs in $\dot{m}$
as at $T=25,000M$.

Figure 27: Time series of fiducial MAD $a=-0.5$ simulation horizon threading
flux $\phi$ (Yellow) and mass accretion rate $\dot{m}$ (Blue) for
$10000M<T<30000M$. The vertical lines at $T=17,730M,25,000M$ and $27,110M$
correspond to sharp peaks in $\phi$, troughs in $\dot{m}$ and flux eruption
loops.

##  7 Discussion and Conclusions

There have been remarkable advances in imaging and simulating AGN jets over
the past couple of decades. Despite this progress there are potentially vital
components — the jet composition and relativistic particle acceleration
remains which remain controversial. Our methodology to address these is to
focus on one well-studied source, M87, and one region of the electromagnetic
spectrum, radio, millimeter and submillimeter, and to incorporate different
phenomenological prescriptions to bridge this divide into the simulations and
then “observe” them. The actual observations, especially those from the Event
Horizon Telescope, can then be used to discern empirically some of the rules
that govern jet formation, collimation, polarization and dissipation. This
approach can be extended using more sources, frequencies and simulations and
statistical comparisons can also be conducted. These extensions will be
discussed in future publications including the completion of this series on
Sgr A*, M87 and 3C 279.

The GRMHD model that we have used to develop a more generally applicable set
of techniques is quite specific in terms of spin ($a/M_{H}=-0.5,0.94$) and
disposition of the surrounding gas (dense orbiting torus with non-relativistic
wind at high latitude outside the jet). The magnetic flux density and polarity
were consequences of the conditions of the simulations. Given these boundary
conditions, the concentration of horizon-crossing magnetic flux and the
formation of an electromagnetic outflow or jet are inevitable.

Within the Bondi radius ($\sim 10^{5}M$), the jet profile is roughly
parabolic, consistent with other simulations, e.g., Penna et al. (2013). The
form of the flux and velocity variation across the jet should also be
reasonably generic, though the stability properties and entrainment at the jet
surface is probably sensitive to the numerical details. In conclusion, we
should have a pretty representative suite of simulations of the flow and the
field to link to the highest resolution mm-observations.

The Observing JAB simulations methodology reproduces a surprising number of
observed signatures of M87 morphology and dynamics. Starting with turbulent
heating models including that used by the EHT, R-$\beta$, and the Critical
$\beta$ model of Anantua et al. (2020)a, we have the expected ring-like global
structure for intensity and EVPAs strongest in local maxima of intensity.
Adding equipartition-inspired models, the jet magnetic substructure for
Constant Electron Beta models characterized by constant electron-gas-to-
magnetic pressure along the jet gives a more broadly distributed profile.

In the case of M87, the radio emissivity is not simply a function of the gas
and/or the magnetic pressure. So the rule for particle acceleration must
depend upon other factors (e.g. $\beta$ and $\Gamma{\bf v}$). We have
implemented models where the emissivity is governed by total plasma $\beta$ in
the turbulent plasma and by conversion of magnetic-to-particle energy
(parametrized by the contribution $\beta_{e}$ of radiating electrons and
positrons) in the relativistic jet.

Our models also go beyond what is currently directly observable in simulating
the effects of incrementally changing positron fraction; however, SANE and MAD
produce a sharp enough dichotomy to currently be distinguished. The key
finding is that polarization is a sharp cleaver distinguishing SANE and MAD
accretion flows. Particularly, we find distinct polarized emission signatures
that depend on the positrons content in radically different ways for SANE and
MAD simulations.

In summary, the primary findings of the observing simulations methodology
applied to M87 include:

  * •

Both $R-\beta$ and Critical $\beta$ turbulent heating models produce ring-like
intensity profiles with some MAD cases satisfying linear polarization
constraints and all satisfying preliminary circular polarization upper bounds.

  * •

The piecewise addition of a Constant $\beta_{e}$ jet tends to broader annular
emission profiles.

  * •

MAD and SANE images with polarization at constant overall flux have markedly
different morphological properties. The MAD can exhibit a prominent flux
eruption in intensity and linear polazrization.

  * •

The Faraday depths of the SANE are 2-3 orders greater than for MAD. The SANE
linear polarization is more disordered and circular polarization structure is
completely scrambled.

  * •

The circular polarization degree for MAD maps dominated by intrinsic $V/I$
exhibit a linear vanishing of $V/I$ in the fraction of paired emitters.

The AGN environment is certainly a messy and chaotic one– replete with winds,
gas, dust and molecular clouds to name a few. The task of emission modeling
jet/accretion flow/black hole systems in such an uncertain setting, on the
other hand, is a fertile wonderland for the creation of theoretical models and
the discovery of new phenomenology. With few constraints on black hole spin or
jet composition, vast libraries of GRMHD libraries remain viable for even the
most well studied sources like M87. The “Observing” JAB Simulations approach
embraces this uncertainty by using piecewise models and generic plasma
compositions to allow for complex interactions leading to unexpected results
such as the positron-mediated Faraday effects leading to the sharp SANE-MAD
dichotomy in polarization signatures illustrated in this work. The present
application leaves us not only closer characterizing M87 as polarized MAD flow
near horizon scales, but also to narrow the possible plasma descriptions for
other JAB systems such as the jetted AGN 3C 279 that will be the third work of
this series– and the vast universe of future horizons to be discovered.

##  8 Future Directions

With our suite of turbulent and sub-equipartition heating models with
positrons, we have taken a key step in bridging rapidly advancing GRMHD
simulations and observations. The stark SANE-MAD dichotomy found in polarized
intensity spatial distribution and time evolution presents a key opportunity
to rule out SANE models of M87 by comparing variability, e.g., EVPA rotation
rate, with the results of M87 2017 combining later observing campaigns.

It has been demonstrated that prescriptions involving dissipation as a
function of effective magnetic field $B_{e}=\mathcal{D}|\vec{n}\times\vec{B}|$
exhibit violation of bilateral symmetry across the jet axis both in the
stationary, axisymmetric, self-similar semi-analytic model (Blandford2017;
Anantua et al., 2020)a (with general relativistic ray tracing in Emami et al.
(2021)), and in the time-dependent 3D GRMHD simulation in Anantua et al.
(2018). Though barely visible in M87 observations, e.g., at 86 GHz Kim et al.
(2018a), this is predicted to be a robust – albeit, generic – observation for
EHT, with details depending on whether EHT sees a jet or disk-jet in the inner
few gravitational radii from the hole. Signs of this bilateral asymmetry from
“Observing” JAB Simulations have appeared in 230 GHz EHT observations of 3C
279 (Kim et al., 2020). We may implement prescriptions in $B_{e}$ in future
emission modeling work to reproduce bilateral asymmetry– particularly for 3C
279.

FRI disk wind momentum carries jet, while FRII jet momentum carries the wind.
Another way jets exchange momentum with their surroundings is through
currents. We can apply the current density model (Anantua et al., 2018) to
investigate whether current sheets account for limb-brightening past 100$M$.
The $B$-field alone struggles to remain toroidal past $100M$ unless it’s
replenished by the disk. In addition to currents, we may systematically
associate the dissipation in JAB systems with a number of plausible physical
mechanisms, such as Shakura-Sunyaev momentum transport and Newtonian shear as
outlined in the Appendix. These phenomenological models give firm theoretical
intuition behind the physical mechanisms powering jets.

In future work, we will also incorporate positrons in a broader range of
emission models. We may use positron production rates from Wong et al. (2021)
to evolve the local positron fraction– a key advance over the single-positron-
ratio maps used here. The computational expense of a three-fluid
($e^{-},e^{+},p$) simulation may be mitigated by spatial symmetry and temporal
stationarity of some simulated flows.

A key feature of the Observing JAB simulations exercise presented here is its
generality: a simulation of a general relativistic magnetohydrodynamic flow
onto a compact object is flexible enough to model jets from neutron stars,
black hole/x-ray binaries and AGN alike. In this work, we started with a suite
of simulations fairly representative of an AGN in that it exhibited the
commonly occurring combination of a thick ion torus confining electromagnetic
flux from a polar outflow from a rotating black hole, then fine-tuned it to
M87 to emulate its observed JAB system polarized substructure. Disk emission
has been emphasized in other  work, starting with Sgr A* at our Galactic
Center, replete with new near horizon observations of photon rings courtesy of
EHT. Our models will also be applied to near-horizon emission in future EHT
observational targets such as the highly variable quasar 3C 279 in the last
work of this trilogy.

## Acknowledgments

Richard Jude Anantua was supported by the California Alliance at the outset of
this investigation and the Oak Ridge Associated Universities Powe Award for
Junior Faculty Enhancement and Simons Collaboration on Extreme Electrodynamics
of Compact Sources towards the end. Roman Shcherbakov and Alexander
Tchekhovskoy have provided excellent guidance and mentorship at the beginning
of this investigation. UTSA undergraduate Noah Heredia has been helpful
through graphic-related activities. Daniel Palumbo provided observational
guidance. Angelo Ricarte was supported by the Black Hole Initiative at Harvard
University, made possible through the support of grants from the Gordon and
Betty Moore Foundation and the John Templeton Foundation. The opinions
expressed in this publication are those of the author(s) and do not
necessarily reflect the views of the Moore or Templeton Foundations. Razieh
Emami acknowledges the support by the Institute for Theory and Computation at
the Center for Astrophysics as well as grant numbers 21-atp21-0077, NSF
AST-1816420 and HST-GO-16173.001-A.

##  9 Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abdo et al. (2009) Abdo A. A. e. a., 2009, The Astrophysical Journal, 707, 55
  * Aharonian (2006) Aharonian F. e. a., 2006, Science, 314, 1424
  * Akiyama et al. (2015) Akiyama K., Lu R.-S., Fish V. L., Doeleman S. S., Broderick A. E., Dexter J., Hada K., Kino M., Nagai H., Honma M., Johnson M. D., e. a., 2015, The Astrophysical Journal, 807, 150
  * Aloy et al. (2015) Aloy M. A., Borovikov S., Pogorelov N., 2015, in Pogorelov N. V., Audit E., Zank G. P., eds, Numerical Modeling of Space Plasma Flows ASTRONUM-2014 Vol. 498 of Astronomical Society of the Pacific Conference Series, A New RMHD Code for Astrophysical Applications. p. 185
  * Anantua et al. (2018) Anantua R., Blandford R., Tchekhovskoy A., 2018, Galaxies, 6, 31
  * Anantua et al. (2023) Anantua R., Dúran J., Ngata N., Oramas L., Röder J., Emami R., Ricarte A., Curd B., Broderick A. E., Wayland J., Wong G. N., Ressler S., Nigam N., Durodola E., 2023, Galaxies, 11, 4
  * Anantua et al. (2020) Anantua R., Emami R., Loeb A., Chael A., 2020, The Astrophysical Journal, 896, 30
  * Anantua et al. (2020) Anantua R., Ressler S., Quataert E., 2020, Monthly Notices of the Royal Astronomical Society, 493, 1404
  * Anantua (2016) Anantua R. J., 2016, PhD thesis, Stanford University, California
  * Asada & Nakamura (2012) Asada K., Nakamura M., 2012, The Astrophysical Journal Letters, 745, L28
  * Beskin et al. (1992) Beskin V. S., Istomin Y. N., Parev V. I., 1992, Soviet Astronomy, 36, 642
  * Biretta et al. (1999) Biretta J. A., Sparks W. B., Macchetto F., 1999, The Astrophysical Journal, 520, 621
  * Biretta et al. (1995) Biretta J. A., Zhou F., Owen F. N., 1995, The Astrophysical Journal, 447, 582
  * Blakeslee et al. (2009) Blakeslee J. P., Jordán A., Mei S., Côté P., Ferrarese L., Infante L., Peng E. W., Tonry J. L., West M. J., 2009, The Astrophysical Journal, 694, 556
  * Blandford et al. (2017) Blandford R., Yuan Y., Hoshino M., Sironi L., 2017, Space Science Reviews, 207, 291
  * Blandford & Begelman (1999) Blandford R. D., Begelman M. C., 1999, Monthly Notices of the Royal Astronomical Society, 303, L1
  * Blandford & Levinson (1995) Blandford R. D., Levinson A., 1995, The Astrophysical Journal, 441, 79
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, Monthly Notices of the Royal Astronomical Society, 199, 883
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, Monthly Notices of the Royal Astronomical Society, 179, 433
  * Bolton et al. (1949) Bolton J. G., Stanley G. J., Slee O. B., 1949, Nature, 164, 101
  * Broderick & Loeb (2009) Broderick A. E., Loeb A., 2009, The Astrophysical Journal, 697, 1164
  * Broderick et al. (2015) Broderick A. E., Narayan R., Kormendy J., Perlman E. S., Rieke M. J., Doeleman S. S., 2015, The Astrophysical Journal Letters, 805, 179
  * Broderick & Tchekhovskoy (2015) Broderick A. E., Tchekhovskoy A., 2015, The Astrophysical Journal, 809, 97
  * Christian & Loeb (2015) Christian P., Loeb A., 2015, Phys. Rev. D, 91, 101301
  * Churazov et al. (2001) Churazov E., Brüggen M., Kaiser C. R., Böhringer H., Forman W., 2001, The Astrophysical Journal, 554, 261
  * Curtis (1918) Curtis H. D., 1918, Publications of Lick Observatory, 13, 9
  * de Gasperin et al. (2012) de Gasperin F., Orrú E., Murgia M., Merloni A., Falcke H., e. a., 2012, Astronomy and Astrophysics, 547, A56
  * Dexter et al. (2012) Dexter J., Agol E., Fragile P. C., McKinney J. C., 2012, Journal of Physics Conference Series, 372, 012023
  * Dexter et al. (2012) Dexter J., McKinney J. C., Agol E., 2012, Monthly Notices of the Royal Astronomical Society, 421, 1517
  * Doeleman et al. (2012) Doeleman S. S. e. a., 2012, Science, 338, 355
  * Emami et al. (2021) Emami R., Anantua R., Chael A. A., Loeb A., 2021, Astrophys. J., 923, 272
  * Event Horizon Telescope Collaboration et al. (2022) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Algaba J. C., Anantua R., Asada K., Azulay R., Bach U., Baczko A.-K., Ball D., Baloković M., et al. 2022, The Astrophysical Journal Letters, 930, L12
  * Event Horizon Telescope Collaboration et al. (2022) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Algaba J. C., et al. 2022, The Astrophysical Journal Letters, 930, L16
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., Azulay R., Baczko A.-K., Ball D., Baloković M., Barrett J., et al. 2019a, The Astrophysical Journal Letters, 875, L1
  * Event Horizon Telescope Collaboration et al. (2019b) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., Azulay R., Baczko A.-K., Ball D., Baloković M., Barrett J., et al. 2019b, The Astrophysical Journal Letters, 875, L2
  * Event Horizon Telescope Collaboration et al. (2019c) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., Azulay R., Baczko A.-K., Ball D., Baloković M., Barrett J., et al. 2019c, The Astrophysical Journal Letters, 875, L3
  * Event Horizon Telescope Collaboration et al. (2019d) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., Azulay R., Baczko A.-K., Ball D., Baloković M., Barrett J., et al. 2019d, The Astrophysical Journal Letters, 875, L4
  * Event Horizon Telescope Collaboration et al. (2019e) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., Azulay R., Baczko A.-K., Ball D., Baloković M., Barrett J., et al. 2019e, The Astrophysical Journal Letters, 875, L5
  * Event Horizon Telescope Collaboration et al. (2019f) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., Azulay R., Baczko A.-K., Ball D., Baloković M., Barrett J., et al. 2019f, The Astrophysical Journal Letters, 875, L6
  * Event Horizon Telescope Collaboration et al. (2023) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., Azulay R., Baczko A.-K., Ball D., Baloković M., Barrett J., et al. 2023, submitted to The Astrophysical Journal Letters, p. (Paper VII)
  * Event Horizon Telescope Collaboration et al. (2021a) Event Horizon Telescope Collaboration Akiyama K., Algaba J. C., Alberdi A., Alef W., Anantua R., Asada K., Azulay R., Baczko A.-K., et al. 2021a, The Astrophysical Journal Letters, 910, L12
  * Event Horizon Telescope Collaboration et al. (2021b) Event Horizon Telescope Collaboration Akiyama K., Algaba J. C., Alberdi A., Alef W., Anantua R., Asada K., Azulay R., Baczko A.-K., et al. 2021b, The Astrophysical Journal Letters, 910, L13
  * Farris et al. (2012) Farris B. D., Gold R., Paschalidis V., Etienne Z. B., Shapiro S. L., 2012, Physical Review Letters, 109, 221102
  * Feng & Wu (2017) Feng J., Wu Q., 2017, Monthly Notices of the Royal Astronomical Society, 470, 612
  * Fish et al. (2016) Fish V., Akiyama K., Bouman K., Chael A., Johnson M., Doeleman S., Blackburn L., Wardle J., Freeman W., 2016, Galaxies, 4, 54
  * Forman et al. (2017) Forman W., Churazov E., Jones C., Heinz S., Kraft R., Vikhlinin A., 2017, The Astrophysical Journal, 844, 122
  * Fragile & Meier (2009) Fragile P. C., Meier D. L., 2009, The Astrophysical Journal, 693, 771–783
  * Gammie et al. (2003) Gammie C. F., McKinney J. C., Tóth G., 2003, The Astrophysical Journal, 589, 444
  * Gebhardt et al. (2011) Gebhardt K., Adams J., Richstone D., Lauer T. R., Faber S. M., Gültekin K., Murphy J., Tremaine S., 2011, The Astrophysical Journal, 729, 119
  * Ghisellini (2012) Ghisellini G., 2012, Monthly Notices of the Royal Astronomical Society: Letters, 424, L26–L30
  * Globus & Levinson (2016) Globus N., Levinson A., 2016, Monthly Notices of the Royal Academy of Sciences, 461, 2605
  * Goddi et al. (2021) Goddi C., Martí-Vidal I., Messias H., Bower G. C., Broderick A. E., Dexter J., Marrone D. P., Moscibrodzka M., Nagai H., Algaba J. C., Asada K., Crew G. B., Gómez J. L., Impellizzeri C. M. V., Janssen M., Kadler M., Krichbaum T. P., Lico R., Matthews L. D., Nathanail A., Ricarte A., Ros E., Younsi Z., et al. 2021, The Astrophysical Journal Letters, 910, L14
  * Gold et al. (2016) Gold R., McKinney J. C., Johnson M. D., Doeleman S. S., 2016, ArXiv e-prints
  * Goldreich & Julian (1969) Goldreich P., Julian W. H., 1969, The Astrophysical Journal, 157, 869
  * Hada et al. (2016) Hada K., Kino M., Doi A., Nagai H., Honma M., Akiyama K., Tazaki F., Lico R., Giroletti M., Giovannini G., Orienti M., Hagiwara Y., 2016, The Astrophysical Journal, 817, 131
  * Hada et al. (2017) Hada K., Park J. H., Kino M., Niinuma K., et al 2017, ArXiv e-prints
  * Howes (2010) Howes G. G., 2010, Monthly Notices of the Royal Astronomical Society: Letters, 409, L104–L108
  * Junor & Biretta (1995) Junor W., Biretta J. A., 1995, Astronom. J., 109, 500
  * Kim et al. (2020) Kim J.-Y., Krichbaum Event Horizon Telescope Collaboration 2020, Astronomy & Astrophysics, 640, A69
  * Kim et al. (2018a) Kim J.-Y., Krichbaum T. P., Lu R.-S., Ros E., Bach U., Bremer M., de Vicente P., Lindqvist M., Zensus J. A., 2018a, ArXiv e-prints
  * Kim et al. (2018b) Kim J. Y., Krichbaum T. P., Lu R. S., Ros E., Bach U., Bremer M., de Vicente P., Lindqvist M., Zensus J. A., 2018b, Astronomy & Astrophysics, 616, A188
  * Kimura & Toma (2020) Kimura S. S., Toma K., 2020, The Astrophysical Journal, 905, 178
  * Koide et al. (2000) Koide S., Meier D. L., Shibata K., Kudoh T., 2000, The Astrophysical Journal, 536, 668
  * Kovalev et al. (2007) Kovalev Y. Y., Lister M. L., Homan D. C., Kellermann K. I., 2007, The Astrophysical Journal Letters, 668, L27
  * Laurent & Titarchuk (2018) Laurent P., Titarchuk L., 2018, The Astrophysical Journal, 859, 89
  * Li et al. (2009) Li Y.-R., Yuan Y.-F., Wang J.-M., Wang J.-C., Zhang S., 2009, The Astrophysical Journal, 699, 513
  * Longair (2011) Longair M. S., 2011, High Energy Astrophysics
  * Lu et al. (2023) Lu R.-S., Asada K., Krichbaum T. P., Park J., Tazaki F., Pu H.-Y., Nakamura M., Lobanov A., Hada K., Akiyama K., Kim J.-Y., Marti-Vidal I., Gómez J. L., Kawashima T., Yuan F., Ros E., et al. 2023, arXiv e-prints, p. arXiv:2304.13252
  * Lu et al. (2014) Lu R.-S., Broderick A. E., Baron F., Monnier J. D., Fish V. L., Doeleman S. S., Pankratius V., 2014, The Astrophysical Journal, 788, 120
  * Marscher & Jorstad (2010) Marscher A. P., Jorstad S. G., 2010, ArXiv e-prints
  * McKinney & Blandford (2009) McKinney J. C., Blandford R. D., 2009, Monthly Notices of the Royal Astronomical Society, 394, L126
  * McKinney et al. (2012) McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, Monthly Notices of the Royal Astronomical Society, 423, 3083
  * Mertens et al. (2016) Mertens F., Lobanov A. P., Walker R. C., Hardee P. E., 2016, Astronomy & Astrophysics, 595, A54
  * Mishra et al. (2016) Mishra B., Fragile P. C., Johnson L. C., Kluźniak W., 2016, Monthly Notices of the Royal Academy of Sciences, 463, 3437
  * Moscibrodzka et al. (2016) Moscibrodzka M., Falcke H., Shiokawa H., 2016, Astronomy & Astrophysics, 586, A38
  * Mościbrodzka & Gammie (2018) Mościbrodzka M., Gammie C. F., 2018, Monthly Notices of the Royal Astronomical Society, 475, 43
  * Mościbrodzka et al. (2011) Mościbrodzka M., Gammie C. F., Dolence J. C., Shiokawa H., 2011, The Astrophysical Journal, 735, 9
  * Narayan & Yi (1994) Narayan R., Yi I., 1994, The Astrophysical Journal Letters, 428, L13
  * Owen et al. (2000) Owen F. N., Eilek J. A., Kassim N. E., 2000, The Astrophysical Journal, 543, 611
  * Owen et al. (1989) Owen F. N., Hardee P. E., Cornwell T. J., 1989, The Astrophysical Journal, 340, 698
  * Parfrey et al. (2019) Parfrey K., Philippov A., Cerutti B., 2019, Phys. Rev. Lett., 122, 035101
  * Pasetto et al. (2021) Pasetto A., Carrasco-González C., Gómez J. L., Martí J.-M., Perucho M., O’Sullivan S. P., Anderson C., Díaz-González D. J., Fuentes A., Wardle J., 2021, The Astrophysical Journal Letters, 923, L5
  * Penna et al. (2013) Penna R. F., Narayan R., Sa̧dowski A., 2013, Monthly Notices of the Royal Academy of Sciences, 436, 3741
  * Perlman et al. (2001) Perlman E. S., Biretta J. A., Sparks W. B., Macchetto F. D., Leahy J. P., 2001, The Astrophysical Journal, 551, 206
  * Porth et al. (2016) Porth O., Olivares H., Mizuno Y., Younsi Z., Rezzolla L., Moscibrodzka M., Falcke H., Kramer M., 2016, ArXiv e-prints
  * Prather et al. (2021) Prather B., Wong G., Dhruv V., Ryan B., Dolence J., Ressler S., Gammie C., 2021, The Journal of Open Source Software, 6, 3336
  * Prieto et al. (2016) Prieto M. A., Fernández-Ontiveros J. A., Markoff S., Espada D., González-Martín O., 2016, Monthly Notices of the Royal Academy of Sciences, 457, 3801
  * Quataert & Gruzinov (1999) Quataert E., Gruzinov A., 1999, The Astrophysical Journal, 520, 248
  * Ricarte et al. (2021) Ricarte A., Qiu R., Narayan R., 2021, Monthly Notices of the Royal Astronomical Society, 505, 523
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, Astronomy and Astrophysics, 24, 337
  * Soltan (1982) Soltan A., 1982, Monthly Notices of the Royal Astronomical Society, 200, 115
  * Walker et al. (2016) Walker R. C., Hardee P. E., Davies F., Ly C., Junor W., Mertens F., Lobanov A., 2016, Galaxies, 4, 46
  * Walker et al. (2008) Walker R. C., Ly C., Junor W., Hardee P. J., 2008, Journal of Physics Conference Series, 131, 012053
  * Walsh et al. (2013) Walsh J. L., Barth A. J., Ho L. C., Sarzi M., 2013, The Astrophysical Journal, 770, 86
  * Wang & Zhou (2009) Wang C.-C., Zhou H.-Y., 2009, Monthly Notices of the Royal Astronomical Society, 395, 301
  * Wardle & Homan (2003) Wardle J. F. C., Homan D. C., 2003, Astrophysics and Space Science, 288, 143
  * Wielgus & et al. (2020) Wielgus M., et al. 2020, The Astrophysical Journal, 901, 67
  * Wilson & Yang (2002) Wilson A. S., Yang Y., 2002, The Astrophysical Journal, 568, 133
  * Wong et al. (2022) Wong G. N., Prather B. S., Dhruv V., Ryan B. R., Mościbrodzka M., kwan Chan C., Joshi A. V., Yarza R., Ricarte A., Shiokawa H., Dolence J. C., Noble S. C., McKinney J. C., Gammie C. F., 2022, The Astrophysical Journal Supplement Series, 259, 64
  * Wong et al. (2021) Wong G. N., Ryan B. R., Gammie C. F., 2021, The Astrophysical Journal, 907, 73
  * Yao et al. (2021) Yao P. Z., Dexter J., Chen A. Y., Ryan B. R., Wong G. N., 2021, Monthly Notices of the Royal Astronomical Society, 507, 4864
  * Yoon et al. (2020) Yoon D., Chatterjee K., Markoff S. B., van Eijnatten D., Younsi Z., Liska M., Tchekhovskoy A., 2020, Monthly Notices of the Royal Astronomical Society, 499, 3178

## Appendix A. Alternative Phenomenological Emission Models

To extend the range of phenomenology covered by emission models in this work,
we add some of the models introduced in Anantua (2016).These models are in
terms of the partial pressure due to electrons that contribute to emission
near the observed frequency.

### A.1 Magnetic Bias Model

The calculated intensity in simple constant electron $\beta$ models typically
declines with radius far faster than the observations. This can be addressed
by stipulating that $\beta\propto P^{-N}$, where the exponent $N$ is typically
in the range $\sim 0.25-0.5$ and the constant of proportionality can be set at
a fiducial radius.

We can incorporate these considerations in a Magnetic Bias Model by scaling
$P_{e}\propto P_{\rm mag}^{n}\propto b^{2n}$, which near the black hole
reduces to $P_{e}\propto P^{n}\propto P^{1-N}$. In this regime, Bias Models
for $n=0$, i.e., $N=1$, have constant pressure along the jet.

### A.2 $\alpha$ Model

A more detailed approach to setting the partial pressure is to suppose that
the free energy in the jet is dissipated at a specific rate $W^{\prime}$ in
the comoving frame and that the associated energy density is proportional to
$u_{e}$ $\sim 0.1W^{\prime}t_{\rm loss}^{\prime}$ where $t_{\rm loss}$ is the
smaller of the radiative cooling time and the expansion time. There are many
possible models for the shear stress but, as with accretion disks ((Shakura &
Sunyaev, 1973)), it suffices to assume that the shear stress $\tau$ is a fixed
fraction ($\alpha$) of the total pressure, $P=P_{\rm mag}+P_{g}+...$. Close to
the hole $P\sim P_{\rm mag}$. The dominant component of the rate of velocity
shear evaluated in the comoving frame is
$S^{\prime}\sim(1-v_{z}^{2})^{-1}|dv_{z}/ds|$ so that
$W^{\prime}=\frac{1}{2}\tau^{\prime}S^{\prime}$.

### A.3 Shear or Velocity Gradient ($S^{2}$) Model

The shear, or velocity gradient model, which we shall also term $S^{2}$-model
due to its dependence on $S$ both via shear stress and velocity shear, is
based on the Newtonian model for viscosity using a characteristic length scale
$L_{S}$ and setting  kinematic viscosity $\nu^{\prime}=cL/3$. The dynamic
viscosity is then obtained by multiplying $\nu^{\prime}$ by density (or energy
density divided by $c^{2}$)

$\mu^{\prime}=\frac{L_{S}}{3c}\sqrt{\left(\rho
c^{2}+\frac{b_{\mu}b^{\mu}}{2\mu_{0}}\right)\left(\frac{u_{g}}{3}+\frac{b_{\mu}b^{\mu}}{2\mu_{0}}\right)}$
(22)

Since the shear stress $\tau^{\prime}=\mu^{\prime}S^{\prime}$, the dissipation
rate in this model is then quadratic in the shear stress:
$W^{\prime}=\frac{1}{2}S^{\prime}\tau^{\prime}=\frac{1}{2}\mu^{\prime}S^{\prime
2}$, and
$\tilde{P}_{e}=W^{\prime}\min[t_{\mathrm{exp}}^{\prime},t_{\mathrm{rad}}^{\prime}]$
as above.

$\mathbf{A.4\ Summary\ of\ Pressure\mathrm{-}to\mathrm{-}Dissipation\
Prescriptions}$

We now write complete formulas for the partial pressures in our models with
explicit dissipation. Putting the parts of the $\alpha$, $S^{2}$ and $j^{2}$
model all together

$\begin{split}P_{\mathrm{e,\ \alpha\
model}}=\frac{1}{2}\alpha\left(\frac{b_{\mu}b^{\mu}}{2\mu_{0}}\right)\left(\gamma^{2}|\frac{dv_{z}}{ds}|\right)t^{\prime}\hskip
113.81102pt\\\ P_{\mathrm{e,\ S^{2}\ model}}=\frac{L}{6c}\sqrt{\left(\rho
c^{2}+\frac{b_{\mu}b^{\mu}}{2\mu_{0}}\right)\left(\frac{u_{g}}{3}+\frac{b_{\mu}b^{\mu}}{2\mu_{0}}\right)}\left(\gamma^{2}\frac{dv_{z}}{ds}\right)^{2}t^{\prime}\hskip
56.9055pt\\\ P_{\mathrm{e,\ j^{2}\
model}}=\mu_{0}cL|j_{\mu}j^{\mu}|t^{\prime}\mathrm{,\ \ where\ }\hskip
142.26378pt\\\
t^{\prime}=\mathrm{min}\left\\{\frac{1}{\gamma|\gamma\vec{\triangledown}\cdot\vec{v}+\vec{v}\cdot\vec{\triangledown}\gamma|},\frac{\mu_{0}c}{\sigma_{T}}\sqrt{\frac{3em_{e}}{2\pi}}\mathcal{D}^{-1}B_{e}^{-2/3}\nu^{-1/2}\right\\}\hskip
56.9055pt\end{split}$ (23)

(note the $\alpha$ model is assumed to be dominated by magnetic pressure).

Numerical implemetation of these modes is facilitated by expressing the
$s-$derivative in velocity shear as

$\left(\frac{dv_{z}}{ds}\right)^{2}=\left(\frac{dv_{z}}{dx}\right)^{2}+\left(\frac{dv_{z}}{dy}\right)^{2}$
(24)

avoiding division by 0 at the axis  and expressing

$\vec{\triangledown}\gamma=\begin{pmatrix}\frac{\partial}{\partial x}\\\
\frac{\partial}{\partial y}\\\ \frac{\partial}{\partial
z}\end{pmatrix}\left(1-\frac{v^{2}}{c^{2}}\right)^{-\frac{1}{2}}=\frac{1}{c^{2}}\gamma^{3}v\begin{pmatrix}\frac{\partial
v}{\partial x}\\\ \frac{\partial v}{\partial y}\\\ \frac{\partial v}{\partial
z}\end{pmatrix}$ (25)

where the derivatives along the axial directions is now

$\frac{\partial v}{\partial
x_{i}}\sim\frac{v\begin{pmatrix}\vec{x}+x_{\mathrm{i\
Step}}\hat{e}_{i}\end{pmatrix}-v\begin{pmatrix}\vec{x}\end{pmatrix}}{x_{\mathrm{i\
Step}}}$ (26)

## Appendix B. Alternative Positron Distribution Models

### B.1 Positron Proxies by Plasma Variables

One strategy for estimating the local distribution of positrons is to relate
positron density to plasma variables that positrons may trace. The radiation
energy density, magnetic pressure, electron temperature and functions thereof
may play this role. Inspired by successful parametrizations of electron
temperature, we may take as a proxy for positron fraction
$\frac{n_{e^{+}}}{n_{e^{-}}}$

$f_{e^{+}}(\beta)=f_{\mathrm{max}}e^{-\frac{\beta}{\beta_{c}}}$ (27)

or

$R_{e^{+}}(\beta)=\frac{\beta^{2}}{1+\beta^{2}}f_{\mathrm{min}}+\frac{1}{1+\beta^{2}}f_{\mathrm{max}}$
(28)

so that

$P_{e}=P_{e^{-}}+P_{e^{+}}=\beta_{e0}(1+(f_{e^{+}}(\beta)\mathrm{\ OR\
}R_{e^{+}}(\beta)))P_{B}.$ (29)

### B.2 Positron Distribution from Geometry

We also propose a geometric model separating jet inner and outer (or funnel)
regions by a critical parabolic streamline such that

  * •

Inner jet is supplied by pairs from the stagnation surface

  * •

Outer funnel wall jet is supplied by the analytic prescription for pair
production given above

$\mathit{Inner-Outer\ Jet\ Model\ Precriptions}$

We specify inner- and outer-jet emission prescriptions here:

Broderick & Tchekhovskoy (2015) argue models in which pairs are produced near
the horizon and $u_{e\pm}\propto B^{2}$ overproduce core emission relative to
jet emission. We may use the spark gap estimate for lepton number density from
Broderick & Tchekhovskoy (2015) and $u_{e\pm}=\gamma_{e\pm}n_{e\pm}m_{e}c^{2}$
to construct inner jet emissivity

$j_{\nu,\mathrm{in}}=\begin{cases}u_{e\pm}P_{B}^{\frac{3+\alpha}{2}}\nu^{-\alpha},\
z<z_{\mathrm{max}},\ 0<\xi<\frac{1}{2}\xi_{\mathrm{max}}\\\
\beta_{\mathrm{e0}}P_{B}^{\frac{3+\alpha}{2}}\nu^{-\alpha},\ z\geq
z_{\mathrm{max}},\ 0<\xi<\frac{1}{2}\xi_{\mathrm{max}}\end{cases}$

where $z_{\mathrm{max}}\approx 10r_{g}$, the altitude of the stagnation
surface.

For the outer jet, we use the model of Mościbrodzka et al. (2011), where
synchrotron photons from a radiatively inefficient accretion flow collide with
the jet funnel wall, where they undergo photon annihilation via Breit-Wheeler
process. We may adopt this as a model for the outer jet source of emitting jet
electrons and positrons.

$j_{\nu}=\begin{cases}j_{\nu,\mathrm{in}},\
0<\xi<\frac{1}{2}\xi_{\mathrm{max}}\\\ \frac{\sqrt{2}\pi
e^{2}n_{e}\nu_{s}}{3cK_{2}(\Theta_{e}^{-1})}(X^{1/2}+2^{11/12}X^{1/6})^{2}e^{-X^{1/3}},\
\frac{1}{2}\xi_{\mathrm{max}}<\xi<\xi_{\mathrm{max}}\end{cases}$

, where $X\equiv\nu/\nu_{s},\nu_{s}=\frac{2}{9}\frac{eB}{2\pi
m_{e}c}\Theta_{e}^{2}\sin\theta$, $\xi=s^{2}/z$ and $\alpha=\frac{p-1}{2}$ .