Decomposing the internal faraday rotation of black hole accretion flows Free

ABSTRACT

1 INTRODUCTION

Using a network of millimetre telescopes around the world, the Event Horizon Telescope (EHT) has recently produced the first images of a black hole (BH) accretion flow (Event Horizon Telescope Collaboration 2019a,b,c,d,e,f). These images resolve the ‘shadow’ of the supermassive BH M87*, corresponding to rays that begin on its event horizon, providing new constraints on the properties of the BH and its accretion disc.

In Event Horizon Telescope Collaboration (2019e), henceforth EHT5, a library of general relativistic magnetohydrodynamic (GRMHD) simulations was produced to compare to the EHT images, exploring three fundamental quantities. The first is the strength of the magnetic field: models that accumulate strong magnetic flux around the BH are able to counteract the ram pressure of in-falling gas with magnetic pressure, resulting in what is termed a magnetically arrested disc (MAD; Igumenshchev, Narayan & Abramowicz 2003; Narayan, Igumenshchev & Abramowicz 2003; Chael, Narayan & Johnson 2019). In contrast, the weaker and more turbulent magnetic fields in Standard And Normal Evolution (SANE) models have a weaker effect on the gas dynamics of the disc (Narayan et al. 2012; Sądowski et al. 2013; Ryan et al. 2018). The second quantity is the BH’s angular momentum, described by the dimensionless spin parameter |$a \equiv Jc/(GM_\bullet ^2) \in (-1,1)$|⁠, where negative values correspond to counter-rotating accretion discs. Mass and spin are the only properties intrinsic to an astrophysical BH, but BH spins are constrained much more loosely than their masses. Most of our understanding of supermassive BH spin evolution originates from theory (e.g. Bardeen & Wagoner 1969; Thorne 1974; Gammie, Shapiro & McKinney 2004; King, Pringle & Hofmann 2008; Volonteri et al. 2013). The third quantity explored in this work is R_high, one parametrization of the relative temperatures of electrons and ions in the plasma. Such a prescription is necessary because the mean free path near the event horizon is so large, the two populations depart from thermal equilibrium and divide into a two-temperature plasma (Shapiro, Lightman & Eardley 1976; Rees et al. 1982; Narayan & Yi 1995a; Sądowski et al. 2017; Ryan et al. 2018). By combining EHT imaging with other multiwavelength constraints such as the X-ray flux and jet power, EHT5 rule out all but 19/60 possible models that span these three variables. In particular, all a = 0 models are excluded.

EHT constraints on M87* thus far have only considered total intensity (Stokes I), when in fact polarized visibilities have been obtained (Stokes I, Q, U, and V), but have not been published. Hence, current constraints have only utilized one fourth of the measured information. Synchrotron emission, the emission mechanism at millimetre wavelengths, has a near-unity intrinsic polarization fraction (Le Roux 1961; Bromley, Melia & Liu 2001; Broderick & Blandford 2003). Polarimetric imaging is consequently predicted to tightly constrain accretion models, which differ substantially in their linear polarization fractions as well as their morphologies (Palumbo, Wong & Prather 2020). Although an image has not yet been constructed for Sgr A*, strong and variable linear polarization has been observed for this source. Previous very long baseline interferometry measurements have revealed a partially ordered magnetic field structure at its centre (Johnson et al. 2015).

As polarized radiation travels through a magnetized plasma, it is transformed by effects sensitive to the local density, temperature, and magnetic field. Faraday rotation turns the electric vector position angle (EVPA) of linearly polarized emission, and Faraday conversion exchanges linearly and circularly polarized radiation (Sazonov 1969; Rybicki & Lightman 1986; Melrose 1997; Shcherbakov 2008). Future analyses with the EHT will be able to further distinguish models in the time and frequency domains (e.g. Broderick & Loeb 2006a,b; Roelofs et al. 2017; Medeiros et al. 2018). It is therefore timely for us to better understand time and frequency-dependent effects that may help us constrain accretion models.

As seen from equation (2), the RM directly traces the electron temperature, number density, and magnetic field along the line of sight. Among the models consistent with the EHT data of M87*, all of these quantities can vary by orders of magnitude. Note that while RM scales as RM ∝ n_eB, the power emitted by synchrotron emission scales as P ∝ n_eB² (Rybicki & Lightman 1986). In principle, this could allow RM to break a degeneracy that exists between n_e and B based on total intensity alone.

RMs have been measured for the two main EHT targets Sgr A* and M87*, as well as a handful of other low-luminosity AGN. For Sgr A*, RM = −5 × 10⁵ rad m⁻² (Bower et al. 2003; Marrone et al. 2007; Bower et al. 2018), while −7 × 10⁴ rad m⁻² has been measured a few arcseconds away (Eatough et al. 2013). For M87*, an upper limit of |RM| < 7 × 10⁵ rad m⁻² was measured using the Submillimeter Array (Kuo et al. 2014). 3C 84 has an RM of 10⁵⁻⁶ rad m⁻² (Kim et al. 2019) at 43 GHz and 9 × 10⁵ rad m⁻² at 230 GHz (Plambeck et al. 2014). Similarly, 3C 273 has an RM of 5 × 10⁵ rad m⁻² at 230 GHz (Hovatta et al. 2019). Neither linear polarization nor RM could be measured for the low-luminosity AGNs M81 and M84, which might imply significant scrambling (Brunthaler et al. 2001; Bower et al. 2017).

These measurements have been used to constrain the accretion rates of EHT targets Sgr A* and M87* by assuming simple analytic models describing the accretion flow that can be input into equation (2; Marrone et al. 2006). The accretion rates of Sgr A* and M87* have thus been constrained to <10⁻⁶ M_⊙ yr (although it may be much lower; Agol 2000; Quataert & Gruzinov 2000b; Marrone et al. 2006) and <9.2 × 10⁻⁴ M_⊙ yr (Kuo et al. 2014; Li, Yuan & Xie 2016), respectively. However, there are important effects that complicate the simple scenario implicitly assumed by equation (2), where a Faraday rotator sits entirely between a source and our line of sight (see also Broderick & McKinney 2010). Most importantly, for BH accretion flows, Faraday rotation and emission occur co-spatially, such that along a given geodesic, not all photons are Faraday rotated by the same material. In addition, general relativity (GR) complicates a photon’s trajectory and can modify its polarization properties by parallel transport alone. This is especially true of emission near the photon ring, where null geodesics can pass through the accretion flow multiple times in different directions, leading to interesting polarization signatures (Johnson et al. 2019; Himwich et al. 2020). Finally, neither the emission nor the Faraday rotation can be assumed to behave in a spatially uniform way, and the magnetic field may switch sign, especially in turbulent SANE discs.

Mościbrodzka et al. (2017) produced the first polarized model images of M87* at millimetre wavelengths from ray traced GRMHD simulations. Using a set of SANE models, they determined that Faraday rotation can be strong enough to spatially scramble the polarization from the counter-jet, which must pass through a larger Faraday depth than the forward-jet on the way to the observer. Their analysis also revealed a significant inclination dependence of the RM, such that even very high accretion rate models could satisfy RM constraints when viewed face-on. Jiménez-Rosales & Dexter (2018) determined that high-accretion rate models are disfavoured, since this scrambling too strongly depolarizes the emission. Tsunetoe et al. (2020) began to explore the spin dependence and favoured a = 0.9 models, although an expanded parameter survey is warranted.

