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

Throughout this Letter we use the following definitions and conventions for polarimetric quantities, following the International Astronomical Union (IAU) definitions of the Stokes parameters $({ \mathcal I },{ \mathcal Q },{ \mathcal U },{ \mathcal V })$ (Hamaker & Bregman 1996; Smirnov 2011). The complex linear polarization field ${ \mathcal P }$ is

$$\begin{eqnarray}&&{ \mathcal P }={ \mathcal Q }+i{ \mathcal U }.\end{eqnarray} \tag{ 1 }$$

Then, the electric-vector position angle (EVPA) is defined as

$$\begin{eqnarray}&&\mathrm{EVPA}\equiv \displaystyle \frac{1}{2}\arg ({ \mathcal P }).\end{eqnarray} \tag{ 2 }$$

The EVPA is measured east of north on the sky. Therefore, positive ${ \mathcal Q }$ is aligned with the north–south direction and negative ${ \mathcal Q }$ with the east–west direction. Positive ${ \mathcal U }$ is at a + 45 deg angle with respect to the positive ${ \mathcal Q }$ axis (in the northeast–southwest direction). Positive Stokes ${ \mathcal V }$ 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 ${ \mathcal V }$ 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 ${ \mathcal V }$ corresponds to θ_B ∈ [0.5π, π].

The linear and circular polarization fractions at a point in the image are defined as

$$\begin{eqnarray}&&| m| \equiv \displaystyle \frac{| { \mathcal P }| }{{ \mathcal I }},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&| v| \equiv \displaystyle \frac{| { \mathcal V }| }{{ \mathcal I }}.\end{eqnarray} \tag{ 4 }$$

We also define the rotation measure between two wavelengths λ₁ and λ₂

$$\begin{eqnarray}&&\mathrm{RM}\equiv \displaystyle \frac{\mathrm{EVPA}({\lambda }_{1})-\mathrm{EVPA}({\lambda }_{2})}{{\lambda }_{1}^{2}-{\lambda }_{2}^{2}}.\end{eqnarray} \tag{ 5 }$$

Unresolved observations measure the net (image-integrated) polarization fractions

