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

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, plasma magnetizations and turbulent heating and equipartition-based emission prescriptions (and piecewise combinations thereof) is chosen to represent distinct possibilities for the M87 jet/accretion flow/black hole (JAB) system. Simulation jet morphology, polarization and variation are then"observed"and compared with real observations to infer the rules that govern the polarized emissivity. Our models support several possible spin/emission model/plasma composition combinations supplying the jet in M87, whose black hole shadow has been observed down to the photon ring at 230 GHz by the Event Horizon Telescope (EHT). Net linear polarization and circular polarization constraints favor magnetically arrested disk (MAD) models whereas resolved linear polarization favors standard and normal evolution (SANE) in our parameter space. We also show that some MAD cases 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.


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 ⋆ E-mail: richard.anantua@utsa.edu 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 ≡ GMH/c 2 , where MH is the mass of the hole).Structure has been identified down to ∼ 10M (Fish et al. 2016)) and beyond by the Event Horizon Telescope (EHT) project, which has made linearly polarized images with resolution limit ∼ 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 ∼ 0.1", ∼ 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 γ-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 gammaray 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 sam-ple 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 γ-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) 1 .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 timedependent 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.

OBSERVATIONS OF M87
Located at the heart of the Virgo Cluster d = 16.7±0.6Mpc 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.

Black Hole
We adopt a black hole mass of (6.5 ± 0.7) × 10 9 M⊙ (Event Horizon Telescope Collaboration et al. 2019a), corresponding to length, time, angular and energy scales 10 13 m, 9 hr, 4 µas and 1.2×10 57 J, respectively.The associated Eddington luminosity is ∼ 8 × 10 40 W. A lower bound of a > 0.2MM87 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.98MH.We first adopt an intermediate angular frequency of ΩH = 0.35/(GMH/c 3 ) of 10 −5 s −1 for M87's black hole.Using the relations JH = GMHa/c and ΩH = a r 2 + +a 2 , where r+ = 1 2 (rS + r 2 S − 4(JH /(MHc)) 2 ) is the radius of the outer Kerr horizon, the corresponding dimensionless black hole spin is a/MH = 0.94.We also consider lower spin prograde a/M = 0.5 (ΩH = 0.13/(GMM87/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 θ = 20 •  (Wang & Zhou 2009; Prieto et al. 2016).The extractable rotational energy and angular momentum are ∼ 2 × 10 56 J and ∼ 4 × 10 61 kg m 2 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., GMM87/c 2 , (GMM87/c 2 )/d, GMM87/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.

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.
• Event Horizon Telescope (EHT) observational data made with effective beam of size ∼ 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 ∼ 30M .They are also broadly consistent with there being a selfabsorbed innermost jet called the core with constant brightness temperature ∼ 2 × 10 10 K and flux density Sν ∼ 1 Jy for 3 GHz ≲ ν ≲ 300 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 ≲ z ≲ 7000M with the separation of the brightened edges roughly given by ∼ 6z 1/2 .At larger radii, the jet expansion is closer to linear.

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 where ξ = X 2 /Y .The observed radio-mm-submm isotropic power is dominated by mm observations and is ∼ 10 34 W (Prieto et al. 2016).A beaming corrected guess might be as high as ∼ 10 35 W.

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×10 9 M⊙ 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.

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 ∼ 80 pc from the nucleus, produces blobs that appear to move up to 6c on the observer plane (Biretta et al. 1999).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 bol ∼ 3 × 10 35 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 disk ⩽ 3.4•10 41 erg/s.At 10% efficiency, L = η ṁc 2 implies an upper bound to the mass accretion rate of 3.8 • 10 21 g/s = 6 × 10 −5 M⊙/yr.The twofold change in M87's observed bolometric luminosity from its quiescent state value to (L = 5.4×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.

X-ray Observations
M87 has been observed at X-ray wavelengths by Chandra (Wilson & Yang 2002).The equipartition magnetic field value for the knot HST-1 was found to be ∼ 3 × 10 −4 G (Owen et al. 1989).
The steady jet isotropic luminosity in 2-10 keV X-rays is ∼ 3 × 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 ∼ 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 ∼ 10 15 M⊙ 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 Ljet = 5 × 10 36 W, noting that M87 could be in a relatively dormant state right now.

Gamma-ray Observations
The 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).

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 ∼ 2 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 ∼ 10 37 W.