In this work, we perform a more comprehensive survey of RM for seven models consistent with EHT5, chosen to bracket the allowed parameter space. We consider variations as a function of time and inclination, and develop new techniques to model Faraday effects in some models. In Section 2, we describe the GRMHD simulations we use as a starting point, the radiative transfer calculations in post-processing, and a novel Taylor expansion model for treating internal Faraday rotation. In Section 3, we describe our results. This includes an exploration of strong spatial variations in RM, RM distribution functions, RM as a measure of accretion rate, inclination dependence, the degree of non-λ² behaviour among models, and a case study of time variability. We discuss our results in Section 4, and end with a summary and conclusion in Section 5.

2 METHODOLOGY

We begin with a set of GRMHD simulations that are consistent with EHT5. We then use ipole¹ (Mościbrodzka & Gammie 2018) to perform polarized radiative transfer calculations, specifying electron properties in this step. Finally, using a first-order Taylor expansion, we create a model of the polarized image in order to capture frequency-dependent effects and compute RMs.

2.1 GRMHD Simulations and Radiative Transfer

In this work, we study seven models for M87* in the EHT GRMHD simulation library, each of which is consistent with all of the observational constraints considered in EHT5. The properties of these simulations, all performed with iharm (Gammie, McKinney & Tóth 2003), are listed in Table 1. Both SANE and MAD models are considered, while values of a and R_high are chosen to bracket the allowed values in EHT5. These models are described in more detail in EHT5, and the SANE a = +0.94 model is included in the GRMHD code comparison project of Porth et al. (2019).

Note that among these seven models, there are only four unique GRMHD simulations, since R_high only affects the radiative transfer in post-processing. In each of these simulations, the angular momentum of the disc is either perfectly aligned (denoted by positive spin) or antialigned (denoted by negative spin) with that of the BH, although misaligned discs remain an active area of research (e.g. Fragile et al. 2007; Chatterjee et al. 2020; Liska et al. 2020). MAD simulations are run with a 384 × 192 × 192 grid with a maximum radius of 10³ GM_•/c², while SANE simulations are run with a 288 × 128 × 128 grid and a maximum radius of 50 GM_•/c². However, we find that these models only exhibit inflow equilibrium within a radius of approximately 20 GM_•/c². Throughout this work, we restrict the integration of our radiative transfer equations to within this radius. Fortunately, this limitation actually has negligible effect on our results for face-on inclinations (i ≲ 30°), as we explore in Appendix  C. This is because the funnel region is largely evacuated in these simulations and does not contribute to Faraday rotation. However, restricting our calculations to ≤20 GM_•/c² may lead us to underestimate the total Faraday rotation at inclinations ≳30°, due to material that may exist in more distant, unequilibrated regions of the simulations.

For each model, we create images for 11 snapshots spanning the last quarter of the corresponding GRMHD run, corresponding to times t/(GM_•/c³) ∈ [7500, 10000], a duration of 880 d for M87*. For each snapshot, we study five inclinations, i ∈ {5°, 17°, 30°, 60°, 90°}.² Then, for each snapshot and inclination, we create six polarized images. We sample three frequencies, 226.999, 227.000, and 227.001 GHz, to construct a frequency-dependent model of the image, detailed in the following section. We find that we must adopt extremely small differences in frequency in order to resolve the RM on geodesics that have high Faraday depth. Otherwise, the RM in some geodesics can be underestimated due to the ‘nπ degeneracy:’ the EVPA may rotate so rapidly within that multiple rotations are missed. Then, for each of these three frequencies, we create two separate images that only include emission from either the positive or negative z domains, where the z-axis is oriented parallel to the BH spin (and perpendicular to the disc in these models).

Separate images including only the top and bottom halves of the emission are helpful for modelling Faraday effects that are often significantly different between the front and back sides of the emitting region. In these models, emission behind the mid-plane often experiences orders of magnitude more Faraday rotation than emission anterior to it, since it must pass through the relatively cold mid-plane (Mościbrodzka et al. 2017). This split remains helpful at high inclinations, since some emission is lensed to the opposite side of the image. A complete image is made by summing together the individual images made with only the near- and far-side emission. Note that when emission from only one side is included, absorption and Faraday effects from both sides remain included.

For all images, we adopt a BH mass of 6.2 × 10⁹ M_⊙ and a distance of 16.9 Mpc (Gebhardt et al. 2011), for consistency with EHT5. Gas densities are scaled such that the average image produces an intensity of 0.5 Jy at an inclination of 17°, which is the most likely inclination at which we are viewing M87* based on its larger scale jet (Walker et al. 2018). The same gas density scalings are used for models at different inclinations, although their average fluxes can depart from 0.5 Jy. We find that the total intensities of images created at an inclination of 90° tend to be approximately a factor of 2 larger than those created with an inclination of 17°. Each image is created at 0.5 |$\mu$|as pixel resolution, with a 160 |$\mu$|as field of view, a factor of 2 finer angular resolution than employed in EHT5.

To summarize, we consider seven models that are consistent with EHT observations of M87*. For each model, we create images for 11 snapshots, 5 inclinations, 3 frequencies, and the 2 sides of the accretion flow. We also create additional images to better resolve frequency-dependent (Section 2.2) and time-dependent (Section 3.7) phenomena. In Fig. 1, we plot total intensity images of the seven models during their final snapshot at 17° inclination. To better visualize low surface brightness features, we intentionally saturate 0.3 per cent of the pixels in this and all subsequent visualizations. In all of the images presented in this work, the forward-jet points straight up, and the material on the right side of the image is moving towards the observer.

Figure 1.

Total intensity images of the seven models considered in this work, taken at the final snapshot of each GRMHD simulation, which occurs at time t = 10000 GM_•/c³. Models are tilted by 17° from face-on towards the top of these images and in subsequent figures. To help visualize low surface brightness features, intensity is scaled with respect to the intensity of the pixel in the 99.7th percentile, I_99.7, saturating 0.3 per cent of the pixels.

Open in new tabDownload slide

2.2 Modelling the polarized image

Here, we introduce a modelling formalism for the polarized image as a function of frequency. For a given image snapshot, we create a two-zone (front and back half) model of the radiation in each pixel. Our model is constructed as follows:

• Treating the front and back halves of the emitting region separately, we first convert the Stokes parameters {I, Q, U, V} into their counterparts on the Poincare sphere, {I, N, ϕ, ψ} (following the notation of Shcherbakov, Penna & McKinney 2012). These variables have well-behaved Taylor series expansions, and are further discussed in Appendix  A.

• Using images at three frequencies, we then compute the first and second derivatives of the spherical Stokes parameters. Each derivative is computed in the variable that best approximates expected physical dependencies: I and N are differentiated in ν, ϕ is differentiated in λ², and ψ is differentiated in λ³.

