First M87 Event Horizon Telescope Results. VIII. Magnetic Field Structure near The Event Horizon

Event Horizon Telescope (EHT) observations at 230 GHz have now imaged polarized emission around the supermassive black hole in M87 on event-horizon scales. This polarized synchrotron radiation probes the structure of magnetic fields and the plasma properties near the black hole. Here we compare the resolved polarization structure observed by the EHT, along with simultaneous unresolved observations with the Atacama Large Millimeter/submillimeter Array, to expectations from theoretical models. The low fractional linear polarization in the resolved image suggests that the polarization is scrambled on scales smaller than the EHT beam, which we attribute to Faraday rotation internal to the emission region. We estimate the average density n_e of order 10^4-7 cm-3, magnetic field strength B of order 1-30 G, and electron temperature Te of order (1-12) x 10^10 K of the radiating plasma in a simple one-zone emission model. We show that the net azimuthal linear polarization pattern may result from organized, poloidal magnetic fields in the emission region. In a quantitative comparison with a large library of simulated polarimetric images from general relativistic magnetohydrodynamic (GRMHD) simulations, we identify a subset of physical models that can explain critical features of the polarimetric EHT observations while producing a relativistic jet of sufficient power. The consistent GRMHD models are all of magnetically arrested accretion disks, where near-horizon magnetic fields are dynamically important. We use the models to infer a mass accretion rate onto the black hole in M87 of (3-20) x 10^-4 Msun yr-1.


Introduction
The Event Horizon Telescope (EHT) Collaboration has recently published total intensity images of event-horizon-scale emission around the supermassive black hole in the core of the M87 galaxy (M87 * ; Event Horizon Telescope Collaboration et al. 2019aCollaboration et al. , 2019bCollaboration et al. , 2019cCollaboration et al. , 2019d, hereafter EHTC I, EHTC II, EHTC III, EHTC IV). The data reveal a 42 ± 3 μas diameter ringlike structure that is broadly consistent with the shadow of a black hole as predicted by Einstein's Theory of General Relativity (Event Horizon Telescope Collaboration et al. 2019e, 2019f; hereafter EHTC V, EHTC VI). The brightness temperature of the ring at 230 GHz (10 10 K) is naturally explained by synchrotron emission from relativistic electrons gyrating around magnetic field lines. The ring brightness asymmetry results from light bending and Doppler beaming due to relativistic rotation of the matter around the black hole. M87 * is best known for launching a kpc-scale FR-I type relativistic jet, whose kinetic power is estimated to be ∼10 42-44 erg s −1 (e.g., Stawarz et al. 2006;de Gasperin et al. 2012). The structure of the relativistic jet has been resolved and studied at radio to X-ray wavelengths (e.g., Di Matteo et al. 2003;Harris et al. 2009;Kim et al. 2018;Walker et al. 2018).
The published EHT image of M87 * together with multiwavelength observations are consistent with the picture that the supermassive black hole in M87 is surrounded by a relativistically hot, magnetized plasma (Rees et al. 1982;Yuan & Narayan 2014;Reynolds et al. 1996;Yuan et al. 2002;Di Matteo et al. 2003). However, it is not clear whether the compact ring emission is produced by plasma that is inflowing (in a thick accretion flow), outflowing (at the jet base or in a wind), or both. Furthermore, the total intensity EHT observations also could not constrain the structure of magnetic fields in the observed emission region. In order to find out which physical scenario is realized in M87 * , additional information is necessary.
Event Horizon Telescope Collaboration et al. (2021, hereafter EHTC VII) reports new results from the polarimetric EHT 2017 observations of M87 * . The polarimetric images of M87 * are reproduced in Figure 1. These images reveal that a significant fraction of the ring emission is linearly polarized, as expected for synchrotron radiation. The EHT polarimetric measurements are consistent with unresolved observations of the radio core at the same frequency with the Submillimeter Array (SMA; Kuo et al. 2014) and the Atacama Large Millimeter/submillimeter Array (ALMA; Goddi et al. 2021). They also provide a detailed view of the polarized emission region on event-horizon scales near the black hole. Polarized synchrotron radiation traces the underlying magnetic field configuration and magnetized plasma properties along the line of sight (Bromley et al. 2001;Broderick & Loeb 2009;Mościbrodzka et al. 2017). These polarimetric measurements allow us to carry out new quantitative tests of horizon-scale scenarios for accretion and jet launching around the M87 * black hole. In this Letter we present our interpretation of the EHTC VII resolved polarimetric images of the ring in M87 * .
Our analysis is presented as follows. In Section 2 we report polarimetric constraints from M87 * EHT 2017 and supplemental observations, and argue that they can be used for scientific interpretation, focusing on several key diagnostics of the degree of order and spatial pattern of the polarization map. In Section 3 we present one-zone estimates of the properties of the synchrotron-emitting plasma. In the transrelativistic parameter regime relevant for the M87 core, a full calculation of polarized radiative transfer using a realistic description of the plasma properties is essential. To that end, in Section 4 we describe a set of numerical simulations of magnetized accretion flows that we use to compare with our set of observables. In Section 5 we show that only a small subset of the simulation parameter space is consistent with the observables. The favored simulations feature dynamically important magnetic fields. We discuss limitations of our models in Section 6, and discuss how future EHT observations can further constrain the magnetic field structure and plasma properties near the supermassive black hole's event horizon in Section 7.

Conventions in Observations and Models
Throughout this Letter we use the following definitions and conventions for polarimetric quantities, following the International Astronomical Union (IAU) definitions of the Stokes parameters ( )     , , , (Hamaker & Bregman 1996;Smirnov 2011). The complex linear polarization field  is Then, the electric-vector position angle (EVPA) is defined as The EVPA is measured east of north on the sky. Therefore, positive  is aligned with the north-south direction and negative  with the east-west direction. Positive  is at a + 45 deg angle with respect to the positive  axis (in the northeast-southwest direction). Positive Stokes  indicates right-handed circular polarization, meaning in our convention that the electric field vector of the incoming electromagnetic wave is rotating counter-clockwise as seen by the observer. In the synchrotron radiation models that we consider, a positive value of emitted Stokes  is associated with an angle θ B between the wavevector k μ and magnetic field b μ as measured in the frame of the emitting plasma in the range θ B ä [0, 0.5π].
Negative  corresponds to θ B ä [0.5π, π]. The linear and circular polarization fractions at a point in the image are defined as We also define the rotation measure between two wavelengths λ 1 and λ 2 where the sums are over all pixels i in the resolved image. In addition to the signed circular polarization fraction v net , we also frequently consider the absolute value |v net |, as circular polarization measurements of the M87 * core at 230 GHz do not constrain its sign (Goddi et al. 2021).
In describing the resolved linear polarization in EHT images, we define the image-average linear polarization fraction, weighted by the total intensity of each image pixel, as Note that 〈|m|〉 depends on the imaging resolution (beam size), while |m| net is the usual unresolved linear polarization fraction and does not depend on resolution.

Unresolved Polarization and Rotation Measure
Measurements toward M87's Core from ALMA As part of the EHT 2017 observation campaign, we obtained ALMA array measurements of the unresolved, net, near 230 GHz, polarimetric properties of M87's core and jet on 2017 April 5, 6, 10, and 11 (hereafter these observations are referred to as ALMA-only observations). ALMA-only measurements are given at ν = 221 GHz, a central frequency of ALMA Band 6, which has four spectral windows, each centered at 213, 215, 227, and 229 GHz. These results, along with details on the observations and data calibration, are presented in Goddi et al. (2021); we summarize them here in Table 1. From the ALMA-only data, the net linear polarization fraction (Equation (6)) of the core is |m| net ; 2.7%. The data also provide an upper limit on the net circular polarization fraction (Equation (7)) of the core of |v| net  0.3%, with a magnitude and sign that vary over the four observed epochs. Goddi et al. (2021) also measured an EVPA that varies with wavelength across the ALMA band; the slope of EVPA with wavelength is consistent with EVPA ∝ λ 2 , as expected for Faraday rotation. The inferred rotation measure (Equation (5), for frequencies/wavelengths in ALMA Band 6) is also time variable and changes sign between 2017 April 5 and 11, with a maximum magnitude |RM| ; 1.5 × 10 5 rad m −2 .
The ALMA-only measurements include extended ∼arcsecond scale structures that are entirely resolved out of the EHT maps of M87's core region. As a result, the total 221 GHz flux density of M87 * measured by ALMA alone is a factor of ;2 larger than that captured by the resolved EHT images (see also EHTC IV). For that reason, we adopt a more conservative upper limit of |v| net < 0.8%, which would account for the case where the large-scale emission is not circularly polarized.

Spatially Resolved Linear Polarization of M87's Core in EHT 2017 Data
The resolved polarimetric images of the M87 core reported in EHTC VII display robust features between different imagereconstruction algorithms and across four days of observations (2017 April 5, 6, 10, and 11). At 20 μas resolution, those images consistently show a region of the highest linear polarized intensity in the southwest portion of the ring, with a fractional linear polarization |m|  30% at its maximum. The image-average linear polarization fraction takes on values 5.7% 〈|m|〉 10.7% across the different observation days and image-reconstruction techniques. The range of the imageintegrated net polarization fraction is 1.0% |m| net 3.7% (see EHTC VII, their Table 2). Because polarized emission outside the EHT field of view but inside the ALMA-only core is unconstrained, we adopt the EHT |m| net range when evaluating models.
The EHT images on all four observing days reveal a characteristic azimuthal pattern of the EVPA angle around the emission ring. To quantify this pattern across the image, we use a decomposition of the complex linear polarization = +    i into azimuthal modes with complex coefficients β m (Palumbo et al. 2020, hereafter PWP). For a polarization field in the image plane given in polar coordinates (ρ, j), the β m coefficients are where I ann is the Stokes  flux density contained inside the annulus set by the limiting radii r min and r max . We take r = 0 min and r max to be large enough to include the entire EHT image.
Within the library of polarized images from general relativistic magnetohydrodynamic (GRMHD) simulations presented in EHTC V, PWP found that the m = 2 coefficient, β 2 , was the most discriminating in identifying the underlying magnetized accretion model. The phase of β 2 maps well onto the qualitative behavior expected of polarization maps with idealized magnetic field configurations. In our convention, radial EVPA patterns have positive real β 2 (∠β 2 = 0 deg), azimuthal EVPA patterns have negative real β 2 (∠β 2 = 180 deg), and left-(right-) handed spiral patterns have positive (negative) pure imaginary β 2 (∠β 2 = 90 deg and −90 deg, respectively). These idealized EVPA pattern configurations and their corresponding β 2 coefficients are summarized in Appendix A and Figure 1 of PWP. The measured range of the complex β 2 coefficient across the different imagereconstruction methods and observing days reported in EHTC VII, their for the phase. Appendix A demonstrates that the information in the β 2 coefficient can be obtained in the visibility domain using measurements of the linear polarization E − (gradient) and B − (curl) modes of the polarization field (e.g., Kamionkowski & Kovetz 2016). Trends in β 2 metric computed across the GRMHD image library (Section 4) can be obtained in the visibility domain using only E-and B-mode measurements taken on EHT 2017 baselines, as long as the visibilities are accurately phase calibrated. Because accurate phase calibration of EHT data is non-trivial and requires fully modeling the polarized source structure, we use image-domain comparisons to the reconstructions presented in EHTC VII for the constraints in this Letter. As in total intensity, both the unresolved and resolved polarimetric properties of the 230 GHz M87 * image changed over the week between 2017 April 5 and April 11. Notably, the integrated EVPA in the EHT image changes by ≈30-40 deg (while the ALMA-only EVPA changes by 10 deg). We will not interpret this variability further in this work; however, Section 6 discusses expectations for time variability from viable simulation models. The observational ranges of the key parameters that we use in comparing theoretical models to data in Section 5-namely |m| net , |v| net , 〈|m|〉, and β 2 amplitude and phase-are summarized in Table 2.