Summary
M87 is fairly extreme in many respects.It has the most massive black hole we can study in detail.Its disk luminosity is ≲ 3 × 10 −7 L Edd while it creates jets with total power ≳ 3 × 10 −4 L Edd which have been shown to be collimated on scales ≲ 10M .Most importantly, we have recently begun to learn much more about the region within ∼ 100M from EHT observations.This makes it an excellent source to model.

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.75rS < x 1 < 20rS (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 Tp and electron temperature Te.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.

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 shockcapturing 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 along with a no-monopoles constraint ∂i √ −gB i = 0. Here, the rest-mass density of the fluid is ρ, u µ is the fluid four-velocity, b µ is the magnetic induction four-vector, and  the magnetohydrodynamic stress-energy tensor is 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 γ via P = (γ − 1) u.The effects of spcaetime are accounted for in the usual way, with g = −detgµν , the determinant of the covariant metric and Γ 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, 000GMH/c 2 respectively.Each simulation was run from t = 0 GMH/c 3 until 30, 000 GMH/c 3 in order to provide a converged characterization of the source, although we use snapshots from the latter 25, 000 GMH/c 3 of each simulation.In particular, we consider two MAD simulations with dimensionless black hole spins a/MH = −0.5 and +0.94.We also consider one SANE simulation with spin a/MH = −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  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 β in Fig. 6.

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.

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. (2018Anantua et al. ( , 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 , where Pe is the (presumed isotropic) partial pressure of the electrons emitting at the frequency ν ′ .The choice to scale nearhorizon jet emissivity in terms of magnetic pressure is motivated by the observation that the jet becomes increasingly simple and electromagnetic -high σ -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 Pe.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 v ⊥ = E × B/B 2 in the simulation frame.In general, the component of the fluid velocity resolved along the field v ∥ 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 γ-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 ∼ 1 − 300 EV, many orders of magnitude greater than the ∼ 1 MV minimum needed to create positrons or the ∼ 1 GV needed to accelerate electrons to the ≲ 1 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, ∼ 1000M and we shall suppose that the associated radio emission is a locally accelerated pair plasma where γ-ray production is balanced by annihilation and other quantum electrodynamical processes.

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.

Constant Electron β Model
The simplest prescription is the Constant βe Model, that the pressure of synchrotron radiating electrons is a constant fraction of the magnetic pressure Close to the hole, the total pressure P is dominated by the magnetic pressure PB = b 2 /2µ0 but at large radius, there is gas pressure contributed by the entrained gas.If we adopt β ≲ 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.

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 ′ -and introducing a resistivity η = µ0Lj where Lj is a length scale which we choose to be a fixed fraction of the jet width.The dissipation rate is then W ′ = ηj ′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.

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 featuresnotable 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 magnetorotational 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.

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

Critical β Model
The Critical β Model is an alternative turbulent heating model with an exponential parameter βc controlling the transition between electron-and ion-dominated heating 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 β, the electron-to-ion temperature ratios in the R − β and the Critical Beta Models have similar asymptotic behavior.However, at high β, the Critical β Model Te/Ti 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 β Model also has transition behavior at intermediate β's controlled by an exponential parameter βc, leading to greater variety of intermediate β emission regions probed between the same range of electron-to-ion temperature ratios compared to the R − β Model.

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 β and σ.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 β-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 σ of 1/2.Armed with these models of particle thermodynamics, we turn to the emitted radiation.

Radiation
An electromagnetically-dominated jet outflow must be continually converting electromagnetic energy into high energy pairs through a −E • 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 νSA =

22
LEM 10 44 erg s −1 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 If a power law is accelerated, the local spectrum should break by ∆α = 0.5 at this frequency.However, note the extreme sensitivity to the bulk Lorentz factor Γ. This probably controls what is actually observed.
The entire nuclear spectrum within r ∼ 3 × 10 5 M has been carefully determined by Prieto et al. (2016).They find a sharp break at ν = ν b ∼ 150 GHz.Presumably the flatter α ∼ 0.2 spectrum at ν ≲ ν 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, ∼ 3 GHz, should originate at r ∼ 1000M .

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 γ − γ 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 γ-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.

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 ≈ 1 MeV ≈ 2 × 1.2 × 10 20 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 stellarmass 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 ≫ 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 ne0 and a jet where all sources of pair production result in a plasma with an equal number density of ionic and pair plasma (fpos ≡ npairs/ne0 = 1).
There are benchmarks for positron content in the Literature, such as the estimate by Ghisellini (2012) of the frac-tion of a jet (opening angle ψ, distance from black hole R0) converted to positrons, as f = 0.1 min{1, ℓ 0 60 } where the compactness is ℓ = σ T L 0 ψR 0 D 4 mec 3 .More possibilities for painting positrons on jets can be found in Appendix.

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 β are relatively azimuthally symmetric.
We now describe the process of ray tracing the resulting emission.

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 xN pixels over a 160 µ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 σ ≡ b 2 /ρ and artificially inject mass and energy in these regions.Since the plasma density (and temperature) are therefore typically unphysically high in regions of large σ, ray-tracing codes like ipole normally introduce a

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 β throughout the inflow/outflow system Event Horizon Telescope Collaboration et al. (2019e, 2021b).We start with models using this approach with the R − β and Critical β turbulent heating prescriptions, and then we refine the models by prescribing jet funnel emission in a region between σmin = 1/2 and σ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.5MMAD and SANE snapshots from Anantua et al. (2023) to include the highly prograde spin a = 0.94, temporal evolution, extreme positron fractions npairs/ne0 = 50, 100, varying frequency from 230 GHz to 86 GHz and varying inclination up to 40 • 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.

SANE R-β
Fig. 9 shows intensity with electric vector polarization angle (EVPA) and circular polarization maps for the SANE a = −0.5 simulation in the R − β model.This model has asymptotic ion-to-electron temperature ratios R low = 1 and R 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µas x 40µ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 munit 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.

SANE Critical β
In Fig. 10 we show our other turbulent heating modelthe Critical Beta Model.This model controls the transition from preferential electron heating to preferential ion heating through the exponential parameter β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 βc = 0.5, total flux is concentrated in a ring at ∼ 20µas and regions along lines of sight close to the polar axis, though polarization morphology trends remain similar to the R − β case.Note our R − β models are more linearly depolarized than Critical β Models even with lower contributions to intrinsic emission at high β in the latter models.This is one of several examples of Faraday effects we will see in this Section.

SANE R-β with Constant β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−β 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 σ separating the constant βe jet from the turbulently heated plasma.

SANE Critical β with Constant β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.

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).

MAD R-β
Starting in Fig. 13 with the R-β 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.

MAD Critical β
In Fig. 14, the Critical β image and V /I map mirror the global structure in the R − β 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.

MAD R-β with Constant β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 βe jets slightly reduced the circular polarization degree both with and without positrons.

MAD Critical β with Constant βE Jet
In Fig. 16 the Critical Beta model with jet exhibits the same trends as its R − β counterpart above.

Comparison of R-Beta and Critical Beta Models
In Fig. 17 we compare Critical Beta Model with R−β 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 β's may contribute given our shallow exponential parameter βc = 0.5.

Extreme Positron Fractions
As mentioned in Sec.6 increasing the pair fraction causes intrinsically emitted circular polarization to decrease, Fara- day 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 npairs/n0 = 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-β 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.

Comparison of Models with Polarization Constraints
In  U and I across the image plane, and its average local (resolved) magnitude 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/resolutiondependent 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|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 Vnet constraint in Goddi et al. (2021).For the structure parameter β2, however, Table 5 shows only MAD a = 0.94M pure turbulent heating models pass.
We may check the robustness of these tendencies by reanalyzing over the temporal evolution of the simulations (cf.Sect.6.8).Though little changes for the comparison against the more reliable unresolved linear polarization (|mnet|) observations, the beam-dependent (< |m| >) comparison changes significantly with temporal evolution.

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.

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 3. Linear polarization |m|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 ⩽ |m|net ⩽ 0.037 and 0.057 << |m| >< 0.107.Note that the bold values refer to fiducial models which satisfy the net linear polarization constraints.

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).

Frequency and Inclination Dependence
In Figs.21 and 22, we search for extended structure at 86 GHz in the R-β with β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 distinc- Table 6.Faraday rotation depth at min and max fpos for fiducial models at T = 25, 000M : .
tive image morphological features, as shown in the 86 GHz maps in Fig. 23 with inclination angle 40 • (instead of the default M87 inclination 17 • used throughout this work).In Fig. 23, the R−β model with β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.

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 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 (ne) 0 (ne) 0 +2n pairs = 1/3 where the plasma is an even mix of electrons, positrons and protons).
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 ṁ and horizon-threading flux ϕ from 10, 000M < T < 30, 000M reveals at T = 25, 000M a sharp rise in ϕ accompanied by a sharp decrease in ṁ, 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 ϕ and troughs in ṁ as at T = 25, 000M .
The MAD advantage found at the fiducial time T = 25, 000M for our emission models with respect to our polarization constraints is fairly consistent in time, with the cavaet that we continue to rely on the net linear polarization |m|net measurement less subject to uncertainties related to the EHT beam resolution than < |m| >).For evenly spaced samples in the range 20, 000M ⩽ T ⩽ 30, 000M ,  MADs pass a higher percentage of |m|net constraints for 82% of times while the corresponding spin SANEs outperform MADs for 18% of times (note for the resolved case, however, SANEs pass more frequently than MADs 73% of times).The SANE/MAD dichotomy in |m|net favoring MAD models while < |m| > favors SANEs manifests in Table 8 comparing the median values in a sample of evenly spaced times from 20, 000M ⩽ T ⩽ 30, 000M .New circular polarization observations by the EHT (?) more are more one-sided (favoring MADs) in our model parameter space.Performing a similar analysis for circular polarization, we find 55% of times MAD has a larger fraction of models passing the |v|net constraint, 55% of times MAD has  laboration et al. 2021b) is the jet power, which heavily excludes SANEs, is observationally uncertain (ranging between 10 42 to and 10 45 erg/s) and which we do not consider here with our limited number of simulations.The overwhelming MAD advantage found by EHT depends on both the choices for simulations/emission models to be included in the library and the particular constraints chosen to make the comparison, underscoring the need for alternative investigations such as this exploring new regions of parameter space.

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, frequen- cies 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/MH = −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 (∼ 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   models, the jet magnetic substructure for Constant Electron Beta models characterized by constant electron-gas-tomagnetic 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.β and Γv).We have implemented models where the emissivity is governed by total plasma β in the turbulent plasma and by conversion of magnetic-to-particle energy (parametrized by the contribution β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 − β and Critical β 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 βe jet tends to broader annular emission profiles.
• MAD and SANE images with polarization at constant overall flux have markedly different morphological proper- ties.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.

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 Be = D|⃗ n × ⃗ B| exhibit violation of bilateral symmetry across the jet axis both in the stationary, axisymmetric, self-similar semi-analytic model (?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 robustalbeit, 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 Be in future emission mod- 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 limbbrightening past 100M .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-positronratio maps used here.The computational expense of a threefluid (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 En-    tion at the Center for Astrophysics as well as grant numbers 21-atp21-0077, NSF AST-1816420 and HST-GO-16173.001-A.

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

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 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 ue± ∝ 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 ue± = γe±ne±mec 2 to construct inner jet emissivity ν −α , z ⩾ zmax, 0 < ξ < 1 2 ξmax where zmax ≈ 10rg, 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.

)Figure 5 .
Figure 5. Vertical slices for electron number density Ne, 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).
7) 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).
Figure 7.Comparison of R-β (solid lines) and Critical Beta (dashed lines) models for reasonable parameter values.