• Using the first derivatives, we construct Taylor expansions to first order in each of the spherical Stokes parameters. These can be expressed

  $$\begin{eqnarray} I(\nu) \approx I(\nu _0) + \frac{\mathrm{ d}I}{\mathrm{ d}\nu }(\nu -\nu _0), \end{eqnarray}$$

  (5)

  $$\begin{eqnarray} N(\nu) \approx N(\nu _0) + \frac{\mathrm{ d}N}{\mathrm{ d}\nu }(\nu -\nu _0), \end{eqnarray}$$

  (6)

  $$\begin{eqnarray} \phi (\nu) \approx \phi (\nu _0) + \frac{\mathrm{ d}\phi }{\mathrm{ d}\lambda ^2}(\lambda ^2-\lambda _0^2) \approx \phi (\nu _0) - \frac{2c^2}{\nu _0^3}\frac{\mathrm{ d}\phi }{\mathrm{ d}\lambda ^2}(\nu -\nu _0), \end{eqnarray}$$

  (7)

  $$\begin{eqnarray} \psi (\nu) \approx \psi (\nu _0) + \frac{\mathrm{ d}\psi }{\mathrm{ d}\lambda ^3}(\lambda ^3-\lambda _0^3) \approx \psi (\nu _0) - \frac{3c^3}{\nu _0^4}\frac{\mathrm{ d}\psi }{\mathrm{ d}\lambda ^3}(\nu -\nu _0). \end{eqnarray}$$

  (8)
• In order to approximate bandwidth integration, we integrate analytically the equations for the Stokes parameters {I, Q, U, V}. For a band spanning the frequencies ν₁ and ν₂, the bandwidth-averaged Stokes parameters are given by

  $$\begin{eqnarray} I_\mathrm{BW} = \frac{1}{\nu _2-\nu _1}\int _{\nu _1}^{\nu _2} I(\nu) \mathrm{ d}\nu , \end{eqnarray}$$

  (9)

  $$\begin{eqnarray} Q_\mathrm{BW} = \frac{1}{\nu _2-\nu _1}\int _{\nu _1}^{\nu _2} N(\nu)\cos [\phi (\nu)]\sin [\psi (\nu)] \mathrm{ d}\nu , \end{eqnarray}$$

  (10)

  $$\begin{eqnarray} U_\mathrm{BW} = \frac{1}{\nu _2-\nu _1}\int _{\nu _1}^{\nu _2} N(\nu)\sin [\phi (\nu)]\sin [\psi (\nu)] \mathrm{ d}\nu , \end{eqnarray}$$

  (11)

  $$\begin{eqnarray} V_\mathrm{BW} = \frac{1}{\nu _2-\nu _1}\int _{\nu _1}^{\nu _2} N(\nu)\cos [\psi (\nu)] \mathrm{ d}\nu . \end{eqnarray}$$

  (12)

  Analytic solutions to these integrals are provided in Appendix  B.

• Finally, the front and back emission halves are summed to produce the complete image.

In practice, dϕ/dλ², which encapsulates Faraday rotation, is the most important frequency-dependent effect to take into account. Faraday conversion, which is encapsulated in dψ/dλ³, is not as strong as Faraday rotation in these models. In some of our models, our split into the front and back halves is necessary for capturing the Faraday depolarization of back half of the emission (due to the Faraday thick mid-plane) while preserving the polarization of the front half. As we will later explore in more detail, some pixels may simultaneously have emission from the forward-jet with an RM of ∼10³ rad m⁻², and emission from the counter-jet with an RM of ∼10⁹ rad m⁻². Consequently, ϕ(λ²) for the pixel is highly non-linear, but can be approximated as the sum of two emitting regions with distinct RMs.

In some pixels, especially within the photon ring, we compute unreasonably high values of |dN/dν| within our Taylor expansion. This is due to multiple emission components in the same pixel Faraday rotating at different rates, which can periodically increase and decrease the total linearly polarized intensity in the pixel depending on the relative phases of these components. We use |$\frac{\mathrm{ d}^2N}{\mathrm{ d}\nu ^2}$| to help identify such problematic pixels. Let Δν = ν − ν₀ be the difference in frequency space between some frequency ν at which Stokes parameters are being evaluated and the frequency at which derivatives have been computed, ν₀ = 228 GHz. Treating |$\frac{1}{2}\frac{\mathrm{ d}^2N}{\mathrm{ d}\nu ^2} \Delta \nu ^2$| as an upper limit on the error of N(ν), we freeze the value of N in pixels where |$|\frac{\mathrm{ d}N}{\mathrm{ d}\nu } \Delta \nu | \lt |\frac{1}{2}\frac{\mathrm{ d}^2N}{\mathrm{ d}\nu ^2} \Delta \nu ^2|$| and |$|\frac{\mathrm{ d}\ln (N)}{\mathrm{ d}\nu }| \gt |\frac{1}{\Delta \nu }|$|⁠. That is, we freeze N if its first derivative implies that it would grow or shrink by more than a factor of e across Δν and the absolute error on the change of N may be larger than the change of N itself. We apply an identical condition on the derivative of I, which almost always affects pixels also affected by the condition on Stokes N.

We examine to what extent this criterion is applied in our final snapshot images, and find that it affects only a minority of pixels. Among the pixels which amount to 99 per cent of the image integrated I, we find that at most 5 per cent of the pixels are affected by this correction, typically in the photon ring. This largest fraction occurs in the SANE, a = −0.94, R_high = 80 model, which as we shall show contains the emission with the largest Faraday depths. In some of the other models, none of the pixels are affected by this criterion.

In Fig. 2, we demonstrate the efficacy of our method by plotting the linearly polarized intensity (⁠|$\sqrt{Q^2+U^2}$|⁠) for the SANE, a = +0.94, R_high = 160 model. The top left-hand panel shows the result of a single-frequency calculation at 228 GHz, while the top right image shows the result of averaging 255 images across a 4 GHZ bandwidth between 226 and 230 GHz. In the bandwidth-averaged image, much of the large-scale emission has been suppressed and two streaks of emission in the upper right of the image have become more prominent. The total spatially unresolved linear polarization fraction has dropped by a factor of 2.4. Although the counter-jet dominates the linearly polarized emission in the single-frequency image due to lensing (Dexter, McKinney & Agol 2012), most of that emission is scrambled away due to the large Faraday depth in the mid-plane (Mościbrodzka et al. 2017).

Figure 2.

Images of total linearly polarized intensity (⁠|$\sqrt{Q^2+U^2}$|⁠) for a case where bandwidth depolarization is important, the SANE, a = +0.94, R_high = 160 model, centred at 228 GHz. In the top left, we plot an image at a single frequency, while in the top right, we plot the average of 255 images across a 2 GHz bandwidth. As further explored in Fig. 3, bandwidth depolarization suppresses the contribution of the counter-jet. In the bottom left, we plot the result of our Taylor expansion model based on six images. Its residual with respect to the 255 averaged images is shown in the bottom right. The Taylor expansion model successfully depolarizes the counter-jet without depolarizing the forward-jet. Since we model the image with two zones, the largest discrepancies occur in and within the photon ring, where geodesics cross the mid-plane multiple times.

Open in new tabDownload slide