External and Internal Faraday Rotation
Faraday rotation in a uniform plasma with rotation measure (RM) rotates the EVPA away from its intrinsic value EVPA 0 according to Equation ( Polarized light rays passing through a uniform medium are subject to the same RM. The net source polarization angle is then coherently rotated away from its intrinsic value without any depolarization. This scenario of "external" Faraday rotation has been used to infer the mass accretion rate for sources where an RM is measured or constrained (e.g., Bower et al. 2003;Marrone et al. 2006Marrone et al. , 2007Kuo et al. 2014), by assuming that the observed radiation passes through the bulk of the accretion flow. Because relativistic electrons suppress the Faraday rotation coefficient as ∝ T 1 e 2 (e.g., Jones & Hardee 1979), these models assume that the RM is produced outside the emission region at the radius where Θ e = kT e /m e c 2 = 1, usually r ∼ 100 r g (where r g = GM/c 2 is the gravitational radius).
However, in accreting systems like M87 * , it is unclear whether this external Faraday rotation model is a good approximation. As we estimate below, one-zone emission models of M87 * predict substantial RM within the emission region itself at radii r  5 r g . At its low viewing inclination, the observed polarized radiation emitted near the horizon may not pass through the bulk of the high-density, infalling gas. Therefore, "internal" Faraday rotation, which can depolarize the emission as well as rotate the EVPA (Burn 1966), is also an important effect to consider.
The observed ;10% linear polarization of the ring at the EHT scale of ∼20 μas is much lower than the intrinsic synchrotron polarization fraction 70% expected locally. This could result from synchrotron self-absorption of the emitted radiation, but one-zone estimates and theoretical models (e.g., EHTC V, and references therein) suggest that the 230 GHz emission is mostly optically thin. It is more likely that the observed depolarization of the resolved emission could be the result of polarization structure that is scrambled at resolutions finer than the EHT beam. Turbulent magnetic fields and Faraday rotation internal to the emission region could produce this scrambling. In Section 4.3 we show that turbulence in GRMHD models alone is insufficient to produce this level of depolarization. Significant internal Faraday rotation of polarization vectors on different rays by different amounts can produce a sufficiently scrambled image that is depolarized when spatially averaged over a telescope resolution element (beam,e.g., Burn 1966;Agol 2000;Quataert & Gruzinov 2000;Beckert & Falcke 2002;Ruszkowski & Begelman 2002;Ballantyne et al. 2007).
From the simultaneous ALMA-only M87 * observations, the RM implied by changes in the EVPA across the ALMA band is |RM|  1.5 × 10 5 rad m −2 . These values are consistent with, but much more constraining than, the range determined from past SMA observations (−3.4-7 × 10 5 rad m −2 , Kuo et al. 2014). The ALMA-only EVPA difference varies by order unity in magnitude and sign over the observing campaign, and includes a large flux contribution from extended emission not captured by EHT 2017 imaging (EHTC IV). Using a twocomponent model, Goddi et al. (2021) show that the RM toward the core emission in the EHT field of view could exceed the RM computed from the ALMA-only data, with allowed values as large as |RM|  10 6 rad m −2 . Because of this uncertainty, we do not use the observed RM as an observational constraint in our analysis. We account for uncertainty related to the observed time variability by using reconstructed polarized EHT images from both 2017 April 5 and 11 to define the acceptable ranges (see Table 2) of the observational parameters used to score theoretical models in Section 5.

Estimates and Phenomenological Models
In this Section, we take a first look at the importance of internal Faraday rotation and magnetic field structure in determining the characteristics of the 230 GHz EHT image. In Section 3.1 we obtain order-of-magnitude estimates of the plasma properties in M87 * by interpreting the observed depolarization as entirely due to the effect of internal Faraday rotation on small scales. In Section 3.2 we explore the effects of different idealized magnetic field configurations on the observed polarization pattern from plasma orbiting a black hole in the absence of Faraday effects.

Parameter Estimates from One-zone Models
Based on a one-zone isothermal sphere model, EHTC V derived order-of-magnitude estimates of the plasma number density n e and magnetic field strength B in the emitting region around M87 * as constrained by the Stokes  image's brightness, size, and total flux density: Note. The ranges for |m| net , 〈|m|〉, and β 2 were taken from EHTC VII Table 2. These ranges represent the minimum −1σ error bound and maximum + 1σ error bound across five different image-reconstruction methods, and incorporate both statistical uncertainty in the polarimetric image reconstruction and systematic uncertainty in the assumptions made by different reconstruction algorithms. The upper limit on |v| net was taken as ;2× the maximum value found by Goddi et al. (2021).
In this model, the emission radius was assumed to be r ; 5r g , and the electron temperature was assumed to be T e = 6.25 × 10 10 K, based on the observed brightness temperature of the EHT image. This temperature corresponds to Θ e = kT e /m e c 2 = 10.5, so the emitting electrons have moderately relativistic mean Lorentz factorsḡ » Q » 3 30 e . The angle between the magnetic field and line of sight is set at θ = π/3. This model ignores several physical effects that are included in more sophisticated models and simulations and which are necessary to fully describe the emission from M87 * . The plasma is considered to be at rest and so there is no Doppler (de)boosting of the emitted intensity from relativistic flow velocities. Redshift from the gravitational potential of the black hole is also not included.
Given n e , B, and T e , we can estimate the strength of the Faraday rotation effect at 230 GHz quantified by the optical depth to Faraday rotation t r V : where ρ V is the Faraday rotation coefficient (e.g., Jones & Hardee 1979). For emission entirely behind an external Faraday screen, t r V is related to the rotation measure RM via t l = r 2RM 2 V , which follows from the radiative transfer equations for spherical Stokes parameters in the absence of other effects (see e.g., Appendix A of Mościbrodzka et al. 2017) and the fact that ρ V ∝ λ 2 .
Our estimated t r V indicates that Faraday rotation internal to the emission region is an important effect and could thus explain the depolarization observed in M87 * . Faraday effects are even more important for the case of polarized light emitted by relativistic electrons that travel through a dense, colder accretion flow (e.g., Mościbrodzka et al. 2017;Ricarte et al. 2020). In addition, for the same parameters, Faraday conversion of linear to circular polarization may also be important (  t r 0.5 Q ), while self-absorption is weak (τ I ; 0.05). Requiring an internal Faraday optical depth t p > r 2 V (large enough to produce significant depolarization) provides an additional constraint on one-zone models independent of those used in EHTC V, which fixed the electron temperature at an assumed value. Assuming t p > r 2 V allows us to break the degeneracy between magnetic field strength, electron temperature, and plasma number density.
Hence, we consider the same model as in EHTC V at several different values of β e = 8πn e kT e /B 2 , constrained by the requirement that the Faraday optical depth t p > r 2 V . To be consistent with the 230 GHz EHT data, we also require that the observed image have a total flux F ν between 0.2 and 1.2 Jy, and that the model has a maximum total intensity optical depth τ I < 1. Figure 2 shows what constraints these requirements put on the electron number density n e and the dimensionless electron temperature Θ e at three different values of β e . For values of 0.01 < β e < 100, in this simple model the electron temperature is constrained to lie in a mildly relativistic regime 2  Θ e  20 (10 10 < T e < 1.2 × 10 11 K), and the magnetic field strength is 1  B  30 G. The number density of the emitting electrons depends more sensitively on the assumed value of β e , taking on values between 10 4 cm −3 and 10 7 cm −3 .
The one-zone model estimates suggest that both the total intensity and polarized emission can be produced in a mildly relativistic plasma in a magnetic field of relatively low strength B  30 G. For higher values of B, the electron temperature would be too small to explain the observed maximum brightness temperature (;10 10 K) in the M87 * EHT image (EHTC IV). Very high values of B are independently disfavored by the small degree of circular polarization |v| net  1% seen in M87 * . For B ; 100 G, the ratio of the Stokes  emissivity to the Stokes  emissivity j V /j I ; 1%. For B ; 10 3 G, j V /j I ; 10%, for all temperatures >10 10 K. We also note that magnetic fields of B  5 G are sufficient to produce jet powers of P jet  10 42 erg s −1 (e.g., EHTC V) via the Blandford & Znajek (1977) process.

EVPA Pattern and Field Geometry
To demonstrate how the intrinsic magnetic field structure in the emission region influences the observed polarization pattern, in this section we present the polarization configurations from three idealized magnetic field geometries around a black hole-a purely toroidal field, a purely radial field, and a purely vertical fieldas seen by a distant observer. In Figure 3 Figure 2. Allowed parameter space in number density and dimensionless electron temperature (n e , Θ e ) (red region) for the one-zone model described in Section 3.1 for three constant values of β e = 8πn e m e c 2 Θ e /B 2 . We require that the optical depth τ I < 1 (green region), the Faraday optical depth t p > r 2 V (blue region), and the total flux density 0.2 < F ν < 1.2 Jy (black region). Contours of constant magnetic field strength are denoted by labeled dashed lines.
we show polarimetric images from these simple field configurations computed with two methods: a numerical model of an optically thin emission region around the black hole (top row of Figure 3), and an analytic treatment of the parallel transport of the polarization vector that is originally perpendicular to the magnetic field (R. Narayan et al. 2021, in preparation, middle row of Figure 3). We show the polarization maps from both methods for the three idealized magnetic field configurations viewed at an inclination angle of i = 163 deg. Both the analytical and numerical calculations assume a zero-spin black hole (Schwarzschild metric), though we have found that the main features of these polarization patterns are insensitive to spin. configuration generated by an orbiting emission region in the shape of a torus at 8r g in three imposed magnetic field geometries and viewed at i = 163 deg (with material orbiting clockwise on the sky). The orbital angular momentum vector is pointing away from the observer and to the east (to the left). Total intensity is shown in the background with higher brightness temperature regions shown as being lighter in color. In the foreground, the observed EVPA direction is shown with white ticks, with the tick length proportional to the polarized flux. (b) Analytic calculations of the polarization configuration from a thin ring of magnetized fluid at 8r g inclined by 163 deg to the observer in the same magnetic field geometries as in (a). While the distribution of emitting material is different in the two models, both the sense of asymmetry in the brightness distributions and the polarization patterns match those from the numerical calculations. (c) Schematic cartoons showing the emitting frame wavevector k, magnetic field direction B, and polarization vector =Ṕ k B^^for each case. In the bottom-right panel,ˆ¢ k denotes the approximate light bending contribution to the wavevector.
In the top row of Figure 3 we show the result of numerical calculations performed with the general relativistic ray tracing code grtrans (Dexter & Agol 2009;Dexter 2016) of polarized emission from an optically and Faraday-thin compact emission region, or "hotspot", in Keplerian orbit around a black hole in the equatorial plane. The hotspot has a radial extent of 3r g and moves in an imposed and idealized magnetic field geometry in a circular orbit at a radius of 8r g (following Gravity Collaboration et al. 2018Collaboration et al. , 2020. We construct a phenomenological model of a torus of emitting, rotating plasma by studying the time-averaged polarized emission images from one revolution of this hotspot around the black hole. We have verified that a semi-analytic implementation (Broderick & Loeb 2006) of a hot accretion flow model (Yuan et al. 2003) produces consistent polarization patterns when using the same field geometry.
In the second row of Figure 3, we compare these numerical results to results from an analytic calculation of the observed polarization pattern generated by the emission of polarized light on a thin ring of radius 8r g in the equatorial plane. In this model (R. Narayan et al. 2021, in preparation) the polarization vectors are emitted perpendicular to the imposed magnetic field geometry in the fluid rest frame; they are transformed on their way to the observer using an approximate, analytic treatment of the effects of light bending, parallel transport, and Doppler beaming. This calculation includes radial inflow as well as rotation in the velocity field; the models shown use purely toroidal motion (clockwise on the sky) with the same idealized magnetic field geometries as in the numerical case. The models match the asymmetric brightness distributions and polarization patterns of the numerical calculations. In particular, both models produce consistent helical EVPA pattern in the case of a vertical magnetic field.
The linear polarization directionP of synchrotron radiation in the emitted frame is perpendicular to the wavevectork and the magnetic field vectorB. We define the toroidal magnetic field as consisting only of magnetic field components in the azimuthal direction, while the poloidal magnetic field consists of the remainder, including both radial and vertical components. In a purely toroidal field case, the EVPA shows a radial pattern (left column in Figure 3). Purely radial magnetic fields (middle column) give a complementary result; the polarization has a toroidal configuration, similar to a 90 deg rotation of the linear polarization ticks from the toroidal case.
In a vertical magnetic field (right column in Figure 3), we might expect thatP should be vertical (north−south) everywhere since a verticalB is tilted east−west for this viewing geometry. We might also expect that¯ P 0 when the black hole is viewed face on, becauseˆ||k B. Instead, the linearly polarized emission from a purely vertical field shows a twisting pattern that wraps around the black hole. The twist results from a combination of light bending and relativistic aberration. Light bending in the emitting region near the black hole contributes a radial componentˆ¢ k to the emitted wavevectork that initially points away from the black hole (see the schematic cartoon in the bottom-right panel of Figure 3). As a result, close to the black hole, the total wavevectorˆˆ= + ¢ k k k emit and the magnetic fieldB are no longer parallel, the polarization is non-zero, and the resulting EVPA pattern is north−south symmetric. Relativistic motion of the emitting material (aberration) breaks the symmetry and gives the twisting pattern a handedness corresponding to the orbital direction. For the pure vertical field considered here, the handedness depends on the rotation direction and the observed pattern is consistent with clockwise rotation. The dependence on direction of motion and magnetic field configuration are discussed in more detail in a forthcoming paper (R. Narayan et al. 2021, in preparation). The EVPA patterns in these images do not show a strong dependence on the black hole spin.
In a rotating flow, weak magnetic fields are sheared into a predominantly toroidal configuration (e.g., Hirose et al. 2004). In the absence of other effects (e.g., external Faraday rotation), the observed azimuthal EVPA pattern suggests the presence of dynamically important magnetic fields in the emission region, which can retain a significant poloidal component in the presence of rotation. In the next sections, we compare numerical simulations of the accretion flow and jet-launching region in M87 * with different field configurations to the EHT2017 data to better constrain the magnetic field structure.

M87 * Model Images from GRMHD Simulations
The low resolved fractional linear polarization observed by the EHT contradicts the results from an idealized magnetic field structure with no disorder. For typical parameters of the 230 GHz emission region, Faraday rotation and conversion are expected to be important. Magnetic field structure, plasma dynamics and turbulence, and radiative transfer effects including Faraday rotation can be realized in images from three-dimensional general relativistic magnetohydrodynamic (3D GRMHD) simulations of magnetized accretion flows. We use 3D GRMHD simulations (described in Section 4.1) in combination with polarized general relativistic radiative transfer (GRRT) models (described in Section 4.2) to model polarized images of M87 * . In Section 4.3, we describe trends of the key observables (|m| net , |v| net , 〈|m|〉, and β 2 ) in our GRMHD polarimetric image library.

GRMHD Model Description
The simulation library generated for the analysis of the EHT 2017 total intensity data in EHTC V consists of a set of 3D GRMHD simulations that were postprocessed to generate simulated black hole images via GRRT. For simulations using black holes with non-zero angular momentum, we only considered accretion flows in which the angular momentum of the flow and the hole were aligned (parallel or anti-parallel). Because the equations of non-radiating 129 GRMHD are scale invariant, each fluid simulation was thus fully parameterized by two values describing the angular momentum of the black hole and the relative importance of the magnetic flux near the horizon of the accretion system. A comparison of several contemporary GRMHD codes, including those used to generate the simulation library, can be found in Porth et al. (2019) and in H. Olivares et al. (2021, in preparation).
The black hole angular momentum J is expressed in terms of the dimensionless black hole spin parameter a * ≡ Jc/GM 2 . In this Letter, we consider simulations run with the iharm code (Gammie et al. 2003;Noble et al. 2006) with a * = −0.94, −0.5, 0, 0.5, and 0.94, where positive (negative) spin implies alignment (anti-alignment) between the accretion disk and the black hole angular momentum. Several studies of "tilted" disks have been conducted (e.g., Fragile et al. 2007;McKinney et al. 2013;Morales Teixeira et al. 2014;Liska et al. 2018;White et al. 2019;Chatterjee et al. 2020). As there does not yet exist a full library of tilted disk simulations spanning a range of spins, we limit our analysis to the aligned and anti-aligned simulations considered in EHTC V.
The strength of the magnetic flux near the horizon qualitatively divides accretion flow solutions into two categories: the Magnetically Arrested Disk (MAD) state (e.g., Bisnovatyi-Kogan & Ruzmaikin 1974;Igumenshchev et al. 2003;Narayan et al. 2003) in which the magnetic flux near the horizon saturates and significantly affects the dynamics of the flow, and the contrasting Standard and Normal Evolution (SANE) state (e.g., De Villiers et al. 2003;Gammie et al. 2003;Narayan et al. 2012). The relative importance of magnetic flux in a simulation is quantitatively described by the dimensionless quantity where Φ BH is the magnitude of the magnetic flux crossing one hemisphere of the event horizon (see Tchekhovskoy et al. 2011;Porth et al. 2019) and  M is the mass accretion rate through the event horizon. The flux saturates at values of f  50, and the flow becomes MAD. The SANE simulations that we consider have lower values of f ≈ 5. 130 Accreted material supplied at large scales could, in principle, supply any value of net vertical flux. Here, we do not explore cases with small or zero net vertical flux f  1. We also do not consider values in the relatively narrow intermediate range 5  f  50.
The SANE simulations considered here used a grid resolution of 288 × 128 × 128, a fluid adiabatic index γ = 4/ 3, and an outer simulation domain of r out = 50 r g . The MAD simulations used a grid resolution of 384 × 192 × 192, an adiabatic index γ = 13/9, and an outer simulation domain of r out = 10 3 r g . The simulations were carried out in modified spherical polar Kerr-Schild coordinates, where grid resolution is concentrated toward the midplane to help resolve the magnetorotational instability. All models in the EHT library are evolved for at least t = 10 4 r g /c and their accretion flows reach steady state within r = 10-20 r g .

Ray-traced Polarimetric Images from GRMHD Simulations
Unlike the equations of GRMHD, the equations of radiative transfer are not scale invariant, and so we must introduce two scales to the simulation when we ray-trace images from the numerical fluid data. The length (and time) scale is set by the mass of the black hole, assumed to be M BH = 6.2 × 10 9 M e in accordance with the value used to generate the EHTC V simulation library. For our models, we also adopt the D = 16.9 Mpc distance to M87 * used in EHTC V. The density scale of the accreting plasma (equal to the scale of the magnetic pressure) is chosen so that on average the simulated images reproduce the observed 230 GHz compact flux density, F ν ; 0.5 Jy.
Images were generated from the set of simulations for several values of the polar inclination angle i chosen to be broadly consistent with observational estimates of the inclination angle of the M87 jet (e.g., Walker et al. 2018). The position angle on the sky can be changed after image generation by rotating both the image and the Stokes  and  components appropriately. Each image has a 320 × 320 pixel resolution over a 160 μas field of view, where each pixel contains full Stokes     , , , intensities. In GRMHD simulations, we make the approximation that the plasma is thermal, i.e., that the electrons and ions are described by a Maxwell-Jüttner distribution function (Jüttner 1911). However, the plasma around M87 * and in other hot accretion flows is most likely collisionless, with electrons and protons that are unable to equilibrate their temperatures (e.g., Shapiro et al. 1976;Ichimaru 1977). We mimic collisionless plasma properties in producing images from the GRMHD simulations by allowing the electron temperature T e to deviate from the proton temperature T i . The simulations used in this work only track the total internal energy density u gas , not the distinct electron and ion temperatures. We set T e after running the simulation according to local plasma parameters following the parameterization introduced by Mościbrodzka et al. (2016; see also Mościbrodzka et al. 2017 and EHTC V). The ratio between the ion and electron temperatures R is determined by the local plasma β = p gas /p mag , where p gas = (γ − 1)u gas , and p mag = B 2 /8π. The temperature ratio is then taken to be where R high (R low ) are the free parameters of the model and give the approximately constant temperature ratio at high (low) β. This approach allows us to associate the electron heating with magnetic properties of the plasma. In calculating the electron temperature, we further assume that the plasma is purely ionized hydrogen and that ions are nonrelativistic with an adiabatic index γ p = 5/3 and electrons are relativistic with γ e = 4/3. Then, given u gas from the simulation and R from Equation (15), (EHTC V): We note that this procedure is not entirely self-consistent, as the γ of the combined electron-ion fluid will change depending on the relative pressure contributions of electrons and protons while we assume it is fixed throughout the simulation domain. See Sądowski et al. (2017) for an alternative, self-consistent approach.
In this Letter, we consider a library of 72,000 simulated images composed of sets of 200 realizations of the same accretion system described by a fixed set of heating and observation parameters. Each set of 200 images is drawn from output files spaced by 25-50 r g /c from the set of 10 GRMHD simulations spanning five spin values in both MAD and SANE field configurations. The inclination angle for each image is set to one of either i = 12, 17, 22 deg (retrograde models, a * < 0) or i = 158, 163, 168 deg (prograde models, a * 0), according to which parity is required to orient the brightest portion of the ring in the southern part of the image while ensuring the orientation of the approaching jet is consistent with large-scale observations.
We use electron heating parameters R low = 1, 10 and R high = 1, 10, 20, 40, 80, or 160 in Equation (15) correspond to lower electron-to-proton temperature ratios in the low β regions (e.g., the jet funnel). This choice is physically motivated for M87 * , where radiative cooling of the electrons may keep T e < T i even in magnetized regions where electron heating is efficient (e.g., Mościbrodzka et al. 2011;Ryan et al. 2018;Chael et al. 2019). Lower electron temperatures in R low = 10 models increase the Faraday rotation depth and can result in increased depolarization in parts of the image.
GRMHD simulations produce a highly magnetized jet funnel above the black hole's poles, away from the accretion disk. In the funnel, where the plasma magnetization parameter σ ≡ B 2 /4πρc 2 ? 1, our numerical methods typically fail to accurately evolve the plasma internal energy. In the image library, we cut off all emission in regions where σ > 1 to ensure that we limit the emitting region to plasma whose internal energy is safely evolved without numerical artifacts (as in EHTC V). We tested the importance of a σ > 1 electron population by generating a supplementary set of images from all models with a cut at σ = 10 and found that it did not change the overall distribution of the derived metrics we use for model scoring in Section 5.
Each set of 200 model images with the same parameters in the image library requires a unique density scaling factor that is determined by matching the average flux density from the model to the observed compact flux density of M87 * measured by the EHT. Hence, the mass accretion rates, radiative efficiencies, and jet powers will differ between two models even if they are derived from the same underlying simulation (e.g., if R high , R low , or i are changed). The additional models discussed in Section 6, which explore the effects of different σ cutoff values and the inclusion nonthermal electrons, also require unique mass-scaling factors.
All of the polarimetric images from GRMHD simulations that we analyze in this Letter were generated using the ipole code (Mościbrodzka & Gammie 2018), which has been tested against analytic solutions and numerical ones produced by other numerical GRRT codes (Dexter 2016;Mościbrodzka 2020). A more comprehensive comparison of various GRRT codes that perform total intensity transport and fully polarized GRRT can be found in Gold et al. (2020) andB. Prather et al. (2021, in preparation), respectively. Preliminary results from B. Prather et al. (2021, in preparation) show that the tested codes are consistent at the fraction of 1% in all Stokes parameters. All calculated images in this work ignore light travel time delays through the emission region (the so-called "fast light" approach), and are calculated at a single frequency of ν = 230 GHz, neglecting the finite observing bandwidth of the EHT. We confirm that neither of those effects are important for models of interest for M87 * .

Sample GRMHD Model Images and Polarization Maps
In Figures 4 and 5 we show images and polarization maps for a subset of library models. In general, because the horizonscale magnetic fields in MAD models are strong enough not to be advected with the accretion flow, they are more likely to have a significant poloidal component and produce azimuthal EVPA patterns (Figure 3). In contrast, SANE models tend to show more radial EVPA patterns. Some MAD a * = 0.94 and SANE a * = 0 images are notable exceptions to this trend. These trends are also apparent in the distributions of the β 2 phase across the full image library that we consider later in Figure 9.
The GRMHD models at their native resolution include notable disorder in the EVPA structure, resulting from both magnetic turbulence and Faraday rotation. Models with larger R high have lower electron temperatures and higher Faraday rotation depths, resulting in the most disordered polarization maps. Many of the EVPA patterns seen in the images blurred with a 20 μas Gaussian kernel to simulate the limited EHT resolution resemble those from the idealized magnetic field models in Figure 3, indicating that the net EVPA pattern after blurring may trace the intrinsic magnetic field structure.
In Figure 6 we show a sample polarization map at full resolution compared to the same map blurred with circular Gaussian kernels of 10 μas and 20 μas FWHM. From tests with synthetic data, blurring (convolving with a circular Gaussian kernel) provides a reasonable approximation to image reconstruction from the EHT data at a comparable resolution (EHTC VII). The resolved average fractional polarization in the blurred images 〈|m|〉 traces the degree of order in the intrinsic polarization map. In the blurred images, disordered polarized structure on small scales produces beam depolarization. The degree of depolarization decreases with increasing spatial resolution (decreasing beam size).
The bottom row of Figure 6 shows the same unblurred and blurred polarization maps, but calculated without the effect of Faraday rotation (ρ V = 0). Those images show more coherent EVPA structure, with much larger |m| net and, particularly when blurred, much larger 〈|m|〉. Evidently, for this particular model, the depolarization visible in the corresponding top panels is due to Faraday rotation internal to the emission region. In addition, the net EVPA pattern shifts by a significant amount. The change in β 2 by ;80 deg would correspond to an apparent RM of ;−4 × 10 5 rad m −2 . Our GRRT calculations include all Faraday rotation occurring inside the GRMHD simulation domain (r out = 50-100 r g ), both external and internal to the 230 GHz emission region. The observables considered here, for the low viewing inclination of M87 * , do not depend strongly on that outer radius, as long as it is at r  40 r g . We cannot rule out the presence of additional Faraday rotating material at larger radii 100 r g , and its effects are not included in our models. Appendix B discusses the origin of the RM in our models in more detail.

GRMHD Model Theory Metrics
We compute the polarimetric observables (|m| net , |v| net , 〈|m|〉, β 2 ) described in Section 2.3 from model images blurred with a circular Gaussian kernel with a FWHM of 20 μas in order to compare them to the ranges measured from EHT and ALMA-only data. Both 〈|m|〉 and β 2 depend on the resolution and hence the size of the Gaussian blurring kernel. The value of β 2 also depends on the choice of the image center. We do not shift the library images before computing β m coefficients for comparison with the range inferred from the EHT image reconstructions, which have been centered by aligning them to the centered, fiducial total intensity images in EHTC IV. As discussed in Palumbo et al. (2020), a centering offset u expressed as a fraction of the diameter of a PWP m = 2 ring causes a quadratic falloff in β 2 power as δβ 2 /|β 2 | ≈ 4u 2 . Effects on the β 2 phase enter at similar order. In the case of the EHT image, u is likely less than one-fifth, meaning that centering errors in β 2 will be sub-dominant to other uncertainties, such as the choice of the blurring kernel or the variation across methods and days. Figure 7 (right panel) shows the resolved average polarization fraction 〈|m|〉 as a function of their image-averaged Faraday rotation depth, t á ñ r V . At small t á ñ r V , the average polarization fraction is 〈|m|〉 ; 20%-50%. Intrinsic disorder in the magnetic field structure due to turbulence is generally insufficient to produce the low observed image-average polarization fraction in EHT 2017 M87 * data (5.7% 〈|m|〉 10.7%). This is especially evident for the SANE models with prograde black hole spin, which have the highest resolved polarization fractions. At large t á ñ r V , strong scrambling from internal Faraday rotation typically results in small predicted polarization fractions of <5% at the scale of the EHT beam.
The clear exceptions to this trend are some SANE retrograde models (a * = −0.9375 for large R high ), which show 〈|m|〉 ; 10%-20% despite their large ⟨ ⟩ t r  10 3 V . In these models, most of the observed polarized flux originates in the forward . Sample snapshot false-color images and polarization maps for a subset of the models in the EHT M87 * simulation image library at their native resolution (top three rows) and blurred with a 20 μas circular Gaussian beam (bottom three rows). The inclination angle for all images is either 17 deg (for negative a * models) or 163 deg (for positive a * model), with the black hole spin vector pointing to the left and away from the observer. The tick length is proportional to the polarized flux, saturated at 0.5 of the maximum value in each panel. Here models with R low = 1 are shown. In general, the EVPA pattern is predominantly azimuthal for MAD models (e.g., MAD a * = 0 R high = 1) and radial for SANE models (e.g., SANE a * = 0.94 R high = 1), although the SANE a * = 0 models in particular are exceptions to this trend. All models show scrambling in the polarization structure on small scales from internal Faraday rotation, with more pronounced scrambling in models with cooler electrons (larger R high parameter).
jet, while most of the computed Faraday depth is accumulated near the midplane. Photons that travel from the forward jet to the observer do not encounter the large Faraday depth. For similar reasons, the inferred RM can be much lower than implied by their large values of integrated t r V .
Distributions of all observables are shown in Figure 7 (〈|m|〉, left panel), Figure 8 (|m| net and |v| net ), and Figure 9 (|β 2 | and ∠β 2 ). SANE models tend to have a lower integrated polarization fraction and larger circular polarization fraction than M87 * at 230 GHz. In many cases this is a result of very large Faraday rotation internal to the emission region. MAD models tend to have larger net linear polarization fraction than observed in M87 * . The resolved average fractional polarization produces similar trends. Most SANE models with prograde spin are too scrambled and most MAD models are too ordered compared to the reconstructed polarization maps of M87 * . Full distributions for all models, including their R high , R low , and a * dependence, are further discussed in Appendix C.

Model Constraints from Polarimetry
To evaluate whether a given GRMHD model is consistent with the EHT observations reported in EHTC VII, we require  Figure 4 but for R low = 10. We find similar trends, but with more scrambling from larger Faraday depths due to lower electron temperatures. images from the model to satisfy constraints on the four parameters derived from the reconstructed EHT images and ALMA-only measurements presented in Table 2 and summarized again here.
1. The image-integrated net linear polarization |m| net is in the measured range from the EHT image reconstructions: 1% |m| net 3.7%.
In comparing models to observables, the β 2 metric is the most constraining. Only 790 snapshot images out of the 72,000 considered fall in the range of those reconstructed in both β 2 amplitude and phase, compared to 11,526 snapshots for both |m| net and |v| net and 7,727 for the resolved image-average linear polarization fraction 〈|m|〉.
Below we explore two quantitative methods for scoring models, either by requiring that at least one single snapshot image from a model simultaneously passes all constraints ("simultaneous scoring;" see Section 5.2) or that each observational constraint is satisfied by at least one snapshot image from a given model ("joint scoring;" see Section 5.3).

Simultaneous Snapshot Model Scoring
In the simultaneous scoring procedure, we rule out models where none of the 600 snapshot images (200 time samples at three inclination angles) can simultaneously satisfy the constraints on all of the polarimetric observables. Only 73 out of 72,000 snapshot images across 15 of 120 models simultaneously pass all of the constraints. Of those, all but two viable snapshot images come from a MAD model. The only models with more than five passing images are MAD a * = 0 R low = 1 R high = 160 and MAD a * = −0.5 R low = 1 R high = 80, 160. Figure 10 shows three viable snapshot images from both SANE and MAD models as well as three snapshot images from models ruled out by simultaneous scoring (i.e., with no snapshots in the entire sample from the model simultaneously satisfying all constraints). These images are representative of the snapshots that simultaneously satisfy all constraints on the image-integrated metrics; they have not been selected based on detailed matching of the resolved polarization structure to the EHT images. Nonetheless, the top row of images show good qualitative agreement with the primary features of the EHT image in Figure 1. In contrast, the snapshots from the ruled-out models tend to be too polarized, too depolarized, or too radial  in their EVPA pattern. Figure 11 shows the distributions of |β 2 | for all 600 snapshots from these three passing and three failing models. Although variability is present, the systematic differences over the five observables considered allow us to constrain the models. The left panel of Figure 12 shows the total number of images that pass simultaneous scoring as a function of model, summing over the six R high values.

Joint Distribution Model Scoring
In the joint scoring procedure, we use the measured distributions of the data metrics to ask whether the observed value of each metric for M87 * is consistent with being drawn from the distribution seen in the GRMHD simulations. To do this, we measure χ 2 values for the five metrics for all snapshots k from a given model as 2 where x j,k are the values of a scoring metric x j for each of the 600 snapshots k from a given model,x j is the mean of those values for the model, and σ j is taken as one half of the observed data range from . The joint likelihood of each model is the product = P   j j of those for the five metrics x j . To produce a non-zero likelihood  in this method, at least one snapshot from a model must lie further from its mean than the data value does. That can be a different snapshot for each metric, which makes this method more lenient than the simultaneous scoring method. We also note that snapshots are allowed to have the wrong sign of the difference with their mean, due to the definition of χ 2 and our use of the mean of the model snapshots themselves. In practice, this makes little difference in the results.
In this method, we consider models viable whose joint likelihood is >1% of the maximum found from any model. The right panel of Figure 12 shows the resulting joint likelihoods summed over R high .

Comparison of Scoring Results
The results of both scoring procedures are summarized in Figure 12, summed over R high . Both scoring methods prefer MAD models to SANE models, with most of the passing models coming from the MAD a * = 0 and a * =±0.5 simulations.
The main difference between the two scoring procedures is that joint scoring prefers R low = 10 models, while R low = 1 is preferred by simultaneous scoring. SANE models with a * = 0.94, R low = 1, 10, and R high = 10 are ruled out by simultaneous scoring, but score fairly well in joint scoring. For the favored MAD models, when R low = 1, there are more Figure 10. Sample 230 GHz image library polarization maps shown in the same style as the EHT image in Figure 1. All maps are blurred with a 20 μas circular Gaussian beam. In all images, the inclination angle is either 17 deg (for negative a * models) or 163 deg (for positive a * model), with the black hole spin vector pointing to the left and away from the observer. The top row displays SANE (a * = 0, R high = 80) and two MAD snapshots (both a * = −0.5 and R high = 160) from left to right. All of the top row images satisfy simultaneous constraints on image-integrated metrics (|m| net , |v| net , 〈|m|〉, |β 2 |, ∠β 2 ) derived from the EHT2017 results. The bottom row displays snapshots from models that fail to produce any images that simultaneously satisfy the observational constraints. These snapshots are from two SANE (a * = 0.5 and R high = 1 and 160) and one MAD (a * = 0.94, R high = 160) model, from left to right. The failing images are either more polarized than the data (left and right panels) or too depolarized (middle panel). They also fail to match the observed EVPA pattern (β 2 phase).
images that simultaneously satisfy all constraints, but when R low = 10, the distributions generally stay closer to the observed data ranges and are thus favored by the joint scoring method. Due to differences between simultaneous and joint scoring results, as well as other methods we have tried, we consider the inferred parameters of R low , R high , and a * from passing models to be less robust than the overall trend that MAD models are favored.
The simultaneous scoring method has the advantage of conceptual simplicity, and of applying each constraint simultaneously per image (thus accounting for correlations between the scoring metrics). Simultaneous scoring is more strict and rules out more models than joint scoring, but it may be more limited by the finite number of images generated per model. The joint scoring procedure has the advantage of being more conservative in disfavoring models, but assumes the observational constraints are independent in calculating a joint likelihood. Instead, they are correlated (in particular |m| net , 〈|m|〉, and |β 2 |).
The number of images in each model that pass each constraint individually (used in joint scoring) and that simultaneously pass all constraints (used in simultaneous scoring) are presented in Appendix D.

Combined EHTC V and Current Polarimetric Constraints
EHTC V presented constraints on the GRMHD simulation models based on fits to the EHT total intensity data, model Figure 12. Results of the simultaneous (left panel) and joint (right panel) scoring methods for comparing GRMHD models to M87 * observables. The simultaneous scoring method shows the total number of viable images for each image library model after summing over R high . Out of a total of 73 passing images, only two are from a SANE model. The right panel shows the joint likelihood of each library model after summing over R high . In this method, R low = 10 MAD models are preferred and SANE a * = +0.94 R high = 10 models are also allowed. Figure 11. Distributions of |β 2 | for the sample passing and failing models in Figure 10. The dashed black lines mark the measured values for the snapshot images in Figure 10, and the blue bands show the range inferred from EHT M87 * data. The models can be constrained using EHT observables even in the presence of significant scatter due to time variability. self-consistency (requiring a radiative efficiency less than that of a thin accretion disk at the same black hole spin), and M87's measured jet power (requiring a simulation to produce a jet power consistent with a conservative lower limit of that from M87 * , >10 42 erg s −1 ). Those constraints ruled out MAD a * = −0.94 models (from failing to satisfy the EHT image morphology), SANE models with a * = −0.5, and all models with a * = 0 (from failing to produce enough jet power). Here we retain only the jet power lower limit, which is the most constraining and straightforward to apply to the expanded image library considered in this work.
Relativistic jets launched in GRMHD simulations (defined here as in EHTC V, with a cutoff of βγ > 1) are fully consistent with being produced via the Blandford−Znajek process (e.g., McKinney & Gammie 2004;McKinney 2006). As a result, a * = 0 models have small or zero jet power, P jet , and are rejected by this constraint. These models can still produce significant total outflow powers (P out in EHTC V) in a mildly relativistic jet or wind. Many other models with low values of R high or moderate black hole spin, including those of SANEs with a * = + 0.94, which are allowed by the joint scoring procedure, are also ruled out by the jet power constraint (see Table 3 in Appendix D and EHTC V). Combining the simultaneous scoring polarimetric constraints with the jet power constraint results in 15 remaining viable models: all MADs, and spanning the full range of non-zero a * explored. This conclusion does not depend on the choice of the simultaneous or joint model-scoring procedure.

Discussion
The resolved EHT 2017 linear polarization map of M87 * shows a predominantly azimuthal linear polarization (EVPA) pattern, and relatively low fractional polarization of 20% on 20 μas scales. We interpret the low fractional polarization as the result of Faraday rotation internal to the emission region, which acts to rotate, scramble, and depolarize the resolved polarized emission. Adopting this constraint in a one-zone model, we estimate typical values of particle density n e , magnetic field strength B, and electron temperature T e . In semi-analytic emission models with externally imposed, idealized magnetic field configurations, azimuthally dominated EVPA patterns are produced by poloidal (radial and/or vertical) magnetic field components. To fully capture the complicated combined effects from magnetic field structure, turbulence, relativity, and Faraday rotation on polarimetric images of M87 * , we turn to radiative transfer calculations from GRMHD simulations.
We compared a large image library of emission models from GRMHD simulations with metrics designed to capture these salient features of the data. The combined constraints of a predominantly azimuthal EVPA pattern and a low but non-zero fractional polarization are inconsistent with most SANE GRMHD models with weaker horizon-scale magnetic fields. Some MAD models with relatively cold electrons, realized in our library by larger values of R high and/or R low , remain consistent with the data. Here we discuss the implications of our results, and limitations in our set of theoretical models that may impact our interpretation.

Near-horizon Plasma and Magnetic Field Properties in Passing Models
Both our one-zone and GRMHD models find similar plasma conditions in the 230 GHz emission region, driven by the requirements of weak 230 GHz absorption and strong 230 GHz Faraday rotation. In viable GRMHD models, we find average, intensity-weighted plasma properties in the emission region of n e ∼ 10 4-5 cm −3 , B ; 7-30 G, and θ e ∼ 8-60. These are in good agreement with our one-zone estimates (Section 3.1). We have also calculated the intensity-weighted values of the absorption and Faraday optical depth, τ I and t r V , over snapshots that simultaneously satisfy all our observational constraints. The median values are τ I ; 0.1 and  t r 50 V . All of our viable images have t p > r 2 V , while 2 out of 73 have τ I  1, consistent with our assumptions in Section 3.1 that the plasma Faraday depth is large while the Stokes  optical depth is small.
By quantitatively evaluating a large library of images based on GRMHD models (Section 5), we identify 25 out of 120 models that remain viable after applying constraints based only on EHT and ALMA-only polarimetric observations. Additionally applying a cut on jet power of P jet > 10 42 erg s −1 (EHTC V) rules out the five viable SANE models and all a * = 0 models. The precise number and identity of the viable models depends mildly on the chosen scoring procedure and on the Gaussian blurring kernel size applied to the EHT image reconstructions and library simulated images. The overall preference for MAD over SANE models is found from both the simultaneous and joint scoring procedures, as well as other variants. After applying the jet power constraint, no viable SANE models remain for any of the scoring methods that we explored.
MAD models are associated with dynamically important magnetic fields. The significant poloidal components of those fields can produce a predominantly azimuthal polarization pattern (Figure 4), similar to those seen in idealized models with prescribed poloidal magnetic fields (Figure 3). Strong Faraday effects complicate a direct interpretation of the observed EHT polarization map in terms of those idealized models. Still, our more detailed comparison favoring MADs suggests the presence of dynamically important magnetic fields in the emission region on event-horizon scales.
In Figure 13 we present mass accretion rate and jet power distributions both for the viable models identified in EHTC V and when adopting the new constraints from polarimetry. 131 Polarimetric constraints break degeneracies present in the single epoch total intensity data, allowing us to estimate a mass accretion rate onto the black hole of ( -)  Ḿ 3 20 Edd -10 6 , where  M Edd is the Eddington accretion rate. 132 The measured radiative efficiency  =  L Mc 2 (where L is the bolometric luminosity) of the passing models is relatively high for a hot accretion flow model: ò  1%. These models have jet powers of P jet ; 10 42-43 erg s −1 .
The mass accretion rate found here is much lower than the Bondi rate calculated from Chandra observations 131 Note differences in some  M and P jet values compared to EHTC V. We have corrected minor tabulation errors from that work, and have used a slightly different time range for averaging the MAD a * = −0.94 simulation. 132 The Eddington rate is defined as  =  M L c Edd Edd Edd 2 , where L Edd = 4πGMm p c/σ T is the Eddington luminosity and we adopt an efficiency factor ò Edd = 0.1. Note that this assumed efficiency factor ò Edd is distinct from the reported radiative efficiency  =  L Mc 2 measured from the simulations. Matteo et al. 2003; see also Russell et al. 2015), and higher than that found from hybrid disk+jet models of the M87 * spectral energy distribution (SED; Prieto et al. 2016). Our inferred jet powers of 10 43 erg s −1 are toward the lower end of the observed range. In particular, the jet power measured at the location of Hubble Space Telescope (HST)-1 is ∼10 43-44 erg s −1 (Stawarz et al. 2006), and Low-Frequency Array (LOFAR) observations suggest that a jet power of ∼10 44 erg s −1 was necessary within the last ∼million years to inflate the observed radio lobes on scales of ∼80 kpc (de Gasperin et al. 2012).

(Di
Measurements of the accretion rate and the radiative efficiency can begin to constrain the microphysical plasma processes that heat electrons in M87 * , for example by inferring the fraction of the dissipated energy in the system that heats electrons, δ e . In axisymmetric, self-similar, hot accretion flow models, a system with   -M M 10 5 Edd and a radiative efficiency ò  1% has a value of δ e in the range 0.1-0.5 (see Figure 2 of Yuan & Narayan 2014). This range is consistent with that produced by simulations of turbulence and reconnection in the β ∼ 1 regime (e.g., Rowan et al. 2017;Werner et al. 2018;Kawazura et al. 2019). Future studies using simulations with self-consistent electron heating and radiative cooling (Section 6.3) can better constrain δ e and its dependence on local plasma parameters throughout the accretion flow and jet-launching region.
We have assumed that all effects responsible for the appearance of the EHT polarized image of M87 * are captured within the relatively small GRMHD simulation spatial domain, 10 2-3 r g . Goddi et al. (2021) developed a two-component model for the ALMA and image-integrated EHT data where each component is Faraday rotated by a different screen. The model demonstrates that the rotation measure of the compact component is unconstrained by the ALMA measurements alone, as the ALMA measurements are also sensitive to the Faraday rotation properties of the larger-scale component. In addition, the observed time variability in ALMA data (e.g., the RM sign change) can be explained by the observed EVPA variation of the compact core seen by the EHT. To produce the observed variability requires an RM of ≈−6 × 10 5 rad m −2 .
The ALMA data do not constrain the location or nature of this Faraday screen, except that it must be relatively close to the compact core, r  10 5 r g .
For our favored plasma parameters for M87 * , we expect substantial Faraday RM internal to the emission region itself, t p r  2 V , consistent with that measured from viable GRMHD images. In a model of uniform, external Faraday rotation this Faraday depth at 230 GHz would correspond to an RM of 10 6 rad m −2 . Figure 14 shows that the apparent RMs measured from our GRMHD images span a wide range, often comparable to or larger than that inferred from the Goddi et al. (2021) two-component model (10 6 rad m −2 ). For the low inclination angle of M87 * , the apparent RM measured from GRMHD images is not a good tracer of the mass accretion rate (Mościbrodzka et al. 2017), and originates close to the emission region and well within the simulation domain (Ricarte et al. 2020, and Appendix B). The RM inferred from lowinclination GRMHD models of M87 * can also vary rapidly and change signs (Ricarte et al. 2020), as seen in the ALMA-only data. As a result, the RM inferred from the two-component model in Goddi et al. (2021) is apparently consistent with the intrinsic properties of the GRMHD models studied here, without invoking an additional, external Faraday screen. At the same time, we cannot rule out that such an external screen could be present. Future EHT observations with wider frequency spacing can directly measure the resolved RM of the core and address this uncertainty.

Electron Acceleration
Magnetic reconnection, magnetohydrodynamic (MHD) turbulence and collective plasma modes in collisionless hot accretion flows likely result in nonthermal particle acceleration. While pioneering attempts have been made (e.g., Ball et al. 2016;Chael et al. 2017;Davelaar et al. 2019), it is not yet known how to properly incorporate electron acceleration in global GRMHD simulations of hot accretion flows.
We adopt an empirical approach to investigate the impact of nonthermal (accelerated) electrons on 230 GHz polarimetric images of M87 * to quantify whether and how neglecting particle acceleration in our models affects our conclusions. We use a specific, but widely explored (e.g., Özel et al. 2000;Markoff et al. 2001;Yuan et al. 2003), description for electron acceleration, namely that accelerated electrons add a power-law tail to the thermal distribution function. The power-law tail is described by the fraction of the thermal energy density in the power-law tail, η, the power-law slope, p, and the maximum Lorentz factor of the accelerated electrons, g max . The minimum Lorentz factor, g min , is calculated self-consistently by assuming a continuous transition between the thermal and power-law distribution functions (e.g., Yuan et al. 2003). In this model, the parameters p, η, and g max are constant across the accretion flow. We assume that g  ¥ max and we explore values of p = 3.5, 4.5 and η = 0.01, 0.1.
However, fixing the accretion rate to that used in the thermal model results in an increased total intensity flux density when we add high-energy nonthermal electrons to our models. If instead we compare the models at fixed flux density, we need to reduce the mass accretion rate of the hybrid model. Therefore, generalizing the distribution function introduces order unity uncertainties in the inferred mass accretion rate, radiative efficiency, and jet power. The changes in the polarimetric observables in a given snapshot are also larger at fixed flux density. For example, in the MAD a * = −0.5 model, the |m| net increases from 4.7% to 6% when adding nonthermal electrons. As a result, in principle, polarimetric observables constrain the distribution function as well.
Such constraints presumably depend on the details of the assumed particle acceleration scenario. Viable scenarios include hybrid electron distribution functions, or models where particle acceleration is a function of local gas conditions or magnetic energy density rather than fixed throughout the flow (e.g., Ball et al. 2016;Davelaar et al. 2019). More realistic particle acceleration scenarios could be considered using resistive GRMHD simulations (e.g., Ripperda et al. 2020).

Coherently Polarized Forward Jet Emission
As discussed above, some SANE retrograde model images in the library show coherently polarized features even when the Faraday depth through the entire emission region is large. The observed polarized flux in those cases originates on the near side of the midplane and is not scrambled from Faraday rotation along the line of sight. A similar effect might be possible if nonthermal electrons could be accelerated efficiently in the low-density, strongly magnetized funnel region in front of the black hole.
It is beyond the scope of this Letter to evaluate whether or how such a model might be realized physically, e.g., whether any process could fill the funnel with high-energy electrons efficiently enough to produce the observed 230 GHz luminosity from the funnel alone. Instead we carry out one sample calculation of polarized emission from the funnel of a prograde a * = 0.94 SANE library snapshot. We assign a nonthermal energy density u nth = αu mag wherever the magnetization σ > 1, with α = 0.02 the fraction of the magnetic energy density u mag that is put into nonthermal particles. We calculate synchrotron Figure 14. Absolute value of RM vs. net linear polarization |m| net for a subset of our EHT GRMHD library models explored in more detail in Ricarte et al. (2020). Closed symbols represent positive RM while open symbols represent negative RM, revealing significant time variability across the 2500 r g /c spanned by these snapshots. In gray, we plot our allowed region of |m| net and bracket the range of core RM inferred from contemporaneous ALMA-only observations, 2-100 × 10 4 rad m −2 (Goddi et al. 2021). The dashed horizontal line demarcates the RM at which an EVPA rotation by π radians would have been observed between the 212 and 230 GHz frequency range used in the ALMA-only measurements, 1.05 × 10 7 rad m −2 . Despite large Faraday depths, a large fraction of these snapshots exhibit RMs consistent with simultaneous ALMA-only constraints. RM and |m| net are anti-correlated, as larger Faraday depths lead to greater scrambling of the intrinsic polarization. Figure 15. Sample polarization maps with varying electron distribution function. Columns display single snapshots from three selected models. Row 1 shows images with a thermal electron distribution function, as assumed in the standard EHT image library. Rows 2 through 5 are the same models but with emission from a hybrid distribution of electrons. Row 6 shows a hybrid model but the mass accretion rate of the model is adjusted to reproduce the same total intensity flux as the purely thermal snapshot. All maps are blurred with a 20 μas circular Gaussian. In all images, i = 17 deg (for negative a * models) or i = 163 deg (for positive a * model), with the black hole spin vector pointing to the left and away from the observer. radiation from a pure power-law distribution of electrons with g = 100 min and p = 3. Figure 16 compares the original thermal snapshot with two realizations of this hybrid thermal+nonthermal funnel emission models. In the purely thermal case, Faraday rotation depolarizes the emission at the EHT beam scale, producing low fractional polarization across the image that is inconsistent with EHT observations of M87 * . Adding power-law electrons in the funnel produces coherent linearly polarized emission. When we assume u nth ∝ u mag (middle panel), the power-law emission is concentrated close to the black hole and lensed into a ring (Dexter et al. 2012). The weak forward jet component is strongly polarized but lies inside the observed ring, and is thus potentially inconsistent with the EHT total intensity and polarimetric image. In the right panel, we exclude nonthermal emission from inside a radius r cut = 6r g . Both nonthermal maps are consistent with our cuts on net and image-average linear polarization fraction, |m| net and 〈|m|〉. However, both are inconsistent with the observed EVPA pattern of M87 * (i.e., the β 2 phase).
For this example, we assume a plasma of protons and electrons rather than e + /e − pairs. The latter are presumably more likely to form in the funnel (Mościbrodzka et al. 2011;Chen et al. 2018;Levinson & Cerutti 2018;Anantua et al. 2020;Crinquand et al. 2020;Wong et al. 2021), and have different circular polarization properties. Future observations that constrain the resolved circular polarization structure could potentially discriminate between pair and electron-ion plasmas in the emitting region. At longer wavelengths and larger scales, the limb-brightened jet structure of M87 (e.g., Walker et al. 2018) also suggests that the radiating electrons are not concentrated inside the funnel as modeled here.

Radiative Models
Our GRMHD images use the parameterization of Mościbrodzka et al. (2016) to model the electron and ion temperatures given the total gas temperature from a simulation. In this prescription, the electron-to-ion temperature ratio is a function entirely of the local plasma β. This functional form (Equation (15)) captures the general behavior seen in many simulations of electron heating in turbulent or reconnecting collisonless plasmas; namely, the electron heating is more efficient (and thus the temperature ratio is closer to unity) when β < 1 (e.g., Howes 2010; Rowan et al. 2017). However, the actual distribution of T e in a hot accretion flow reflects the balance of heating, cooling, and advection of hot electrons throughout the system. Furthermore, the GRMHD simulations in the library considered here do not include radiative cooling. Our passing models for M87 * favor a radiative efficiency of ò ∼ 1% (Section 6.1), however, and we may begin to worry if cooling is dynamically important in M87 * .
To assess these uncertainties, it will be useful to compare the results in this work with results from simulations performed with radiative GRMHD codes. These codes typically use either the M1 closure method (e.g., Sądowski et al. 2013;McKinney et al. 2014;Sądowski et al. 2017) or a Monte Carlo approach (e.g., Ryan et al. 2015) to track radiation and its interactions with the plasma near the black hole. In addition to the effects of cooling on the gas temperature, these codes can also evolve the separate electron and ion temperatures under the influence of cooling and different subgrid prescriptions for the electron heating efficiency (e.g Ressler et al. 2015Ressler et al. , 2017Chael et al. 2018Chael et al. , 2019Ryan et al. 2018;Dexter et al. 2020).
In Figure 17, we show a comparison of the temperature ratio T i /T e obtained directly from the two MAD radiative simulations of M87 * in Chael et al. (2019;right column), with the temperature ratio obtained from the same simulations using Equation (15) with R low = 1, R high = 20 (left column). The two rows show simulations using different underlying models for electron heating: the top row (H10) uses the turbulent heating prescription of Howes (2010), and the bottom row (R17) uses the reconnection heating prescription of Rowan et al. (2017). The simulation data in both rows is averaged in time and azimuth. Figure 17 shows that the temperature ratios obtained from the electron-ion evolution and from the postprocessing prescription both transition from smaller to larger values when moving from the jet/funnel region (low β) to the disk (high β). In simulation H10 (top row), T i /T e ≈ 1 in the funnel, which well matches the result from Equation (15) with R low = 1; in simulation R17, electrons are cooler in the funnel (T i /T e ≈ 5), and even colder in the disk-jet interface directly outside the σ = 1 contour (T i /T e ≈ 15). The cooler electrons in the funnel Figure 16. Sample library SANE a * = 0.94 R high = 160 snapshot (left panel) and with power-law emission from nonthermal electrons added in σ > 1 regions (in the funnel near the pole, middle panel) and in the funnel only outside of a radius r cut = 6 r g (right panel). The library model is heavily depolarized due to Faraday rotation. Nonthermal radiation from the forward jet is coherently polarized; these images look qualitatively similar to polarimetric images of the forward-jet dominated semianalytic models of Broderick & Loeb (2009). Even a small total flux contribution can increase the net and image-averaged linear polarization fractions to lie within the observed range. However, in this example the EVPA patterns (β 2 phase) of the TH+PL images remain inconsistent with M87 * data. regions of R17 reflect the decreased efficiency of low-β electron heating in the Rowan et al. (2019) model than in Howes (2010). In Section 5.4, we note that the results of joint scoring (but not simultaneous scoring) favor R low = 10 models over R low = 1 models. Model R17 shows that some radiative models can naturally produce T i /T e > 1 in low-β regions. However, the funnel value of T i /T e ≈ 5 in this simulation is in between the two cases R low = 1, R low = 10 that we considered in the image library.
While the disk electrons are cooler than the jet electrons in both radiative simulations shown in Figure 17, T i /T e in the disk of model H10 increases at small radii. In contrast, model R17 better matches the postprocessing model prediction of a roughly constant disk temperature ratio. (Note that the 230 GHz emission is entirely produced at radii r  10r g .) In addition to the MAD simulations from Chael et al. (2019) in Figure 17, we have also checked the temperature ratios in the radiative SANE simulations of Ryan et al. (2018). In these simulations, the average temperature ratio in the EHT 230 GHz emission region can also be roughly approximated by the Mościbrodzka et al. (2016) prescription, with R low = 1 and R high in the range R high ∼ 10-20.
Our preliminary results indicate that while some important features of the temperature-ratio distributions produced in radiative simulations can be described by the R low , R high model (Equation (15)), the current postprocessing model cannot capture all of the behavior produced in radiative simulations. A more detailed comparison is left for future work.

Predictions
We have identified a subset of a large parameter space of GRMHD models that is consistent with constraints derived from current EHT total intensity and polarimetric observations of M87 * . The models that pass our constraints on the polarimetric structure and jet power from M87 * are all magnetically arrested (MAD) accretion flows. Here we make predictions for testing our interpretation with future observations.

Repeated Observations
Repeated EHT observations of M87 * at 230 GHz (both in total intensity, e.g., Wielgus et al. 2020, and in linear polarization) will continue to constrain the model parameter space. Figure 18 shows the time evolution of β 2 amplitude and phase for 200 snapshots of three viable library models: MAD a * = −0.5, R low = 10, R high = 20; MAD a * = +0.5, R low = 10, R high = 80; and MAD a * = +0.94, R low = 10, R high = 80. The observer inclination was 17 deg and 163 deg for the retrograde and prograde models, respectively.
Both quantities show variations on timescales from days to months. The phase and amplitude of β 2 should change over the course of a week of observations. In EHTC VII, we observe changes in the the β 2 amplitude and phase over the week of observations in 2017, and use the results from two epochs to define our acceptable parameter ranges. Figure 18 suggests that occasionally the observed changes in β 2 on ∼week timescales can be much more dramatic than we observe in 2017, with variations in β 2 phase of 90 deg for some models on short timescales.
The scatter in both quantities on longer ∼month timescales is much larger than the uncertainty range derived from the EHT 2017 measurements. If our passing GRMHD models accurately describe the 230 GHz emitting region in M87 * , future EHT observations should detect variability in the polarization structure. According to current models, the time-averaged β 2 amplitude and 〈|m|〉 should remain similar to the current values for prograde spin models, and tend toward larger values for retrograde spin models. For high-prograde spin (or many SANE models), the β 2 phase should on average be closer to zero than we observe in 2017.

Future Observations at 260 and 345 GHz
In selecting models, we have focused on metrics corresponding to salient features of the data. We have not attempted to compare models in detail to specific features of the reconstructed polarimetric images, most notably the apparently depolarized bright patch in the eastern part of the image (Figure 1). We do note that such depolarized features occur in many of our library images, particularly in MAD models with R low = 10. If the eastern patch in the 2017 image is depolarized due to Faraday rotation, it may be possible to tell with future higher-frequency observations. Figure 19 shows images of two of the favored MAD models at the current observing frequency of 230 GHz and at two additional frequencies planned for the future EHT observing campaigns. In addition to internal Figure 17. Comparison of the temperature ratio obtained directly from two radiative GRMHD simulations with the values obtained from the postprocessing model used in this work. The top and bottom rows show, respectively, time-and azimuth-averaged data from the two radiative MAD simulations of M87 * presented in Chael et al. (2019). The top row shows data from simulation H10, where electrons are heated by the turbulent heating prescription of Howes (2010), and the bottom row shows the R17 simulation, where electrons are heated by the magnetic reconnection prescription of Rowan et al. (2017). The left column shows the spatial distribution of the ion-to-electron temperature ratio obtained from these simulations by applying our β-dependent postprocessing equation to the total gas temperature (Equation (15) with R low = 1, R high = 20). The right column shows the temperature ratio obtained directly from the independently evolved electron and ion species in the radiative simulations. The white contour shows the σ = 1 surface.
Faraday rotation, the sense of the EVPA pattern may also be subject to a net, coherent rotation due to external Faraday rotation. At higher frequency, Faraday rotation is suppressed and EHT observations will see the intrinsic magnetic field pattern more clearly.
For a snapshot from the MAD a * = −0.5 (R high = 160 and R low = 1) model, Figure 19 shows that the |m| net and 〈|m|〉 values are predicted to increase with frequency. The 230, 260, and 345 GHz net EVPAs are −77, −70, and −82 deg, respectively, corresponding to an (apparent) rotation measure RM ∼ 1 × 10 5 rad m −2 between 230 and 345 GHz. The net circular polarization |v| net remains small and nearly constant with frequency; it is 0.42%, 0.35%, 0.32% for 230, 260, and 345 GHz, respectively. A similar trend is observed in a MAD Figure 18. Amplitude (left panel) and phase (right panel) of β 2 as a function of time for three viable GRMHD library models identified here (points, all with R low = 10) compared to ranges measured from EHT 2017 M87 * data (gray shaded region). The dashed lines show the median values for each model. The retrograde spin model predicts higher β 2 amplitude in future observations. In the high-prograde spin model, the median β 2 phase is closer to zero than the observed range in 2017. Changes in both quantities occur on timescales of weeks to months, and should be apparent in future EHT data sets. Figure 19. Random snapshot of MAD models with a * = −0.5 (top row) and with a * = 0.5 (bottom row) at the current EHT frequency of 230 GHz and two higher frequencies of 260 and 345 GHz, which are planned for the future EHT observations. All images are convolved with a 20 μas Gaussian. In all images the black hole spin vector is pointing to the left and away from the observer. In all cases, the ring fractional polarization increases slightly with frequency. The EVPA pattern, as measured by the β 2 , is similar at all three frequencies.
a * = 0.5 (R high = 80 and R low = 10) model. The image |m| net and average polarization 〈|m|〉 are again expected to increase with frequency. The corresponding net EVPAs are −38, −40, and −30 deg, corresponding to an apparent rotation measure of RM ∼ −1 × 10 5 rad m −2 . The net circular polarization fraction |v| net remains roughly constant and close to zero, 0.33, 0.06, and 0.2% from low to high frequency.
Both of these models display similar EVPA structure at all three frequencies, indicating that in this example the net EVPA pattern is due to magnetic field structure rather than coherent Faraday rotation. Future multi-frequency observations will be able to infer the core RM and intrinsic EVPA pattern set by the near-horizon magnetic fields.

Conclusions
The EHT has produced resolved polarized intensity maps in the near-horizon region around the supermassive black hole in M87. Taken together with image-integrated data from simultaneous observations with ALMA, these images constrain the space of accretion flow and jet models used to interpret the EHT total intensity image with broad implications for jet launching near a black hole event horizon. Here we summarize the main results of that analysis.
1. We interpret the depolarization seen in EHT images as the result of beam depolarization due to Faraday rotation internal to the emission region. In the context of one-zone models and combined with the size and brightness temperature of the total intensity image, we estimate an average emission region plasma density of n e ∼ 10 4-7 cm −3 , magnetic field strength of B ∼ (1-30) G, and T e = (1-12) × 10 10 K. 2. The net EVPA pattern of the M87 * polarization maps is predominantly azimuthal. In the context of semi-analytic models with imposed, idealized magnetic field geometry, such a pattern can be reproduced using a significant component of poloidal (radial and/or vertical) magnetic field. The presence of such magnetic fields in a rotating fluid would imply that the magnetic fields are dynamically important. However, significant Faraday rotation may be present and it is not clear whether the observed EVPA pattern can be interpreted in terms of magnetic field structure alone. 3. To capture the effects of realistic magnetic field structure, plasma conditions, and Faraday rotation and conversion, we have compared salient observables to a large library of images from the GRMHD simulation library of EHTC V. The observables are the net circular polarization fraction constrained by ALMA (|v| net ), the net and imageaveraged linear polarization fraction measured by the EHT (|m| net and 〈|m|〉), and the m = 2 coefficient of a Fourier expansion of the azimuthal EVPA pattern (β 2 ). Of these, β 2 is the most constraining metric. 4. The model-scoring procedures disfavor most models from the GRMHD image library from polarimetric observations alone. Many weakly magnetized (SANE) models are too depolarized, or show an EVPA pattern that is too radial. Many strongly magnetized (MAD) models are too coherently polarized. The polarization fraction is generally set by the Faraday rotation depth close to the emission region. MAD models more frequently produce azimuthal EVPA patterns, as expected for magnetic field structures that include a significant poloidal field component. Combined with a conservative lower limit on the jet power of M87, only strongly magnetized (MAD) models remain viable. We use those remaining models to estimate the mass accretion rate onto the central supermassive black hole as  = M ( -) 10 yr 4 1 . The average plasma parameters found from GRMHD images are in good agreement with those inferred from one-zone models. 5. The model space considered in this Letter is incomplete, and systematic uncertainties remain a challenge. While the radiative efficiency that we find is relatively high, we consider only non-radiative GRMHD models. We do not consider GRMHD models with misalignment between the disk and the black hole angular momentum. We also only consider one parameterization for determining the electron distribution function from the simulation data. Of these three major areas of uncertainty, we have explored a small sample of alternative models for determining the electron distribution function, including both alternative prescriptions for electron heating in strongly magnetized regions and including a nonthermal component. The quantitative estimates of mass accretion rate and jet power found here depend on the assumed electron distribution function and are uncertain at the order unity level. The alternative electron distribution functions considered here do not change the main finding that MAD models with dynamically important near-horizon magnetic fields appear more viable for explaining the first polarimetric EHT observations of M87 * . 6. Our favored models show time variability in the polarization metrics used here. The median values found at several epochs should be sufficiently well measured to distinguish between the current retrograde and prograde spin models. At higher frequencies of 260 and 345 GHz, weaker Faraday effects should result in an increased degree of polarization. Measurements of the EVPA pattern at higher frequencies can distinguish between Faraday rotation along the line of sight and the imprint of the underlying magnetic field structure. Continued imaging with the EHT and advances in radiative and nonthermal theoretical models will further constrain the electron distribution and magnetic field structure in the jet-launching region near the supermassive black hole event horizon in M87.
Delaney  AST-1903847, AST-1935980, AST-2034306); the Natural Science Foundation of China (grants 11573051, 11633006, 11650110427, 10625314, 11721303, 11725312, 11933007, 11991052, 11991053 This work used the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562, and CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442. XSEDE Stam-pede2 resource at TACC was allocated through TG-AST170024 and TG-AST080026N. XSEDE JetStream resource at PTI and TACC was allocated through AST170028. The simulations were performed in part on the SuperMUC cluster at the LRZ in Garching, on the LOEWE cluster in CSC in Frankfurt, and on the HazelHen cluster at the HLRS in Stuttgart. This research was enabled in part by support provided by Compute Ontario (http://computeontario. ca), Calcul Quebec (http://www.calculquebec.ca) and Compute Canada (http://www.computecanada.ca). We thank the staff at the participating observatories, correlation centers, and institutions for their enthusiastic support.
This Letter makes use of the following ALMA data: ADS/JAO. ALMA#2016.1.01154.V. ALMA is a partnership of the European Southern Observatory (ESO; Europe, representing its member states), NSF, and National Institutes of Natural Sciences of Japan, together with National Research Council (Canada), Ministry of Science and Technology (MOST; Taiwan), Academia Sinica Institute of Astronomy and Astrophysics (ASIAA; Taiwan), and Korea Astronomy and Space Science Institute (KASI; Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, Associated Universities, Inc. (AUI)/NRAO, and the National Astronomical Observatory of Japan (NAOJ). The NRAO is a facility of the NSF operated under cooperative agreement by AUI. This paper has made use of the following APEX data: Project ID T-091.F-0006-2013. APEX is a collaboration between the Max-Planck-Institut für Radioastronomie (Germany), ESO, and the Onsala Space Observatory (Sweden). The SMA is a joint project between the SAO and ASIAA and is funded by the Smithsonian Institution and the Academia Sinica. The JCMT is operated by the East Asian Observatory on behalf of the NAOJ, ASIAA, and KASI, as well as the Ministry of Finance of China, Chinese Academy of Sciences, and the National Key R&D Program (No. 2017YFA0402700) of China. Additional funding support for the JCMT is provided by the Science and Technologies Facility Council (UK) and participating universities in the UK and Canada. The LMT is a project operated by the Instituto Nacional de Astrófisica, Óptica, y Electrónica (Mexico) and the University of Massachusetts at Amherst (USA), with financial support from the Consejo Nacional de Ciencia y Tecnología and the National Science Foundation. The IRAM 30-m telescope on Pico Veleta, Spain is operated by IRAM and supported by CNRS (Centre National de la Recherche Scientifique, France), MPG (Max-Planck-Gesellschaft, Germany) and IGN (Instituto Geográfico Nacional, Spain). The SMT is operated by the Arizona Radio Observatory, a part of the Steward Observatory of the University of Arizona, with financial support of operations from the State of Arizona and financial support for instrumentation development from the NSF. The SPT is supported by the National Science Foundation through grant PLR-1248097. Partial support is also provided by the NSF Physics Frontier Center grant PHY-1125897 to the Kavli Institute of Cosmological Physics at the University of Chicago, the Kavli Foundation and the Gordon and Betty Moore Foundation grant GBMF 947. The SPT hydrogen maser was provided on loan from the GLT, courtesy of ASIAA. The EHTC has received generous donations of FPGA chips from Xilinx Inc., under the Xilinx University Program. The EHTC has benefited from technology shared under open-source license by the Collaboration for Astronomy Signal Processing and Electronics Research (CASPER). The EHT project is grateful to T4Science and Microsemi for their assistance with Hydrogen Masers. This research has made use of NASAʼs Astrophysics Data System. We gratefully acknowledge the support provided by the extended staff of the ALMA, both from the inception of the ALMA Phasing Project through the observational campaigns of 2017 and 2018. We would like to thank A. Deller and W. Brisken for EHT-specific support with the use of DiFX. We acknowledge the significance that Maunakea, where the SMA and JCMT EHT stations are located, has for the indigenous Hawaiian people.

Appendix A Relationship between β 2 Coefficient and Eand B-modes
The β 2 coefficients of the azimuthal decomposition of the complex linear polarization = +    i used in Section 5 are directly related to the decomposition of polarization fields into E-and B-modes familiar from cosmology. In this Appendix, we illustrate that relationship and demonstrate that the information in the image-space decomposition of GRMHD library images in Section 5 can also be accessed directly in calibrated visibility domain data sampled on EHT 2017 baselines, provided the data are accurately phase calibrated.

A.1. Definitions
In defining flat sky E-and B-modes, we follow the conventions of Kamionkowski & Kovetz (2016), Section 4.1 (up to factors of 2 ). E-and B-modes are naturally defined in the visibility space sampled by an interferometer. For a baseline vector u with a magnitude u and PA θ, theẼ andB visibilities are related to the Stokes visibilities and by a rotation of 2θ in the Fourier plane cos 2 sin 2 sin 2 cos 2 , , .
A 1 In real space, the E-and B-mode fields are analogous to the gradient and curl of the polarization tensor: where ò ac is the 2D Levi-Cevita symbol and the polarimetric tensor is A 3 ab P transforms as a trace-free tensor under rotations; for a rotation matrix R(α) that rotates the coordinate axes by an angle α, P → R(α)PR T (α) (equivalently, the complex field  a  e i 2 ). While the values of the  and  images depend on the choices of coordinate axes, the real space E-and B-mode images are coordinate-independent scalars.

A.2. Relationship between (E, B) and β 2 Coefficients
Consider a linearly polarized image = +    i in 2D image-domain polar coordinates (ρ, f). We can expand the image in a multipole series: A 4 m m m im 0 In the decomposition of Equation (A4), I 0 is the total flux density of the Stokes  image, the β m coefficients are complex, and the radial envelope function f m (ρ) is normalized so that The β m coefficients defined in this way then correspond to those defined in Equation (9): In particular, β 0 is the image-integrated complex fractional polarization, and β 2 encodes the same information on the gradient and curl of the polarization field that is available in the E-and B-modes. Note that because  is a complex image, in  Figure 20 illustrates the values of the β 2 and the corresponding signs of theẼ andB visibilities for these four azimuthally symmetric cases.

A.3.Ẽ andB Distributions of GRMHD Library Images
Starting with Equation (A10), and integrating over the baseline angle θ, only the β 2 mode survives.
Thus, if we have calibrated visibility measurements ofQ and U (and thusẼ andB), we can measure the β 2 mode by averagingẼ andB in visibility space.
To illustrate this connection between E-and B-mode visibilities and the β 2 coefficient, we compute the following averages on images in the GRMHD library  where we take the average only over (u, v) points sampled by the EHT in 2017, including conjugate baselines. As atmospheric phase errors and D-term miscalibration will affect the phase of theẼ andB visibility and thus the results for the averages, we perform this test on synthetic data from the image library using perfectly calibrated data with no noise.
Because we include conjugates in the average over all data points, the average of the imaginary part is zero. We perform the averaging only over a u, v range [1, 10]Gλ in order to remove any effects from large-scale structure, which may have a different net sense of polarization than the resolved emission ring.
We finally combine  R and  R in a complex quantity: where the negative signs are chosen to match the angle convention for β 2 in Equation (9). In Figure 21 we show histograms of the magnitude and angle of the b  for comparison with the β 2 histograms in Figure 9. The distributions broken down by model type (MAD/SANE, prograde/retrograde) reproduce the general behavior of β 2 amplitude and phase from the image-domain calculations in Figure 9, although the normalization of the b  amplitude is different. Comparing Figure 9 with Figure 21, it is apparent that all of the essential information on the EVPA structure used in this Letter can, in principle, be extracted from EHT visibilities without image reconstruction. However, because phase and amplitude calibration of the EHT visibilities is necessary for extracting the E-and B-modes from the visibilities, modeling the source structure in the image domain would remain a necessary part of the analysis even if we were to use  R and  R instead of β 2 .
Appendix B Faraday Rotation in GRMHD Models of M87 * As linear polarization travels through magnetized plasma, Faraday rotation shifts its EVPA by t r 2 V radians. If  t r 1 V , as is the case for most of our models (see Figure 7), Faraday rotation can, in principle, scramble otherwise observable polarimetric signals. In this section, we explore in more detail the sources of Faraday rotation in our models, and demonstrate that observable linear polarization signals can, in fact, exist in models with  t r 1 V . This is because Faraday rotation occurs co-spatially with the emission and should not be conceived of as a purely external screen. Ricarte et al. (2020) studied the resolved Faraday rotation properties of a subset of the same models used in this work. Figure 14 shows their inferred |RM| versus |m| net for those images. For each model, 11 snapshots spaced between 7500 and 10000 r g /c are included. Each of these models was found to pass the constraints of EHTC V. Snapshots with positive RMs are plotted with filled symbols, while those with negative RMs are plotted with open symbols. In gray, we overplot the allowed range of |m| net as well as the range of RM for the core region inferred from simultaneous ALMA-only observations.
Despite the large Faraday depths of these models, many of them are capable of producing RMs that are consistent with the observed data. RM and |m| net are anti-correlated, as expected, because a greater amount of Faraday rotation should both increase the RM and cause a greater amount of scrambling of the polarized emission. Note that the RM varies by orders of magnitude and even flips sign over time in these models. This is due to the summation of time-variable regions with significantly different and even oppositely signed Faraday rotation depths, which also contributes to highly non-λ 2 evolution of the EVPA with wavelength.
One potential source of uncertainty is the limited volume of our GRMHD simulations. In our image library, the radiative transfer equation is solved within a radius of only 50 or 100 r g depending on the model, while in principle significant Faraday rotation may occur at much larger radius. Figure 22 demonstrates that for the gas distributions studied, most of the Faraday rotation occurs at radii much smaller than the outer domain. Five example models are visualized here: (a) MAD a * = +0.94 R high = 20, (b) SANE a * = +0.5 R high = 1, (c) MAD a * = −0.5 R high = 160, (d) SANE a * = +0.5 R high = 160, and (e) SANE a * = 0 R high = 80. The brightness of each pixel scales with the total intensity (intentionally saturating 0.3% of the pixels), while the color of each pixel scales with the intensity-weighted Faraday depth t r I , V , as shown in the colorbar. t r I , V is distinct from t r V because it is intensity weighted along each ray, such that in each pixel As also shown in Figure 7, these models span a wide range of Faraday depths. Typically, SANE models and models with a larger R high have larger Faraday depths than MAD models and those with smaller R high . SANE models require a larger accretion rate to reproduce the total intensity of M87 * , which increases the amount of Faraday rotating material. Meanwhile, increasing R high lowers the temperature of the midplane by construction. This makes Faraday rotation more efficient, and also requires a larger accretion rate to compensate for the lower electron temperatures (see Mościbrodzka et al. 2017, for an extended discussion).
In Figure 23, we confirm that the linear polarization parameters used in this study are not strongly evolving at the outer simulation domain. Here, we provide polarization maps and the linear polarization parameters for these models blurred with a Gaussian beam with a FWHM of 20 μas. The rows display models with outer integration radii of 10, 20, 40, and 50 r g . For the three leftmost models, there is little difference between images constructed with = r r 10 g max and those with = r r 50 g max , echoing our previous findings that there are small fractional differences in the Faraday depth between these scales. Faraday rotation thick models (d) and (e) show modest differences between images calculated with outer boundaries of 10 r g and 50 r g . Those images also appear to be converged by 40 r g . Some of our models produce observable polarimetric signatures despite t r V being large enough to potentially depolarize all of the emission. This apparent contradiction is resolved by the fact that not all emission is Faraday rotated by the same amount. Because Faraday rotation occurs co-spatially with the emission, instead of as an external screen, there can exist emission traveling on Faraday-thin paths to the camera even in models where  t r 1 V when integrated along the entire geodesic. In Figure 24, we illustrate this phenomenon by splitting these images into the emission originating in front of or behind the midplane. Here, models are plotted as in Figure  Emission from in front of the midplane that is observed within the photon ring has a much larger t r I , V than the rest of the image because those geodesics pass through the midplane and around the black hole through Faraday-thick material. Figure 22. Intensity-weighted Faraday depths visualized with five example models: (a) MAD a * = +0.94 R high = 20, (b) SANE a * = +0.5 R high = 1, (c) MAD a * = −0.5 R high = 160, (d) SANE a * = +0.5 R high = 160, and (e) SANE a * = 0 R high = 80. The brightness of each pixel scales with the total intensity (intentionally saturating 0.3% of the pixels), while the color indicates the intensity-weighted Faraday depth, t r I , V . In the top row, the maximum integration radius is set to 20 r g , while in the bottom row, the maximum integration radius is set to 50 r g . We find very little difference, confirming that our results should be insensitive to the outer radius of the simulation domain. Prograde SANE models show rapidly decreasing |m| net with increasing R high ( Figure 25) and are significantly depolarized when R high > 10. This behavior was previously demonstrated by Mościbrodzka et al. (2017). The accretion flow electron temperature decreases with increasing R high , increasing the strength of Faraday rotation while also concentrating the emission at high latitudes behind the black hole (see also EHTC V). The emission is then depolarized when traveling through the Faraday-thick midplane plasma.
Retrograde SANE models, however, show nearly the opposite behavior, with depolarization maximized for R high = 1. At larger values of R high , linearly polarized emission appears on the near side of the midplane, producing coherent linear polarization structure that is not Faraday depolarized.
MAD models at all spins show a mild degree of depolarization with increasing R high . The accretion flow electron temperature remains high even for large values of R high , as much of the plasma has β ; 1.
Similar qualitative behavior is seen in 〈|m|〉 ( Figure 26) and the amplitude of β 2 (Figure 27). However, those quantities show less time variability (narrower distributions) than is seen in |m| net . As a result, observed ranges of those values are more constraining. In particular, the MAD models show consistent offsets where 〈|m|〉 and |β 2 | are lower for R low = 10 than R low = 1 models. Some spin dependence is also apparent, with high-prograde spin usually corresponding to the highest degrees of ordered polarization.
When the β 2 amplitude is not strongly suppressed (e.g., by Faraday rotation), the β 2 phase distributions are related to intrinsic magnetic field structure (e.g., Figure 3 and PWP). Prograde spin, R high = 1 SANE models and retrograde spin, large R high SANE models both show radial EVPA patterns, resulting in β 2 phase distributions near zero. MAD models show spin-dependent β 2 phase distributions for low values of R high , ranging from spiral patterns (∠β 2 ; −90 deg) for retrograde spin to more radial patterns at high-prograde spin. The patterns are relatively constant functions of R high and R low , although with some shift of MAD prograde distributions to twistier EVPA patterns, particularly for R low = 10.
Most models show distributions of v net centered on zero, near the observed range ( Figure 29). MAD models generally show low circular polarization fractions, while heavily depolarized SANE models (retrograde low R high , prograde large R high ) tend to show larger |v net | than observed, which can be explained by stronger Faraday conversion in the emission region.    Figure 25, but for the image-averaged polarization fraction 〈|m|〉.

Appendix D Detailed Model Scoring Results
In Table 3 we provide a summary of the number of images for each model that fall within the observed range of each individual theory metric (used in the joint scoring procedure) and within the observed ranges of all metrics simultaneously (used in the simultaneous scoring). Boldfaced type is used for models that are deemed viable by one of the scoring systems.
For simultaneous scoring, a viable model contains at least one image that simultaneously satisfies all constraints. For joint scoring, a viable model has a joint likelihood >1% that of the maximum found across all models. We also provide a summary score-"pass" indicates a model that satisfies the polarimetric constraints according to either scoring procedure, as well as the jet-power cut of P jet > 10 42 erg s −1 (EHTC V).   Note. The number of images passing each polarimetric constraint are given along with the number N sim simultaneously passing all of them. The accretion rate  -M 3 is in units of 10 −3 M e yr −1 , and P jet,42 is the jet power in units of 10 42 erg s −1 . Models that pass according to either the simultaneous or joint scoring method (boldface) and have P jet,42 1 are given a summary score of pass.