Figure 9 .
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.

Figure 10 .
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.

Figure 11 .
Figure 11.For the a = −0.5 SANE at T = 25, 000M : Top panel: R-Beta with β 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.

Figure 12 .
Figure 12.For the a = −0.5 SANE at T = 25, 000M : Top panel: Critical Beta with β 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.

Figure 13 .
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.

Figure 14 .
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.

Figure 15 .
Figure 15.For the a = −0.5 MAD at T = 25, 000M : Top panel: R-Beta with β 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.

Figure 16 .
Figure 16.For the a = −0.5 MAD at T = 25, 000M : Top panel: Critical Beta with β e0 = 0.01 jet at 230 GHz without positrons Bottom panel: Critical Beta with β e0 = 0.01 jet at 230 GHz for an even mix of pair and ionic plasma.
morphology and dynamics.Starting with turbulent heating models including that used by the EHT, R-β, and the Critical β 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 Figure 20.MAD a = +0.94R-Beta at 230 GHz.Top panel: fPos = 0 Bottom panel: fPos = 100
hancement towards the end.This work was supported by a grant from the Simons Foundation (00001470, RA, LO, JD and BC).Roman Shcherbakov and Alexander Tchekhovskoy have provided excellent guidance and mentorship at the beginning of this investigation.UTSA undergraduate Noah Heridia has been helpful through graphic-related activities.BASIS Shavano San Antonio student Luke Fehlis provided valuable input on data analysis.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 Computa-