The bottom left-hand panel is the result of our Taylor approximated model based on six images (three frequencies, two sides), which successfully captures the depolarization of the counter-jet, but not the forward-jet. In the bottom right-hand panel, we subtract this model from the averaged images and plot the residual. Weighting by the linear polarized intensity of the properly averaged image, total linear polarization is recovered on average with an absolute error of 0.32 |$\mu$|Jy and a relative error of 0.081. In contrast, the single-frequency image has an average absolute error 5.7 |$\mu$|Jy and relative error of 4.0. In this case, the single-frequency image also overestimates the total linear polarization somewhat, even if the source is spatially unresolved. When Stokes parameters are summed across the entire image, the properly averaged image has a linear polarization fraction of 8.3 × 10⁻³, the Taylor approximated image has a linear polarization fraction of 7.3 × 10⁻³, and the single-frequency image has a linear polarization fraction of 1.9 × 10⁻². A model containing at least two separate zones is necessary in order to reproduce the depolarization of the counter-jet without also depolarizing the forward-jet. The remaining residuals require more than two regions to fully capture the complexity of the Faraday rotating structure along the line of sight. Notice that the most significant errors occur within the photon ring or its interior, where geodesics pass through the mid-plane multiple times.

Figure 3.

Linearly polarized intensity (⁠|$\sqrt{Q^2+U^2}$|⁠) for the SANE, a = +0.94, R_high = 160 model at 227 GHz (0 bandwidth), broken into forward-jet (Top Left) and counter-jet (Top Right) emission components, which appear very different due to lensing. For the geodesic marked by a blue circle, we plot in the lower panel the frame-invariant versions of the radiative transfer coefficients that correspond to polarized emissivity |$(j_Q \rm {, red})$| and Faraday rotation |$(\rho _V \rm {, blue})$|⁠. For ρ_V, the dotted lines correspond to negative values, and there are many sign flips due to the turbulent nature of a SANE disc. Emission from the counter-jet passes through the cold mid-plane on the way to the observer, which produces a very large RM. Thus, when properly averaging over bandwidth, the counter-jet emission is suppressed by bandwidth depolarization. (See Fig. 2, top right-hand panel.).

Open in new tabDownload slide

2.3 Rotation measure

However, when discussing the RM for spatially unresolved measurements in subsequent sections, it is important to recognize the complicated evolution of χ(λ²) that results from the complex RM structure we uncover in Section 3.1. When assigning a single value of the RM to an entire spatially unresolved image, we use the Taylor expansion methodology developed in Section 2.2 to approximate 16 polarized images, each integrated over a bandwidth of 0.25 GHz, equally spaced in λ² space between 226 and 230 GHz to emulate EHT observations. These images are used to compute χ(λ²), and Δχ is computed from the endpoints of the band. We correct for phase wraps by adding or subtracting π to χ(λ²) as necessary to obtain the correct sign of dχ/dλ² based on equation (15). Using 16 bands of width 0.25 GHz, the maximum measurable RM is |$|\mathrm{RM}|_\mathrm{max} = \pi /\Delta \lambda ^2 \approx \pi \nu _0^3/(2c^2\Delta \nu) = 8.3\times 10^8 \ \mathrm{rad} \,\, \mathrm{m}^{-2}$|⁠.

3 RESULTS

3.1 Structure of spatially resolved RM

In this section, we study the RM of each of our models within single simulation snapshots. We find significant spatial variation across the image due to inhomogeneities in the accretion flow on event horizon scales. At different locations, RM can vary by orders of magnitude and even flip sign. As a result, for spatially unresolved polarized measurements, some of these models may exhibit highly non-linear χ(λ²), making them difficult to characterize with a single RM.

In Fig. 4, we visualize the spatially resolved RM structure of the images presented in Fig. 1. In this figure, the brightness of each pixel scales with the linearly polarized specific intensity |$(\sqrt{Q^2+U^2})$|⁠, while the colouration encodes the RM. As in Fig. 1, 0.3 per cent of the pixels are saturated to more clearly display low surface brightness structures. The colour scale is normalized separately for each model, spanning ± the 90th percentile of |RM| of the pixels that have at least 50 per cent of the maximum linear polarized intensity plotted. The red regions have positive RM, while the blue regions have negative RM. The top two rows depict the RM with Δν = 2 MHz, where our Taylor expansion model is constructed. The bottom two rows plot the RM across a 4 GHz band, a more realistic bandwidth, within which bandwidth depolarization becomes important for some models. Images of 4 GHz are plotted with the same brightness and RM scale as their 2 MHz counterparts.

Figure 4.

Visualization of the RM structure for the snapshots plotted in Fig. 1. The brightness of each pixel is proportional to the linearly polarized intensity in the pixel (again saturating 0.3 per cent of the pixels), while the colouration encodes the RM. The colour scale is normalized separately for each model, spanning ± the 90th percentile of |RM| of the pixels which have at least 50 per cent of the maximum linear polarized intensity. The red regions have positive RM, while the blue regions have negative RM. The top two rows are not bandwidth corrected, while the bottom two are. This can dramatically affect the SANE models, but has little effect on the MAD models. Complex structure and sign flips are a generic prediction of these simulations.

Open in new tabDownload slide

As shown in this figure, complex spatial variation and frequent sign-flips are a generic feature of these RM maps. This behaviour is not surprising in SANE models, which are characterized by weaker, disordered magnetic fields, but is less expected in MAD models, which are characterized by strong poloidal fields. In one suggestive snapshot, we confirm that these RM sign-flips are due to sign-flips in the magnetic field with respect to the geodesic. Fig. 5 plots the intensity-weighted Faraday depth in each pixel, τ_F = ∫ρ_Vds, for a snapshot of the MAD, a = 0.94, R_high = 20 model. Here, ρ_V is the (frame-invariant) radiative transfer coefficient responsible for Faraday rotation, and s is the affine parameter describing the geodesic. The sign of this quantity, shown to exhibit both positive to negative values, directly encodes the direction of the magnetic field with respect to the photon trajectory. RM sign flips have been predicted by earlier MHD simulations without GR, but only at large inclinations (Sharma et al. 2007).

Figure 5.

Intensity-weighted Faraday depth in one snapshot of the MAD, a = 0.94, R_high = 20 model, revealing clear sign flips in the magnetic field parallel to the line of sight. The colour of each pixel encodes the total Faraday depth, while the brightness is proportional to the linear polarized intensity. The sign of the Faraday depth directly encodes the direction of the magnetic field parallel to the geodesic. Subsequent 3D visualization of this snapshot reveals that these sign flips occur in the tangential magnetic field in the plane of the disc.

Open in new tabDownload slide

By performing a 3D visualization of this snapshot with visit (Childs et al. 2012), we find that sign flips in the magnetic field occur in its tangential and radial components when crossing the disc mid-plane. Near the event horizon, field lines on the Northern hemisphere point west, while field lines on the Southern hemisphere point east, although they point north overall (in the positive z-direction). This is a natural consequence of the tangential stretching of vertical field lines as they are dragged into the BH by accreting material, as well as frame dragging (see e.g. Contopoulos et al. 2009; Gabuzda 2018, for helpful schematics). Since this sign flip occurs when crossing the mid-plane, accreting streams of gas that straddle the mid-plane can exhibit streaks of positive RM adjacent to streaks of negative RM.

