Decomposing the Internal Faraday Rotation of Black Hole Accretion Flows
Angelo Ricarte, Ben S. Prather, George N. Wong, Ramesh Narayan, Charles Gammie, Michael Johnson
MMNRAS , ?? – ?? (2020) Preprint 8 September 2020 Compiled using MNRAS L A TEX style file v3.0
Decomposing the Internal Faraday Rotation of Black HoleAccretion Flows
Angelo Ricarte , , Ben S. Prather , George N. Wong , Ramesh Narayan , ,Charles Gammie , , and Michael D. Johnson , Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA Black Hole Initiative, 20 Garden Street, Cambridge, MA 02138, USA Department of Physics, University of Illinois at UrbanaâĂŞChampaign, 1110 West Green Street, Urbana, IL 61801, USA Department of Astronomy, University of Illinois at UrbanaâĂŞChampaign, 1002 West Green Street, Urbana, IL 61801, USA
ABSTRACT
Faraday rotation has been seen at millimeter wavelengths in several low luminosityactive galactic nuclei, including Event Horizon Telescope (EHT) targets M87* andSgr A*. The observed rotation measure (RM) probes the density, magnetic field, andtemperature of material integrated along the line of sight. To better understand howaccretion disc conditions are reflected in the RM, we perform polarized radiative trans-fer calculations using a set of general relativistic magneto-hydrodynamic (GRMHD)simulations appropriate for M87*. We find that in spatially resolved millimetre wave-length images on event horizon scales, the RM can vary by orders of magnitude andeven flip sign. The observational consequences of this spatial structure include signif-icant time-variability, sign-flips, and non- λ evolution of the polarization plane. Forsome models, we find that internal rotation measure can cause significant bandwidthdepolarization even across the relatively narrow fractional bandwidths observed by theEHT. We decompose the linearly polarized emission in these models based on theirRM and find that emission in front of the mid-plane can exhibit orders of magnitudeless Faraday rotation than emission originating from behind the mid-plane or withinthe photon ring. We confirm that the spatially unresolved (i.e., image integrated)RM is a poor predictor of the accretion rate, with substantial scatter stemming fromtime variability and inclination effects. Models can be constrained with repeated ob-servations to characterise time variability and the degree of non- λ evolution of thepolarization plane. Key words: accretion, accretion discs — black hole physics — galaxies: individual(M87) — magnetohydrodynamics (MHD) — polarization — techniques: polarimetric
Using a network of millimetre telescopes around the world,the Event Horizon Telescope (EHT) has recently producedthe first images of a black hole (BH) accretion flow (EventHorizon Telescope Collaboration et al. 2019a,b,c,d,e,f).These images resolve the “shadow” of the supermassive BHM87*, corresponding to rays that begin on its event horizon,providing new constraints on the properties of the BH andits accretion disc.In Event Horizon Telescope Collaboration et al. (2019e),henceforth EHT5, a library of general relativistic magneto-hydrodynamic (GRMHD) simulations was produced to com-pare to the EHT images, exploring three fundamental quan- tities. The first is the strength of the magnetic field: mod-els that accumulate strong magnetic flux around the BHare able to counteract the ram pressure of in-falling gaswith magnetic pressure, resulting in what is termed a Mag-netically Arrested Disc (MAD) (Igumenshchev et al. 2003;Narayan et al. 2003; Chael et al. 2019). In contrast, theweaker and more turbulent magnetic fields in Standard AndNormal Evolution (SANE) models have a weaker effect onthe gas dynamics of the disc (Narayan et al. 2012; Sądowskiet al. 2013; Ryan et al. 2018). The second quantity is theBH’s angular momentum, described by the dimensionlessspin parameter a ≡ Jc/ ( GM • ) ∈ ( − , , where negativevalues correspond to counter-rotating accretion discs. Massand spin are the only properties intrinsic to an astrophys-ical BH, but BH spins are constrained much more looselythan their masses. Most of our understanding of supermas- c (cid:13) a r X i v : . [ a s t r o - ph . H E ] S e p Ricarte et al. sive BH spin evolution originates from theory (e.g., Bardeen& Wagoner 1969; Thorne 1974; Gammie et al. 2004; Kinget al. 2008; Volonteri et al. 2013). The third quantity ex-plored in this work is R high , one parameterisation of therelative temperatures of electrons and ions in the plasma.Such a prescription is necessary because the mean free pathnear the event horizon is so large, the two populations departfrom thermal equilibrium and divide into a two-temperatureplasma (Shapiro et al. 1976; Rees et al. 1982; Narayan & Yi1995a; Sądowski et al. 2017; Ryan et al. 2018). By combin-ing EHT imaging with other multi-wavelength constraintssuch as the X-ray flux and jet power, EHT5 rule out all but19/60 possible models which span these three variables. Inparticular, all a = 0 models are excluded.EHT constraints on M87* thus far have only consideredtotal intensity (Stokes I ), when in fact polarized visibilitieshave been obtained (Stokes I , Q , U , and V ), but have notbeen published. Hence, current constraints have only utilized1/4 of the measured information. Synchrotron emission, theemission mechanism at millimetre wavelengths, has a near-unity intrinsic polarization fraction (Le Roux 1961; Bromleyet al. 2001; Broderick & Blandford 2003). Polarimetric imag-ing is consequently predicted to tightly constrain accretionmodels, which differ substantially in their linear polariza-tion fractions as well as their morphologies (Palumbo et al.2020). Although an image has not yet been constructed forSgr A*, strong and variable linear polarization has been ob-served for this source. Previous very long baseline interfer-ometry (VLBI) measurements have revealed a partially or-dered magnetic field structure at its centre (Johnson et al.2015).As polarized radiation travels through a magnetizedplasma, it is transformed by effects sensitive to the localdensity, temperature, and magnetic field. Faraday rotationturns the electric vector position angle (EVPA) of linearlypolarized emission, and Faraday conversion exchanges lin-early and circularly polarized radiation (Sazonov 1969; Ry-bicki & Lightman 1986; Melrose 1997; Shcherbakov 2008).Future analyses with the EHT will be able to further dis-tinguish models in the time and frequency domains (e.g.,Broderick & Loeb 2006a,b; Roelofs et al. 2017; Medeiroset al. 2018). It is therefore timely for us to better under-stand time and frequency dependent effects which may helpus constrain accretion models.One such polarimetric observable is the rotation mea-sure (RM), defined by the change in the EVPA, χ , as afunction of the change of wavelength squared, RM = χ − χ λ − λ , (1)where the subscripts 1 and 2 denote two different wave-lengths. RM probes Faraday rotation, and has been usedin a variety of contexts to infer magnetic field properties(e.g., Zavala & Taylor 2004; Brentjens & de Bruyn 2005;Frick et al. 2011; Pasetto et al. 2018; Agudo et al. 2018).For a source of polarized emission that is entirely behinda Faraday rotating medium, the RM can be written as anintegral of plasma properties along the line-of-sight via RM = 8 . × rad m − (cid:90) observersource f rel (Θ e ) n e − B || G ds pc , (2)where n e is the electron number density, B || is the parallelcomponent of the magnetic field, and f rel is a correction termsuppressing Faraday rotation at relativistic temperatures(Gardner & Whiteoak 1966). At relativistically hot tempera-tures, f rel (Θ e ) ≈ log(Θ e ) / (2Θ e ) , whereas at sub-relativistictemperatures, f rel asymptotes to 1. Here Θ e ≡ k B T e /m e c , k B is the Boltzmann constant, T e is the electron tempera-ture, m e is the electron rest mass, and c is the speed of light(Jones & Odell 1977).As seen from equation 2, the RM directly traces theelectron temperature, number density, and magnetic fieldalong the line of sight. Among the models consistent with theEHT data of M87*, all of these quantities can vary by ordersof magnitude. Note that while RM scales as RM ∝ n e B , thepower emitted by synchrotron emission scales as P ∝ n e B (Rybicki & Lightman 1986). In principle, this could allowRM 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 tar-gets Sgr A* and M87*, as well as a handful of other low-luminosity AGN. For Sgr A*, RM = − × rad m − (Bower et al. 2003; Marrone et al. 2007; Bower et al. 2018)while − × rad m − has been measured a few arc-seconds away (Eatough et al. 2013). For M87*, an upperlimit of | RM | < × rad m − was measured usingthe Submillimeter Array (Kuo et al. 2014). 3C 84 has anRM of − rad m − (Kim et al. 2019) at 43 GHz and × rad m − at 230 GHz (Plambeck et al. 2014). Sim-ilarly, 3C 273 has an RM of × rad m − at 230 GHz(Hovatta et al. 2019). Neither linear polarization nor RMcould be measured for the low-luminosity AGNs M81 andM84, which might imply significant scrambling (Brunthaleret al. 2001; Bower et al. 2017).These measurements have been used to constrain theaccretion rates of EHT targets Sgr A* and M87* by assum-ing simple analytic models describing the accretion flow thatcan be input into equation 2 (Marrone et al. 2006). The ac-cretion rates of Sgr A* and M87* have thus been constrainedto < − M (cid:12) yr (although it may be much lower) (Agol2000; Quataert & Gruzinov 2000b; Marrone et al. 2006) and < . × − M (cid:12) yr (Kuo et al. 2014; Li et al. 2016) respec-tively. However, there are important effects that complicatethe simple scenario implicitly assumed by equation 2, wherea Faraday rotator sits entirely between a source and ourline-of-sight (Broderick & McKinney 2010, see also). Mostimportantly, for BH accretion flows, Faraday rotation andemission 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’strajectory and can modify its polarization properties by par-allel transport alone. This is especially true of emission nearthe photon ring, where null geodesics can pass through theaccretion flow multiple times in different directions, leadingto interesting polarization signatures (Johnson et al. 2019;Himwich et al. 2020). Finally, neither the emission nor theFaraday rotation can be assumed to behave in a spatiallyuniform way, and the magnetic field may switch sign, espe-cially in turbulent SANE discs.Mościbrodzka et al. (2017) produced the first polarizedmodel images of M87* at millimetre wavelengths from ray MNRAS , ?? – ????
Using a network of millimetre telescopes around the world,the Event Horizon Telescope (EHT) has recently producedthe first images of a black hole (BH) accretion flow (EventHorizon Telescope Collaboration et al. 2019a,b,c,d,e,f).These images resolve the “shadow” of the supermassive BHM87*, corresponding to rays that begin on its event horizon,providing new constraints on the properties of the BH andits accretion disc.In Event Horizon Telescope Collaboration et al. (2019e),henceforth EHT5, a library of general relativistic magneto-hydrodynamic (GRMHD) simulations was produced to com-pare to the EHT images, exploring three fundamental quan- tities. The first is the strength of the magnetic field: mod-els that accumulate strong magnetic flux around the BHare able to counteract the ram pressure of in-falling gaswith magnetic pressure, resulting in what is termed a Mag-netically Arrested Disc (MAD) (Igumenshchev et al. 2003;Narayan et al. 2003; Chael et al. 2019). In contrast, theweaker and more turbulent magnetic fields in Standard AndNormal Evolution (SANE) models have a weaker effect onthe gas dynamics of the disc (Narayan et al. 2012; Sądowskiet al. 2013; Ryan et al. 2018). The second quantity is theBH’s angular momentum, described by the dimensionlessspin parameter a ≡ Jc/ ( GM • ) ∈ ( − , , where negativevalues correspond to counter-rotating accretion discs. Massand spin are the only properties intrinsic to an astrophys-ical BH, but BH spins are constrained much more looselythan their masses. Most of our understanding of supermas- c (cid:13) a r X i v : . [ a s t r o - ph . H E ] S e p Ricarte et al. sive BH spin evolution originates from theory (e.g., Bardeen& Wagoner 1969; Thorne 1974; Gammie et al. 2004; Kinget al. 2008; Volonteri et al. 2013). The third quantity ex-plored in this work is R high , one parameterisation of therelative temperatures of electrons and ions in the plasma.Such a prescription is necessary because the mean free pathnear the event horizon is so large, the two populations departfrom thermal equilibrium and divide into a two-temperatureplasma (Shapiro et al. 1976; Rees et al. 1982; Narayan & Yi1995a; Sądowski et al. 2017; Ryan et al. 2018). By combin-ing EHT imaging with other multi-wavelength constraintssuch as the X-ray flux and jet power, EHT5 rule out all but19/60 possible models which span these three variables. Inparticular, all a = 0 models are excluded.EHT constraints on M87* thus far have only consideredtotal intensity (Stokes I ), when in fact polarized visibilitieshave been obtained (Stokes I , Q , U , and V ), but have notbeen published. Hence, current constraints have only utilized1/4 of the measured information. Synchrotron emission, theemission mechanism at millimetre wavelengths, has a near-unity intrinsic polarization fraction (Le Roux 1961; Bromleyet al. 2001; Broderick & Blandford 2003). Polarimetric imag-ing is consequently predicted to tightly constrain accretionmodels, which differ substantially in their linear polariza-tion fractions as well as their morphologies (Palumbo et al.2020). Although an image has not yet been constructed forSgr A*, strong and variable linear polarization has been ob-served for this source. Previous very long baseline interfer-ometry (VLBI) measurements have revealed a partially or-dered magnetic field structure at its centre (Johnson et al.2015).As polarized radiation travels through a magnetizedplasma, it is transformed by effects sensitive to the localdensity, temperature, and magnetic field. Faraday rotationturns the electric vector position angle (EVPA) of linearlypolarized emission, and Faraday conversion exchanges lin-early and circularly polarized radiation (Sazonov 1969; Ry-bicki & Lightman 1986; Melrose 1997; Shcherbakov 2008).Future analyses with the EHT will be able to further dis-tinguish models in the time and frequency domains (e.g.,Broderick & Loeb 2006a,b; Roelofs et al. 2017; Medeiroset al. 2018). It is therefore timely for us to better under-stand time and frequency dependent effects which may helpus constrain accretion models.One such polarimetric observable is the rotation mea-sure (RM), defined by the change in the EVPA, χ , as afunction of the change of wavelength squared, RM = χ − χ λ − λ , (1)where the subscripts 1 and 2 denote two different wave-lengths. RM probes Faraday rotation, and has been usedin a variety of contexts to infer magnetic field properties(e.g., Zavala & Taylor 2004; Brentjens & de Bruyn 2005;Frick et al. 2011; Pasetto et al. 2018; Agudo et al. 2018).For a source of polarized emission that is entirely behinda Faraday rotating medium, the RM can be written as anintegral of plasma properties along the line-of-sight via RM = 8 . × rad m − (cid:90) observersource f rel (Θ e ) n e − B || G ds pc , (2)where n e is the electron number density, B || is the parallelcomponent of the magnetic field, and f rel is a correction termsuppressing Faraday rotation at relativistic temperatures(Gardner & Whiteoak 1966). At relativistically hot tempera-tures, f rel (Θ e ) ≈ log(Θ e ) / (2Θ e ) , whereas at sub-relativistictemperatures, f rel asymptotes to 1. Here Θ e ≡ k B T e /m e c , k B is the Boltzmann constant, T e is the electron tempera-ture, m e is the electron rest mass, and c is the speed of light(Jones & Odell 1977).As seen from equation 2, the RM directly traces theelectron temperature, number density, and magnetic fieldalong the line of sight. Among the models consistent with theEHT data of M87*, all of these quantities can vary by ordersof magnitude. Note that while RM scales as RM ∝ n e B , thepower emitted by synchrotron emission scales as P ∝ n e B (Rybicki & Lightman 1986). In principle, this could allowRM 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 tar-gets Sgr A* and M87*, as well as a handful of other low-luminosity AGN. For Sgr A*, RM = − × rad m − (Bower et al. 2003; Marrone et al. 2007; Bower et al. 2018)while − × rad m − has been measured a few arc-seconds away (Eatough et al. 2013). For M87*, an upperlimit of | RM | < × rad m − was measured usingthe Submillimeter Array (Kuo et al. 2014). 3C 84 has anRM of − rad m − (Kim et al. 2019) at 43 GHz and × rad m − at 230 GHz (Plambeck et al. 2014). Sim-ilarly, 3C 273 has an RM of × rad m − at 230 GHz(Hovatta et al. 2019). Neither linear polarization nor RMcould be measured for the low-luminosity AGNs M81 andM84, which might imply significant scrambling (Brunthaleret al. 2001; Bower et al. 2017).These measurements have been used to constrain theaccretion rates of EHT targets Sgr A* and M87* by assum-ing simple analytic models describing the accretion flow thatcan be input into equation 2 (Marrone et al. 2006). The ac-cretion rates of Sgr A* and M87* have thus been constrainedto < − M (cid:12) yr (although it may be much lower) (Agol2000; Quataert & Gruzinov 2000b; Marrone et al. 2006) and < . × − M (cid:12) yr (Kuo et al. 2014; Li et al. 2016) respec-tively. However, there are important effects that complicatethe simple scenario implicitly assumed by equation 2, wherea Faraday rotator sits entirely between a source and ourline-of-sight (Broderick & McKinney 2010, see also). Mostimportantly, for BH accretion flows, Faraday rotation andemission 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’strajectory and can modify its polarization properties by par-allel transport alone. This is especially true of emission nearthe photon ring, where null geodesics can pass through theaccretion flow multiple times in different directions, leadingto interesting polarization signatures (Johnson et al. 2019;Himwich et al. 2020). Finally, neither the emission nor theFaraday rotation can be assumed to behave in a spatiallyuniform way, and the magnetic field may switch sign, espe-cially in turbulent SANE discs.Mościbrodzka et al. (2017) produced the first polarizedmodel images of M87* at millimetre wavelengths from ray MNRAS , ?? – ???? (2020) H Internal Faraday Rotation traced GRMHD simulations. Using a set of SANE models,they determined that Faraday rotation can be strong enoughto spatially scramble the polarization from the counter-jet,which must pass through a larger Faraday depth than theforward-jet on the way to the observer. Their analysis alsorevealed a significant inclination dependence of the RM, suchthat even very high accretion rate models could satisfy RMconstraints when viewed face-on. Jiménez-Rosales & Dex-ter (2018) determined that high-accretion rate models aredisfavoured, since this scrambling too strongly depolarizesthe emission. Tsunetoe et al. (2020) began to explore thespin dependence and favoured a = 0 . models, although anexpanded parameter survey is warranted.In this work, we perform a more comprehensive sur-vey of RM for 7 models consistent with EHT5, chosen tobracket the allowed parameter space. We consider variationsas a function of time and inclination, and develop new tech-niques to model Faraday effects in some models. In §2, wedescribe the GRMHD simulations we use as a starting point,the radiative transfer calculations in post-processing, and anovel Taylor expansion model for treating internal Faradayrotation. In §3, we describe our results. This includes anexploration of strong spatial variations in RM, RM distri-bution functions, RM as a measure of accretion rate, incli-nation dependence, the degree of non- λ behaviour amongmodels, and a case study of time variability. We discuss ourresults in §4, and end with a summary and conclusion in §5. We begin with a set of GRMHD simulations that are con-sistent with EHT5. We then use ipole (Mościbrodzka &Gammie 2018) to perform polarized radiative transfer cal-culations, specifying electron properties in this step. Finally,using a first order Taylor expansion, we create a model ofthe polarized image in order to capture frequency dependenteffects and compute rotation measures. In this work, we study 7 models for M87* in the EHTGRMHD simulation library, each of which is consistent withall of the observational constraints considered in EHT5. Theproperties of these simulations, all performed with iharm (Gammie et al. 2003), are listed in Table 1. Both SANE andMAD models are considered, while values of a and R high arechosen to bracket the allowed values in EHT5. These mod-els are described in more detail in EHT5, and the SANE a = +0 . model is included in the GRMHD code compari-son project of Porth et al. (2019).As defined by Mościbrodzka et al. (2016), R high pre-scribes the electron temperature via T i T e = R high β p β p + 11 + β p (3)where T i and T e are the ion and electron temperatures re-spectively, and β p is the ratio of gas to magnetic pressure. T i https://github.com/moscibrodzka/ipole/ Magnetic Field State a R high MAD +0.94 160MAD +0.94 20MAD -0.5 160MAD -0.5 20SANE +0.94 160SANE -0.94 80SANE -0.94 10
Table 1.
Parameters of the 7 models considered in this paper.Each of these models passes all metrics considered in EHT5 andare chosen to bracket the allowed parameter space. is determined by the GRMHD simulation. In the mid-plane,where β p is high, T e → T i /R high , since turbulent plasmamodels reveal that heating preferentially affects ions, whichthen cannot efficiently transfer energy to electrons (Reeset al. 1982; Narayan & Yi 1995b; Quataert & Gruzinov 1999;Howes 2010; Kawazura et al. 2019). Consequently, increas-ing R high has the effect of decreasing the emission from andincreasing the Faraday rotation within the mid-plane. As aresult of decreasing mid-plane emission, the accretion ratemust also be scaled upwards to obtain the correct total in-tensity for M87*. As we shall show, both of these effectshave important implications for the RM.Note that among these 7 models, there are only 4 uniqueGRMHD simulations, since R high only affects the radiativetransfer in post-processing. In each of these simulations, theangular momentum of the disc is either perfectly aligned (de-noted by positive spin) or anti-aligned (denoted by negativespin) with that of the BH, although misaligned discs remainan active area of research (e.g., Fragile et al. 2007; Liska et al.2020; Chatterjee et al. 2020). MAD simulations are run witha × × grid with a maximum radius of GM • /c ,while SANE simulations are run with a × × gridand a maximum radius of GM • /c . However, we findthat these models only exhibit inflow equilibrium within aradius of approximately GM • /c . Throughout this work,we restrict the integration of our radiative transfer equationsto within this radius. Fortunately, this limitation actuallyhas negligible effect on our results for face-on inclinations( i (cid:46) ◦ ), as we explore in Appendix C. This is because thefunnel region is largely evacuated in these simulations anddoes not contribute to Faraday rotation. However, restrict-ing our calculations to (cid:54) GM • /c may lead us to underes-timate the total Faraday rotation at inclinations (cid:38) ◦ , dueto material that may exist in more distant, unequilibratedregions of the simulations.We create polarized ray-traced images using ipole (Mościbrodzka & Gammie 2018), which first solves the nullgeodesic equation backwards from the image plane, then in-tegrates forward the radiative transfer equations for the 4Stokes parameters, { I, Q, U, V } . Here, I is the total inten-sity, (cid:112) Q + U and arg( Q + iU ) are the linearly polarizedintensity and EVPA respectively, and V is the circularly po-larized intensity. The radiative transfer equations accountfor synchrotron emission, synchrotron self-absorption, Fara-day rotation, and Faraday conversion. Radiative transfer co-efficients follow Dexter (2016) for a thermal electron distri-bution function, with a slight modification to ρ ν,V , the coef-ficient responsible for Faraday rotation. As also discussed in MNRAS , ?? – ?? (2020) Ricarte et al. y [ a s ] MAD, a=+0.94, R high =160 MAD, a=+0.94, R high =20 MAD, a=-0.5, R high =160 MAD, a=-0.5, R high =20
30 0 30 x [ as] y [ a s ] SANE, a=+0.94, R high =160
30 0 30 x [ as] SANE, a=-0.94, R high =10
30 0 30 x [ as] SANE, a=-0.94, R high =80 I / I . Figure 1.
Total intensity images of the 7 models considered in this work, taken at the final snapshot of each GRMHD simulation, whichoccurs at time t = 10000 GM • /c . Models are tilted by 17 ◦ from face-on towards the top of these images and in subsequent figures. Tohelp visualise low surface brightness features, intensity is scaled with respect to the intensity of the pixel in the 99.7th percentile, I . ,saturating 0.3 per cent of the pixels. Dexter et al. (2020), minor modifications are needed to en-sure continuous and accurate behaviour at low temperatureand frequency. Following Shcherbakov (2008), we set ρ ν,V = 2 ne ν B m e cν K (Θ − e ) K (Θ − e ) cos(Θ e ) g ( X ) , (4)where K and K are modified Bessel functions of the sec-ond kind, e is the electron charge, m e is the electron mass, ν B = eB/ πm e c , X = (cid:104) √ − νν c (cid:105) − / , and g ( X ) =1 − .
11 ln(1 + 0 . X ) for cyclotron frequency ν c .For each model, we create images for 11 snapshots span-ning the last quarter of the corresponding GRMHD run,corresponding to times t/ ( GM • /c ) ∈ [7500 , , a du-ration of 880 days for M87*. For each snapshot, we study5 inclinations, i ∈ { ◦ , ◦ , ◦ , ◦ , ◦ } . Then, for eachsnapshot and inclination, we create 6 polarized images. Wesample 3 frequencies, 226.999, 227.000 and 227.001 GHz, toconstruct a frequency-dependent model of the image, de-tailed in the following section. We find that we must adoptextremely small differences in frequency in order to resolvethe RM on geodesics that have high Faraday depth. Other-wise, the RM in some geodesics can be underestimated dueto the “n π degeneracy:” the EVPA may rotate so rapidlywithin that multiple rotations are missed. Then, for each of Throughout this paper, for models with positive spin, we actu-ally compute images for i = 180 ◦ − i written , as in EHT5. these 3 frequencies, we create 2 separate images that onlyinclude emission from either the positive or negative z do-mains, where the z -axis is oriented parallel to the black holespin (and perpendicular to the disc in these models).Separate images including only the top and bottomhalves of the emission are helpful for modelling Faraday ef-fects that are often significantly different between the frontand back sides of the emitting region. In these models, emis-sion behind the mid-plane often experiences orders of mag-nitude 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 athigh inclinations, since some emission is lensed to the oppo-site side of the image. A complete image is made by summingtogether the individual images made with only the near- andfar-side emission. Note that when emission from only oneside is included, absorption and Faraday effects from bothsides remain included.For all images, we adopt a black hole mass of . × M (cid:12) and a distance of 16.9 Mpc (Gebhardt et al. 2011),for consistency with EHT5. Gas densities are scaled suchthat the average image produces an intensity of 0.5 Jy atan inclination of 17 ◦ , which is the most likely inclinationat which we are viewing M87* based on its larger scale jet(Walker et al. 2018). The same gas density scalings are usedfor models at different inclinations, although their averagefluxes can depart from 0.5 Jy. We find that the total inten-sities of images created at an inclination of ◦ tend to beapproximately a factor of 2 larger than those created with MNRAS , ?? – ????
11 ln(1 + 0 . X ) for cyclotron frequency ν c .For each model, we create images for 11 snapshots span-ning the last quarter of the corresponding GRMHD run,corresponding to times t/ ( GM • /c ) ∈ [7500 , , a du-ration of 880 days for M87*. For each snapshot, we study5 inclinations, i ∈ { ◦ , ◦ , ◦ , ◦ , ◦ } . Then, for eachsnapshot and inclination, we create 6 polarized images. Wesample 3 frequencies, 226.999, 227.000 and 227.001 GHz, toconstruct a frequency-dependent model of the image, de-tailed in the following section. We find that we must adoptextremely small differences in frequency in order to resolvethe RM on geodesics that have high Faraday depth. Other-wise, the RM in some geodesics can be underestimated dueto the “n π degeneracy:” the EVPA may rotate so rapidlywithin that multiple rotations are missed. Then, for each of Throughout this paper, for models with positive spin, we actu-ally compute images for i = 180 ◦ − i written , as in EHT5. these 3 frequencies, we create 2 separate images that onlyinclude emission from either the positive or negative z do-mains, where the z -axis is oriented parallel to the black holespin (and perpendicular to the disc in these models).Separate images including only the top and bottomhalves of the emission are helpful for modelling Faraday ef-fects that are often significantly different between the frontand back sides of the emitting region. In these models, emis-sion behind the mid-plane often experiences orders of mag-nitude 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 athigh inclinations, since some emission is lensed to the oppo-site side of the image. A complete image is made by summingtogether the individual images made with only the near- andfar-side emission. Note that when emission from only oneside is included, absorption and Faraday effects from bothsides remain included.For all images, we adopt a black hole mass of . × M (cid:12) and a distance of 16.9 Mpc (Gebhardt et al. 2011),for consistency with EHT5. Gas densities are scaled suchthat the average image produces an intensity of 0.5 Jy atan inclination of 17 ◦ , which is the most likely inclinationat which we are viewing M87* based on its larger scale jet(Walker et al. 2018). The same gas density scalings are usedfor models at different inclinations, although their averagefluxes can depart from 0.5 Jy. We find that the total inten-sities of images created at an inclination of ◦ tend to beapproximately a factor of 2 larger than those created with MNRAS , ?? – ???? (2020) H Internal Faraday Rotation an inclination of ◦ . Each image is created at . µ as pixelresolution, with a 160 µ as field of view, a factor of 2 finerangular resolution than employed in EHT5.To summarise, we consider 7 models that are consis-tent with EHT observations of M87*. For each model, wecreate images for 11 snapshots, 5 inclinations, 3 frequencies,and the 2 sides of the accretion flow. We also create addi-tional images to better resolve frequency-dependent (§2.2)and time-dependent (§3.7) phenomena. In Figure 1, we plottotal intensity images of the 7 models during their finalsnapshot at ◦ inclination. To better visualise low surfacebrightness features, we intentionally saturate 0.3 per centof the pixels in this and all subsequent visualisations. In allof the images presented in this work, the forward-jet pointsstraight up, and the material on the right side of the imageis moving towards the observer. Here, we introduce a modelling formalism for the polarizedimage as a function of frequency. For a given image snapshot,we create a two-zone (front and back half) model of theradiation in each pixel. Our model is constructed as follows: • Treating the front and back halves of the emittingregion 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 et al.2012). These variables have well-behaved Taylor series ex-pansions, and are further discussed in Appendix A. • Using images at three frequencies, we then compute thefirst and second derivatives of the spherical Stokes parame-ters. Each derivative is computed in the variable that bestapproximates expected physical dependencies: I and N aredifferentiated in ν , φ is differentiated in λ , and ψ is differ-entiated in λ . • Using the first derivatives, we construct Taylor expan-sions to first order in each of the spherical Stokes parameters.These can be expressed I ( ν ) ≈ I ( ν ) + dIdν ( ν − ν ) , (5) N ( ν ) ≈ N ( ν ) + dNdν ( ν − ν ) , (6) φ ( ν ) ≈ φ ( ν ) + dφdλ ( λ − λ ) ≈ φ ( ν ) − c ν dφdλ ( ν − ν ) , (7) ψ ( ν ) ≈ ψ ( ν ) + dψdλ ( λ − λ ) ≈ ψ ( ν ) − c ν dψdλ ( ν − ν ) . (8) • In order to approximate bandwidth-integration, we in-tegrate 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 I BW = 1 ν − ν (cid:90) ν ν I ( ν ) dν, (9) Q BW = 1 ν − ν (cid:90) ν ν N ( ν ) cos[ φ ( ν )] sin[ ψ ( ν )] dν, (10) U BW = 1 ν − ν (cid:90) ν ν N ( ν ) sin[ φ ( ν )] sin[ ψ ( ν )] dν, (11) V BW = 1 ν − ν (cid:90) ν ν N ( ν ) cos[ ψ ( ν )] dν. (12)Analytic solutions to these integrals are provided in Ap-pendix B. • Finally, the front and back emission halves are summedto produce the complete image.In practice, dφ/dλ , which encapsulates Faraday rotation,is the most important frequency-dependent effect to takeinto account. Faraday conversion, which is encapsulated in dψ/dλ , is not as strong as Faraday rotation in these mod-els. In some of our models, our split into the front and backhalves is necessary for capturing the Faraday depolarizationof back half of the emission (due to the Faraday thick mid-plane) while preserving the polarization of the front half. Aswe will later explore in Figure 6, some pixels may simulta-neously have emission from the forward-jet with an RM of ∼ rad m − , and emission from the counter-jet with anRM of ∼ rad m − . Consequently, φ ( λ ) for the pixel ishighly non-linear, but can be approximated as the sum oftwo emitting regions with distinct RMs.In some pixels, especially within the photon ring, wecompute unreasonably high values of | dN/dν | within ourTaylor expansion. This is due to multiple emission compo-nents in the same pixel Faraday rotating at different rates,which can periodically increase and decrease the total lin-early polarized intensity in the pixel depending on the rela-tive phases of these components. We use d Ndν to help identifysuch problematic pixels. Let ∆ ν = ν − ν be the difference infrequency space between some frequency ν at which Stokesparameters are being evaluated and the frequency at whichderivatives have been computed, ν = 228 GHz . Treating d Ndν ∆ ν as an upper limit on the error of N ( ν ) , we freezethe value of N in pixels where | dNdν ∆ ν | < | d Ndν ∆ ν | and | d ln( N ) dν | > | ν | . That is, we freeze N if its first derivativeimplies that it would grow or shrink by more than a factor of e across ∆ ν and the absolute error on the change of N maybe larger than the change of N itself. We apply an identicalcondition on the derivative of I , which almost always affectspixels also affected by the condition on Stokes N .We examine to what extent this criterion is applied inour final snapshot images, and find that it affects only a mi-nority of pixels. Among the pixels which amount to 99 percent of the image integrated I , we find that at most 5 percent of the pixels are affected by this correction, typically inthe photon ring. This largest fraction occurs in the SANE, a = − . , R high = 80 model, which as we shall show con-tains the emission with the largest Faraday depths. In someof the other models, none of the pixels are affected by thiscriterion.In Figure 2, we demonstrate the efficacy of our methodby plotting the linearly polarized intensity ( (cid:112) Q + U ) forthe SANE, a = +0 . , R high = 160 model. The top left MNRAS , ?? – ?? (2020) Ricarte et al. y [ a s ] One Frequency Averaged Images
40 20 0 20 40 x [ as] y [ a s ] Taylor Approximated
40 20 0 20 40 x [ as]
Average - Taylor Sp e c i f i c I n t e n s i t y [ J y p x ] Sp e c i f i c I n t e n s i t y [ J y p x ] Sp e c i f i c I n t e n s i t y [ J y p x ] Sp e c i f i c I n t e n s i t y [ J y p x ] Figure 2.
Images of total linearly polarized intensity ( (cid:112) Q + U ) for a case where bandwidth depolarization is important, the SANE, a = +0 . , 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, weplot the average of 255 images across a 2 GHz bandwidth. As further explored in Figure 3, bandwidth depolarization suppresses thecontribution of the counter-jet. In the bottom left, we plot the result of our Taylor expansion model based on 6 images. Its residual withrespect to the 255 averaged images is shown in the bottom right. The Taylor expansion model successfully depolarizes the counter-jetwithout depolarizing the forward-jet. Since we model the image with two zones, the largest discrepancies occur in and within the photonring, where geodesics cross the mid-plane multiple times. panel shows the result of a single-frequency calculation at228 GHz, while the top right image shows the result of av-eraging 255 images across a 4 GHZ bandwidth between 226and 230 GHz. In the bandwidth-averaged image, much of thelarge-scale emission has been suppressed and two streaks ofemission in the upper right of the image have become moreprominent. The total spatially unresolved linear polariza-tion fraction has dropped by a factor of 2.4. Although thecounter-jet dominates the linearly polarized emission in thesingle-frequency image due to lensing (Dexter et al. 2012),most of that emission is scrambled away due to the largeFaraday depth in the mid-plane (Mościbrodzka et al. 2017).The bottom left panel is the result of our Taylor ap-proximated model based on 6 images (3 frequencies, 2 sides), which successfully captures the depolarization of thecounter-jet, but not the forward-jet. In the bottom rightpanel, we subtract this model from the averaged images andplot the residual. Weighting by the linear polarized intensityof the properly averaged image, total linear polarization isrecovered on average with an absolute error of . µ Jy anda relative error of . . In contrast, the single-frequency im-age has an average absolute error . µ Jy and relative errorof . . In this case, the single-frequency image also over-estimates the total linear polarization somewhat, even if thesource is spatially unresolved. When Stokes parameters aresummed across the entire image, the properly averaged im-age has a linear polarization fraction of . × − , the Tay-lor approximated image has a linear polarization fraction of MNRAS , ?? – ????
Images of total linearly polarized intensity ( (cid:112) Q + U ) for a case where bandwidth depolarization is important, the SANE, a = +0 . , 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, weplot the average of 255 images across a 2 GHz bandwidth. As further explored in Figure 3, bandwidth depolarization suppresses thecontribution of the counter-jet. In the bottom left, we plot the result of our Taylor expansion model based on 6 images. Its residual withrespect to the 255 averaged images is shown in the bottom right. The Taylor expansion model successfully depolarizes the counter-jetwithout depolarizing the forward-jet. Since we model the image with two zones, the largest discrepancies occur in and within the photonring, where geodesics cross the mid-plane multiple times. panel shows the result of a single-frequency calculation at228 GHz, while the top right image shows the result of av-eraging 255 images across a 4 GHZ bandwidth between 226and 230 GHz. In the bandwidth-averaged image, much of thelarge-scale emission has been suppressed and two streaks ofemission in the upper right of the image have become moreprominent. The total spatially unresolved linear polariza-tion fraction has dropped by a factor of 2.4. Although thecounter-jet dominates the linearly polarized emission in thesingle-frequency image due to lensing (Dexter et al. 2012),most of that emission is scrambled away due to the largeFaraday depth in the mid-plane (Mościbrodzka et al. 2017).The bottom left panel is the result of our Taylor ap-proximated model based on 6 images (3 frequencies, 2 sides), which successfully captures the depolarization of thecounter-jet, but not the forward-jet. In the bottom rightpanel, we subtract this model from the averaged images andplot the residual. Weighting by the linear polarized intensityof the properly averaged image, total linear polarization isrecovered on average with an absolute error of . µ Jy anda relative error of . . In contrast, the single-frequency im-age has an average absolute error . µ Jy and relative errorof . . In this case, the single-frequency image also over-estimates the total linear polarization somewhat, even if thesource is spatially unresolved. When Stokes parameters aresummed across the entire image, the properly averaged im-age has a linear polarization fraction of . × − , the Tay-lor approximated image has a linear polarization fraction of MNRAS , ?? – ???? (2020) H Internal Faraday Rotation . × − , and the single-frequency image has a linear po-larization fraction of . × − . A model containing at leasttwo separate zones is necessary in order to reproduce the de-polarization of the counter-jet without also depolarizing theforward-jet. The remaining residuals require more than tworegions to fully capture the complexity of the Faraday ro-tating structure along the line-of-sight. Notice that the mostsignificant errors occur within the photon ring or its interior,where geodesics pass through the mid-plane multiple times.In Figure 3, we further decompose the single-frequency228 GHz image into emission from its forward-jet and fromcounter-jet components. In the lower panel, we examineframe-invariant radiative transfer coefficients in the pixelmarked by a blue circle in these images. In this pixel, po-larized emission is emitted roughly equally by forward-jetand counter-jet components. However, the counter-jet emis-sion must pass through the enormous Faraday depth of thecold mid-plane, with RM > rad m − in some regions.The polarized intensity of emission passing through materialwith a rotation measure of RM and a bandwidth of ∆ ν issuppressed by a factor f BW given by f BW = sinc (cid:18) c ν RM∆ ν (cid:19) , (13)where sinc( x ) = sin( x ) /x , c is the speed of light, and ν isthe central frequency, assuming narrow fractional bandwidthand uniform sampling in ν . We define the critical rotationmeasure RM crit via
12 = sinc (cid:18) c ν RM crit BW (cid:19) . (14)That is, RM crit is the minimum RM required to suppress lin-ear polarization by a factor of 2. So far, EHT observationshave been performed using a central frequency of 228 GHzand a bandwidth of 4 GHz. Using these EHT values, we findthat RM crit = 3 × rad m − . Since rad m − (cid:29) RM crit ,the counter-jet’s polarization is significantly suppressed. No-tice that the photon ring emission from the forward-jet isalso suppressed, since this emission must also pass throughthe cold mid-plane. Jiménez-Rosales & Dexter (2018) deter-mine that strong Faraday effects also scramble the imageon a pixel-by-pixel basis, resulting in beam depolarization.This spatial decoherence helps compensate for bandwidthdepolarization when blurred images are constructed. In §2.2, we compute Stokes parameters and their derivativesat a central frequency of 228 GHz and a small bandwidthof 2 MHz. From these, we can directly compute the RM ineach pixel at 228 GHz via RM ≡ dχdλ = ddλ
12 arctan (cid:18) UQ (cid:19) = 12 U (cid:48) Q − Q (cid:48) UQ + U , (15)which follows from a straightforward application of the chain Since ν varies along the geodesic, it is useful to define the frame-invariant quantities ρ V = νρ ν,V and j Q = j ν,Q /ν (for additionaldiscussion, see Mościbrodzka & Gammie 2018). rule; here the (cid:48) symbols denote d/dλ . When we discuss theRM of individual pixels, this is how RM is computed. Usinga small band of 2 MHz, the maximum measurable RM inan individual pixel is | RM | max = π/ ∆ λ ≈ πν / (2 c ∆ ν ) =1 . × rad m − .However, when discussing the RM for spatially unre-solved measurements in subsequent sections, it is importantto recognise the complicated evolution of χ ( λ ) that resultsfrom the complex RM structure we uncover in §3.1. Whenassigning a single value of the RM to an entire spatially un-resolved image, we use the Taylor expansion methodologydeveloped in §2.2 to approximate 16 polarized images, eachintegrated over a bandwidth of . GHz, equally spaced in λ space between 226 and 230 GHz to emulate EHT obser-vations. These images are used to compute χ ( λ ) , and ∆ χ is computed from the endpoints of the band. We correct forphase wraps by adding or subtracting π to χ ( λ ) as neces-sary to obtain the correct sign of dχ/dλ based on Equa-tion 15. Using 16 bands of width 0.25 GHz, the maximummeasurable RM is | RM | max = π/ ∆ λ ≈ πν / (2 c ∆ ν ) =8 . × rad m − . In this section, we study the RM of each of our models withinsingle simulation snapshots. We find significant spatial vari-ation across the image due to inhomogeneities in the accre-tion flow on event horizon scales. At different locations, RMcan vary by orders of magnitude and even flip sign. As a re-sult, for spatially unresolved polarized measurements, someof these models may exhibit highly non-linear χ ( λ ) , makingthem difficult to characterise with a single RM.In Figure 4, we visualise the spatially resolved RMstructure of the images presented in Figure 1. In this fig-ure, the brightness of each pixel scales with the linearly po-larized specific intensity ( (cid:112) Q + U ) , while the colorationencodes the RM. As in Figure 1, 0.3 per cent of the pixelsare saturated to more clearly display low surface brightnessstructures. The colour-scale is normalised separately for eachmodel, spanning ± the 90th-percentile of | RM | of the pixelswhich have at least 50 per cent of the maximum linear polar-ized intensity plotted. Red regions have positive RM, whileblue regions have negative RM. The top two rows depict theRM with ∆ ν = 2 MHz, where our Taylor expansion modelis constructed. The bottom two rows plot the RM across a 4GHz band, a more realistic bandwidth, within which band-width depolarization becomes important for some models. 4GHz images are plotted with the same brightness and RMscale as their 2 MHz counterparts.As shown in this figure, complex spatial variation andfrequent sign-flips are a generic feature of these RM maps.This behaviour is not surprising in SANE models, whichare characterised by weaker, disordered magnetic fields, butis less expected in MAD models, which are characterisedby strong poloidal fields. In one suggestive snapshot, weconfirm that these RM sign-flips are due to sign-flips inthe magnetic field with respect to the geodesic. Figure 5plots the intensity-weighted Faraday depth in each pixel, τ F = (cid:82) ρ V ds , for a snapshot of the MAD, a = 0 . , MNRAS , ?? – ?? (2020) Ricarte et al.
30 20 10 0 10 20 30 x [ as] y [ a s ] Q + U Forward-jet
30 20 10 0 10 20 30 x [ as]Counter-jet Sp e c i f i c I n t e n s i t y [ J y p x ] z [ GM / c ] j Q [ e r g s c m H z s t e r ] ( P o l a r i z e d E m i ss i o n ) To ObserverForward-Jet Counter-JetCold Midplane V [ r a d m s ] ( F a r a d a y R o t a t i o n ) Figure 3.
Linearly polarized intensity ( (cid:112) Q + U ) for the SANE, a = +0 . , R high = 160 model at 227 GHz (0 bandwidth), brokeninto forward-jet (Top Left) and counter-jet (Top Right) emission components, which appear very different due to lensing. For the geodesicmarked by a blue circle, we plot in the lower panel the frame-invariant versions of the radiative transfer coefficients that correspond topolarized emissivity ( j Q , red) and Faraday rotation ( ρ V , blue) . For ρ V , dotted lines correspond to negative values, and there are manysign flips due to the turbulent nature of a SANE disc. Emission from the counter-jet passes through the cold mid-plane on the way tothe observer, which produces a very large RM. Thus, when properly averaging over bandwidth, the counter-jet emission is suppressedby bandwidth depolarization. (See Fig. 2, top right panel.) R high = 20 model. Here, ρ V is the (frame-invariant) radia-tive transfer coefficient responsible for Faraday rotation, and s is the affine parameter describing the geodesic. The signof this quantity, shown to exhibit both positive to negativevalues, directly encodes the direction of the magnetic fieldwith respect to the photon trajectory. RM sign flips havebeen predicted by earlier MHD simulations without GR, butonly at large inclinations (Sharma et al. 2007).By performing a 3D visualisation of this snapshot with VisIt (Childs et al. 2012), we find that sign flips in themagnetic field occur in its tangential and radial componentswhen crossing the disc mid-plane. Near the event horizon, field lines on the northern hemisphere point west, while fieldlines on the southern hemisphere point east, although theypoint north overall (in the positive z-direction). This is anatural consequence of the tangential stretching of verticalfield lines as they are dragged into the BH by accreting ma-terial, as well as frame dragging (see e.g., Contopoulos et al.2009; Gabuzda 2018, for helpful schematics). Since this signflip occurs when crossing the mid-plane, accreting streamsof gas that straddle the mid-plane can exhibit streaks ofpositive RM adjacent to streaks of negative RM.Returning to Figure 4 and comparing the 2 MHz band-width visualisations to the 4 GHz bandwidth visualisations,
MNRAS , ?? – ????
MNRAS , ?? – ???? (2020) H Internal Faraday Rotation y [ a s ] MAD, a=+0.94, R high =160
MAD, a=+0.94, R high =20 MAD, a=-0.5, R high =160 MAD, a=-0.5, R high =20
30 0 30 x [ as] y [ a s ] SANE, a=+0.94, R high =160
30 0 30 x [ as] SANE, a=-0.94, R high =10
30 0 30 x [ as] SANE, a=-0.94, R high =80 y [ a s ] MAD, a=+0.94, R high =160
MAD, a=+0.94, R high =20 MAD, a=-0.5, R high =160 MAD, a=-0.5, R high =20
30 0 30 x [ as] y [ a s ] SANE, a=+0.94, R high =160
30 0 30 x [ as] SANE, a=-0.94, R high =10
30 0 30 x [ as] SANE, a=-0.94, R high =80 -90th percentile+90th percentile R M [ r a d m ] Figure 4.
Visualisation of the RM structure for the snapshots plotted in Figure 1. The brightness of each pixel is proportional to thelinearly polarized intensity in the pixel (again saturating 0.3 per cent of the pixels), while the coloration encodes the RM. The colourscale is normalised separately for each model, spanning ± the 90th-percentile of | RM | of the pixels which have at least 50 per cent ofthe maximum linear polarized intensity. Red regions have positive RM, while blue regions have negative RM. The top two rows are notbandwidth 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. SANE models are more strongly affected by bandwidth de-polarization than MAD models. In each of the SANE mod-els, counter-jet polarized emission is especially suppressed.This more strongly changes the morphology of the progradeSANE than the retrograde SANEs, because the image mor-phologies of the counter-jet and forward-jet in the retrogradeSANEs are more similar.Closely examining the 2 MHz counter-rotating SANE models, one may notice that the linearly polarized inten-sity appears to vary strongly among adjacent pixels, caus-ing the appearance of “static” across the image. This effectis an artefact of single-frequency radiative transfer calcula-tions, an instance of Faraday rotation randomizing not onlythe 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
MNRAS , ?? – ?? (2020) Ricarte et al.
40 20 0 20 40 x [ as] y [ a s ] F [ r a d ] Figure 5.
Intensity-weighted Faraday depth in one snapshot ofthe MAD, a = 0 . , R high = 20 model, revealing clear sign flipsin the magnetic field parallel to the line-of-sight. The colour ofeach pixel encodes the total Faraday depth, while the brightnessis proportional to the linear polarized intensity. The sign of theFaraday depth directly encodes the direction of the magnetic fieldparallel to the geodesic. Subsequent 3D visualisation of this snap-shot reveals that these sign flips occur in the tangential magneticfield in the plane of the disc. components, neither of which exhibit this “static” if plot-ted individually. The phase of the counter-jet emission iseffectively randomized by the enormous Faraday depth inthe mid-plane. Depending on the relative phase, the rotatedcounter-jet polarized emission may sometimes cancel withthe forward-jet polarized emission. This effect would not oc-cur in real observations integrated over a finite bandwidth. In these models, polarized emission exhibits a wide rangeof Faraday depths, and the front and back halves of theemitting regions can differ by many orders of magnitude. InFigure 6, we plot the distribution functions of log | RM | among the pixels of each model during their final snap-shot, weighted by each pixel’s linear polarized intensity. Thisquantity, which we denote as d | p | /d log | RM | , is closely re-lated to F ( φ ) in rotation measure synthesis theory, the com-plex polarized surface brightness per unit Faraday depth(Burn 1966; Brentjens & de Bruyn 2005). Unlike F ( φ ) , wetake absolute values, normalise with respect to the totalemitted polarized surface brightness, and adopt logarithmi-cally spaced bins.In each pixel, | RM | is computed directly from the gradi-ent at 228 GHz across ∆ ν = 2 MHz . These distributions arecomputed for each of the 11 snapshots that we study, and thefilled regions span the range permitted by all of these snap-shots. Blue regions only include emission originating fromthe front half of the emitting region, while red regions onlyinclude 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 band-width depolarized. The distributions of front half of theemission region can be significantly displaced from that ofthe back half, especially in SANE models. Again, this is dueto the large Faraday depth of the sub-relativistic mid-planein these models. In some models, including all of the SANEs,the forward-jet distributions exhibit two distinct peaks. Thepeak at higher | RM | is due to photon ring orbits, which passthrough this mid-plane. Note the extreme difference in | RM | between front and back components in the jet-dominatedretrograde SANE models. These models may exhibit muchlower spatially unresolved | RM | than would be expected byintegrating their Faraday depths across geodesics, since thisFaraday depth mainly only affects the counter-jet. In §3.1, we demonstrate that GRMHD models exhibit richspatial structure, whereby the RM can vary by many or-ders of magnitude and flip sign. As a consequence, χ ( λ ) rotates at very different rates at different locations withinthe image, which may result in clear departures from a λ law if these structures are not spatially resolved. In Fig-ure 7, we investigate this non-linearity by plotting χ ( λ ) forthe final snapshots of our 7 models. In blue, we laboriouslycompute 255 individual images across the 4 GHz bandwidth,then average Stokes parameters within 16 smaller 0.25 GHzbands to estimate χ ( λ ) . In red, we instead use the Taylorexpansion model based on 6 images and analytic integralsdescribed in §2.2 to estimate χ ( λ ) .All models exhibit significantly non- λ behaviour evenwithin 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 16MHz. This is because as shown in Figure 6, a signifi-cant amount of the intensity has an associated | RM | of ≈ rad m − . The MAD R high = 160 models exhibit rel-atively mild non-linearity, since | RM | only just approaches | RM crit | in Figure 6.Recall that we use ∆ χ/ ∆ λ across the bandwidth fromour Taylor expansion model to assign spatially unresolvedRMs to images. This figure reveals some of this model’s lim-itations. The model poorly reproduces χ ( λ ) for the MAD, a = +0 . , R high = 160 model, possibly due to 228 GHzbeing a local extremum of χ ( λ ) where the first derivativeis small. Interestingly, the retrograde SANE Taylor expan-sion models appears to broadly follow the structure of thetrue χ ( λ ) , but with a vertical offset. This is due to incor-rect evolution of emission superposed on top of the photonring. 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. Inthe SANE, a = − . , R high = 80 model, this superposedcomponent is immediately bandwidth depolarized and sub-tracted, leading to the offset. This effect is more delayed as afunction of ∆ λ in the SANE, a = − . , R high = 10 model,which by construction has a warmer mid-plane and thereforeless Faraday rotation. Fortunately, this effect is symmetricabout the Taylor expansion point at 228 GHz and we can MNRAS , ?? – ????
Intensity-weighted Faraday depth in one snapshot ofthe MAD, a = 0 . , R high = 20 model, revealing clear sign flipsin the magnetic field parallel to the line-of-sight. The colour ofeach pixel encodes the total Faraday depth, while the brightnessis proportional to the linear polarized intensity. The sign of theFaraday depth directly encodes the direction of the magnetic fieldparallel to the geodesic. Subsequent 3D visualisation of this snap-shot reveals that these sign flips occur in the tangential magneticfield in the plane of the disc. components, neither of which exhibit this “static” if plot-ted individually. The phase of the counter-jet emission iseffectively randomized by the enormous Faraday depth inthe mid-plane. Depending on the relative phase, the rotatedcounter-jet polarized emission may sometimes cancel withthe forward-jet polarized emission. This effect would not oc-cur in real observations integrated over a finite bandwidth. In these models, polarized emission exhibits a wide rangeof Faraday depths, and the front and back halves of theemitting regions can differ by many orders of magnitude. InFigure 6, we plot the distribution functions of log | RM | among the pixels of each model during their final snap-shot, weighted by each pixel’s linear polarized intensity. Thisquantity, which we denote as d | p | /d log | RM | , is closely re-lated to F ( φ ) in rotation measure synthesis theory, the com-plex polarized surface brightness per unit Faraday depth(Burn 1966; Brentjens & de Bruyn 2005). Unlike F ( φ ) , wetake absolute values, normalise with respect to the totalemitted polarized surface brightness, and adopt logarithmi-cally spaced bins.In each pixel, | RM | is computed directly from the gradi-ent at 228 GHz across ∆ ν = 2 MHz . These distributions arecomputed for each of the 11 snapshots that we study, and thefilled regions span the range permitted by all of these snap-shots. Blue regions only include emission originating fromthe front half of the emitting region, while red regions onlyinclude 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 band-width depolarized. The distributions of front half of theemission region can be significantly displaced from that ofthe back half, especially in SANE models. Again, this is dueto the large Faraday depth of the sub-relativistic mid-planein these models. In some models, including all of the SANEs,the forward-jet distributions exhibit two distinct peaks. Thepeak at higher | RM | is due to photon ring orbits, which passthrough this mid-plane. Note the extreme difference in | RM | between front and back components in the jet-dominatedretrograde SANE models. These models may exhibit muchlower spatially unresolved | RM | than would be expected byintegrating their Faraday depths across geodesics, since thisFaraday depth mainly only affects the counter-jet. In §3.1, we demonstrate that GRMHD models exhibit richspatial structure, whereby the RM can vary by many or-ders of magnitude and flip sign. As a consequence, χ ( λ ) rotates at very different rates at different locations withinthe image, which may result in clear departures from a λ law if these structures are not spatially resolved. In Fig-ure 7, we investigate this non-linearity by plotting χ ( λ ) forthe final snapshots of our 7 models. In blue, we laboriouslycompute 255 individual images across the 4 GHz bandwidth,then average Stokes parameters within 16 smaller 0.25 GHzbands to estimate χ ( λ ) . In red, we instead use the Taylorexpansion model based on 6 images and analytic integralsdescribed in §2.2 to estimate χ ( λ ) .All models exhibit significantly non- λ behaviour evenwithin 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 16MHz. This is because as shown in Figure 6, a signifi-cant amount of the intensity has an associated | RM | of ≈ rad m − . The MAD R high = 160 models exhibit rel-atively mild non-linearity, since | RM | only just approaches | RM crit | in Figure 6.Recall that we use ∆ χ/ ∆ λ across the bandwidth fromour Taylor expansion model to assign spatially unresolvedRMs to images. This figure reveals some of this model’s lim-itations. The model poorly reproduces χ ( λ ) for the MAD, a = +0 . , R high = 160 model, possibly due to 228 GHzbeing a local extremum of χ ( λ ) where the first derivativeis small. Interestingly, the retrograde SANE Taylor expan-sion models appears to broadly follow the structure of thetrue χ ( λ ) , but with a vertical offset. This is due to incor-rect evolution of emission superposed on top of the photonring. 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. Inthe SANE, a = − . , R high = 80 model, this superposedcomponent is immediately bandwidth depolarized and sub-tracted, leading to the offset. This effect is more delayed as afunction of ∆ λ in the SANE, a = − . , R high = 10 model,which by construction has a warmer mid-plane and thereforeless Faraday rotation. Fortunately, this effect is symmetricabout the Taylor expansion point at 228 GHz and we can MNRAS , ?? – ???? (2020) H Internal Faraday Rotation d | p |/ d l o g | R M | MAD, a=+0.94, R=160 MAD, a=+0.94, R=20 MAD, a=-0.5, R=160 MAD, a=-0.5, R=203 5 7 9 11log |RM| [rad m ]0.00.51.0 d | p |/ d l o g | R M | SANE, a=+0.94, R=160 3 5 7 9 11log |RM| [rad m ]SANE, a=-0.94, R=10 3 5 7 9 11log |RM| [rad m ]SANE, a=-0.94, R=80 log RM crit (4 GHz)Front HalfBack Half Figure 6.
Distribution functions of log | RM | for our models, weighted by the linear polarized intensity of each pixel. Emission originatingfrom front side of the emitting region is coloured blue, while emission originating from the back side is coloured red. The relatively coldmid-plane is the dominant source of Faraday rotation in these models, which can result in significant offsets in these distributions betweenthe two halves, especially in SANE models. Forward-jet emission in SANE models exhibits a second peak at high log | RM | due tophoton ring geodesics, which do pass through the mid-plane. still recover the spatially unresolved RM from χ at the endpoints of the band. RM is often used to approximate the accretion rate ˙ M • ,based on simple analytic models (Marrone et al. 2006).These models are based on advection or convection domi-nated accretion flows (Narayan & Yi 1994; Narayan et al.2000; Quataert & Gruzinov 2000a) and make many simplify-ing assumptions. These include spherical symmetry, equipar-tition of energy, and rather arbitrary inner and outer radiito truncate the model. Adapted from Marrone et al. (2006), RM = (cid:0) . × rad m − (cid:1) (cid:104) − ( r out /r in ) − (3 β − / ) (cid:105) × (cid:18) M • . × M (cid:12) (cid:19) − (cid:18) β − (cid:19) r − / ˙ M / • , (16)where we have corrected the exponent of r in , as noted byMacquart et al. (2006). Here, r out and r in are radii usedto truncate the model in units of Schwarzschild radii, and β ∈ [0 . , . describes the slope of the density profile. Notethat the model is insensitive to the choice of r out unless r out /r in ≈ . ˙ M • carries units of M (cid:12) yr − .Here, we test the relationship between RM and accre-tion rate in our suite of images for M87*. In Figure 8, weplot spatially unresolved RM versus accretion rate for the 7 models considered in this work. Each symbol of the samecolour represents a different snapshot of the same model.Positive RMs are plotted with filled symbols, while negativeRMs are plotted with open symbols. Notice the ubiquitousflips in the sign of the RM that occur in all models. Thefilled blue region demarcates the relation in Marrone et al.(2006), spanning variations of the slope of the density pro-file, β ∈ [0 . , . . We set r in = 3 and r out = ∞ . The dashedline shows the upper limit from Kuo et al. (2014). A ◦ inclination appropriate for M87* is shown on the left, whilefor comparison a ◦ inclination is shown on the right.Overall, we find that a spatially unresolved RM is apoor predictor of the accretion rate, especially if the cor-rect model is not known a priori. RM and the accretionrate differ by orders of magnitude both within and amongthe different models. Even within a single model, there isno correlation between RM and accretion rate, which weexplore in more detail for one model in §3.7. Since theirhigher accretion rates imply higher number densities, SANEmodels typically have larger spatially unresolved RMs thanMAD models. However, these models exhibit such strongtime variability that they cannot be distinguished solely bythe upper limit of Kuo et al. (2014). Rather, repeat observa-tions on timescales of months to years will be necessary tocharacterise the distribution of | RM | over time and detectpotential sign flips.For a given accretion rate, the GRMHD models in thiswork produce much lower RM than the analytic model of MNRAS , ?? – ?? (2020) Ricarte et al. [] MAD, a=-0.5, R high =160 61.3061.2561.20 [] MAD, a=-0.5, R high =2039.2539.3039.35 [] MAD, a=+0.94, R high =160 36.0436.0636.08 [] MAD, a=+0.94, R high =2082.580.077.575.0 [] SANE, a=-0.94, R high =10 1.76e-061.73e-061.70e-06 [m ]80.082.585.087.5 [] SANE, a=-0.94, R high =801.76e-061.73e-061.70e-06 [m ]0100200300 [] SANE, a=+0.94, R high =160 Taylor Expanded ModelAveraged Images
Figure 7.
Non- λ behaviour of the 7 models we consider during their final simulation snapshot. Within 16 bands that are each 0.25GHz 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 6images around 228 GHz. Most of these models exhibit non-linear behaviour even within this narrow fractional bandwidth. Marrone et al. (2006) at an inclination of ◦ . As we furtherexplore in §3.5, this is because these simulations are viewedthrough an evacuated funnel region at low inclination. Withan inclination of ◦ , the | RM | more closely matches thatpredicted by Marrone et al. (2006), although they still re-main systematically offset. A similar inclination dependenceis found for SANE models in Mościbrodzka et al. (2017).The retrograde SANE models remain the most offset fromthe analytic model. Even at ◦ , the large Faraday depthoccurs in an area with little emission, since the electrons areassigned low temperatures in the mid-plane. Here, we study the dependence of RM on inclination ingreater detail. In Figure 9, we plot the distributions of RMfor all 11 snapshots of all 7 models at 5 different inclinations.In this plot, boxes contain the 25th to 75th quantiles, thehorizontal black or yellow line marks the median, and theerror bars span the full range of the 11 snapshots studied.This plot omits the sign flips observed in Figure 8, but wecomment that they remain ubiquitous at all inclinations forthese calculations which terminate at a radius of GM • /c .As in Mościbrodzka et al. (2017), we find that the ab-solute value of the RM depends on the inclination angle,but it is not large compared to the substantial scatter be-tween snapshots. This weaker dependence is likely due to MNRAS , ?? – ????
Non- λ behaviour of the 7 models we consider during their final simulation snapshot. Within 16 bands that are each 0.25GHz 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 6images around 228 GHz. Most of these models exhibit non-linear behaviour even within this narrow fractional bandwidth. Marrone et al. (2006) at an inclination of ◦ . As we furtherexplore in §3.5, this is because these simulations are viewedthrough an evacuated funnel region at low inclination. Withan inclination of ◦ , the | RM | more closely matches thatpredicted by Marrone et al. (2006), although they still re-main systematically offset. A similar inclination dependenceis found for SANE models in Mościbrodzka et al. (2017).The retrograde SANE models remain the most offset fromthe analytic model. Even at ◦ , the large Faraday depthoccurs in an area with little emission, since the electrons areassigned low temperatures in the mid-plane. Here, we study the dependence of RM on inclination ingreater detail. In Figure 9, we plot the distributions of RMfor all 11 snapshots of all 7 models at 5 different inclinations.In this plot, boxes contain the 25th to 75th quantiles, thehorizontal black or yellow line marks the median, and theerror bars span the full range of the 11 snapshots studied.This plot omits the sign flips observed in Figure 8, but wecomment that they remain ubiquitous at all inclinations forthese calculations which terminate at a radius of GM • /c .As in Mościbrodzka et al. (2017), we find that the ab-solute value of the RM depends on the inclination angle,but it is not large compared to the substantial scatter be-tween snapshots. This weaker dependence is likely due to MNRAS , ?? – ???? (2020) H Internal Faraday Rotation Accretion Rate [M yr ] R M [ r a d m ] K14 Upper LimitMarrone+2006 10 Accretion Rate [M yr ] MAD, a=+0.94, R=20MAD, a=+0.94, R=160MAD, a=-0.5, R=20MAD, a=-0.5, R=160 SANE, a=+0.94, R=160SANE, a=-0.94, R=10SANE, a=-0.94, R=80
Figure 8.
RM as a function of accretion rate for the 7 models considered in this paper. Filled symbols have positive RM, while opensymbols have negative RM. In the left panel, the observer is oriented at a 17 ◦ inclination, while in the right panel, the observer isoriented at 90 ◦ . At an inclination of ◦ , we find that accretion rates are systematically higher than those that would be inferred bysimple analytic models (Marrone et al. 2006). At 90 ◦ , we find RMs in better agreement with analytic models, although there remainssubstantial scatter. Retrograde SANE models are outliers, since their forward-jet emission does not intercept the large Faraday depth inthe mid-plane. the small radius at which we truncate our calculations. Wenotice no differences between our calculations at low incli-nation: ◦ , ◦ , and ◦ . This is fortunate for our study ofM87*, as this indicates that we need not be concerned withsmall deviations from our fiducial inclination of ◦ .At low inclinations, we view the accretion flow throughan evacuated funnel region with less Faraday rotating ma-terial than at high inclinations. We demonstrate this by cal-culating the characteristic distance of Faraday rotating ma-terial as a function of inclination. By modifying the ipole source code, we compute for each geodesic (cid:104) R FR (cid:105) ≡ (cid:90) R BL | ρ V | LP ds (cid:30) (cid:90) | ρ V | LP ds, (17)where ρ V is the frame-invariant radiative transfer coefficientresponsible for Faraday rotation, R BL is the radius of thematerial in Boyer-Lindquist (or equivalently Kerr-Schild) co-ordinates, LP = (cid:112) Q + U is the total amount of linearlypolarized emission that has been emitted along the geodesicso far (on the way to the camera), and s is the affine param-eter describing the geodesic. In other words, this is the char-acteristic distance of Faraday rotating material, weighted bythe fraction of the final linearly polarized emission that hasalready been added to the pixel on the way to the cam-era. Once R FR is computed for each pixel, a single valueis calculated for the model by computing an average across the image, weighted by total final linear polarization of eachpixel.In Figure 10, we plot the characteristic distance of Fara-day rotating material of these models during their final snap-shot as a function of inclination. For inclinations < ◦ , mostof the Faraday rotation occurs at < GM • /c , while (cid:104) R FR (cid:105) increases at higher inclinations. The innermost stable circu-lar orbit exists at smaller radius for prograde models thanfor retrograde models, which leads to a noticeable differencein (cid:104) R FR (cid:105) between these two classes of models at low incli-nation.Recall that we restrict the domain of our calculations towithin GM • /c , the radius within which the simulationsare in inflow equilibrium. As we further explore in AppendixC, if this radius is increased to M • , we find consistentresults for i (cid:54) ◦ , but substantially larger (cid:104) R FR (cid:105) for i (cid:62) ◦ . Therefore, we believe that our | RM | values throughoutthe paper should be considered lower limits for i > ◦ , asmaterial from beyond the converged region may contributeto Faraday rotation at larger inclinations.Since the RM is highly non-uniform across these images,spatially resolved RM distribution functions provide greaterinsight into the inclination dependence. In Figure 11, we plotthe RM distribution functions as in Figure 6, for inclinations i ∈ { ◦ , ◦ , ◦ , ◦ } . Unlike in Figure 6, we do not splitthe distributions into their front and back halves. At lowinclination, these models exhibit significant emission with MNRAS , ?? – ?? (2020) Ricarte et al. | R o t a t i o n M e a s u r e | [ r a d m ] MAD, a=+0.94, R high =20MAD, a=+0.94, R high =160MAD, a=-0.5, R high =20MAD, a=-0.5, R high =160 SANE, a=+0.94, R high =160SANE, a=-0.94, R high =10SANE, a=-0.94, R high =80K14 Upper Limit
Figure 9.
RM as a function of inclination. 11 snapshots are shown for each model. For M87*, 17 degrees is considered the most likelyinclination 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 25thto 75th quantiles, 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. low | RM | . As the inclination approaches ◦ , this fractionof emission with low | RM | diminishes, and the distributionis skewed towards higher values. Note that the distributionsat ◦ are indistinguishable from those at ◦ . χ ( λ ) as a Model Discriminant Since emission and Faraday rotation occur co-spatially andnon-uniformly throughout these models, χ ( λ ) need not belinear. We find that the degree of non-linearity varies sig-nificantly among models, due to the complex RM structuredescribed in §3.1. As a metric of non-linearity, we fit linesto χ ( λ ) from the 16 small bands spanning 226 to 230 GHzand obtain the coefficient of determination, R , defined via R = 1 − SS res SS tot . (18)Here, SS res is the regression sum of squares, and SS tot is the total sum of squares. R describes the fraction of vari-ation within the data that can be ascribed to the simplelinear dependence.Our results are shown in Figure 12, following the sameformatting as Figure 9. Based on these results, SANE mod-els should exhibit non-linear behaviour most of the time. Incontrast, MAD models are almost always well described bya linear χ ( λ ) law, especially those with R high = 20 . Thiscan be understood by returning to Figure 6. Pixels with | RM | > | RM | crit have individual EVPAs that rotate sub- stantially across the 4 GHz bandwidth. Images consistingof a substantial fraction of such pixels will therefore exhibitstructure in χ ( λ ) within the bandwidth. Among the modelsconsidered in this study, this behaviour appears much morelikely among SANEs. We study one model with higher time resolution in order toquantify the variability of its RM. For the MAD, a = 0 . , R high = 20 model, we create images for every available snap-shot within t/ ( GM • /c ) ∈ [7500 , , which are each sep-arated by GM • /c . This model has the smallest accretionrate and | RM | of the models explored in this work. Its | RM | is sufficiently smaller than RM crit that χ ( λ ) remains lin-ear within the bandwidth (see Figure 7), and thus we createonly 2 polarized images (two frequencies) at each snapshotinstead of the usual 6 (three frequencies, each separately fortwo sides of the disc).The time variability of this model at an inclination of ◦ is visualised in Figure 13. In the top row, panels areseparated by about half a year, while in the bottom rowpanels are separated by about 5 days. As in previous fig-ures, pixel brightness encodes the linear polarized intensity,while the colour encodes the RM. Both positive and neg-ative RM regions can be found in a typical snapshot. Thespatially unresolved RM is written at the bottom of eachpanel. The bottom row illustrates how the dynamics of animage with both positive and negative RM regions can re- MNRAS , ?? – ????
RM as a function of inclination. 11 snapshots are shown for each model. For M87*, 17 degrees is considered the most likelyinclination 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 25thto 75th quantiles, 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. low | RM | . As the inclination approaches ◦ , this fractionof emission with low | RM | diminishes, and the distributionis skewed towards higher values. Note that the distributionsat ◦ are indistinguishable from those at ◦ . χ ( λ ) as a Model Discriminant Since emission and Faraday rotation occur co-spatially andnon-uniformly throughout these models, χ ( λ ) need not belinear. We find that the degree of non-linearity varies sig-nificantly among models, due to the complex RM structuredescribed in §3.1. As a metric of non-linearity, we fit linesto χ ( λ ) from the 16 small bands spanning 226 to 230 GHzand obtain the coefficient of determination, R , defined via R = 1 − SS res SS tot . (18)Here, SS res is the regression sum of squares, and SS tot is the total sum of squares. R describes the fraction of vari-ation within the data that can be ascribed to the simplelinear dependence.Our results are shown in Figure 12, following the sameformatting as Figure 9. Based on these results, SANE mod-els should exhibit non-linear behaviour most of the time. Incontrast, MAD models are almost always well described bya linear χ ( λ ) law, especially those with R high = 20 . Thiscan be understood by returning to Figure 6. Pixels with | RM | > | RM | crit have individual EVPAs that rotate sub- stantially across the 4 GHz bandwidth. Images consistingof a substantial fraction of such pixels will therefore exhibitstructure in χ ( λ ) within the bandwidth. Among the modelsconsidered in this study, this behaviour appears much morelikely among SANEs. We study one model with higher time resolution in order toquantify the variability of its RM. For the MAD, a = 0 . , R high = 20 model, we create images for every available snap-shot within t/ ( GM • /c ) ∈ [7500 , , which are each sep-arated by GM • /c . This model has the smallest accretionrate and | RM | of the models explored in this work. Its | RM | is sufficiently smaller than RM crit that χ ( λ ) remains lin-ear within the bandwidth (see Figure 7), and thus we createonly 2 polarized images (two frequencies) at each snapshotinstead of the usual 6 (three frequencies, each separately fortwo sides of the disc).The time variability of this model at an inclination of ◦ is visualised in Figure 13. In the top row, panels areseparated by about half a year, while in the bottom rowpanels are separated by about 5 days. As in previous fig-ures, pixel brightness encodes the linear polarized intensity,while the colour encodes the RM. Both positive and neg-ative RM regions can be found in a typical snapshot. Thespatially unresolved RM is written at the bottom of eachpanel. The bottom row illustrates how the dynamics of animage with both positive and negative RM regions can re- MNRAS , ?? – ???? (2020) H Internal Faraday Rotation Inclination [deg] R F R [ R g ] MAD, a=+0.94, R=20MAD, a=+0.94, R=160MAD, a=-0.5, R=20MAD, a=-0.5, R=160SANE, a=+0.94, R=160SANE, a=-0.94, R=10SANE, a=-0.94, R=80
Figure 10.
Characteristic distance of Faraday rotating materialin these models during their final snapshot as a function of incli-nation. Due to the evacuated jet region in these simulations, theFaraday rotating material is confined to low radius at low incli-nation, but extends to larger radius at higher inclination. Recallthat our calculations terminate at R = 20 GM • /c , within whichthese simulations are in inflow equilibrium. sult in variability and RM sign flips on a timescale of a fewdays. The RM sign flip does not require a drastic change inglobal source structure; rather, the balance between positiveand negative RM regions is shifted.In Figure 14, we plot the RM as a function of time forthis model, as well as its auto-correlation function. The ge-ometrized time unit is converted to days via t • = GM • /c ,which for M87* is 8.5 hours. The grey band encloses the16th to 84th ( σ ) percentiles. At ◦ , these include bothpositive and negative values, such that RM = − . +1 . − . × rad m − . Examining the images, we do not notice anyobvious special behaviour in the accretion flow during peri-ods of large | RM | . Since the RM is determined by the mo-tion of material on event horizon scales, the auto-correlationfunction of this time series drops rapidly, falling below 0.5in less than the separation between snapshots.In Figure 15, we plot the joint probability distributionsof log | RM | with accretion rate, linear polarized intensity,and circular polarized intensity for the model at ◦ . One,two, and three sigma contours are overlaid in white. We donot find any correlation between | RM | and ˙ M • , indicatingthat within a single model, a change in RM does not im-ply a change in the accretion rate, as might be suggestedby analytic models. Rather, as we have discussed, | RM | andits sign appears to result from a complicated and stochas-tic cancellation of positive and negative regions. For SgrA*, Bower et al. (2018) found an anti-correlation betweenlinear polarized intensity and RM, but no such correlationwith circular polarized intensity. For this particular model of M87*, we recover qualitatively similar results: a linearregression yields log | RM | = − . LP + 2 . ± . ,where LP = (cid:112) Q + U in Jy and RM is in units of rad m − ,with a moderate r-value of -0.57. No statistically significantcorrelation is found between | RM | and circular polarization.An anti-correlation between | RM | and LP is not surprising,since greater Faraday rotation implies greater scrambling ofthe polarization vector field. Here, we summarize the qualitative commonalities and dif-ferences between the different models we have considered. • Prograde MAD : These models require the lowestaccretion rate to generate the appropriate total intensity,and consequently exhibit the lowest | RM | . Compared to theother models, there is not too much difference in | RM | be-tween the two halves of the emitting region, since both com-ponents occur close to the mid-plane. These models usuallyexhibit linear χ ( λ ) within 4 GHz. • Retrograde MAD : Retrograde MAD models requirelarger accretion rates than their prograde counterparts, butexhibit similar values of | RM | . Some areas of the R high =160 models have large enough | RM | to weakly bandwidthdepolarize or scramble portions of the image. • Prograde SANE : Due to lensing, the counter-jetspans a larger angular scale and contributes more to the totalintensity than the forward-jet. Bandwidth depolarization ef-fects are severe, and emission from the counter-jet is entirelydepolarized, assuming a 4 GHz bandwidth. Consequently,the total intensity image (dominated by the counter-jet) ap-pears morphologically different compared to the linearly po-larized image (dominated by the forward-jet). χ ( λ ) exhibitsstrongly non-linear evolution. • Retrograde SANE : The counter-jet experiences 6 or-ders of magnitude more Faraday rotation than the forward-jet. As with their prograde counterparts, χ ( λ ) is highlynon-linear. However, the difference in morphology betweenthe total intensity image and the linearly polarized image isless significant than the prograde case, since the two emis-sion components subtend more similar angular scales.By construction, increasing R high decreases the temperatureof electrons in the mid-plane, which therefore increases theFaraday depth for emission that passes through it. In allmodels, larger R high results in greater non-linearity of χ ( λ ) . If anywhere in an image, RM > RM crit , then that region’slinear polarization should be suppressed as described byequation 13. We find that this is more likely to occur inSANE models, but may also affect some parts of MAD mod-els. Bandwidth depolarized regions manifest as areas withlower than average linear polarization fraction.In some of the worst cases, like that presented in Figure2, we find that basic image properties are affected by band-width depolarization, such as the linear polarization fractionas well as the morphology of the linear polarization vector MNRAS , ?? – ?? (2020) Ricarte et al. d | p |/ d l o g | R M | MAD, a=+0.94, R=160 MAD, a=+0.94, R=20 MAD, a=-0.5, R=160 MAD, a=-0.5, R=203 5 7 9 11log |RM| [rad m ]0.00.51.0 d | p |/ d l o g | R M | SANE, a=+0.94, R=160 3 5 7 9 11log |RM| [rad m ]SANE, a=-0.94, R=10 3 5 7 9 11log |RM| [rad m ]SANE, a=-0.94, R=80 log RM crit (4 GHz)5306090 Figure 11.
RM distribution functions as in Figure 6, now shown as a function of inclination. Both halves of the emitting region arecombined in this figure. At higher inclinations, we find that these distributions are skewed towards higher values, as the population ofphotons experiencing comparatively little RM diminishes. field, even after images are blurred. This may be importantfor studies which use images computed at a single frequencyto compare to observations taken over a finite bandwidth.For future studies, the methodology introduced in §2.2 canbe generalized by applying the Taylor expansion to eachpoint along the ray-traced geodesic instead of just two sep-arate emission regions. This would allow the appropriatebandwidth integrations to occur within the ray-tracing codeitself.
At present, the electron distribution function is poorly con-strained. For creating these images, a thermal distributionfunction is assumed, along with the R high prescription de-veloped by Mościbrodzka et al. (2016). Mao et al. (2017)studied the effects of adding a power-law component to thedistribution function and found that even if a few percent ofthe total energy is put into a non-thermal power law compo-nent, a diffuse halo of emission can be produced. The effectsof non-thermal electron distribution functions on polarizedimages remain to be studied.Recall that we truncate our radiative transfer calcula-tions at a radius of GM • /c , only within which theseGRMHD simulations exhibit inflow equilibrium. However,Faraday rotating material can plausibly exist at larger ra-dius, especially at higher inclinations. Using very long du-ration simulations, Dexter et al. (2020) find that Faradayrotation can peak at radii R ∼ − GM • /c , depend-ing on the electron prescriptions. A more distant Faraday screen would be expected to uniformly rotate all EVPAs by afixed amount, which may leave signatures in the EVPA vec-tor field (Palumbo et al. 2020). Additionally, such a screenshould maintain a consistent RM for longer timescales thanthese models. Such distant screens may in fact be requiredto explain the consistent RMs of Sgr A* (Bower et al. 2018)and 3C 84 (Plambeck et al. 2014) over timescales of years. We have investigated rotation measure (RM) within a sub-set of the EHT simulation library that is consistent with theobservational constraints on M87* considered in EHT5. Wefind more information in the RM structure of GRMHD sim-ulations than can be described by a single scalar, RM . Wesummarise our results below: • In a single snapshot, we find extreme variations in theRM between different regions. The RM may vary by ordersof magnitude and even flip sign across the image. The RMinferred from a spatially unresolved measurement is there-fore the result of the complicated interplay of these differentregions. • Emission originating from in front of the disk mid-planemay be orders of magnitude less Faraday rotated than emis-sion from the back. In the high accretion rate SANE models,the RM is large enough to completely depolarize emissionfrom the counter-jet, which may in fact dominate the to-tal intensity. The sub-relativistic mid-plane is the dominantsource of Faraday rotation in these models.
MNRAS , ?? – ????
MNRAS , ?? – ???? (2020) H Internal Faraday Rotation Inclination [ ] R , L i n e a r F i t MAD, a=+0.94, R high =20MAD, a=+0.94, R high =160MAD, a=-0.5, R high =20MAD, a=-0.5, R high =160 SANE, a=+0.94, R high =160SANE, a=-0.94, R high =10SANE, a=-0.94, R high =80
Figure 12.
Non-linearity of χ ( λ ) among the models considered in this study. R of a linear fit to χ ( λ ) describes the fraction ofvariation in the data that can be explained by a simple linear model. SANE models exhibit more non-linearity than MAD models in thisstudy, since their images contain a significant amount of polarized intensity with | RM | > | RM | crit . • Many models exhibit clear departures from a λ laweven across a narrow fractional bandwidth of 4 GHz. Non-linearity is a more common feature among the SANE mod-els, and increases with R high . • The RM structure changes as material moves on eventhorizon scales. These models all exhibit strong time variabil-ity, causing the spatially unresolved RM to vary and evenflip sign on a time-scale of days. • RM is a poor predictor of the accretion rate. Thesemodels predict several orders of magnitude spread in RMfor a given accretion rate, and within a single model, thesequantities are not correlated as a function of time. In addi-tion, analytic models used to infer the accretion rate basedon the RM (Marrone et al. 2006) systematically underesti-mate the accretion rate onto M87*, since the source shouldbe viewed through an evacuated funnel region.In future work, a more thorough investigation of theEHT simulation library is merited, including models for SgrA*. Alternative models for the electron distribution functionshould also be considered. Repeated observations of both SgrA* and M87* will be useful to probe the time variabilitypredicted by these models.
We thank Jason Dexter for valuable feedback which greatlyimproved the content of this paper. We also thank Lau-rent Loinard and the anonymous referee for their careful and thorough reviews. This material is based upon worksupported by the National Science Foundation under GrantNo. OISE 1743747 and NSF AST 1716327. Computationswere performed using the resources of the Black Hole Ini-tiative (BHI). Computations at the BHI were made pos-sible through the support of grants from the Gordon andBetty Moore Foundation and the John Templeton Founda-tion. This work used the Extreme Science and Engineer-ing Discovery Environment (XSEDE) resource stampede2at TACC through allocation TG-AST170024. The opinionsexpressed in this publication are those of the author(s) anddo not necessarily reflect the views of the Moore or Temple-ton Foundations.
The data underlying this article will be shared on reasonablerequest to the corresponding author.
REFERENCES
Agol E., 2000, ApJ, 538, L121Agudo I., Thum C., Ramakrishnan V., Molina S. N., Casadio C.,Gómez J. L., 2018, MNRAS, 473, 1850Bardeen J. M., Wagoner R. V., 1969, ApJ, 158, L65Bower G. C., Wright M. C. H., Falcke H., Backer D. C., 2003,ApJ, 588, 331Bower G. C., Dexter J., Markoff S., Rao R., Plambeck R. L., 2017,ApJ, 843, L31MNRAS , ?? – ?? (2020) Ricarte et al. y [ a s ] RM = 1.20 × 10 rad m t = 0 days RM = 0.07 × 10 rad m t = 177 days RM = 0.28 × 10 rad m t = 354 days RM = 4.92 × 10 rad m t = 530 days
30 0 30 x [ as] y [ a s ] RM = 1.27 × 10 rad m t = 629 days
30 0 30 x [ as] RM = 0.25 × 10 rad m t = 633 days
30 0 30 x [ as] RM = 0.21 × 10 rad m t = 636 days
30 0 30 x [ as] RM = 0.24 × 10 rad m t = 640 days R M [ r a d m ] Figure 13.
Visualisation of the time variability of the RM structure in the MAD, a = 0 . , R high = 20 model. The brightness of eachpixel scales with its linear polarized intensity, while the colour represents its rotation measure. The RM that would be inferred from aspatially unresolved measurement is written at the bottom of each panel. The spatially unresolved RM can change on time-scales of daysas the Faraday rotating gas moves on event horizon scales. The sign flip does not require a dramatic change in the source structure. Time [days] U n r e s o l v e d R M [ r a d m ] MAD, a = 0.94, R high =20 Time [days] A u t o c o rr e l a t i o n Figure 14.
Left:
RM as a function of time in the MAD, a = 0 . , R high = 20 model that we study with greater time resolution. Thegrey band encloses the σ percentile region over this time, which includes both positive and negative values. Right:
Auto-correlationfunction of this time series. The auto-correlation drops below 50 per cent in less than the separation between snapshots.MNRAS , ?? – ????
Auto-correlationfunction of this time series. The auto-correlation drops below 50 per cent in less than the separation between snapshots.MNRAS , ?? – ???? (2020) H Internal Faraday Rotation log M [M yr ] l o g | R M | [ r a d m ] log Q + U [Jy] log| V | [Jy] Figure 15.
Joint probability distributions of log | RM | with accretion rate, linear polarized intensity, and circular polarized intensityin the MAD, a = 0 . , R high = 20 model with an inclination of 17 ◦ . One, two, and three sigma contours are overlaid in white, usingsolid, dashed, and dotted lines respectively. We find no correlation between | RM | and ˙ M • , implying that a change in RM does not implya change in the accretion rate. As observed by Bower et al. (2018) for Sgr A*, we recover an anti-correlation between | RM | and linearpolarization, but not circular polarization.Bower G. C., et al., 2018, ApJ, 868, 101Brentjens M. A., de Bruyn A. G., 2005, A&A, 441, 1217Broderick A., Blandford R., 2003, MNRAS, 342, 1280Broderick A. E., Loeb A., 2006a, MNRAS, 367, 905Broderick A. E., Loeb A., 2006b, ApJ, 636, L109Broderick A. E., McKinney J. C., 2010, ApJ, 725, 750Bromley B. C., Melia F., Liu S., 2001, ApJ, 555, L83Brunthaler A., Bower G. C., Falcke H., Mellon R. R., 2001, ApJ,560, L123Burn B. J., 1966, MNRAS, 133, 67Chael A., Narayan R., Johnson M. D., 2019, MNRAS, 486, 2873Chatterjee K., et al., 2020, arXiv e-prints, p. arXiv:2002.08386Childs H., et al., 2012, in , High Performance Visualization–Enabling Extreme-Scale Scientific Insight. pp 357–372Contopoulos I., Christodoulou D. M., Kazanas D., GabuzdaD. C., 2009, ApJ, 702, L148Dexter J., 2016, MNRAS, 462, 115Dexter J., McKinney J. C., Agol E., 2012, MNRAS, 421, 1517Dexter J., et al., 2020, arXiv e-prints, p. arXiv:2004.00019Eatough R. P., et al., 2013, Nature, 501, 391Event Horizon Telescope Collaboration et al., 2019a, ApJ, 875,L1Event Horizon Telescope Collaboration et al., 2019b, ApJ, 875,L2Event Horizon Telescope Collaboration et al., 2019c, ApJ, 875,L3Event Horizon Telescope Collaboration et al., 2019d, ApJ, 875,L4Event Horizon Telescope Collaboration et al., 2019e, ApJ, 875,L5Event Horizon Telescope Collaboration et al., 2019f, ApJ, 875, L6Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007,ApJ, 668, 417Frick P., Sokoloff D., Stepanov R., Beck R., 2011, MNRAS, 414,2540Gabuzda D., 2018, Galaxies, 7, 5Gammie C. F., McKinney J. C., Tóth G., 2003, ApJ, 589, 444 Gammie C. F., Shapiro S. L., McKinney J. C., 2004, ApJ, 602,312Gardner F. F., Whiteoak J. B., 1966, ARA&A, 4, 245Gebhardt K., Adams J., Richstone D., Lauer T. R., Faber S. M.,Gültekin K., Murphy J., Tremaine S., 2011, ApJ, 729, 119Himwich E., Johnson M. D., Lupsasca A. r., Strominger A., 2020,arXiv e-prints, p. arXiv:2001.08750Hovatta T., O’Sullivan S., Martí-Vidal I., Savolainen T.,Tchekhovskoy A., 2019, A&A, 623, A111Howes G. G., 2010, MNRAS, 409, L104Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ,592, 1042Jiménez-Rosales A., Dexter J., 2018, MNRAS, 478, 1875Johnson M. D., et al., 2015, Science, 350, 1242Johnson M. D., et al., 2019, arXiv e-prints, p. arXiv:1907.04329Jones T. W., Odell S. L., 1977, ApJ, 214, 522Kawazura Y., Barnes M., Schekochihin A. A., 2019, Proceedingsof the National Academy of Science, 116, 771Kim J. Y., et al., 2019, A&A, 622, A196King A. R., Pringle J. E., Hofmann J. A., 2008, MNRAS, 385,1621Kuo C. Y., et al., 2014, ApJ, 783, L33Le Roux E., 1961, Annales d’Astrophysique, 24, 71Li Y.-P., Yuan F., Xie F.-G., 2016, ApJ, 830, 78Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M.,Markoff S. B., Van Moer M., 2020, MNRAS,Macquart J.-P., Bower G. C., Wright M. C. H., Backer D. C.,Falcke H., 2006, ApJ, 646, L111Mao S. A., Dexter J., Quataert E., 2017, MNRAS, 466, 4307Marrone D. P., Moran J. M., Zhao J.-H., Rao R., 2006, ApJ, 640,308Marrone D. P., Moran J. M., Zhao J.-H., Rao R., 2007, ApJ, 654,L57Medeiros L., Chan C.-k., Özel F., Psaltis D., Kim J., MarroneD. P., Sadowski A., 2018, ApJ, 856, 163Melrose D. B., 1997, Journal of Plasma Physics, 58, 735Mościbrodzka M., Gammie C. F., 2018, MNRAS, 475, 43MNRAS , ?? – ?? (2020) Ricarte et al.
Mościbrodzka M., Falcke H., Shiokawa H., 2016, A&A, 586, A38Mościbrodzka M., Dexter J., Davelaar J., Falcke H., 2017, MN-RAS, 468, 2214Narayan R., Yi I., 1994, ApJ, 428, L13Narayan R., Yi I., 1995a, ApJ, 452, 710Narayan R., Yi I., 1995b, ApJ, 452, 710Narayan R., Igumenshchev I. V., Abramowicz M. A., 2000, ApJ,539, 798Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ,55, L69Narayan R., SÄ dowski A., Penna R. F., Kulkarni A. K., 2012,MNRAS, 426, 3241Palumbo D. C. M., Wong G. N., Prather B. S., 2020, arXiv e-prints, p. arXiv:2004.01751Pasetto A., Carrasco-González C., O’Sullivan S., Basu A., BruniG., Kraus A., Curiel S., Mack K.-H., 2018, A&A, 613, A74Plambeck R. L., et al., 2014, ApJ, 797, 66Porth O., et al., 2019, ApJS, 243, 26Quataert E., Gruzinov A., 1999, ApJ, 520, 248Quataert E., Gruzinov A., 2000a, ApJ, 539, 809Quataert E., Gruzinov A., 2000b, ApJ, 545, 842Rees M. J., Begelman M. C., Blandford R. D., Phinney E. S.,1982, Nature, 295, 17Roelofs F., Johnson M. D., Shiokawa H., Doeleman S. S., FalckeH., 2017, ApJ, 847, 55Ryan B. R., Ressler S. M., Dolence J. C., Gammie C., QuataertE., 2018, ApJ, 864, 126Rybicki G. B., Lightman A. P., 1986, Radiative Processes in As-trophysicsSazonov V. N., 1969, Soviet Ast., 13, 396Shapiro S. L., Lightman A. P., Eardley D. M., 1976, ApJ, 204,187Sharma P., Quataert E., Hammett G. W., Stone J. M., 2007, ApJ,667, 714Shcherbakov R. V., 2008, ApJ, 688, 695Shcherbakov R. V., Penna R. F., McKinney J. C., 2012, ApJ, 755,133Sądowski A., Narayan R., Penna R., Zhu Y., 2013, MNRAS, 436,3856Sądowski A., Wielgus M., Narayan R., Abarca D., McKinneyJ. C., Chael A., 2017, MNRAS, 466, 705Thorne K. S., 1974, ApJ, 191, 507Tsunetoe Y., Mineshige S., Ohsuga K., Kawashima T., AkiyamaK., 2020, PASJ,Volonteri M., Sikora M., Lasota J. P., Merloni A., 2013, ApJ, 775,94Walker R. C., Hardee P. E., Davies F. B., Ly C., Junor W., 2018,ApJ, 855, 128Zavala R. T., Taylor G. B., 2004, ApJ, 612, 749
APPENDIX A: SPHERICAL STOKESPARAMETERS
In §2.2, we convert the standard Stokes parameters { Q, U, V } to the spherical Stokes parameters { N, φ, ψ } asin Shcherbakov et al. (2012) for the purposes of a more sta-ble Taylor expansion. This transformation is defined by: Q = N cos φ sin ψ (A1) U = N sin φ sin ψ (A2) V = N cos ψ (A3) while its inversion can be derived from trigonometric iden-tities as N = (cid:112) Q + U + V (A4) φ = arctan( U/Q ) (A5) ψ = arctan( (cid:112) Q + U /V ) . (A6)From these equations, we can see that N is the total amountof both linear and circular polarization, φ describes thephase of the linear polarization, and ψ describes the lin-ear to circular polarization ratio. Note the lack of a factorof 1/2 in φ , such that φ = 2 χ , where χ is the EVPA.In a pixel with a large | RM | , Faraday rotation causes Q and U to oscillate rapidly with frequency. In sphericalStokes parameters, N remains stable, while φ changes lin-early with wavelength squared, with dφ/dλ = 2RM . Thismakes the Spherical stokes parameters more stable to Taylorexpansion.When we calculate derivatives of the { N, φ, ψ } , we com-pute them in terms of { Q, U, V } and their derivatives. Thisallows us to avoid mistakes due to phase wrapping, as dis-cussed in §2.2. By simply differentiating according to thechain rule, the derivatives are given by N (cid:48) = QQ (cid:48) + UU (cid:48) + V V (cid:48) (cid:112) Q + U + V (A7) φ (cid:48) = U (cid:48) Q − Q (cid:48) UQ + U (A8) ψ (cid:48) = V ( QQ (cid:48) + UU (cid:48) ) − V (cid:48) ( Q + U )( Q + U + V ) (cid:112) Q + U (A9)where (cid:48) denotes differentiation with respect to frequency (oranother physically similar quantity such as wavelength orwavelength-squared). APPENDIX B: ANALYTIC BANDWIDTHINTEGRALS
Here, we provide analytic solutions to the integrals describedin §2.2. Spherical Stokes parameters { I , N , φ , ψ } andtheir derivatives { dI/dν, dN/dν, dφ/dλ , dψ/dλ } are esti-mated at frequency ν . Let ν c be the central frequency of aband extending between ν = ν c − ∆ ν/ and ν = ν c +∆ ν/ .We then define a dimensionless frequency x = ( ν − ν ) / ∆ ν such that x = x ( ν ) and x = x ( ν ) . We can then write I BW = (cid:90) x x ( I + I x ) dx (B1) Q BW = (cid:90) x x ( N + N x ) cos( φ + φ x ) sin( ψ + ψ x ) dx (B2) U BW = (cid:90) x x ( N + N x ) sin( φ + φ x ) sin( ψ + ψ x ) dx (B3) V BW = (cid:90) x x ( N + N x ) cos( ψ + ψ x ) dx (B4)where x = ( ν c − ν − ∆ ν/ / ∆ ν , x = ( ν c − ν +∆ ν/ / ∆ ν ,and we define the following quantities from the derivativesof the spherical Stokes parameters: MNRAS , ?? – ?? (2020) H Internal Faraday Rotation I = dIdν ∆ ν (B5) N = dNdν ∆ ν (B6) φ = − c ∆ νν dφdλ (B7) ψ = − c ∆ νν dψdλ (B8)These integrals have analytic solutions given by I BW = I ( x − x ) + I (cid:16) x − x (cid:17) (B9) Q BW = 12 (cid:104) N (cid:16) cos( φ + φ x − ψ − ψ x ) φ − ψ − cos( φ + φ x − ψ − ψ x ) φ − ψ +cos( φ + x ( φ + ψ ) + ψ φ + ψ ) − cos( φ + x ( φ + ψ ) + ψ ) φ + ψ (cid:17) + N (cid:16) sin( φ + φ x − ψ − ψ x )( φ − ψ ) − x ( φ − ψ ) cos( φ + φ x − ψ − ψ x )( φ − ψ ) + x ( φ + ψ ) cos( φ + x ( φ + ψ ) + ψ )( φ + ψ ) − sin( φ + x ( φ + ψ ) + ψ )( φ + ψ ) + x ( φ − ψ ) cos( φ + φ x − ψ − ψ x )( φ − ψ ) − sin( φ + φ x − ψ − ψ x )( φ − ψ ) +sin( φ + x ( φ + ψ ) + ψ )( φ + ψ ) − x ( φ + ψ ) cos( φ + x ( φ + ψ ) + ψ )( φ + ψ ) (cid:17)(cid:105) (B10) U BW = 12 N (cid:104) − sin( φ + φ x − ψ − ψ x ) φ − ψ +sin( φ + x ( φ + ψ ) + ψ ) φ + ψ +sin( φ + φ x − ψ − ψ x ) φ − ψ − sin( φ + x ( φ + ψ ) + ψ ) φ + ψ (cid:105) +12 N (cid:104) − x ( φ − ψ ) sin( φ + φ x − ψ − ψ x )( φ − ψ ) +cos( φ + φ x − ψ − ψ x )( φ − ψ ) + x ( φ + ψ ) sin( φ + x ( φ + ψ ) + ψ )( φ + ψ ) +cos( φ + x ( φ + ψ ) + ψ )( φ + ψ ) + x ( φ − ψ ) sin( φ + φ x − ψ − ψ x )( φ − ψ ) +cos( φ + φ x − ψ − ψ x )( φ − ψ ) − x ( φ + ψ ) sin( φ + x ( φ + ψ ) + ψ )( φ + ψ ) +cos( φ + x ( φ + ψ ) + ψ )( φ + ψ ) (cid:105) (B11) V BW = 1 ψ (cid:104) − ψ ( N + N x ) sin( ψ + ψ x )+ ψ ( N + N x ) sin( ψ + ψ x ) − N cos( ψ + ψ x ) + N cos( ψ + ψ x ) (cid:105) (B12) APPENDIX C: EFFECTS OF CHANGING THEMAXIMUM INTEGRATION RADIUS In ipole , the parameter rmax_geo (henceforth R out ) setsthe radius in Boyer-Lindquist (or Kerr-Schild) coordinateswithin which radiative transfer coefficients are calculated.Although the MAD and SANE simulations have outerboundaries of GM • /c and GM • /c respectively, wefind that these simulations exhibit inflow equilibrium onlywithin GM • /c . Hence, for the results throughout thispaper, R out is therefore set to GM • /c . Here, we explorethe effects of changing this outer radius to GM • /c . Al-though this choice now includes material from unconvergedregions, this allows us to gain some insight into how gas inmore distant regions might affect our predictions.In Figure C1, we compare the characteristic distance ofFaraday rotating material, (cid:104) R FR (cid:105) , as in Figure 10. The boldsolid lines originate from the R out = 20 GM • /c models (asshown in Figure 10), while the faint thin lines originate fromthe R out = 50 GM • /c models. Interestingly, there is littledifference in (cid:104) R FR (cid:105) for inclinations i (cid:54) ◦ , an evacuatedfunnel region in these simulations. At larger inclinations,distant material from unconverged regions could potentiallycontribute the majority of the Faraday rotation.In Figure C2, we plot the RM distribution functionsas in Figure 11, where our R out = 20 GM • /c results are MNRAS , ?? – ?? (2020) Ricarte et al.
Inclination [deg] R F R [ R g ] MAD, a=+0.94, R=20MAD, a=+0.94, R=160MAD, a=-0.5, R=20MAD, a=-0.5, R=160SANE, a=+0.94, R=160SANE, a=-0.94, R=10SANE, a=-0.94, R=80
Figure C1.
Characteristic distance of Faraday rotating material, as in Figure 10, for R out = 20 GM • /c models as solid lines, and R out = 50 GM • /c models as faint lines. There is little difference for inclinations i (cid:54) ◦ , where the BH is viewed through an evacuatedfunnel region. In contrast, material from more distant, unconverged regions can dominate the Faraday rotation at higher inclinations ifcalculations are allowed to proceed into this area. shown as solid lines and alternative R out = 50 GM • /c re-sults are plotted as dotted lines. For clarity, we only plotthe median values at a given log | RM | instead of the fullrange plotted in Figure 11. We find that there is negligibledifference in our results for inclinations i (cid:54) ◦ , the rangerelevant for M87*. At higher inclinations, the MAD distri-butions are skewed towards higher values for i (cid:62) ◦ , whilethe retrograde SANE models only differ at i = 90 ◦ . Theprograde SANE model shows negligible difference between R out = 20 GM • /c and R out = 50 GM • /c even at i = 90 ◦ .In Figure C3, we plot the effect this has on the spa-tially unresolved RM observed for these sources. Here, thesolid boxes correspond to the R out = 20 GM • /c mod-els (as in Figure 9), while the faint boxes correspond to R out = 50 GM • /c models. Lines demarcating medians havebeen removed for clarity. As expected, there is no noticeabledifference for inclinations i (cid:54) ◦ . For inclinations of (cid:62) ◦ ,the RMs of R out = 50 GM • /c models can be a factor of afew to orders of magnitude larger, depending on the model.Finally, we notice that at inclinations of ◦ , RM signflips occur less frequently for the R out = 50 GM • /c mod-els than for the R out = 20 GM • /c models. This is to beexpected, since material at larger radii evolves on longertimescales. MNRAS , ?? – ????
Characteristic distance of Faraday rotating material, as in Figure 10, for R out = 20 GM • /c models as solid lines, and R out = 50 GM • /c models as faint lines. There is little difference for inclinations i (cid:54) ◦ , where the BH is viewed through an evacuatedfunnel region. In contrast, material from more distant, unconverged regions can dominate the Faraday rotation at higher inclinations ifcalculations are allowed to proceed into this area. shown as solid lines and alternative R out = 50 GM • /c re-sults are plotted as dotted lines. For clarity, we only plotthe median values at a given log | RM | instead of the fullrange plotted in Figure 11. We find that there is negligibledifference in our results for inclinations i (cid:54) ◦ , the rangerelevant for M87*. At higher inclinations, the MAD distri-butions are skewed towards higher values for i (cid:62) ◦ , whilethe retrograde SANE models only differ at i = 90 ◦ . Theprograde SANE model shows negligible difference between R out = 20 GM • /c and R out = 50 GM • /c even at i = 90 ◦ .In Figure C3, we plot the effect this has on the spa-tially unresolved RM observed for these sources. Here, thesolid boxes correspond to the R out = 20 GM • /c mod-els (as in Figure 9), while the faint boxes correspond to R out = 50 GM • /c models. Lines demarcating medians havebeen removed for clarity. As expected, there is no noticeabledifference for inclinations i (cid:54) ◦ . For inclinations of (cid:62) ◦ ,the RMs of R out = 50 GM • /c models can be a factor of afew to orders of magnitude larger, depending on the model.Finally, we notice that at inclinations of ◦ , RM signflips occur less frequently for the R out = 50 GM • /c mod-els than for the R out = 20 GM • /c models. This is to beexpected, since material at larger radii evolves on longertimescales. MNRAS , ?? – ???? (2020) H Internal Faraday Rotation d | p |/ d l o g | R M | MAD, a=+0.94, R=160 MAD, a=+0.94, R=20 MAD, a=-0.5, R=160 MAD, a=-0.5, R=203 5 7 9 11log |RM| [rad m ]0.00.51.0 d | p |/ d l o g | R M | SANE, a=+0.94, R=160 3 5 7 9 11log |RM| [rad m ]SANE, a=-0.94, R=10 3 5 7 9 11log |RM| [rad m ]SANE, a=-0.94, R=80 R out = 20 GM / c R out = 50 GM / c log RM crit (4 GHz) Figure C2.
Rotation measure distribution functions, as in Figure 6, where R out = 20 GM • /c models are shown as solid lines and R out = 50 GM • /c models are shown as dotted lines. For clarity, only median values at a given log | RM | are plotted. There is negligibledifference for all models when the inclination (cid:54) ◦ . The distributions are skewed towards higher values at ◦ for the retrograde SANEs,and for (cid:62) ◦ for MADs.MNRAS , ?? – ?? (2020) Ricarte et al. | R o t a t i o n M e a s u r e | [ r a d m ] MAD, a=+0.94, R high =20MAD, a=+0.94, R high =160MAD, a=-0.5, R high =20MAD, a=-0.5, R high =160 SANE, a=+0.94, R high =160SANE, a=-0.94, R high =10SANE, a=-0.94, R high =80K14 Upper Limit
Figure C3.
Rotation measure as a function of inclination, as in Figure 9, for R out = 20 GM • /c models as solid boxes, and R out =50 GM • /c models as faint boxes. While including material at R out > GM • /c makes little difference for inclinations i (cid:54) ◦ , it mayincrease the RM by factors of a few to orders of magnitude at larger inclinations, depending on the model. MNRAS , ?? – ????