$$\begin{eqnarray}&&| m{| }_{\mathrm{net}}=\displaystyle \frac{\sqrt{{\left({\displaystyle \sum }_{i}{{ \mathcal Q }}_{i}\right)}^{2}+{\left({\displaystyle \sum }_{i}{{ \mathcal U }}_{i}\right)}^{2}}}{{\displaystyle \sum }_{i}{{ \mathcal I }}_{i}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{v}_{\mathrm{net}}=\displaystyle \frac{{\displaystyle \sum }_{i}{{ \mathcal V }}_{i}}{{\displaystyle \sum }_{i}{{ \mathcal I }}_{i}},\end{eqnarray} \tag{ 7 }$$

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

$$\begin{eqnarray}&&\langle | m| \rangle =\displaystyle \frac{{\displaystyle \sum }_{i}\sqrt{{{ \mathcal Q }}_{i}^{2}+{{ \mathcal U }}_{i}^{2}}}{{\displaystyle \sum }_{i}{{ \mathcal I }}_{i}}.\end{eqnarray} \tag{ 8 }$$

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.

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 ∝ λ², 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⁵ rad m⁻².

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.

The resolved polarimetric images of the M87 core reported in EHTC VII display robust features between different image-reconstruction 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 image-integrated 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 ${ \mathcal P }={ \mathcal Q }+i{ \mathcal U }$ 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 (ρ, φ), the β_m coefficients are

$$\begin{eqnarray}&&{\beta }_{m}=\displaystyle \frac{1}{{I}_{\mathrm{ann}}}{\int }_{{\rho }_{\min }}^{{\rho }_{\max }}{\int }_{0}^{2\pi }{ \mathcal P }(\rho ,\varphi )\,{e}^{-{im}\varphi }\ \rho d\varphi d\rho ,\end{eqnarray} \tag{ 9 }$$

where I_ann is the Stokes ${ \mathcal I }$ flux density contained inside the annulus set by the limiting radii ${\rho }_{\min }$ and ${\rho }_{\max }$. We take ${\rho }_{\min }=0$ and ${\rho }_{\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, β₂, was the most discriminating in identifying the underlying magnetized accretion model. The phase of β₂ maps well onto the qualitative behavior expected of polarization maps with idealized magnetic field configurations. In our convention, radial EVPA patterns have positive real β₂ (∠β₂ = 0 deg), azimuthal EVPA patterns have negative real β₂ (∠β₂ = 180 deg), and left- (right-) handed spiral patterns have positive (negative) pure imaginary β₂ (∠β₂ = 90 deg and −90 deg, respectively). These idealized EVPA pattern configurations and their corresponding β₂ coefficients are summarized in Appendix A and Figure 1 of PWP. The measured range of the complex β₂ coefficient across the different image-reconstruction methods and observing days reported in EHTC VII, their Table 2, is 0.04 ≤ ∣β₂∣ ≤ 0.07 for the amplitude and $-163\,\deg \leqslant \arg [{\beta }_{2}]\leqslant -129\,\deg$ for the phase.

Appendix A demonstrates that the information in the β₂ 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 β₂ 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 β₂ amplitude and phase—are summarized in Table 2.

Faraday rotation in a uniform plasma with rotation measure (RM) rotates the EVPA away from its intrinsic value EVPA₀ according to Equation (5). The change in EVPA from its intrinsic value at 230 GHz (λ ≃ 1.3 mm) is

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{EVPA}\simeq 9.7\left(\displaystyle \frac{\mathrm{RM}}{{10}^{5}\,\mathrm{rad}\,{{\rm{m}}}^{-2}}\right)\,\deg .\end{eqnarray} \tag{ 10 }$$

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. 2006, 2007; Kuo 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 ∝ $1/{T}_{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² = 1, usually r ∼ 100 r_g (where r_g = GM/c² 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⁵ rad m⁻². These values are consistent with, but much more constraining than, the range determined from past SMA observations (−3.4–7 × 10⁵ rad m⁻², 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 two-component 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⁶ rad m⁻². 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.

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 ${ \mathcal I }$ image's brightness, size, and total flux density:

$$\begin{eqnarray}&&{n}_{e}\simeq 2.9\times {10}^{4}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&B\simeq 4.9\,{\rm{G}}.\end{eqnarray} \tag{ 12 }$$

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¹⁰ K, based on the observed brightness temperature of the EHT image. This temperature corresponds to Θ_e =kT_e /m_e c² = 10.5, so the emitting electrons have moderately relativistic mean Lorentz factors $\bar{\gamma }\approx 3{{\rm{\Theta }}}_{e}\approx 30$. 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 ${\tau }_{{\rho }_{V}}$:

$$\begin{eqnarray}&&{\tau }_{{\rho }_{V}}\approx r\times {\rho }_{V}\simeq 5.2\,\left(\displaystyle \frac{r}{5{r}_{g}}\right),\end{eqnarray} \tag{ 13 }$$

where ρ_V is the Faraday rotation coefficient (e.g., Jones & Hardee 1979). For emission entirely behind an external Faraday screen, ${\tau }_{{\rho }_{V}}$ is related to the rotation measure RM via ${\tau }_{{\rho }_{V}}=2\mathrm{RM}{\lambda }^{2}$, 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 ∝ λ².

Our estimated ${\tau }_{{\rho }_{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 (${\tau }_{{\rho }_{Q}}\simeq 0.5$), while self-absorption is weak (τ_I ≃ 0.05). Requiring an internal Faraday optical depth ${\tau }_{{\rho }_{V}}\gt 2\pi$ (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 ${\tau }_{{\rho }_{V}}\gt 2\pi$ 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², constrained by the requirement that the Faraday optical depth ${\tau }_{{\rho }_{V}}\gt 2\pi$. 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¹⁰ < T_e < 1.2 × 10¹¹ 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⁴ cm⁻³ and 10⁷ cm⁻³.

Zoom In Zoom Out Reset image size

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¹⁰ 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 ${ \mathcal V }$ emissivity to the Stokes ${ \mathcal I }$ emissivity j_V /j_I ≃ 1%. For B ≃ 10³ G, j_V /j_I ≃ 10%, for all temperatures >10¹⁰ K. We also note that magnetic fields of B ≳ 5 G are sufficient to produce jet powers of P_jet ≳ 10⁴² erg s⁻¹ (e.g., EHTC V) via the Blandford & Znajek (1977) process.

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 field— as seen by a distant observer. In Figure 3 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.

Zoom In Zoom Out Reset image size

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. 2018, 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 direction $\bar{P}$ of synchrotron radiation in the emitted frame is perpendicular to the wavevector $\hat{k}$ and the magnetic field vector $\bar{B}$. 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 that $\bar{P}$ should be vertical (north−south) everywhere since a vertical $\bar{B}$ is tilted east−west for this viewing geometry. We might also expect that $\bar{P}\simeq 0$ when the black hole is viewed face on, because $\hat{k}| | \bar{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 $\hat{k}^{\prime}$ to the emitted wavevector $\hat{k}$ 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 ${\hat{k}}_{\mathrm{emit}}=\hat{k}+\hat{k}^{\prime}$ and the magnetic field $\bar{B}$ 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.

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 ¹²⁹ 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². 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

$$\begin{eqnarray}&&\phi \equiv {{\rm{\Phi }}}_{{\rm{BH}}}{\left(\dot{M}{r}_{g}^{2}c\right)}^{-1/2},\end{eqnarray} \tag{ 14 }$$

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 $\dot{M}$ is the mass accretion rate through the event horizon. The flux saturates at values of ϕ ≳ 50, and the flow becomes MAD. The SANE simulations that we consider have lower values of ϕ ≈ 5. ¹³⁰ 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 ϕ ≲ 1. We also do not consider values in the relatively narrow intermediate range 5 ≲ ϕ ≲ 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³ 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⁴ r_g /c and their accretion flows reach steady state within r = 10–20 r_g .

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⁹ M_⊙ 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 ${ \mathcal Q }$ and ${ \mathcal U }$ components appropriately. Each image has a 320 × 320 pixel resolution over a 160 μas field of view, where each pixel contains full Stokes ${ \mathcal I },{ \mathcal Q },{ \mathcal U },{ \mathcal V }$ 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²/8π. The temperature ratio is then taken to be

$$\begin{eqnarray}&&R=\displaystyle \frac{{T}_{i}}{{T}_{e}}={R}_{\mathrm{high}}\ \displaystyle \frac{{\beta }^{2}}{1+{\beta }^{2}}+{R}_{\mathrm{low}}\ \displaystyle \frac{1}{1+{\beta }^{2}},\end{eqnarray} \tag{ 15 }$$

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

$$\begin{eqnarray}&&{T}_{e}=\displaystyle \frac{2\,{m}_{p}{u}_{\mathrm{gas}}}{3\rho k(2+R)}.\end{eqnarray} \tag{ 16 }$$

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). EHTC V only considered models with R_low = 1. Larger values of R_low 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²/4π ρ c² ≫ 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) and B. 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*.

In Figures 4 and 5 we show images and polarization maps for a subset of library models. In general, because the horizon-scale 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 β₂ phase across the full image library that we consider later in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image 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 β₂ by ≃80 deg would correspond to an apparent RM of ≃−4 × 10⁵ rad m⁻². 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.

We compute the polarimetric observables (∣m∣_net, ∣v∣_net, 〈∣m∣〉, β₂) 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 β₂ depend on the resolution and hence the size of the Gaussian blurring kernel. The value of β₂ 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 β₂ power as δ β₂/∣β₂∣ ≈ 4u². Effects on the β₂ phase enter at similar order. In the case of the EHT image, u is likely less than one-fifth, meaning that centering errors in β₂ 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, $\langle {\tau }_{{\rho }_{V}}\rangle$. At small $\langle {\tau }_{{\rho }_{V}}\rangle$, 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 $\langle {\tau }_{{\rho }_{V}}\rangle$, strong scrambling from internal Faraday rotation typically results in small predicted polarization fractions of <5% at the scale of the EHT beam.

Zoom In Zoom Out Reset image size

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 $\langle {\tau }_{{\rho }_{V}} \rangle \gtrsim {10}^{3}$. In these models, most of the observed polarized flux originates in the forward 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 ${\tau }_{{\rho }_{V}}$.

Distributions of all observables are shown in Figure 7 (〈∣m∣〉, left panel), Figure 8 (∣m∣_net and ∣v∣_net), and Figure 9 (∣β₂∣ and ∠β₂). 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To evaluate whether a given GRMHD model is consistent with the EHT observations reported in EHTC VII, we require 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%.
• 2.  
  The image-integrated net circular polarization ∣v∣_net satisfies an upper limit from ALMA-only measurements reported in Goddi et al. (2021): ∣v∣_net ≤ 0.8%.
• 3.  
  The image-averaged linear polarization 〈∣m∣〉 is in the measured range from the EHT image reconstructions at 20 μas scale resolution: 5.7% ≤ ∣m∣_net ≤ 10.7%.
• 4.  
  The amplitude and phase of the complex β₂ coefficient quantifying coherent azimuthal structure are in the measured range: 0.04 ≤ ∣β₂∣ ≤ 0.07 and $-163\,\deg \,\leqslant \arg [{\beta }_{2}]\leqslant -129\,\,\deg$.

We use 72,000 library images (from Section 4) with 200 time snapshots per model at three inclination angles, six values of R_high = 1, 10, 20, 40, 80, 160, two values of R_low = 1, 10, five values of a_* = −0.9375, −0.5, 0, +0.5, +0.9375, and realized with both MAD and SANE magnetic field configurations.

In comparing models to observables, the β₂ metric is the most constraining. Only 790 snapshot images out of the 72,000 considered fall in the range of those reconstructed in both β₂ 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).

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 ∣β₂∣ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 χ² values for the five metrics ${x}_{j}\in \left\{| m{| }_{\mathrm{net}},| v{| }_{\mathrm{net}}\right.$, $\left.\langle | m| \rangle ,| {\beta }_{2}| ,\angle {\beta }_{2}\right\}$ for all snapshots k from a given model as

$$\begin{eqnarray}&&{\chi }_{j,k}^{2}=\displaystyle \frac{{\left({x}_{j,k}-{\bar{x}}_{j}\right)}^{2}}{{\sigma }_{j}^{2}},\end{eqnarray} \tag{ 17 }$$

where x_j,k are the values of a scoring metric x_j for each of the 600 snapshots k from a given model, ${\bar{x}}_{j}$ is the mean of those values for the model, and σ_j is taken as one half of the observed data range from Table 2. Note that the scoring results of this method do not depend on the choice of σ_j . We then calculate an analogous ${\chi }_{j,\mathrm{data}}^{2}$ value for the midpoint of the measured range from Table 2. A likelihood value ${{ \mathcal L }}_{j}$ of the data being drawn from the model distribution is defined as the fraction of images with ${\chi }_{j,k}^{2}\gt {\chi }_{j,\mathrm{data}}^{2}$. The joint likelihood of each model is the product ${ \mathcal L }={{\rm{\Pi }}}_{j}{{ \mathcal L }}_{j}$ of those for the five metrics x_j .

To produce a non-zero likelihood ${ \mathcal L }$ 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 χ² 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.

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 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 ∣β₂∣).

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.

EHTC V presented constraints on the GRMHD simulation models based on fits to the EHT total intensity data, model 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⁴² erg s⁻¹). 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.

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⁻³, 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 ${\tau }_{{\rho }_{V}}$, over snapshots that simultaneously satisfy all our observational constraints. The median values are τ_I ≃ 0.1 and ${\tau }_{{\rho }_{V}}\simeq 50$. All of our viable images have ${\tau }_{{\rho }_{V}}\gt 2\pi$, 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 ${ \mathcal I }$ 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⁴² erg s⁻¹ (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. ¹³¹ 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 $\dot{M}\simeq (3{\unicode{x02013}}20)\,\times {10}^{-4}{M}_{\odot }\,{{\rm{yr}}}^{-1}$. This corresponds to $\dot{m}=\dot{M}/{\dot{M}}_{\mathrm{Edd}}\,\simeq (2\mbox{--}15)\,\times {10}^{-6}$, where ${\dot{M}}_{\mathrm{Edd}}$ is the Eddington accretion rate. ¹³² The measured radiative efficiency $\epsilon =L/\dot{M}{c}^{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⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mass accretion rate found here is much lower than the Bondi rate calculated from Chandra observations (Di 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⁴³ erg s⁻¹ 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⁻¹ (Stawarz et al. 2006), and Low-Frequency Array (LOFAR) observations suggest that a jet power of ∼10⁴⁴ erg s⁻¹ was necessary within the last ∼million years to inflate the observed radio lobes on scales of ∼80 kpc (de Gasperin et al. 2012).

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 $\dot{M}\sim {10}^{-5}{\dot{M}}_{{\rm{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⁵ rad m⁻². 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⁵ r_g .

For our favored plasma parameters for M87*, we expect substantial Faraday RM internal to the emission region itself, ${\tau }_{{\rho }_{V}}\gtrsim 2\pi$, 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⁶ rad m⁻². 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⁶ rad m⁻²). 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 low-inclination 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.

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, ${\gamma }_{\max }$. The minimum Lorentz factor, ${\gamma }_{\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 ${\gamma }_{\max }$ are constant across the accretion flow. We assume that ${\gamma }_{\max }\to \infty$ and we explore values of p = 3.5, 4.5 and η = 0.01, 0.1.

In Figure 15 we present linear polarization maps from two MAD models and one SANE model comparing purely thermal and hybrid electron distributions. Using a hybrid distribution function does not affect the structure of the EVPA map (β₂ phase), but it changes the image-integrated and resolved linear polarization fractions. For example, in the MAD a_* = −0.5 (MAD a_* = 0.5) model, with the selected hybrid parameters, the ∣m∣_net, v_net and 〈∣m∣〉 ranges are 4.3%–4.6% (2.5%–3.8%), 0.25%–0.37% ((−0.5) to (−0.12)%), and 10.6%–11.5% (12%–14%), respectively. Slightly larger deviations from the thermal model are measured in the SANE a_* = −0.94 scenario, where the ∣m∣_net, v_net and 〈∣m∣〉 ranges are 2.2%–4.1%, −0.004% to 0.31%, and 14%–20%, respectively.

Zoom In Zoom Out Reset image size

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

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 radiation from a pure power-law distribution of electrons with ${\gamma }_{\min }=100$ 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 β₂ phase).

Zoom In Zoom Out Reset image size

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.

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. 2015, 2017; Chael et al. 2018, 2019; Ryan 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.

Zoom In Zoom Out Reset image size

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

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 β₂ 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.

Zoom In Zoom Out Reset image size

Both quantities show variations on timescales from days to months. The phase and amplitude of β₂ should change over the course of a week of observations. In EHTC VII, we observe changes in the the β₂ 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 β₂ on ∼week timescales can be much more dramatic than we observe in 2017, with variations in β₂ 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 β₂ 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 β₂ phase should on average be closer to zero than we observe in 2017.

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

Zoom In Zoom Out Reset image size

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⁵ rad m⁻² 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 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⁵ rad m⁻². 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.

In defining flat sky E- and B-modes, we follow the conventions of Kamionkowski & Kovetz (2016), Section 4.1 (up to factors of $\sqrt{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 $\tilde{E}$ and $\tilde{B}$ visibilities are related to the Stokes visibilities $\tilde{{ \mathcal Q }}$ and $\tilde{{ \mathcal U }}$ by a rotation of 2θ in the Fourier plane

$$\begin{eqnarray}&&\left(\begin{array}{c}\tilde{E}(u,\theta )\\ \tilde{B}(u,\theta )\end{array}\right)=\left[\begin{array}{cc}\cos 2\theta & \sin 2\theta \\ -\sin 2\theta & \cos 2\theta \end{array}\right]\left(\begin{array}{c}\tilde{{ \mathcal Q }}(u,\theta )\\ \tilde{{ \mathcal U }}(u,\theta )\end{array}\right).\end{eqnarray} \tag{ A1 }$$

In real space, the E- and B-mode fields are analogous to the gradient and curl of the polarization tensor:

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}E={{\rm{\partial }}}_{a}{{\rm{\partial }}}_{b}{{\boldsymbol{P}}}_{ab},\,\,\,{{\rm{\nabla }}}^{2}B={\epsilon }_{ac}{{\rm{\partial }}}_{b}{{\rm{\partial }}}_{c}{{\boldsymbol{P}}}_{ab},\end{eqnarray} \tag{ A2 }$$

where _ac is the 2D Levi-Cevita symbol and the polarimetric tensor is

$$\begin{eqnarray}{{\boldsymbol{P}}}_{{ab}}=\left[\begin{array}{cc}{ \mathcal Q } & { \mathcal U }\\ { \mathcal U } & -{ \mathcal Q }\end{array}\right].\end{eqnarray} \tag{ A3 }$$

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 ${ \mathcal P }\to { \mathcal P }{e}^{2i\alpha }$). While the values of the ${ \mathcal Q }$ and ${ \mathcal U }$ images depend on the choices of coordinate axes, the real space E- and B-mode images are coordinate-independent scalars.

Consider a linearly polarized image ${ \mathcal P }={ \mathcal Q }+i{ \mathcal U }$ in 2D image-domain polar coordinates (ρ, ϕ). We can expand the image in a multipole series:

$$\begin{eqnarray}&&{ \mathcal P }(\rho ,\phi )={I}_{0}\displaystyle \sum _{m=-\infty }^{\infty }{\beta }_{m}{f}_{m}(\rho ){e}^{{im}\phi }.\end{eqnarray} \tag{ A4 }$$

In the decomposition of Equation (A4), I₀ is the total flux density of the Stokes ${ \mathcal I }$ image, the β_m coefficients are complex, and the radial envelope function f_m (ρ) is normalized so that

$$\begin{eqnarray}&&2\pi {\int }_{{\rho }_{\min }}^{{\rho }_{\max }}{f}_{m}(\rho )\rho \,d\rho =1.\end{eqnarray} \tag{ A5 }$$

The β_m coefficients defined in this way then correspond to those defined in Equation (9):

$$\begin{eqnarray}&&{\beta }_{m}=\displaystyle \frac{1}{{I}_{0}}{\displaystyle \int }_{{\rho }_{\min }}^{{\rho }_{\max }}{\displaystyle \int }_{0}^{2\pi }{ \mathcal P }(\rho ,\phi ){e}^{-{im}\phi }\rho \,d\rho \,d\phi .\end{eqnarray} \tag{ A6 }$$

In particular, β₀ is the image-integrated complex fractional polarization, and β₂ 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 ${ \mathcal P }$ is a complex image, in general ${\beta }_{m}\ne {\beta }_{-m}^{* }$.

The Fourier transform of P(ρ, ϕ) is

$$\begin{eqnarray}&&\tilde{P}(u,\theta )=2\pi {I}_{0}\displaystyle \sum _{m=-\infty }^{\infty }{i}^{-m}{\beta }_{m}{e}^{{im}\theta }{F}_{m}(u),\end{eqnarray} \tag{ A7 }$$

where F_m (u) is the mth order Hankel transform of the radial function f_m (ρ):

$$\begin{eqnarray}&&{F}_{m}(u)={\int }_{{\rho }_{\min }}^{{\rho }_{\max }}{f}_{m}(\rho ){J}_{m}(2\pi \rho u)\rho \,d\rho .\end{eqnarray} \tag{ A8 }$$

The Fourier transforms of ${ \mathcal Q }$ and ${ \mathcal U }$ (the linear Stokes visibilities) are then

$$\begin{eqnarray}\begin{array}{lll}\tilde{{ \mathcal Q }}(u,\theta ) & = & 2\pi {I}_{0}\displaystyle \sum _{m=-\infty }^{\infty }{i}^{-m}\mathrm{Re}\left[{\beta }_{m}{e}^{{im}\theta }{F}_{m}(u)\right],\\ \tilde{{ \mathcal U }}(u,\theta ) & = & 2\pi {I}_{0}\displaystyle \sum _{m=-\infty }^{\infty }{i}^{-m}\mathrm{Im}\left[{\beta }_{m}{e}^{{im}\theta }{F}_{m}(u)\right],\end{array}\end{eqnarray} \tag{ A9 }$$

and from the definition in Equation (A1) the E- and B-mode visibilites are

$$\begin{eqnarray}\begin{array}{rcl}\tilde{E}(u,\theta ) & = & 2\pi {I}_{0}\displaystyle \sum _{m=-\infty }^{\infty }{i}^{-m}\,\mathrm{Re}\left[{\beta }_{m}{e}^{i(m-2)\theta }{F}_{m}(u)\right],\\ \tilde{B}(u,\theta ) & = & 2\pi {I}_{0}\displaystyle \sum _{m=-\infty }^{\infty }{i}^{-m}\,\mathrm{Im}\left[{\beta }_{m}{e}^{i(m-2)\theta }{F}_{m}(u)\right].\end{array}\end{eqnarray} \tag{ A10 }$$

From Equation (A10), we can see immediately that an image with real β₂ and all other β_m≠2 = 0 is a pure E-mode; for instance, if β₂ = 1 (radial polarization vectors), $\tilde{E}\lt 0$. If β₂ = −1 (toroidal polarization vectors), $\tilde{E}\gt 0$. Similarly, if an image has purely imaginary β₂ (and all other β_m≠2 = 0) it is a pure B-mode; if e.g., β₂ = i (right-handed helical polarization vectors), then $\tilde{B}\lt 0$, and if β₂ = −i (left-handed helix), $\tilde{B}\gt 0$. Figure 20 illustrates the values of the β₂ and the corresponding signs of the $\tilde{E}$ and $\tilde{B}$ visibilities for these four azimuthally symmetric cases.

Zoom In Zoom Out Reset image size

Starting with Equation (A10), and integrating over the baseline angle θ, only the β₂ mode survives.

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }\tilde{E}(u,\theta )d\theta =-2\pi {I}_{0}\,\mathrm{Re}\left[{\beta }_{2}{F}_{2}(u)\right]\end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }\tilde{B}(u,\theta )d\theta =-2\pi {I}_{0}\,\mathrm{Im}\left[{\beta }_{2}{F}_{2}(u)\right].\end{eqnarray} \tag{ A12 }$$

Thus, if we have calibrated visibility measurements of $\tilde{Q}$ and $\tilde{U}$ (and thus $\tilde{E}$ and $\tilde{B}$), we can measure the β₂ mode by averaging $\tilde{E}$ and $\tilde{B}$ in visibility space.

To illustrate this connection between E- and B-mode visibilities and the β₂ coefficient, we compute the following averages on images in the GRMHD library

$$\begin{eqnarray}&&{{ \mathcal E }}_{{\rm{R}}}=\displaystyle \frac{{\left\langle \mathrm{Re}[\tilde{E}]\right\rangle }_{(u,v)}}{{\left\langle | \tilde{I}| \right\rangle }_{(u,v)}},\quad {{ \mathcal B }}_{{\rm{R}}}=\displaystyle \frac{{\left\langle \mathrm{Re}[\tilde{B}]\right\rangle }_{(u,v)}}{{\left\langle | \tilde{I}| \right\rangle }_{(u,v)}},\end{eqnarray} \tag{ A13 }$$

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 $\tilde{E}$ and $\tilde{B}$ 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 ${{ \mathcal E }}_{{\rm{R}}}$ and ${{ \mathcal B }}_{{\rm{R}}}$ in a complex quantity:

$$\begin{eqnarray}&&{\beta }_{{ \mathcal E }{ \mathcal B }}=-{{ \mathcal E }}_{{\rm{R}}}-i{{ \mathcal B }}_{{\rm{R}}},\end{eqnarray} \tag{ A14 }$$

where the negative signs are chosen to match the angle convention for β₂ in Equation (9). In Figure 21 we show histograms of the magnitude and angle of the ${\beta }_{{ \mathcal E }{ \mathcal B }}$ for comparison with the β₂ histograms in Figure 9. The distributions broken down by model type (MAD/SANE, prograde/retrograde) reproduce the general behavior of β₂ amplitude and phase from the image-domain calculations in Figure 9, although the normalization of the ${\beta }_{{ \mathcal E }{ \mathcal 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 ${{ \mathcal E }}_{{\rm{R}}}$ and ${{ \mathcal B }}_{{\rm{R}}}$ instead of β₂.

Zoom In Zoom Out Reset image size