Returning to Fig. 4 and comparing the 2 MHz bandwidth visualizations to the 4 GHz bandwidth visualizations, SANE models are more strongly affected by bandwidth depolarization than MAD models. In each of the SANE models, counter-jet polarized emission is especially suppressed. This more strongly changes the morphology of the prograde SANE than the retrograde SANEs, because the image morphologies of the counter-jet and forward-jet in the retrograde SANEs are more similar.

Closely examining the 2 MHz counter-rotating SANE models, one may notice that the linearly polarized intensity appears to vary strongly among adjacent pixels, causing the appearance of ‘static’ across the image. This effect is an artefact of single-frequency radiative transfer calculations, an instance of Faraday rotation randomizing not only the phase of the linear polarization, but also its amplitude. In this model, the linear polarization in each pixel is well approximated by the sum of its forward-jet and counter-jet components, neither of which exhibit this ‘static’ if plotted individually. The phase of the counter-jet emission is effectively randomized by the enormous Faraday depth in the mid-plane. Depending on the relative phase, the rotated counter-jet polarized emission may sometimes cancel with the forward-jet polarized emission. This effect would not occur in real observations integrated over a finite bandwidth.

3.2 RM distribution functions for M87*

In these models, polarized emission exhibits a wide range of Faraday depths, and the front and back halves of the emitting regions can differ by many orders of magnitude. In Fig. 6, we plot the distribution functions of log₁₀|RM| among the pixels of each model during their final snapshot, weighted by each pixel’s linear polarized intensity. This quantity, which we denote as d|p|/dlog₁₀|RM|, is closely related to F(ϕ) in RM synthesis theory, the complex polarized surface brightness per unit Faraday depth (Burn 1966; Brentjens & de Bruyn 2005). Unlike F(ϕ), we take absolute values, normalize with respect to the total emitted polarized surface brightness, and adopt logarithmically spaced bins.

Figure 6.

Distribution functions of log₁₀|RM| for our models, weighted by the linear polarized intensity of each pixel. Emission originating from front side of the emitting region is coloured blue, while emission originating from the back side is coloured red. The relatively cold mid-plane is the dominant source of Faraday rotation in these models, which can result in significant offsets in these distributions between the two halves, especially in SANE models. Forward-jet emission in SANE models exhibits a second peak at high log₁₀|RM| due to photon ring geodesics, which do pass through the mid-plane.

Open in new tabDownload slide

In each pixel, |RM| is computed directly from the gradient at 228 GHz across Δν = 2 MHz. These distributions are computed for each of the 11 snapshots that we study, and the filled regions span the range permitted by all of these snapshots. The blue regions only include emission originating from the front half of the emitting region, while the red regions only include emission originating from the back half. The dashed vertical line marks RM_crit (eqn. 14) for a 4 GHz bandwidth. Any emission to the right of this line is significantly bandwidth depolarized. The distributions of front half of the emission region can be significantly displaced from that of the back half, especially in SANE models. Again, this is due to the large Faraday depth of the sub-relativistic mid-plane in these models. In some models, including all of the SANEs, the forward-jet distributions exhibit two distinct peaks. The peak at higher |RM| is due to photon ring orbits, which pass through this mid-plane. Note the extreme difference in |RM| between front and back components in the jet-dominated retrograde SANE models. These models may exhibit much lower spatially unresolved |RM| than would be expected by integrating their Faraday depths across geodesics, since this Faraday depth mainly only affects the counter-jet.

3.3 Non-linear structure of spatially unresolved RM

In Section 3.1, we demonstrate that GRMHD models exhibit rich spatial structure, whereby the RM can vary by many orders of magnitude and flip sign. As a consequence, χ(λ²) rotates at very different rates at different locations within the image, which may result in clear departures from a λ² law if these structures are not spatially resolved. In Fig. 7, we investigate this non-linearity by plotting χ(λ²) for the final snapshots of our seven models. In blue, we laboriously compute 255 individual images across the 4 GHz bandwidth, then average Stokes parameters within 16 smaller 0.25 GHz bands to estimate χ(λ²). In red, we instead use the Taylor expansion model based on six images and analytic integrals described in Section 2.2 to estimate χ(λ²).

Figure 7.

Non-λ² behaviour of the seven models we consider during their final simulation snapshot. Within 16 bands that are each 0.25 GHz wide, we plot χ(λ²) of a band-averaged image. In blue, we plot the result from 255 images evenly spaced between 226 and 230 GHz. Within each band, the appropriate subset of images is averaged. In red, we plot the result from our Taylor expansion model based on six images around 228 GHz. Most of these models exhibit non-linear behaviour even within this narrow fractional bandwidth.

Open in new tabDownload slide

All models exhibit significantly non-λ² behaviour even within this small fractional bandwidth except for the MAD R_high = 20 models. SANE models are especially non-linear, and in fact exhibit spectrally unresolved structure in χ(λ²) even among the 255 images separated in frequency by 16 MHz. This is because as shown in Fig. 6, a significant amount of the intensity has an associated |RM| of ≈10⁹ rad m⁻². The MAD |$\mathrm{R_\mathrm{high}} = 160$| models exhibit relatively mild non-linearity, since |RM| only just approaches |RM_crit| in Fig. 6.

Recall that we use Δχ/Δλ² across the bandwidth from our Taylor expansion model to assign spatially unresolved RMs to images. This figure reveals some of this model’s limitations. The model poorly reproduces χ(λ²) for the MAD, a = +0.94, R_high = 160 model, possibly due to 228 GHz being a local extremum of χ(λ²) where the first derivative is small. Interestingly, the retrograde SANE Taylor expansion models appears to broadly follow the structure of the true χ(λ²), but with a vertical offset. This is due to incorrect evolution of emission superposed on top of the photon ring. Our Taylor expansion model assigns a large dϕ/dλ² to photon ring pixels, but these pixels also contain forward-jet emission that does not pass through the mid-plane. In the SANE, a = −0.94, R_high = 80 model, this superposed component is immediately bandwidth depolarized and subtracted, leading to the offset. This effect is more delayed as a function of Δλ² in the SANE, a = −0.94, R_high = 10 model, which by construction has a warmer mid-plane and therefore less Faraday rotation. Fortunately, this effect is symmetric about the Taylor expansion point at 228 GHz and we can still recover the spatially unresolved RM from χ at the end points of the band.

3.4 RM as a measure of accretion rate

Here, we test the relationship between RM and accretion rate in our suite of images for M87*. In Fig. 8, we plot spatially unresolved RM versus accretion rate for the seven models considered in this work. Each symbol of the same colour represents a different snapshot of the same model. Positive RMs are plotted with the filled symbols, while negative RMs are plotted with the open symbols. Notice the ubiquitous flips in the sign of the RM that occur in all models. The blue-filled region demarcates the relation in Marrone et al. (2006), spanning variations of the slope of the density profile, β ∈ [0.5, 1.5]. We set r_in = 3 and r_out = ∞. The dashed line shows the upper limit from Kuo et al. (2014). A 17° inclination appropriate for M87* is shown on the left, while for comparison a 90° inclination is shown on the right.

Figure 8.