Figure 26 .
Figure 26.Degree of circular polarization as a function of n pairs (ne) 0 (Top Panel) and unpaired synchrotron emitter frac-
Proxies by Plasma VariablesOne 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 e − + P e + = βe0(1 + (f e + (β) OR R e + (β)))PB.(20)

Table 1 .
M87 Geometry Distance from Earth 1 Schwarzschild radius 2 Apparent angular width 3 Jet opening angle 4 Jet viewing angle 5

Table 3
, we compare observations from Event Horizon Telescope Collaboration et al. (2021a) against fiducial model linear polarization: both summed from net (unresolved) Q,

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

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

Table 8 .
Median linear polarization tables for simulation times 20,000M ⩽ T ⩽ 30, 000M .Note that the bold values refer to SANE and MAD fiducial models with a/M = −0.5 which satisfy the linear polarization |m|net and < |m| > constraints from EHT M87 Paper VII, which take the form of the polarization ranges 0.01 ⩽ |m|net ⩽ 0.037 and 0.057 << |m| >< 0.107.We indicate models satisfying polarization constraints only at the number of constraint significant figures in italics.

Table 9 .
Median circular polarization tables for simulation times 20,000M ⩽ T ⩽ 30, 000M .Note that the italic values refer to SANE and MAD fiducial models with a/M = −0.5 that fail to satisfy the constraints |v|net ⩽ 0.008 and < |v| >⩽ 0.037