RM as a function of accretion rate for the seven models considered in this paper. The filled symbols have positive RM, while the open symbols have negative RM. In the left-hand panel, the observer is oriented at a 17° inclination, while in the right-hand panel, the observer is oriented at 90°. At an inclination of 17°, we find that accretion rates are systematically higher than those that would be inferred by simple analytic models (Marrone et al. 2006). At 90°, we find RMs in better agreement with analytic models, although there remains substantial scatter. Retrograde SANE models are outliers, since their forward-jet emission does not intercept the large Faraday depth in the mid-plane.

Open in new tabDownload slide

Overall, we find that a spatially unresolved RM is a poor predictor of the accretion rate, especially if the correct model is not known a priori. RM and the accretion rate differ by orders of magnitude both within and among the different models. Even within a single model, there is no correlation between RM and accretion rate, which we explore in more detail for one model in Section 3.7. Since their higher accretion rates imply higher number densities, SANE models typically have larger spatially unresolved RMs than MAD models. However, these models exhibit such strong time variability that they cannot be distinguished solely by the upper limit of Kuo et al. (2014). Rather, repeat observations on time-scales of months to years will be necessary to characterize the distribution of |RM| over time and detect potential sign flips.

For a given accretion rate, the GRMHD models in this work produce much lower RM than the analytic model of Marrone et al. (2006) at an inclination of 17°. As we further explore in Section 3.5, this is because these simulations are viewed through an evacuated funnel region at low inclination. With an inclination of 90°, the |RM| more closely matches that predicted by Marrone et al. (2006), although they still remain systematically offset. A similar inclination dependence is found for SANE models in Mościbrodzka et al. (2017). The retrograde SANE models remain the most offset from the analytic model. Even at 90°, the large Faraday depth occurs in an area with little emission, since the electrons are assigned low temperatures in the mid-plane.

3.5 Dependence of RM on inclination

Here, we study the dependence of RM on inclination in greater detail. In Fig. 9, we plot the distributions of RM for all 11 snapshots of all seven models at five different inclinations. In this plot, boxes contain the 25th to 75th quantiles, the horizontal black or yellow line marks the median, and the error bars span the full range of the 11 snapshots studied. This plot omits the sign flips observed in Fig. 8, but we comment that they remain ubiquitous at all inclinations for these calculations that terminate at a radius of 20 GM_•/c².

Figure 9.

RM as a function of inclination. A total of 11 snapshots are shown for each model. For M87*, 17 deg is considered the most likely inclination based on its large-scale jet. We plot the Kuo et al. (2014) upper limit on the RM with a dashed line. Boxes contain the 25th to 75th quantiles, the horizontal black or yellow lines mark the median, and the error bars span the full range of the 11 snapshots considered. We report a noticeable inclination dependence, due to the evacuated jet region through which the BH is observed at low inclination.

Open in new tabDownload slide

As in Mościbrodzka et al. (2017), we find that the absolute value of the RM depends on the inclination angle, but it is not large compared to the substantial scatter between snapshots. This weaker dependence is likely due to the small radius at which we truncate our calculations. We notice no differences between our calculations at low inclination: 5°, 17°, and 30°. This is fortunate for our study of M87*, as this indicates that we need not be concerned with small deviations from our fiducial inclination of 17°.

In Fig. 10, we plot the characteristic distance of Faraday rotating material of these models during their final snapshot as a function of inclination. For inclinations <30°, most of the Faraday rotation occurs at <10 GM_•/c², while 〈R_FR〉 increases at higher inclinations. The innermost stable circular orbit exists at smaller radius for prograde models than for retrograde models, which leads to a noticeable difference in 〈R_FR〉 between these two classes of models at low inclination.

Figure 10.

Characteristic distance of Faraday rotating material in these models during their final snapshot as a function of inclination. Due to the evacuated jet region in these simulations, the Faraday rotating material is confined to low radius at low inclination, but extends to larger radius at higher inclination. Recall that our calculations terminate at R = 20 GM_•/c², within which these simulations are in inflow equilibrium.

Open in new tabDownload slide

Recall that we restrict the domain of our calculations to within 20 GM_•/c², the radius within which the simulations are in inflow equilibrium. As we further explore in Appendix  C, if this radius is increased to 50M_•, we find consistent results for i ≤ 17°, but substantially larger 〈R_FR〉 for i ≥ 60°. Therefore, we believe that our |RM| values throughout the paper should be considered lower limits for i > 17°, as material from beyond the converged region may contribute to Faraday rotation at larger inclinations.

Since the RM is highly non-uniform across these images, spatially resolved RM distribution functions provide greater insight into the inclination dependence. In Fig. 11, we plot the RM distribution functions as in Fig. 6, for inclinations i ∈ {5°, 30°, 60°, 90°}. Unlike in Fig. 6, we do not split the distributions into their front and back halves. At low inclination, these models exhibit significant emission with low |RM|. As the inclination approaches 90°, this fraction of emission with low |RM| diminishes, and the distribution is skewed towards higher values. Note that the distributions at 5° are indistinguishable from those at 30°.

Figure 11.

RM distribution functions as in Fig. 6, now shown as a function of inclination. Both halves of the emitting region are combined in this figure. At higher inclinations, we find that these distributions are skewed towards higher values, as the population of photons experiencing comparatively little RM diminishes.

Open in new tabDownload slide

3.6 Non-Linear χ(λ²) as a model discriminant

Here, SS_res is the regression sum of squares, and SS_tot is the total sum of squares. R² describes the fraction of variation within the data that can be ascribed to the simple linear dependence.

Our results are shown in Fig. 12, following the same formatting as Fig. 9. Based on these results, SANE models should exhibit non-linear behaviour most of the time. In contrast, MAD models are almost always well described by a linear χ(λ²) law, especially those with R_high = 20. This can be understood by returning to Fig. 6. Pixels with |RM| > |RM|_crit have individual EVPAs that rotate substantially across the 4 GHz bandwidth. Images consisting of a substantial fraction of such pixels will therefore exhibit structure in χ(λ²) within the bandwidth. Among the models considered in this study, this behaviour appears much more likely among SANEs.

Figure 12.

Non-linearity of χ(λ²) among the models considered in this study. R² of a linear fit to χ(λ²) describes the fraction of variation in the data that can be explained by a simple linear model. SANE models exhibit more non-linearity than MAD models in this study, since their images contain a significant amount of polarized intensity with |RM| > |RM|_crit.

Open in new tabDownload slide

3.7 Case study: RM time variability

We study one model with higher time resolution in order to quantify the variability of its RM. For the MAD, a = 0.94, R_high = 20 model, we create images for every available snapshot within t/(GM_•/c³) ∈ [7500, 10000], which are each separated by 5 GM_•/c³. This model has the smallest accretion rate and |RM| of the models explored in this work. Its |RM| is sufficiently smaller than RM_crit that χ(λ²) remains linear within the bandwidth (see Fig. 7), and thus we create only two polarized images (two frequencies) at each snapshot instead of the usual six (three frequencies, each separately for two sides of the disc).

The time variability of this model at an inclination of 17° is visualized in Fig. 13. In the top row, panels are separated by about half a year, while in the bottom row panels are separated by about 5 d. As in previous figures, pixel brightness encodes the linear polarized intensity, while the colour encodes the RM. Both positive and negative RM regions can be found in a typical snapshot. The spatially unresolved RM is written at the bottom of each panel. The bottom row illustrates how the dynamics of an image with both positive and negative RM regions can result in variability and RM sign flips on a time-scale of a few days. The RM sign flip does not require a drastic change in global source structure; rather, the balance between positive and negative RM regions is shifted.

Figure 13.

Visualization of the time variability of the RM structure in the MAD, a = 0.94, R_high = 20 model. The brightness of each pixel scales with its linear polarized intensity, while the colour represents its rotation measure. The RM that would be inferred from a spatially unresolved measurement is written at the bottom of each panel. The spatially unresolved RM can change on time-scales of days as the Faraday rotating gas moves on event horizon scales. The sign flip does not require a dramatic change in the source structure.

Open in new tabDownload slide

In Fig. 14, we plot the RM as a function of time for this model, as well as its autocorrelation function. The geometrized time unit is converted to days via t_• = GM_•/c³, which for M87* is 8.5 hr. The grey band encloses the 16th–84th (1σ)percentiles. At 17°, these include both positive and negative values, such that |$\mathrm{RM} = {-0.00}^{+1.27}_{-0.83} \times 10^5 \ \mathrm{rad} \,\, \mathrm{m}^{-2}$|⁠. Examining the images, we do not notice any obvious special behaviour in the accretion flow during periods of large |RM|. Since the RM is determined by the motion of material on event horizon scales, the autocorrelation function of this time series drops rapidly, falling below 0.5 in less than the separation between snapshots.

Figure 14.

Left: RM as a function of time in the MAD, a = 0.94, R_high = 20 model that we study with greater time resolution. The grey band encloses the 1σpercentile region over this time, which includes both positive and negative values. Right: Autocorrelation function of this time series. The autocorrelation drops below 50 per cent in less than the separation between snapshots.

Open in new tabDownload slide

In Fig. 15, we plot the joint probability distributions of log₁₀|RM| with accretion rate, linear polarized intensity, and circular polarized intensity for the model at 17°. 1σ, 2σ, and 3σ contours are overlaid in white. We do not find any correlation between |RM| and |$\dot{M}_\bullet$|⁠, indicating that within a single model, a change in RM does not imply a change in the accretion rate, as might be suggested by analytic models. Rather, as we have discussed, |RM| and its sign appears to result from a complicated and stochastic cancellation of positive and negative regions. For Sgr A*, Bower et al. (2018) found an anticorrelation between linear polarized intensity and RM, but no such correlation with circular polarized intensity. For this particular model of M87*, we recover qualitatively similar results: a linear regression yields log₁₀|RM| = −1.2log₁₀LP + 2.7 ± 0.08, where |$LP = \sqrt{Q^2+U^2}$| in Jy and RM is in units of rad m⁻², with a moderate r-value of −0.57. No statistically significant correlation is found between |RM| and circular polarization. An anticorrelation between |RM| and LP is not surprising, since greater Faraday rotation implies greater scrambling of the polarization vector field.

Figure 15.

Joint probability distributions of log₁₀|RM| with accretion rate, linear polarized intensity, and circular polarized intensity in the MAD, a = 0.94, R_high = 20 model with an inclination of 17°. 1σ, 2σ, and 3σ contours are overlaid in white, using the solid, dashed, and dotted lines, respectively. We find no correlation between |RM| and |$\dot{M}_\bullet$|⁠, implying that a change in RM does not imply a change in the accretion rate. As observed by Bower et al. (2018) for Sgr A*, we recover an anticorrelation between |RM| and linear polarization, but not circular polarization.

Open in new tabDownload slide

4 DISCUSSION

4.1 Comparing different models

Here, we summarize the qualitative commonalities and differences between the different models we have considered.

• Prograde MAD: These models require the lowest accretion rate to generate the appropriate total intensity, and consequently exhibit the lowest |RM|. Compared to the other models, there is not too much difference in |RM| between the two halves of the emitting region, since both components occur close to the mid-plane. These models usually exhibit linear χ(λ²) within 4 GHz.

• Retrograde MAD: Retrograde MAD models require larger accretion rates than their prograde counterparts, but exhibit similar values of |RM|. Some areas of the R_high = 160 models have large enough |RM| to weakly bandwidth depolarize or scramble portions of the image.

• Prograde SANE: Due to lensing, the counter-jet spans a larger angular scale and contributes more to the total intensity than the forward-jet. Bandwidth depolarization effects are severe, and emission from the counter-jet is entirely depolarized, assuming a 4 GHz bandwidth. Consequently, the total intensity image (dominated by the counter-jet) appears morphologically different compared to the linearly polarized image (dominated by the forward-jet). χ(λ²) exhibits strongly non-linear evolution.

• Retrograde SANE: The counter-jet experiences six orders of magnitude more Faraday rotation than the forward-jet. As with their prograde counterparts, χ(λ²) is highly non-linear. However, the difference in morphology between the total intensity image and the linearly polarized image is less significant than the prograde case, since the two emission components subtend more similar angular scales.

By construction, increasing R_high decreases the temperature of electrons in the mid-plane, which therefore increases the Faraday depth for emission that passes through it. In all models, larger R_high results in greater non-linearity of χ(λ²).

4.2 RM and bandwidth depolarization

If anywhere in an image, RM > RM_crit, then that region’s linear polarization should be suppressed as described by equation (13). We find that this is more likely to occur in SANE models, but may also affect some parts of MAD models. Bandwidth depolarized regions manifest as areas with lower than average linear polarization fraction.

In some of the worst cases, like that presented in Fig. 2, we find that basic image properties are affected by bandwidth depolarization, such as the linear polarization fraction as well as the morphology of the linear polarization vector field, even after images are blurred. This may be important for studies which use images computed at a single frequency to compare to observations taken over a finite bandwidth. For future studies, the methodology introduced in Section 2.2 can be generalized by applying the Taylor expansion to each point along the ray-traced geodesic instead of just two separate emission regions. This would allow the appropriate bandwidth integrations to occur within the ray-tracing code itself.

4.3 Caveats and limitations

At present, the electron distribution function is poorly constrained. For creating these images, a thermal distribution function is assumed, along with the R_high prescription developed by Mościbrodzka et al. (2016). Mao, Dexter & Quataert (2017) studied the effects of adding a power-law component to the distribution function and found that even if a few per cent of the total energy is put into a non-thermal power-law component, a diffuse halo of emission can be produced. The effects of non-thermal electron distribution functions on polarized images remain to be studied.

Recall that we truncate our radiative transfer calculations at a radius of 20 GM_•/c², only within which these GRMHD simulations exhibit inflow equilibrium. However, Faraday rotating material can plausibly exist at larger radius, especially at higher inclinations. Using very long-duration simulations, Dexter et al. (2020) find that Faraday rotation can peak at radii R ∼ 30−90 GM_•/c², depending on the electron prescriptions. A more distant Faraday screen would be expected to uniformly rotate all EVPAs by a fixed amount, which may leave signatures in the EVPA vector field (Palumbo et al. 2020). Additionally, such a screen should maintain a consistent RM for longer time-scales than these models. Such distant screens may in fact be required to explain the consistent RMs of Sgr A* (Bower et al. 2018) and 3C 84 (Plambeck et al. 2014) over time-scales of years.

5 SUMMARY AND CONCLUSIONS

We have investigated RM within a subset of the EHT simulation library that is consistent with the observational constraints on M87* considered in EHT5. We find more information in the RM structure of GRMHD simulations than can be described by a single scalar, RM. We summarize our results below:

• In a single snapshot, we find extreme variations in the RM between different regions. The RM may vary by orders of magnitude and even flip sign across the image. The RM inferred from a spatially unresolved measurement is therefore the result of the complicated interplay of these different regions.

• Emission originating from in front of the disc mid-plane may be orders of magnitude less Faraday rotated than emission from the back. In the high accretion rate SANE models, the RM is large enough to completely depolarize emission from the counter-jet, which may in fact dominate the total intensity. The subrelativistic mid-plane is the dominant source of Faraday rotation in these models.

• Many models exhibit clear departures from a λ² law even across a narrow fractional bandwidth of 4 GHz. Non-linearity is a more common feature among the SANE models, and increases with R_high.

• The RM structure changes as material moves on event horizon scales. These models all exhibit strong time variability, causing the spatially unresolved RM to vary and even flip sign on a time-scale of days.

• RM is a poor predictor of the accretion rate. These models predict several orders of magnitude spread in RM for a given accretion rate, and within a single model, these quantities are not correlated as a function of time. In addition, analytic models used to infer the accretion rate based on the RM (Marrone et al. 2006) systematically underestimate the accretion rate on to M87*, since the source should be viewed through an evacuated funnel region.

In future work, a more thorough investigation of the EHT simulation library is merited, including models for Sgr A*. Alternative models for the electron distribution function should also be considered. Repeated observations of both Sgr A* and M87* will be useful to probe the time variability predicted by these models.

ACKNOWLEDGEMENTS

We thank Jason Dexter for valuable feedback that greatly improved the content of this paper. We also thank Laurent Loinard and the anonymous referee for their careful and thorough reviews. This material is based upon work supported by the National Science Foundation under Grant No. OISE 1743747 and NSF AST 1716327. Computations were performed using the resources of the Black Hole Initiative (BHI). Computations at the BHI were made possible through the support of grants from the Gordon and Betty Moore Foundation and the John Templeton Foundation. This work used the Extreme Science and Engineering Discovery Environment resource stampede2 at Texas Advanced Computing Center (TACC) through allocation TG-AST170024. The opinions expressed in this publication are those of the author(s) and do not necessarily reflect the views of the Moore or Templeton Foundations.

DATA AVAILABILITY

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

Footnotes

REFERENCES

APPENDIX A: SPHERICAL STOKES PARAMETERS

In a pixel with a large |RM|, Faraday rotation causes Q and U to oscillate rapidly with frequency. In spherical Stokes parameters, N remains stable, while ϕ changes linearly with wavelength squared, with dϕ/dλ² = 2RM. This makes the Spherical stokes parameters more stable to Taylor expansion.

APPENDIX B: ANALYTIC BANDWIDTH INTEGRALS

APPENDIX C: EFFECTS OF CHANGING THE MAXIMUM INTEGRATION RADIUS

In ipole, the parameter rmax_geo (henceforth R_out) sets the radius in Boyer–Lindquist (or Kerr–Schild) coordinates within which radiative transfer coefficients are calculated. Although the MAD and SANE simulations have outer boundaries of 10³ GM_•/c² and 50 GM_•/c², respectively, we find that these simulations exhibit inflow equilibrium only within 20 GM_•/c². Hence, for the results throughout this paper, R_out is therefore set to 20 GM_•/c². Here, we explore the effects of changing this outer radius to 50 GM_•/c². Although this choice now includes material from unconverged regions, this allows us to gain some insight into how gas in more distant regions might affect our predictions.

In Fig. C1, we compare the characteristic distance of Faraday rotating material, 〈R_FR〉, as in Fig. 10. The bold solid lines originate from the R_out = 20 GM_•/c² models (as shown in Fig. 10), while the faint thin lines originate from the R_out = 50 GM_•/c² models. Interestingly, there is little difference in 〈R_FR〉 for inclinations i ≤ 17°, an evacuated funnel region in these simulations. At larger inclinations, distant material from unconverged regions could potentially contribute the majority of the Faraday rotation.

Figure C1.

Characteristic distance of Faraday rotating material, as in Fig. 10, for R_out = 20 GM_•/c² models as the solid lines, and R_out = 50 GM_•/c² models as the faint lines. There is little difference for inclinations i ≤ 17°, where the BH is viewed through an evacuated funnel region. In contrast, material from more distant, unconverged regions can dominate the Faraday rotation at higher inclinations if calculations are allowed to proceed into this area.

Open in new tabDownload slide

In Fig. C2, we plot the RM distribution functions as in Fig. 11, where our R_out = 20 GM_•/c² results are shown as the solid lines and alternative R_out = 50 GM_•/c² results are plotted as the dotted lines. For clarity, we only plot the median values at a given log₁₀|RM| instead of the full range plotted in Fig. 11. We find that there is negligible difference in our results for inclinations i ≤ 30°, the range relevant for M87*. At higher inclinations, the MAD distributions are skewed towards higher values for i ≥ 60°, while the retrograde SANE models only differ at i = 90°. The prograde SANE model shows negligible difference between R_out = 20 GM_•/c² and R_out = 50 GM_•/c² even at i = 90°.

Figure C2.

Rotation measure distribution functions, as in Fig. 6, where R_out = 20 GM_•/c² models are shown as the solid lines and R_out = 50 GM_•/c² models are shown as the dotted lines. For clarity, only median values at a given log₁₀|RM| are plotted. There is negligible difference for all models when the inclination ≤30°. The distributions are skewed towards higher values at 90° for the retrograde SANEs, and for ≥60° for MADs.

Open in new tabDownload slide

In Fig. C3, we plot the effect this has on the spatially unresolved RM observed for these sources. Here, the solid boxes correspond to the R_out = 20 GM_•/c² models (as in Fig. 9), while the faint boxes correspond to R_out = 50 GM_•/c² models. Lines demarcating medians have been removed for clarity. As expected, there is no noticeable difference for inclinations i ≤ 30°. For inclinations of ≥60°, the RMs of R_out = 50 GM_•/c² models can be a factor of a few to orders of magnitude larger, depending on the model.

Figure C3.

Rotation measure as a function of inclination, as in Fig. 9, for R_out = 20 GM_•/c² models as the solid boxes, and R_out = 50 GM_•/c² models as the faint boxes. While including material at R_out > 20 GM_•/c² makes little difference for inclinations i ≤ 30°, it may increase the RM by factors of a few to orders of magnitude at larger inclinations, depending on the model.

Open in new tabDownload slide

Finally, we notice that at inclinations of 60°, RM sign flips occur less frequently for the R_out = 50 GM_•/c² models than for the R_out = 20 GM_•/c² models. This is to be expected, since material at larger radii evolves on longer time-scales.