A publishing partnership

The following article is Open access

PATOKA: Simulating Electromagnetic Observables of Black Hole Accretion

, , , , , , , , ,

Published 2022 April 13 © 2022. The Author(s). Published by the American Astronomical Society.
, , Citation George N. Wong et al 2022 ApJS 259 64DOI 10.3847/1538-4365/ac582e

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0067-0049/259/2/64

Abstract

The Event Horizon Telescope (EHT) has released analyses of reconstructed images of horizon-scale millimeter emission near the supermassive black hole at the center of the M87 galaxy. Parts of the analyses made use of a large library of synthetic black hole images and spectra, which were produced using numerical general relativistic magnetohydrodynamics fluid simulations and polarized ray tracing. In this article, we describe the PATOKA pipeline, which was used to generate the Illinois contribution to the EHT simulation library. We begin by describing the relevant accretion systems and radiative processes. We then describe the details of the three numerical codes we use, iharm, ipole, and igrmonty, paying particular attention to differences between the current generation of the codes and the originally published versions. Finally, we provide a brief overview of simulated data as produced by PATOKA and conclude with a discussion of limitations and future directions.

Export citation and abstractBibTeXRIS

1. Introduction

The Event Horizon Telescope (EHT) is a globe-spanning network of millimeter wavelength observatories that is capable of imaging two supermassive black holes at event-horizon-scale resolutions. In 2019, the EHT published the first total intensity images of horizon-scale emission from the center of the giant elliptical galaxy Messier 87 (Event Horizon Telescope Collaboration et al. 2019a, 2019b, 2019c, 2019d, 2019e, 2019f, hereafter EHTC I–VI), and in 2021, the first linearly polarized images of the same source (Event Horizon Telescope Collaboration et al. 2021a, 2021b, hereafter EHTC VII–VIII). Representative images from the publications are reproduced in Figure 1. Improvements to the EHT including new antennas, increased sensitivities, and support for measurements at different frequencies lie on the horizon.

Figure 1. Refer to the following caption and surrounding text.

Figure 1. Image of the supermassive black hole at the center of the Messier 87 galaxy at 230 GHz produced by the Event Horizon Telescope Collaboration from observations taken on 2017 April 11. Left: total intensity of light shown in units of brightness temperature Tb = S λ2/2kΩ, where S is the (total) flux density, λ is the observing wavelength, k is the Boltzmann constant, and Ω is the solid angle of the observing element. Right: linear polarimetric data from the same observation plotted atop gray-scale total intensity as in the left panel. The orientation of each tick corresponds to the electric-vector position angle projected onto the image plane, and its color encodes the local fractional linear polarization. Tick marks are only shown in regions where the local total intensity is >1/10 its maximum value and fractional linear polarization is >1/5 its maximum. The FWHM of the largest Gaussian blurring kernel used in image reconstruction is represented as a circle with diameter 20 μas in the bottom right of the panel. See EHTC IV and EHTC VII for more details.

Standard image High-resolution image

The EHT operates at millimeter wavelengths, and the first images of the presumed black hole at the center of Messier 87 (hereafter M87) were produced at 230 GHz (1.3 mm; the EHT observes at two 2 GHz-wide bands centered at 227.1 and 229.1 GHz, EHTC I). The reconstructed images exhibit an asymmetric ring-like structure with diameter 42 ± 3 μ, which is consistent with the image of the shadow of a black hole with mass in the background of emission produced by hot, accreting plasma, as predicted by general relativity 17 (EHTC I; EHTC V; EHTC VI). EHT modeling connected the ring to the effects of strong lensing of near-horizon emission and linked the orientation of the brightness asymmetry to the motion of the emitting plasma and the spin of the hole (EHTC V; see also Wong et al. 2021a). The data also showed a peak in linear polarization along the southwestern segment of the ring that is consistent with the presence of organized, poloidal magnetic fields in the emission region. The overall low fractional linear polarization suggests scrambling of the polarimetric signal on small scales; this scrambling is attributed to internal Faraday rotation (EHTC VIII; see also Mościbrodzka et al. 2017; Jiménez-Rosales & Dexter 2018; Ricarte et al. 2020). Combined with data from other sources, the aggregate EHT data were used to infer constraints on the spin, mass, and accretion rate, as well as the magnetic-field strength and structure in the M87 accretion system (EHTC I; EHTC VIII). The theory constraints rely in part on understanding the complicated relationship between the detailed image features and the underlying spacetime geometry and accretion flow.

The EHT analysis relies in large part on a library of synthetic observations. This library comprises a set of numerical simulations of magnetized, relativistic black hole accretion flow models as well as accompanying polarized ray-traced simulated images and spectral energy distributions. These outputs can then be processed into synthetic observations via pipelines that simulate the effects of the observing process, e.g., eht-imaging (Chael et al. 2016), SYMBA (Roelofs et al. 2020), or THEMIS (Broderick et al. 2020). Synthetic observations serve two primary functions: they can be used to validate image reconstruction procedures, and they can be used in the forward modeling pipelines that compare the observational and synthetic data to infer physical parameters of the system.

General relativistic magnetohydrodynamics (GRMHD) simulations are a mainstay in analysis of black hole accretion flows (e.g., Gammie et al. 2003; Del Zanna et al. 2007; Mignone et al. 2007; Narayan et al. 2012; Sadowski et al. 2013b; White et al. 2016; Porth et al. 2017; Liska et al. 2018), and many general relativistic ray-tracing codes have been developed for both imaging (e.g., Zane et al. 1996; Fuerst & Wu 2004; Noble et al. 2007; Psaltis & Johannsen 2012; Younsi et al. 2012; Chan et al. 2013; Dexter 2016; Chan et al. 2018; Mościbrodzka & Gammie 2018; Pihajoki et al. 2018; Kawashima et al. 2021, building on earlier work by Mihalas & Mihalas 1984) and the generation of spectra (e.g., Dolence et al. 2009; Mościbrodzka et al. 2009; Dolence et al. 2012; Schnittman & Krolik 2013; Zhang et al. 2019; Mościbrodzka 2020; Kawashima et al. 2021).

In this article, we describe the PATOKA modeling pipeline, which was used for the Illinois contribution of fluid accretion models, images, and spectra to the EHT library. The details of this pipeline are messy, like the details of the pipelines used for hydrocarbon extraction around downstate Illinois's Patoka oil terminal (see Figure 2), for which our pipeline is named. In this paper, we describe some of the messy details, including differences between the published code descriptions and the versions that were used in generating the library. In Section 2, we provide a summary of the theoretical underpinnings of the accretion model and radiation physics. In Section 3, we review the details of the numerical codes used to perform the simulations and ray tracing. In Section 4, we present some results that were generated with the pipeline. We end with a discussion of future directions in Section 5.

Figure 2. Refer to the following caption and surrounding text.

Figure 2. Photograph of an oil field near Patoka in Marion County, Illinois published in 1940 February. The PATOKA pipeline is named for the Patoka Oil Terminal in Patoka, Illinois, which connects several oil pipelines in the second Petroleum Administration for Defense District. Photo credit: Arthur Rothstein/Farm Security Administration, LC-DIG-fsa-8a13148.

Standard image High-resolution image

2. Theoretical Background

We begin with an overview of the theoretical considerations associated with the library generation. We cover the parameters of the black hole accretion model and review the equations governing radiative transfer in hot astrophysical plasmas.

2.1. Black Hole Accretion

The supermassive black holes at the centers of the M87 galaxy and the Galactic Center (Sgr A*) are associated with compact radio sources (see, e.g., Lynden-Bell 1969; Balick & Brown 1974; Doeleman et al. 2008; Fish et al. 2011) and are the primary targets for the EHT. M87 and Sgr A* are usually described as low-luminosity active galactic nuclei (see Greene & Ho 2007; Yuan & Narayan 2014 for reviews), so the accretion onto their central black holes is expected to proceed at a low rate and therefore be radiatively inefficient (RIAF for radiatively inefficient accretion flow or ADAF for advection dominated accretion flow; see Ichimaru 1977; Narayan & Yi 1994, 1995; Quataert & Narayan 1999; Yuan et al. 2003). We express the mass accretion rate in terms of the Eddington accretion rate , where c is the speed of light, the nominal accretion efficiency is chosen to be η = 0.1, and the Eddington luminosity is LEdd = 4π GMmp c/σT , where G is Newton's gravitational constant, the mass of a proton is mp , the Thomson cross section is σT , and the mass of the central black hole is M.

For orderly accretion to occur, the gravitational binding energy of the system must be dissipated as matter falls from large to small radii. In RIAF models, the gravitational binding energy of the flow is converted into heat and either advected across the event horizon or lost in mechanical outflows. RIAFs are typically modeled as geometrically thick disks of relativistically hot, magnetized plasma (Rees et al. 1982; Narayan & Yi 1995; Narayan et al. 1995; Reynolds et al. 1996; Yuan et al. 2002; Di Matteo et al. 2003; Yuan & Narayan 2014). M87 also supports a relativistic jet that extends to kiloparsec scales and has an estimated jet power of ≈ 1042–1045 erg s−1 (Stawarz et al. 2006; de Gasperin et al. 2012; Prieto et al. 2016). Whether or not Sgr A* also launches a jet is less certain (e.g., Falcke & Markoff 2000; Yusef-Zadeh et al. 2006; Markoff et al. 2007; Falcke et al. 2009; Brinkerink et al. 2015; Issaoun et al. 2019; Brinkerink et al. 2021).

The Coulomb mean free paths of the ions and electrons in typical RIAF models are much larger than the typical length scales in the disks (see Mahadevan & Quataert 1997), so the accreting plasma should be collisionless, suggesting that nonideal processes may be important. Particle-in-cell (PIC) simulations (e.g., Kunz et al. 2014; Riquelme et al. 2015; Sironi & Narayan 2015), measurements of the solar wind, and phenomena like collisionless shocks infer that instabilities may increase the effective particle collision rate, allowing nonideal effects to be treated as a small correction (see Foucart et al. 2016 and references therein). Extended GRMHD simulations that include these small corrections have shown that the low collisionality does not significantly alter the time-averaged structure of the flow (e.g., Foucart et al. 2017).

Kinetic plasma physics phenomena may ultimately drive the system out of equilibrium in ways that our methods do not account for; however, global general relativistic kinetic simulations of ion–electron plasmas in black hole accretion flows at realistic accretion rates are prohibitively expensive and have never been run. Either way, differences between electron and ion heating mechanisms and the absence of an efficient coupling mechanism between the two species implies at least that the electrons and ions in the flows are likely described by different temperatures (see Shapiro et al. 1976; Mahadevan & Quataert 1997; Quataert 1998; Ressler et al. 2015; Zhdankin et al. 2021).

Analytic models are useful in understanding the broad characteristics of accretion flows (e.g., Novikov & Thorne 1973; Shakura & Sunyaev 1973; Abramowicz et al. 1978; Kozlowski et al. 1978; Narayan & Yi 1994; Yuan et al. 2003; Broderick & Loeb 2006; Komissarov 2006; Broderick et al. 2011; Abolmasov 2018), but they provide an incomplete picture of the dynamics. For example, the analytic treatments often propose or rely on thin-disk models, do not explicitly include magnetic fields, and do not account for time dependence in the flows; yet the development of turbulence and shocks make an analytic treatment particularly challenging. 18

Since a large part of the millimeter emission produced in a RIAF flow is likely produced near the horizon (e.g., Mościbrodzka et al. 2009), a full general relativistic treatment of the system is necessary to adequately recover the precise features of the model. GRMHD simulations serve as a practical method for studying the plasma dynamics in detail. GRMHD simulations evolve the distribution of plasma and electromagnetic energy through time and produce a time-series description of the plasma state. Global GRMHD simulations have been widely used to study RIAF accretion and have been shown to reproduce the jet powers of the Blandford & Znajek (1977) model (e.g., Koide et al. 1999; De Villiers & Hawley 2003; Gammie et al. 2003; McKinney & Gammie 2004; Narayan et al. 2012). Numerical accretion models have been shown to be consistent with various observational constraints provided by jet powers, variability, and the submillimeter spectrum (e.g., Noble et al. 2007; Mościbrodzka et al. 2009; Dexter et al. 2010, 2012; Shcherbakov et al. 2012; Dexter & Fragile 2013; Mościbrodzka et al. 2014; Dexter et al. 2020). The consistency between simulations and the recent horizon-scale images of M87 published by the EHT further supports the RIAF/ADAF thick-disk picture (EHTC V; EHTC VIII).

At low accretion rates, radiative cooling has a negligible impact on the fluid dynamics (Dibi et al. 2012; Ryan et al. 2017; and see also Yoon et al. 2020), so the evolution is invariant under rescalings of both the metric length GM/c2 and the mass density. Ideal GRMHD simulations thus only introduce the following physical parameters: the angular momentum of the system and the magnetic flux near the event horizon.

2.1.1. Angular Momenta and Tilted Disks

The angular momentum of the central black hole J is often expressed in terms of the dimensionless black hole spin parameter a*Jc/GM2, with . The angular momentum vectors of the black hole and the accretion flow may be misaligned; the angular separation between the two axes is called tilt. Although there are plausible scenarios that produce accretion flows with zero tilt, there is at present no way of rejecting models with strong or even maximal (180°, i.e., retrograde) tilt. In fact, retrograde accretion is not particularly unexpected; it exists and has been resolved in the context of stellar disks (see Young et al. 2020 and references therein).

Systems with nonzero tilt may have clear observational signatures. Fragile & Blaes (2008), Dexter & Fragile (2013), Shiokawa (2013), White et al. (2019), Chatterjee et al. (2020) found that tilted accretion flows produce a two-arm shock that would significantly alter the morphology of radio images and the shape of the spectrum. Some simulations suggest that the inner regions of tilted disks ultimately align with the central black hole via an analog of the Bardeen & Petterson (1975) effect (McKinney et al. 2013; Liska et al. 2018).

In our galaxy, there is strong evidence that Sgr A* was recently in an active period of higher accretion (Totani 2006; Ponti et al. 2013), and other observational studies and simulations report evidence of clumpiness within a few parsecs of Sgr A* in the circumnuclear disk (Montero-Castaño et al. 2009; Blank et al. 2016). Simulated mass feeding in the Galactic Center via magnetized stellar winds by Ressler et al. (2018, 2020a, 2020b) found that a wide variety of different initial configurations led to the development of similar nondisk-like accretion flows with strongly poloidal horizon-scale magnetic-field configurations. Evidently, the angular momentum of the boundary condition that supplies the inner accretion flow with mass is likely to be variable on galactic timescales (e.g., Cuadra et al. 2006).

More generally, a nontrivial time-dependent boundary feeding condition is not unlikely, as there is a large discrepancy between the long timescale that governs the changes to properties of the black hole and the shorter ones over which the local environment of the black hole (i.e., the environment that governs the feeding) changes. Thus, the bulk angular momentum of the accreting matter (far from the horizon) may exhibit a time-dependent tilt with respect to the spin angular momentum of the central black hole.

Time-dependent tilts have been invoked to explain the wobble observed in some relativistic jets (see Natarajan & Armitage 1999). These arguments find extra support in the analysis of episodic jet variations and spectral analyses of the broad Fe line (e.g., Hjellming & Rupen 1995; Rout et al. 2020). Theoretical studies that consider the initial spin distribution of supermassive black holes (De Luca et al. 2019) and simulate black hole growth via accretion and mergers in hierarchical galaxy formation models also often favor random alignment between the black holes and their disks (see, e.g., Volonteri et al. 2005).

The prograde/retrograde dichotomy has also been used to explain differences between radio-loud and radio-quiet active galactic nuclei sources (Garofalo 2009; Garofalo et al. 2010; but see Tchekhovskoy & McKinney 2012). Some observational studies have claimed detection of retrograde accretion systems: Morningstar et al. (2014), Chen et al. (2016), and Mikhailov et al. (2018) presented criteria for classifying systems as retrograde and singled out several known systems based on emission models and low estimated Eddington ratios.

2.1.2. Magnetic Flux

In ideal MHD, accreting plasma carries magnetic-field lines with it as it falls onto a central black hole. If the magnetic polarity is constant over time, then field lines will accumulate and increase the magnetic pressure near the hole, eventually saturating when the magnetic pressure is large enough to counterbalance the inward ram pressure of the accretion flow. The amount of magnetic flux threading the event horizon thus qualitatively divides accretion flows into two categories: the magnetically arrested disk (MAD) state (Bisnovatyi-Kogan & Ruzmaikin 1974; Igumenshchev et al. 2003; Narayan et al. 2003), in which the magnetic pressure is large, and the large-scale component of the magnetic field is dynamically important; and the alternate standard and normal evolution (SANE) state (Narayan et al. 2012; Sadowski et al. 2013a).

The magnetic flux threading one hemisphere of the black hole horizon ∫dA · B is

where g is the determinant of the covariant metric, and the Hodge dual of the electromagnetic Faraday field tensor Fα β is

Here, we have used the Levi–Civita tensor , which differs from the permutation symbol by a factor of and sign. 19 We compute a proxy for the magnetic flux as

Equation (3) approximates Equation (1) when the field is well-organized on the horizon. The flux ΦBH is conventionally normalized by

where rg GM/c2 is the gravitational radius of the hole. In the rationalized Lorentz–Heaviside units used in GRMHD simulations, ϕ saturates at ϕc ≈ 15 (see Tchekhovskoy et al. 2011; and Porth et al. 2019). 20 It is not at present possible to calculate what value of ϕ will result from an arbitrary field configuration without evolving the flow.

MAD accretion is spatially and temporally variable and tends to proceed in isolated, thin plasma streams that begin far from the hole. MAD flows are often punctuated by violent magnetic bubble eruption events that release excess trapped magnetic flux. Although events where magnetic flux escapes from the black hole are not fully understood, they have been interpreted as a magnetically driven convection between the disk and the hole (Marshall et al. 2018).

Simulations suggest that the MAD versus SANE dichotomy is observationally encoded in signatures of polarization, variability, and the details of the jet–disk connection (e.g., Gold et al. 2017; Palumbo et al. 2020; Wong et al. 2021a). Analysis of observations of radio-loud active galaxies and their jets are consistent with the MAD picture (e.g., Zamaninasab et al. 2014), and recent analysis of EHT data also suggests that M87 accretion is MAD (EHTC V, EHTC VIII).

2.1.3. Electron Temperature and Radiative Effects

GRMHD simulations typically assume that the dynamics of the accreting plasma can be recovered by treating the flow as a single, thermal fluid; they therefore track only the internal energy (or temperature) of the bulk plasma. In the real world, however, the plasma typically has a long Coulomb mean free path, so the electron temperature Te will not necessarily be equal to the fluid temperature (Shapiro et al. 1976; Rees et al. 1982; Narayan & Yi 1995). Conventionally, the electron temperature is assigned after the fact according to a model that assigns Te as a function of the local fluid parameters, such as the internal energy of the fluid and the local ratio of gas pressure to magnetic pressure, plasma βPgas/Pmag (e.g., Goldston et al. 2005; Mościbrodzka et al. 2009; Shcherbakov et al. 2012; Mościbrodzka & Falcke 2013; Mościbrodzka et al. 2014, 2016—see Anantua et al. 2020b for a comparison of several different models). We discuss this approximation in more detail in Section 2.3.3.

GRMHD simulations that recover the turbulent dynamics rely on implicit, large eddy simulations (see Boris 1990; Boris et al. 1992; Miesch et al. 2015; Grete 2017), wherein it is assumed that the numerical dissipation at the grid scale emulates turbulence and dissipation below the grid scale and thus can be used to infer the electron (and ion) heating rates. Ressler et al. (2015), Sadowski et al. (2017) described a method to track numerical dissipation in relativistic magnetohydrodynamics simulations for this purpose. Alternative descriptions of (nonideal) fluid mechanics have been proposed with direct treatments of viscosity and resistivity; such methods explicitly capture the dissipation.

In the systems we consider, plasma is typically energized through a combination of collisionless shocks, magnetic reconnection, and relativistic turbulence. For shocks, interactions between the charged particles and the magnetic field set the length scale (i.e., via the gyroradius), and the strongly nonequilibrium perturbations produce irreversible heating of the electrons, e.g., via first-order Fermi acceleration. In regions where reconnection occurs, the electrodynamics within the resultant current sheets do work on the plasma and transform magnetic energy into heat. Finally, turbulent cascades carry energy from large scales to small ones, ultimately energizing particles in the plasma at the smallest of scales. The energization mechanisms can produce anisotropic and highly nonthermal particle energy distributions that may differ between the different components of the plasma. These different energization mechanisms have been studied in different regimes and for different plasma compositions (e.g., Numata & Loureiro 2015; Sironi & Narayan 2015; Sironi 2015; Boldyrev & Loureiro 2017; Loureiro & Boldyrev 2017; Mallet et al. 2017; Comisso & Sironi 2018; Shay et al. 2018; Comisso & Sironi 2021; Nättilä & Beloborodov 2021; Zhdankin et al. 2021; Ripperda et al. 2022).

Our pipeline is ultimately agnostic to the details of the energization mechanism. Given a model for how the dissipated energy is partitioned into electrons and ions (e.g., Sharma et al. 2007; Howes 2010; Werner et al. 2016; Rowan et al. 2017; Werner et al. 2018; Zhdankin et al. 2018; Kawazura et al. 2019; Rowan et al. 2019; Zhdankin et al. 2019; Kawazura et al. 2020, 2021, 2021) and an assumption about the particle distribution functions (e.g., the electrons are well described by a thermal distribution), one can independently track and evolve the electron temperature.

At low accretion rates , radiative losses are too slow to materially alter the internal energy of the flow (see, e.g., Sharma et al. 2007; Dibi et al. 2012) and nonradiative GRMHD simulations sufficiently recover the dynamics. As approaches 10−5, however, the synchrotron emission and Compton upscattering become important to the flow dynamics, and it is necessary to use a radiation-GRMHD scheme like KORAL (Sadowski et al. 2013b, implementing a gray M1-closure scheme for the radiation), ebhlight (Ryan et al. 2015, 2019, using a full Monte Carlo treatment of the radiation field), or MOCMC (Ryan & Dolence 2020, using an adaptively refined Monte Carlo treatment of the radiation field). We note that M1-closure schemes may be inadequate where the radiation field is not symmetric with respect to the mean flux, for example in regions where there are multiple, distinct radiation sources. The transition regime in which radiation becomes important has been studied by, e.g., Fragile & Meier (2009), Wu et al. (2016), Sadowski & Gaspari (2017), Sadowski et al. (2017), Ryan et al. (2017).

2.2. Radiative Transfer

Simulated 230 GHz images of supermassive black holes at sufficiently low accretion rates are often dominated by a ring-like feature with high brightness temperature Tb . In terms of observational parameters, Tb = S λ2/2kΩ, where S is an observed flux density describing power received per unit area per unit frequency, λ is the observing wavelength, k is the Boltzmann constant, and Ω is the solid angle of the observing element. In physical terms, Tb is the temperature of a blackbody that would reproduce an observed intensity at a given frequency. The location of the ring is broadly consistent with the critical curve boundary that separates the null geodesics that terminate on the black hole event horizon from those that do not.

Geodesics close to but outside of the critical curve are also lensed by the hole, and the set of all geodesics that are lensed enough to complete n half-orbits around the hole defines the nth photon subring (Johnson et al. 2020; see also Ohanian 1987); as n , the subrings approach the critical curve. This half-orbit definition is nearly equivalent to a definition whereby subrings are indexed by the number of times they pass through the plane perpendicular to the symmetry axis of the spacetime, i.e., the disk midplane in aligned accretion flows, though differences may arise depending on whether the orbit is defined with respect to θ, ϕ, or z. Subrings are geometric regions on a black hole image and are therefore, like the black hole shadow and critical curve, formally unrelated to the visual appearance of the image, which depends on the properties of the underlying emission. Each subring corresponds to a relatively delayed, demagnified image of the universe. See Section 4.2 and Figure 7 in particular for a detailed example of the subring decomposition.

In the astrophysical scenarios we consider, the emission drops off sharply with distance from the black hole, so the images are dominated by emission produced near the hole. The emission near the black hole will thus produce a bright region in the direct n = 0 image near the shadow. Since the image displayed by each next subring corresponds to a demagnified image of the previous one, the full image, which is produced by summing over each of the subring images, will display a sharp, logarithmically narrowing feature in intensity whose peak converges to the critical curve in the limit that the optical depth of the system is negligible. The feature seen in the full, composite image is known as the "photon ring" (Bardeen 1973; Luminet 1979; Gralla & Lupsasca 2020; Johnson et al. 2020), and the region it encircles is the black hole shadow.

The location of the critical curve and the shape of the black hole shadow are controlled by the spacetime geometry. Strategies for constraining properties of the spacetime by measuring the critical curve via the shape of the photon ring or due to autocorrelations and subring image delays have been proposed (e.g., Falcke et al. 2000; Takahashi 2004; Bambi & Freese 2009; Hioki & Maeda 2009; Amarilla et al. 2010; Amarilla & Eiroa 2013; Tsukamoto et al. 2014; Younsi et al. 2016; Mizuno et al. 2018; Wielgus et al. 2020a; Johnson et al. 2020; Medeiros et al. 2020; Olivares et al. 2020; Chesler et al. 2021; Hadar et al. 2021; Wong 2021).

Particle trajectories in curved spacetimes are determined by the geodesic equations

where Γ is a Christoffel symbol, λ the affine parameter, and kα the photon wavevector. The frequency ν of the photon as measured in a frame with four-velocity uα is in general given by 2π ν = −kα uα , although the scale factor in this relation depends on the choice of normalization. In PATOKA, we normalize the photon wavevector kα such that the frequency is given by , where h is Planck's constant.

The properties of the radiation field are often quantified in terms of the frequency- and angle-dependent specific intensity of the light, Iν = the amount of energy per area per solid angle per time per frequency. 21 Specific intensity changes along the geodesics according to interactions with the local matter in four ways: emission can increase the intensity of the light, scattering into the geodesic can increase the intensity of light, the intensity can decrease either due to absorption or due to scattering out of the geodesic (collectively called extinction), and the local plasma properties can mix linear and circular polarizations. The extinction term depends on the local intensity of the incoming light, the path length along the geodesic ds, and an extinction coefficient, which itself depends on the local plasma parameters and, e.g., for synchrotron absorption, is a function of the angle the light ray makes with the magnetic field.

2.2.1. Polarization

In the supermassive black hole accretion systems targeted by the EHT, photons are emitted primarily by the synchrotron process (see, e.g., Yuan & Narayan 2014), when electrons in the accreting plasma are accelerated by magnetic fields. The local orientation of the magnetic field sets a preferred direction for the electromagnetic radiation, and thus synchrotron emission is, in general, polarized according to the direction of the magnetic field.

The full description of polarized light is often given in terms of the specific intensities of the Stokes parameters, which can be written in terms of the Stokes vector . Here, Iν is the familiar total specific intensity of the light, Qν and Uν represent the linearly polarized specific intensities, and Vν represents the circularly polarized one. Not all light must be polarized, so . The electric-vector position angle (EVPA), 22 denoted χ, is

Physically, χ describes the angle the electric field oscillations make with respect to a fiducial direction. The orientation of the electric field exhibits a π-fold symmetry, so the representation of linear polarization χ is a pseudovector. The factor of one-half is due to this symmetry and can be seen in, e.g., the quadratic relationship between the Stokes parameters and the electric field vector E.

The polarized intensities are governed by the polarized radiative transfer equation. In order to produce polarized millimeter radio images, we need only treat the emission and absorption processes, since scattering contributes negligibly at the relevant frequencies and plasma parameters. Neglecting scattering, the polarized radiative transfer equation is

where the emissivities jν , absorptivities αν , and rotativities ρν are frame-dependent quantities (Chandrasekhar 1960). We work in a frame determined by the orientations of the wavevector and the magnetic field; this choice leads to αU = ρU = 0.

As polarized light travels through a magnetized plasma, several effects may cause χ to rotate. The first is due to parallel transport of the polarization vector in the curved geometry. The magnitude of this effect is a pure function of the underlying spacetime (see Pihajoki et al. 2018; and also see Gelles et al. 2021a for a quantitative study of the magnitude in the context of EHT-like sources). The others are due to interactions between the propagating light and the local plasma and are due to emission, absorption, and Faraday rotation and conversion. For example, in a magnetized plasma, the dielectric constant is a tensorial quantity, so components of the light with different polarizations propagate at different speeds. This magnetically induced birefringence produces a characteristic Faraday rotation of χ that is a function of the plasma properties along the line of sight

where ρV is the Faraday rotation coefficient. In the second line, we have substituted a high-frequency expression for ρV for the sub-to-mildly relativistic case; here, ne is the electron number density, is the dimensionless electron temperature, B∣∣ is the component of the magnetic field along the line of sight, and for the Θe ≫ 1 case and ≈1 otherwise (e.g., Rybicki & Lightman 1979; Quataert & Gruzinov 2000; see also Shcherbakov 2008 for a discussion). Evidently the observed polarization pattern is related to the structure of the magnetic field in a complicated way.

Broderick & Blandford (2004) described a method to write the transport equation in terms of quantities that account for rotations of the observer frame along the line of sight. To provide a manifestly covariant description, Gammie & Leung (2012) produced an alternative formulation of the polarized transport equation in terms of a Hermitian coherency tensor Nα β , which can be related to the Stokes parameters. The polarized transfer equation can then be written

where kμ μ is the derivative along the geodesic, Jα β is an emissivity tensor, and Hα β γ δ is a tensor that accounts for absorption and generalized Faraday rotation. Using this description to solve the polarized transport equation amounts to evaluating the nonrelativistic emissivities, absorptivities, and rotativities in the plasma frame, constructing Jα β and Hα β γ δ , and solving Equation (11). In particular, Jα β and Hα β γ δ are constructed in the tetrad basis defined by the fluid four-velocity, the photon wavevector, and the orientation of the local magnetic field.

In some cases, we may solve an approximate form of the transfer equation that does not account for polarization. Beginning with Equation (8) and suppressing the polarized Stokes Q, U, V terms produces

This approximation is sometimes called unpolarized transport, although it accounts for the total intensity of the light rather than only the unpolarized component. To recast the unpolarized transfer equation in a manifestly covariant form, notice that particle number dN and phase space volume, d3 x d3 p, are conserved. Suppressing factors of h and c, the phase space volume for a parcel of radiation at a frequency ν can be written d3 x d3 p = d A d t ν2 d ν dΩ. Thus, the quantity

is invariant. Similarly, jν /ν2 and ν αν are invariant, so we can rewrite Equation (12) in the covariant form

We often write the invariant intensity .

2.2.2. Scattering

Photons undergo Compton scattering as they travel through the plasma in the accretion flow. Compton scattering in the regime of interest here typically increases the energy of the scattered photon. Collectively, these scattering events alter the spectrum; this is called Comptonization. The Compton y parameter gauges the significance of Comptonization:

where in the second line we have given an approximate form for the regime relevant to EHT sources. For a system of size l, the Thomson scattering optical depth is approximately

The rate of interaction between a photon with wavevector kμ and massive particles in a distribution function d nm /d3 p is

where in the second line, we have used a "hot cross section" (see Appendix III of Canfield et al. 1987)

where we rewrite the dot product in terms of the particle speed in the plasma frame βm and the cosine of the angle between the particle momentum and the photon momentum μm . For electron scattering, the cross section is the Klein–Nishina cross section, which is most conveniently written

where the energy of the photon in the electron rest frame is epsilon = −pμ kμ /me .

We use a Monte Carlo approach to process Compton scattering events. First, based on the incident photon, we use rejection sampling to draw an electron from the local electron distribution function (eDF), then we use rejection sampling again to sample the differential scattering cross section for that electron and incident photon in order to pick the energy and scattering angle that sets the wavevector for the scattered (outgoing) photon. We use the Thomson differential cross section for low-energy photons; otherwise we use the Klein–Nishina differential cross section.

2.3. Connecting GRMHD and Radiative Transfer

We produce the above-described electromagnetic observables (images and spectra) from the fluid simulations (which output a description of the fluid state) by performing radiative transfer calculations. In order to use the GRMHD output, the simulations must be scaled (assigned cgs units) and oriented, and an eDF must be specified.

2.3.1. Scaling and Orienting GRMHD Output

Unlike in GRMHD, the equations of radiative transfer are not scale invariant, so the numerical fluid data must be translated into physical units in order to perform the ray tracing. General relativistic radiative transfer (GRRT) introduces two scales: a length scale and a density scale. The length scale determines the absolute size of the accretion system and is often written in terms of the mass of the central black hole . The length scale sets the characteristic timescale, .

The density scale provides units to the fluid rest-mass density, internal energy, and magnetic-field strength while respecting the constraint that magnetization σ = b2/ρ and plasma β are independent of the choice of units. The density scale is specified as the ratio of a mass scale to the volume scale determined by . 23 The density scale is chosen so that the simulated images match an observational constraint: the observed flux density must be correct. Since the simulated images have a limited field of view, it is possible for the flows with large optical depths to produce diffuse, extended emission with image-integrated fluxes that are only consistent with the observation because of the limited image size (see Appendix D). The correct value of is typically found via a root-finding procedure (see Appendix D for more detail). When possible, the resulting is compared to the predictions of single-zone models (see EHTC V; EHTC VIII) and the observational estimates based on rotation measure (e.g., Bower et al. 2003; Marrone et al. 2006, 2007; Kuo et al. 2014).

With these choices made, the state of the system is determined. Note that the choice of length and density scales is nearly degenerate with the choice of observer-to-source distance . We discuss this point in Appendix D.

Finally, the accretion flow must be oriented properly with respect to the camera location on Earth. This amounts to choosing two of the three Euler angles: an inclination angle i between the spin axis of the black hole and the elevation of the camera, and a position angle for the camera (the orientation that the projected spin axis makes on the image plane). Since we consider systems with zero tilt, our flows are statistically axisymmetric, so the third angle is negligible because there is no preferred azimuthal orientation of the disk about the spin axis of the system. The inclination and position angle are free parameters, although other observational information about the system, such as the orientation of a jet, may constrain both.

2.3.2. Assigning Electron Temperatures

Although our GRMHD simulations often make the approximation that the plasma is thermal and described by a single temperature, we account for the likely collisionless nature of the flow by allowing electron temperatures to deviate from the ion temperatures. In PATOKA, when not using the entropy tracking procedure of Ressler et al. (2015), we assign electron temperatures according to the prescription of (Mościbrodzka et al. 2016; see also Mościbrodzka et al. 2017; EHTC V), wherein the ion-to-electron temperature ratio is determined by the local plasma β. This prescription is motivated by the idea that electron heating may be relatively stronger in a strongly magnetized plasma (e.g., Quataert 1998; Quataert & Gruzinov 1999), which can produce a sigmoidal dependence of the temperature ratio on β.

In Mościbrodzka et al. (2016), the temperature ratio is parameterized by rlow, rhigh, and βcrit, with

where . We typically adopt 1 ≲ rlow ≲ 10 and 1 ≲ rhigh ≲ 160 with rhighrlow, based on more detailed studies of electron heating.

To recover the electron temperature from the total fluid internal energy, we partition the fluid internal energy into two components, associated with the electrons and the ions. 24 Assuming an ideal gas equation of state, the energies associated with the ions and electrons will be

where and are the adiabatic indices of the ions and electrons respectively. The ion and electron number densities are related to the total mass density by

where y and z are the number of electrons and nucleons per unionized atom, respectively. 25 The electron temperature is thus

The ions are typically nonrelativistic, so , and that the electrons are relativistic, so . Since we assign electron temperatures after the fact, we assume that , and are constant across the simulation domain. This treatment is not entirely self-consistent, since the adiabatic index of the fluid should change depending on local contributions from the ion and electron fluid components (see Mignone & McKinney 2007; Choi & Wiita 2010; Mizuno 2013; Shiokawa 2013; Sadowski et al. 2017 for self-consistent treatments).

2.3.3. Approximations and Pathologies

GRMHD schemes are not robust in regions with high magnetization σ ≫ 1 or where plasma β ≪ 1. In order to ensure numerical stability, σ and β are computed in each cell of the simulation domain for each time step, and mass or internal energy is injected into the simulation zones to ensure that neither σ nor β−1 exceed some preset ceilings (see Section 4.1 for more detail). We have typically varied the values of the ceilings across different trial-run simulations to ensure that they have minimal impact on the evolution of the flow.

In black hole accretion simulations, the ceilings are generally triggered in the highly magnetized funnel regions near the poles, where we expect very low densities . Although truncation errors dominate the evolution of the fluid internal energy in regions where the ceilings are activated, we expect that the plasma density in the jet will be low enough that it has negligible effect on the radiation. We therefore set the particle number density to zero in regions with σ > σcut. We typically set σcut = 1. We also zero the particle number density within a few degrees of the poles, since some treatments of the polar boundary condition during the fluid evolution may cause σ to artificially drop below unity there. The effect of the σ cutoff is considered in Chael et al. (2019), EHTC VIII, and Appendix C. In reality, the jet may be populated by an "injected" electron–positron pair plasma produced through photon–photon interactions in either pair cascades or pair drizzle (see Section 5.2), although the plasma injection rate (and composition) will, in general, differ from the plasma injected due to floors.

In images produced at edge-on inclinations, much of the emission travels through the disk on its way to the camera. At intermediate-to-large radii, the disk plasma has not had time to reach a turbulent quasi-equilibrium state whose average properties evolve only on the accretion timescale, so its state is largely a relic of the initial conditions. Thus, care must be taken that the polarization is not determined by Faraday rotation in the unequilibrated outer disk. For the same reason, care must be taken that the polarization is not dependent on the size of the simulation domain (see, e.g., Ricarte et al. 2020).

Although we have assumed that the plasma is a pure electron–ion plasma, it is likely that electron–positron pairs populate parts of the domain. Pair plasmas produce different emission signatures (Anantua et al. 2020a; Emami et al. 2021), especially in polarization and in the regions of the flow where the background plasma density is low. Furthermore, the presence of a pair plasma in low-density regions—like the jet—may influence the structure and dynamics of those regions.

Finally, we use an ideal GRMHD scheme to evolve the accretion flow, which implicitly assumes that the plasma is in equilibrium and that the conductivity of the plasma is infinite, such that there are no unscreened electric fields. In fact, neither of these assumptions is likely to be valid globally across our simulation domain. It is unclear how strongly the ideal fluid assumption influences our results, but it is likely that the implicit assumption that the eDF is thermal and isotropic would have a significant effect on simulated observables. For example, kinetic simulations of pair plasmas infer strongly anisotropic, nonthermal eDFs (Comisso & Sironi 2021; Nättilä & Beloborodov 2021), which can have significant consequences for polarimetric observables and even the distribution of emission across the domain, since emission is pitch-angle dependent.

There is yet a limited consensus between different kinetic simulations, so it is not clear which model to use for subgrid phenomena and deviations from the ideal fluid assumption. We choose to address these uncertainties in part through our parametric treatment of the electron temperature assignment, which affords us great flexibility in translating simulation quantities into an eDF. Although our prescription is motivated, it cannot account for all of the potential deviations between our fluid treatment and a full kinetic treatment of the plasma distribution function. Understanding the differences between results from the ideal (and extended) GRMHD fluid treatments and kinetic simulations of ion–electron–positron plasmas is important; however, it requires a detailed study that is beyond the scope of this paper.

3. Code Detail

We now describe the PATOKA pipeline, in which simulated observables are generated by ray-tracing snapshots produced from GRMHD simulations of the black hole accretion flows. These two stages, GRMHD and GRRT, are separated for computational efficiency, since GRMHD simulations are costly and multiple radiation models can be applied to a single fluid snapshot without rerunning the fluid simulation. We describe the details of the three codes we use, placing an emphasis on the differences between the versions we use and the code as described upon release. All codes used in PATOKA compute metric derivatives via numerical differences, and therefore nearly all coordinate dependence is encapsulated in a single line of code that specifies the line element.

3.1. General Relativistic Magnetohydrodynamics

We use the iharm code (Gammie et al. 2003; Noble et al. 2006, 2009; Prather et al. 2021) to integrate the equations of ideal GRMHD. Porth et al. (2019) provides a comparison of contemporary GRMHD codes in the context of SANE accretion flows. iharm is a conservative second-order explicit shock-capturing finite-volume code for arbitrary stationary spacetimes. iharm is a descendant of harm2d, and based on the harm scheme of Gammie et al. (2003).

iharm is designed to vectorize efficiently and achieves good performance and scaling on the systems used for this study. Validation and scaling tests are described in Prather et al. (2021) and show second-order convergence on a suite of test problems. The code is publicly available. 26

The governing equations of ideal GRMHD take the form of a hyperbolic system of conservation laws. In conservation form and written in a coordinate basis, the equations are

with the constraint

where the plasma is defined by its rest-mass density ρ0, its four-velocity uμ , and bμ is the magnetic-field four-vector following McKinney & Gammie (2004; see also Lichnerowicz 1967). Here, Γ is a Christoffel symbol, and i and j denote spatial coordinates. In Equations (30) and 31, we use the ideal MHD condition 27 uμ Fμ ν = 0 to express Fμ ν in terms of Bi ≡ ⋆Fit for notational simplicity. Notice that Bi is the magnetic field as measured in the coordinate frame. The stress–energy tensor contains contributions from both the fluid and the electromagnetic field

where u is the internal energy of the fluid, and the fluid pressure P is related to its internal energy through a constant adiabatic index with . The four-velocity of the fluid is encoded in the three numbers , whose definitions were chosen to improve numerical stability and are given in Appendix B.

The equations are solved over a logically Cartesian, 3D grid in arbitrary coordinates. Eight primitive variables, 28 ρ, u, , and Bi , are stored at the center of each grid zone and evolved forward in coordinate time (see Appendix B for more detail). Fluxes are computed using the local Lax–Friedrichs method (Rusanov 1962), and the divergence constraint of Equation (31) is enforced using the Flux-CT scheme described in Tóth (2000; see also Evans & Hawley 1988). The fifth-order WENO5 scheme (Jiang & Shu 1996) is typically used for spatial reconstruction. The code is parallelized across contiguous domain chunks using the Message Passing Interface (Message Passing Interface Forum 1994) and within each chunk using OpenMP (Dagum & Menon 1998). iharm imposes a constant floor on plasma β and a "geometric" floor constraint on ρ that varies with spatial coordinate; it imposes ceilings on σ, Θe , the fluid velocity measured with respect to the coordinate frame, and optionally the logarithm of the fluid entropy (see Section 4.1 for more detail).

Simulations are typically run in the augmented versions of the horizon-penetrating modified Kerr–Schild (MKS) coordinates introduced in Gammie et al. (2003; see Appendix F for more detail). In MKS, the three spatial coordinates x1, x2, x3 are direct functions of radius, latitude, and azimuth respectively. The inner radial boundary of the simulation is chosen to ensure that ≳5 zones lie within the event horizon. The outer edge of the boundary is typically chosen so that the torus lies comfortably within the simulation domain. We use outflow boundary conditions along the two radial boundaries, a periodic boundary condition in the azimuthal direction, and a pseudo-reflecting boundary condition at the two poles that mirrors the latitudinal 2-components of the magnetic field and fluid velocity across the one-dimensional border.

iharm supports the subgrid electron heating scheme of Ressler et al. (2015) to track the local dissipation rate and then apportion some fraction of the dissipated energy to the electrons using any prescription that can be specified in terms of the fluid state. When the electron thermodynamics module is active, the code outputs two extra variables—functions of the total fluid entropy and the electron entropy—in addition to the principal eight primitive variables listed above.

We initialize the fluid sector of our accretion disk simulations with a (Fishbone & Moncrief 1976, hereafter FM) torus. We perturb the initial conditions by sending uu + δ u, with drawn random uniformly per the simulation zone. Typically ujitter ≈ 0.1, although this value has been varied substantially between different simulations. This perturbation helps break the model symmetries and seed instabilities, including the magneto-rotational instability (Balbus & Hawley 1991), and begin accretion onto the hole. 29 The FM torus has two parameters: the midplane radius of the inner edge r = rin and the radius of the maximum pressure at . See Appendix A for more detail.

iharm is configured to run both SANE and MAD simulations and can be extended to support other initial conditions. Although the steady-state magnetic flux is not trivially related to the initial conditions of the simulation, we have implemented configurations that have been identified in previous work and shown to produce either a SANE or a MAD flow. The initial condition for the magnetic field is determined by a prescribed, axisymmetric electromagnetic vector potential Aϕ (r, ϕ), which is computed at the simulation-zone corners. The magnetic field is then calculated from the curl of the vector potential using a finite different operator that is compatible with the constrained transport scheme to ensure that the magnetic field obeys the no-monopoles constraint (see Zilhão & Noble 2014 for more detail).

The initial conditions are as follows: for SANE disks,

where is the maximum initial plasma density; and for MAD disks,

where r0 is chosen to be the inner boundary of the simulation domain (see Hawley 2000; Gammie et al. 2003; Wong et al. 2021a). The MAD prescription concentrates the initial magnetic field toward the inner edge of the disk and forces it to taper at large r. In both the SANE and MAD cases, the magnetic field is normalized so that the ratio between the maximum Pgas and the maximum Pmag = b2/2 is some target value, typically 100. See Appendix A for more details regarding our fluid initial condition.

Note that other initial conditions are possible, including alternative fluid tori solutions (Fishbone 1977; Abramowicz et al. 1978; Kozlowski et al. 1978; Chakrabarti 1985; Penna et al. 2013) as well as implementations of more realistic "boundary conditions" determined by, e.g., simulations of magnetized stellar winds near the putative supermassive black hole in the Galactic Center (Ressler et al. 2018, 2020a, 2020b).

3.2. Ray-tracing Images

ipole (Mościbrodzka & Gammie 2018) is a publicly available GRRT code for covariant polarized radiative transfer and is a descendant of the unpolarized image code ibothros (Noble et al. 2007). ipole produces a rectangular polarimetric image defined by a field of view (width in GM/c2, or translated to μ as in the context of EHT sources), distance from the black hole dsrc, and orientation with respect to the black hole spin axis and midplane (inclination and position angle). Each pixel reports the specific intensities for the Stokes parameters Iν , Qν , Uν , Vν at the pixel center as well as the total optical and Faraday depths along the pixel-centered geodesic.

ipole is an observer-to-emitter or backward-in-time code, meaning that only geodesics that intersect a predefined pinhole camera (observer) are considered. The camera is defined by a particular coordinate and an orthonormal tetrad specified by the normal observer velocity , 30 a radially directed wavevector, and the black hole spin axis. The pinhole camera defines an image with pixels that intersect and are regularly spaced in angle as measured in the tetrad frame. This prescription ensures that the central pixel corresponds to a geodesic with the impact parameter zero. The geodesic for each pixel is integrated backward toward the black hole (i.e., the emitting matter) until it leaves the simulation's user-specified radiation domain, either at large radius or at the event horizon. The radiative transfer equation is then solved forward toward the camera.

ipole uses operator splitting to solve Equation (11) in two independent stages, which separately account for the relativistic parallel transport and nonrelativistic emission, absorption, and mixing radiative transfer updates. In the first substep, ipole numerically integrates and parallel transports the coherency tensor Nα β along the geodesic. In the second substep, it projects Nα β into a Stokes vector in a local fluid tetrad, 31 evaluates transfer coefficients in that tetrad, and uses a nonrelativistic analytic solution for the case of constant coefficients (Landi Degl'Innocenti & Landi Degl'Innocenti 1985) to perform transport in the tetrad frame. The use of an analytic solution ensures numerical stability when the plasma has large optical or Faraday depths. The code is coordinate agnostic, since all geometric factors are computed numerically from the metric. A comparison of contemporary GRRT codes is available in Gold et al. (2020).

The Illinois version of ipole, which is used in the PATOKA pipeline, implements several new features and differs from the originally published version in several ways. First, Illinois ipole implements an optional adaptive ray-tracing scheme that allows resolution to be concentrated in the regions of the image with sharp, detailed features (for details see Gelles et al. 2021b). The other differences are described below. Both the originally published version of ipole and the Illinois version are publicly available. 32 ipole converges at second order.

3.2.1. Tracking the Emission Source

Since we aim to connect observables to the accretion flow structure and the spacetime geometry, it is useful to be able to study where the emission is produced. In order to compute what fraction of the emission produced by a given volume will contribute to an image, it is necessary to specify the location of the observer, since the emission coefficients depend on the photon wavevectors kμ through the local sampling frequency ν = − kμ uμ and the pitch angle = arccos. Furthermore, strong lensing allows multiple geodesics to sample the same region of space; although the photons may have the same frequency at infinity, in general each will be at a different frequency in the rest frame of the plasma since the geodesics are not parallel. Therefore a one-to-one mapping between points in space and emissivities does not exist.

Moreover, not all emitted light escapes to infinity. Light emitted near the event horizon may fall into the hole, and the optical depth of the plasma between emission and the observer will cause some light to be reabsorbed. It is straightforward to compute what fraction of light makes it to the observer at infinity for the Stokes I total intensity by saving the values computed when solving the approximate radiative transfer equation, Equation (14), as follows.

At each step along the geodesic from s1s2, we save the local optical depth and the local contribution to the invariant intensity, . After the full geodesic has been traced, the observed fraction of the photons emitted at any point along the ray that make it to the camera can be computed as , where τ is the optical depth to the camera from that point. For each step along the geodesic, the local contribution to the final intensity is computed from and the observed fraction; these values are then saved in an array representing the simulation domain. The process is validated after the complete image has been generated by summing the flux density contributed by each simulation zone and comparing to the final total image flux density. Since the output of this procedure is the total contribution from each simulation zone, it is important to account for the coordinate transformation Jacobian when plotting, e.g., the apparent density of the observed emitted photons.

3.2.2. Subring Decomposition

Simulated images of EHT sources commonly display a clear ring-like structure. This structure is identified with the photon ring produced as strong lensing allows light paths to wind around the black hole multiple times (see, e.g., Johnson et al. 2020). In aggregate, the composition of the subrings produces a logarithmically divergent brightness temperature that is limited in part by the optical depth of the plasma.

Observable signatures of the photon ring structure have been studied analytically (e.g., Gralla & Lupsasca 2020; Himwich et al. 2020; Chael et al. 2021; Vincent et al. 2021), but the detailed structure of the observable photon ring is heavily influenced by the structure of the emitting plasma. We have included a subring decomposition routine in ipole that allows the code to produce separate images for each subring. The subring structure is particularly evident when the emission exhibits a nontrivial latitudinal structure, e.g., if the emission is concentrated near the midplane and each subring contribution can be easily separated from its neighbors.

Geodesics that pass close to the hole experience strong lensing and may undergo latitudinal oscillations—this effect is most obvious for the unstable, constant-radius bound orbits, which circle the hole indefinitely and undergo periodic oscillations as they sweep through a range of latitudes (see Teo 2003). Each event xμ along a nonequatorial geodesic can be assigned an orbit number n, corresponding to the number of latitudinal turning points between the event and the camera in the region near where the geodesic passes through the midplane. 33 This value can be directly tracked and saved during the backward geodesic integration. Note that the subrings are nested and overlap: the nth subring lies within the , with different segments of the same geodesic that undergoes m latitudinal oscillations contributing to the first m + 1 subrings. The n = 0 component covers the full image, and the higher-order subrings have exponentially decreasing areas as they tighten around the critical curve.

In order to synthesize the subring image, the emission coefficients are zeroed during the forward radiative transfer integration in the regions of the geodesic with nntarget. It is important to include parallel transport as well as absorption and rotation in regions with nntarget, since although the subring image comprises only photons that were emitted along the corresponding geodesic segment, we must account for how those photons interact with matter and spacetime as they propagate to the observer. Thus, although jI,Q,U,V are set to zero, αI,Q,U,V and ρQ,U,V are not.

3.2.3. Differences from the Original Version

We have slightly modified ipole from the originally published version in several ways:

  • 1.  
    ipole has been modified to read fluid snapshot files with arbitrary logical coordinate systems through the use of a new module called simcoords. Using simcoords, ipole performs all ray-tracing in exponential Kerr–Schild (eKS; see Appendix F) coordinates and uses an interpolated grid map between eKS coordinates and the input snapshot coordinates. In order to use the simcoords module, a fluid snapshot must report the location of each grid zone in Kerr–Schild (KS) coordinates, provide the velocity and magnetic-field vectors with KS components, and correspond to a single slice of constant KS time.
  • 2.  
    Several different methods for treating the σ magnetization cutoff have been implemented. The default version treats σ as an interpolated scalar and simply zeros the transfer coefficients in regions with σ > σcutoff by setting the number density of electrons ne = 0. Alternative methods include gradually suppressing ne according to a sigmoidal function of σ or zeroing ne on a per-zone basis according to zone-centered values of σ. The choice of cutoff procedure may introduce image artifacts—see Appendix C.
  • 3.  
    When evaluating fluid parameters along geodesics, we interpolate the scalars ne , u/ρ, magnetic-field strength as well as the six primitive variables that describe the fluid velocity and magnetic-field orientation. This procedure ensures that temperatures, which are derived from u/ρ, remain reasonable and that the interpolation scheme does not lower the magnetic-field strength in regions where the magnetic field oscillates wildly, like near the jet–disk boundary. Then, by constructing uμ and bμ from the primitives, we can ensure uμ uμ = −1 and uμ bμ = 0. We have observed that interpolating the four-vector components directly tends to produce deviations from these criteria within ≈2 GM/c2 of the event horizon. Deviations have been observed to be catastrophic in some cases when the simcoords module is used.
  • 4.  
    During evaluation of the polarized transfer coefficients, ipole enforces positivity constraints on transfer coefficients to ensure that jν,I ≥ 0 and αν,I ≥ 0 and that and . Failure to enforce the latter condition leads to division-by-zero errors in evaluation of the analytic, constant-transfer-coefficients solution that ipole uses to advance the Stokes vector. Failure to enforce the condition on jν,I can lead to polarization fractions > 100%, and failure to enforce positivity on jν,I leads to unphysical negative intensities. Furthermore, ipole reorthogonalizes the fluid-frame tetrads to ensure orthonormality to machine precision. Finally, consistent limiting expressions in areas of low absorption or rotation improve the accuracy of the maps of the polarization fraction.
  • 5.  
    ipole includes implementations of transfer coefficients for various nonthermal eDFs, such as the power law and kappa distributions. The coefficients are computed using symphony (Pandya et al. 2016, 2018; Marszewski et al. 2021), which is now a submodule of the code.
  • 6.  
    Several analytic and test models have been added to check ipole against existing radiative transport schemes (see, e.g., Gold et al. 2020; and Appendix E).
  • 7.  
    The coefficients of the radiative transfer equation, along with the fluid parameters that produced them, can be recorded along the geodesic corresponding to any pixel in the image and output as a trace file. This enables easier post-hoc study of regions that may produce characteristic emission and also allows for easier error diagnosis.

3.3. Producing Spectra

igrmonty (Dolence et al. 2009) is a Monte Carlo radiative transport code designed to generate spectra from GRMHD fluid simulation snapshot files of optically thin ionized plasmas. It accounts for the full angle- and frequency-dependence of emission and absorption, and it treats single Compton scattering exactly. igrmonty is an emitter-to-observer code, also known as forward-in-time, meaning it simulates emission at all frequencies and angles across the entire domain. This procedure is slower than the observer-to-emitter procedure, but it produces full spectra ν Lν of the source as seen from all inclinations and longitudes around the black hole.

igrmonty tracks a Monte Carlo sample of the radiation field in the form of superphotons. A superphoton with weight w is a packet of w ≫ 1 photons, where each photon has the same position xμ and wavevector kμ . Superphotons are created across the computational domain according to the local emissivity of the plasma. In order to emit a target number Ntarget of superphotons over the full simulation domain, the bolometric luminosity due to emission is precomputed for each zone and used to determine what fraction of the Ntarget total superphotons should be emitted per zone. This heuristic can produce a noisy signal when the fluid simulation resolution is increased if there is nontrivial structure in the accretion flow, thus higher fluid resolutions typically require a larger number of superphotons during the spectrum-generation step.

To optimize signal in the spectrum, the weight factor w is conventionally chosen to be frequency-dependent according to the heuristic that each logarithmic bin in energy space in the final output spectrum should contain approximately the same number of superphotons. Since it is impossible to predetermine how many superphotons survive extinction on their journeys to the observer, we estimate how the precomputed weight factors should be set by assuming that all emitted photons escape to infinity—note that we relax this assumption when simultaneously treating both bremsstrahlung and synchrotron, as the two components tend to be peak at significantly different frequencies and thus undergo different amounts of extinction. We also neglect factors like redshift, scattering, and the ultimate angular dependence of the spectrum. When a new superphoton is created, its frequency is chosen according to rejection sampling, and the rest of its wavevector is initialized in a local orthonormal tetrad according to the pitch-angle-dependent emissivity prescription. Each superphoton saves its initial position; this information can be used to infer the properties of the emission region.

After emission, the optical depths to both absorption and scattering are computed as the superphotons propagate along their geodesics. Absorption is accounted for by decreasing the weight w of the superphoton packet, which is directly proportional to the intensity of the radiation packet. igrmonty treats scattering in a probabilistic sense: a superphoton will scatter with some probability at each step along its geodesic; if it scatters, the wavevector of the scattered superphoton is evaluated in a local orthonormal tetrad and determined probabilistically from the differential electron scattering cross section.

In order to boost the signal in the Compton-upscattered component of the spectrum, igrmonty uses a variant of the bias method introduced in Kahn (1950). Although the likelihood p that a full superphoton will scatter may be small in the optically thin plasmas we consider, the bias procedure enables a fraction 1/b of a fiducial superphoton to scatter with probability bp. Here, b is the value of the bias parameter, which can in general vary across the different simulation zones and time steps. If the superphoton scatters, its weight is decreased to 1 − 1/b, and a new superphoton representing the scattered component is created with weight 1/b.

igrmonty is publicly available, 34 converges at second order in geodesic integration, and converges like in the number of Monte Carlo radiation field samples. The PATOKA version of igrmonty only solves the approximate (unpolarized) radiative transfer equation, but recently Mościbrodzka (2020) introduced the RADPOL scheme, which accounts for fully polarized synchrotron emission, absorption, Faraday rotation and conversion, and Compton scattering. Appendix E presents a comparison between igrmonty and ipole in the lower-frequency part of the spectrum where Compton upscattering is unimportant.

3.3.1. Scattering Bias Factor

It is difficult to determine how the bias should be set before running the simulation. If it is too low, too few superphotons will scatter, and the Compton contribution to the spectrum will be unusably noisy. If it is too high, the code could reach a supercritical state, in which superphotons produced through scattering also undergo scattering events, and the total number of superphotons diverges. In all cases, the bias factor must be ≥1.

The bias parameter is a function of the local fluid state and typically scaled by the square of the local electron temperature to improve resolution in the scattering events with higher energies compared to the average. This location-dependent value is multiplied by a global tuning factor btuning, which scales the scattering rate across the entire simulation. igrmonty begins each run with a low-resolution "bias tuning" step, during which btuning is varied until the ratio between the number of superphotons created through scattering, and the number created through emission is approximately unity. The ratio is tracked during the tuning runs, and the evaluation is halted if the ratio exceeds a large number to ensure that the supercritical state can be preempted.

3.3.2. Bremsstrahlung

Yarza et al. (2020; see also Narayan & Yi 1995; Esin et al. 1996; Mahadevan 1997) found that, in radiative simulations of SANE accretion flows, bremsstrahlung emission may dominate the high-energy ( ≳512 keV) component of the spectrum and contribute to the bolometric luminosity as the accretion rate is increased. In relativistic plasmas, both electron–electron and electron–ion bremsstrahlung contribute to the total emissivity; igrmonty supports several contemporary prescriptions for both emissivities (see the Appendix of Yarza et al. 2020 for a comparison of the prescriptions). Since synchrotron and bremsstrahlung emission may be simultaneously nonzero, each emitted superphoton is assigned a bremsstrahlung-fraction value between zero and one. This value is saved when the spectrum is recorded and can be used to determine how much emission is associated with which emission process.

Bremsstrahlung emission dominates the direct synchrotron at higher frequencies where the optical depth of the plasma is low. A naive application of the superphoton weight assignment scheme described above may thus preferentially improve the signal in the bremsstrahlung component of the spectrum. This shortcoming can be addressed by either of the following: modifying the weighting procedure to independently generate weights for the lower-frequency synchrotron emission and the higher-frequency bremsstrahlung emission; or independently generating spectra from the synchrotron- and bremsstrahlung-emission components and adding them together. Compton scattering is typically unimportant for bremsstrahlung emission, so the bremsstrahlung calculation may be performed with scattering turned off.

3.3.3. Arbitrary Electron Distribution Functions

igrmonty relies on the full eDF both when computing the Compton scattering cross section and during the scattering procedure when choosing an electron off which to scatter. The version of the code published in Dolence et al. (2009) implemented a semianalytic treatment of the cross section and scattering sampling routines; here, we instead use an analytic expression only for the distribution function itself and numerically evaluate the cross section. We use the rejection sampling technique to sample the distribution rather than rely on an inversion of the distribution function. Because this procedure relies only on a prescription for the distribution function, it is quite general and supports any isotropic distribution that can be written as a function of the local fluid parameters.

4. Sample Simulation Data Products

We now present example data products from the PATOKA pipeline. The data we show here were generated for the EHT M87 simulation library that was used for validation in EHTC IV, EHTC VII, and for analysis in EHTC V; EHTC VI; EHTC VIII. We also present examples of unpublished data that were produced for EHTC IV; EHTC V; EHTC VI; EHTC VII; EHTC VIII.

4.1. GRMHD Models

The GRMHD library generated by PATOKA for the EHT M87 analysis included two parts. The first part comprised ten "canonical" simulations spanning the MAD and SANE accretion states and five spins a* = −15/16, −1/2, 0, +1/2, +15/16. 35 Hereafter, we write 1/2 as 0.5 and 15/16 as 0.94 to be consistent with the EHT paper sequence. The second part was composed of several dozen ancillary simulations used to test different initial conditions, disk sizes, adiabatic indices, and simulation resolutions. Table 1 lists all simulations and identifies which simulations belong to the canonical set. Each of the models generated for the library was evolved until at least t = 104 GM/c3, during which time its accretion flow reaches a statistical steady state within r ≤ 10–20 GM/c2. GRMHD fluid snapshots were saved every 5 GM/c3, corresponding to ≈43 hours for a black hole with mass ≈6.2 × 109 M.

Table 1. List of GRMHD Simulations

Flux a* rin rout Resolution(s)Notes
SANE0.94612504/3192×192×192, 288288 canonical
SANE0.75612504/3192, 288...
SANE0.5612504/3192×192×192, 288288 canonical
SANE0614504/3288canonical
SANE−0.258175013/9288...
SANE−0.5817504/3288canonical
SANE−0.51020505/3192×128×128, 288×192×192...
SANE−0.7510205013/9288...
SANE−0.941020504/3288canonical
SANE−0.941020505/3192×128×128, 288×192×192...
MAD0.94204110005/3288, 384...
MAD0.942041100013/9192, 288, 384, 448384 canonical
MAD0.752041100013/9384...
MAD0.52041100013/9192, 288, 384, 448384 canonical
MAD0.252041100013/9192, 384...
MAD02041100013/9192, 288, 384, 448384 canonical
MAD−0.252041100013/9192, 384...
MAD−0.52041100013/9192, 288, 384, 448384 canonical
MAD−0.942041100013/9192, 288, 384, 448384 canonical

Note. GRMHD fluid simulation parameters. Canonical simulations are identified in the notes column. Flux labels the relative strength of the magnetic flux at the horizon, a* describes the spin of the black hole, rin and are parameters for the initial Fishbone–Moncrief torus, rout is the outer edge of the simulation domain, is the adiabatic index of the fluid, and resolution gives the Nr × Nθ × Nϕ number of grid zones in the simulation. In the case of resolution, single-number abbreviations mean the following: 192 → 192 × 96 × 96, 288 → 288 × 128 × 128, 384 → 384 × 192 × 192, and 448 → 448 × 224 × 224.

Download table as:  ASCIITypeset image

The following set of limiting conditions was imposed on the fluid evolution to ensure stability:

  • 1.  
    density ρ > 10−6 k for ,
  • 2.  
    internal energy ,
  • 3.  
    ρ and u were increased until σ < 400 and , where is the magnetic pressure,
  • 4.  
    ρ was increased until ,
  • 5.  
    when evolving electron temperatures, u was decremented until , and
  • 6.  
    the velocity components were downscaled until the Lorentz factor with respect to the normal observer was Γ = −uμ nμ < 50.

Global disk simulations inevitably invoke some of these bounds, most frequently the ones on σ and Γ. The former is visible in the top-right panel of Figure 3 as a halo of accretion near the event horizon.

Figure 3. Refer to the following caption and surrounding text.

Figure 3. GRMHD simulation of a MAD accretion flow with a* = 0.5. This simulation was initialized with an rin = 20, torus and run at a resolution of 448 × 224 × 224 zones in the r, θ, and ϕ directions, respectively. Top panels show rest-frame plasma density in arbitrary units at three different times over the course of the simulation; the line integral convolution technique is used to represent the motion of the plasma. The bottom panels show magnetization at the same times as the top panels with overplotted contours of the axisymmetrized vector potential component Aϕ , which give a sense of field strength and disorder. The black hole event horizon is plotted as a dark circle. As the simulation evolves, the flow becomes increasingly turbulent and a high-magnetization jet region opens around the poles.

Standard image High-resolution image

We only produced simulations with purely corotating (aligned) or counterrotating (antialigned) accretion flows, since varying the disk tilt adds another dimension to the parameter space and is thus prohibitively expensive. We only produced SANE simulations with ϕ ≈ 1, although initial conditions that produce ϕ ≪ 1 and 1 ≲ ϕϕc accretion states are also known. The fluid was assumed to have a uniform constant adiabatic index for each simulation , although the value of was varied between different simulations.

Each simulation was run on a 3D regular grid defined in the horizon-penetrating funky MKS (FMKS; see Appendix F) coordinates, which concentrate resolution toward the midplane and away from the jet at small radius. The simulation domain extended from ≥5 simulation zones within the event horizon to an outer radius rout ≥ 50GM/c2, depending on the simulation. We found that evolving a large disk over a long time could cause an initially SANE accretion flow to go MAD, so we chose to initialize the canonical SANEs with smaller accretion disks, allowing us to use smaller simulation domains and lower absolute resolutions.

Figure 3 shows snapshots of the plasma rest-mass density ρ and magnetization σ over the course of evolution of an intermediate-spin a* = 0.5 MAD simulation. The initial condition with a large Fishbone–Moncrief torus and an ordered, looping magnetic field, progressively gives way to a turbulent accretion flow with a low-density high-magnetization funnel containing a strong, ordered magnetic field. The qualitative difference between the MAD and SANE accretion morphologies is shown in Figure 4. The MAD simulation is punctuated by magnetic bubbles corresponding to flux ejection events.

Figure 4. Refer to the following caption and surrounding text.

Figure 4. Comparison of density, plasma β, and magnetization σ in the midplane of MAD and SANE simulations in steady state. Both simulations were run with an intermediate spin a* = 0.5. In the MAD simulation, the accretion proceeds in thin strands, contrasting the steady but turbulent disk-like SANE accretion mode. The evacuated regions in the MAD simulation with low β and high σ are magnetic bubbles produced during the flux ejection events that episodically recur when an excess of magnetic flux has been trapped on the event horizon.

Standard image High-resolution image

4.2. Electromagnetic Observables

The GRMHD simulations were postprocessed to generate images and spectra using ipole and igrmonty. In total, about 120,000 images were generated for the canonical M87 total intensity and polarization analyses, and about 3 million were generated for supporting analyses, including resolution and field-of-view studies, explorations of the analyses robustness to changes in the numerical parameters like the geodesic step size, and machine-learning projects. The full set of the radiative transport model parameters is described below. The canonical images for the total intensity analysis were generated at 160 × 160 pixels resolution, and the images generated for the polarization analysis were at 320 × 320 pixels. All images were produced to have a 160 μ as field of view.

For the images produced for the EHT M87 sequence, we assume that the accreting plasma is composed of pure ionized hydrogen, so that y = z = 1 in Equation (27). We fixed βcrit = 1 and allowed the two parameters in Equation (22) to vary between rlow = 1, 10 and rhigh = 1, 10, 20, 40, 80, 160. See EHTC V and EHTC VIII for a discussion of the motivation behind these choices.

The library discussed here was generated to compare against observations of M87, so it was generated using the physical parameters that would target that system. The inclination angle i was chosen to be consistent with the orientation of the M87 jet at large scales, i ≈ 17° (Hada et al. 2017; Kim et al. 2018; Walker et al. 2018), so we produced a library with inclinations ranging over multiple inclinations i = 7°, 12°, 17°, 22°, 27° relative to the line of sight. We do not know a priori whether the black hole spin axis is directed toward or away from us. An exploratory survey of the library showed that it was necessary to orient the black hole spin vector away from Earth in order to reproduce both the image brightness asymmetry and the position angle of the large scale jet (EHTC V). The position angle of the spin axis can be reoriented during analysis by rotating both the image and the per-pixel EVPA.

The GRMHD simulation mass density was scaled to physical units by requiring that the simulated compact flux density at 230 GHz be consistent with the observed contemporaneous flux density of between 0.5–0.7 Jy—see EHTC IV for more detail on identifying this target flux density. The relationship between the scaling factor and total flux density Ftot is a complicated function of the accretion flow, but it tends to be monotonic near the target value, so identifying the appropriate scaling factor corresponds to a simple root-finding procedure. Since the fitting procedure is expensive, the flux density is typically fit using the approximate total intensity solution over a regular subsample of the snapshots at low resolution. The quality of the fit is substantiated when the high-resolution data are generated. After identifying the value of required to produce the target flux density, every snapshot from each GRMHD model is typically imaged, producing a sequence with a 5 GM/c3 cadence.

The result of running the flux-density-fitting procedure is shown in Figure 5 for the canonical MAD a* = 0.5 GRMHD simulation with rlow = 1 and i = 17°. The unscaled mass accretion rate of the system and the dimensionless magnetic flux parameter ϕ are also plotted for comparison. Here, the light curves have been fit so that the average 230 GHz flux density matches 1 Jy over the last 5000 GM/c3 of the simulation.

Figure 5. Refer to the following caption and surrounding text.

Figure 5. Time series of flux variables from the MAD a* = 0.5 simulation. Top panel: mass accretion rate in arbitrary code units vs. time. Center panel: dimensionless magnetic flux ϕ vs. time. Bottom panel: light curves at ν = 230 GHz for rlow = 1 and different rhigh values. The light curves have been scaled to match a 1 Jy target flux density over the last 5000 GM/c3 of the simulation. The first half of the simulation is dominated by a transient from the torus initial condition. Notice that stability in and ϕ does not necessarily equate to stability in the 1.3 mm light curve.

Standard image High-resolution image

Figure 6 shows an example M87-like synthetic image from the MAD a* = 0.5 model with rlow = 1 and rhigh = 40 at i = 17°. Each pixel contains the full polarimetric Stokes I, Q, U, V specific intensities, which can be processed to provide information about the linear polarization fraction , EVPA χ, and circular polarization fraction CP = V/I. Blurring is performed by convolving the Stokes intensities with a 15 μas Gaussian beam. In the right three columns, the pixel brightness is determined by the total intensity so that the low-intensity regions, which contribute little flux in an observation, appear black.

Figure 6. Refer to the following caption and surrounding text.

Figure 6. Example polarimetric images of MAD a* = 0.5 simulation produced by ipole with (bottom) and without (top) blurring to 15 μ as Gaussian beam. Panels show total intensity, linear polarization fraction, electric-vector position angle (EVPA), and circular polarization fraction. The brightness of each pixel has been scaled to the total intensity; the patterns in linear polarization, EVPA, and circular polarization continue into the dark regions. Blurring decreases the observed linear polarization fraction in regions where EVPA is rapidly varying.

Standard image High-resolution image

Figure 7 shows the result of the ring decomposition procedure used to isolate the different subrings in a ray-traced image (see, e.g., Figure 3 of Johnson et al. 2020 for cross sections of a similar decomposition). As n increases, the extent of each next subring, which corresponds to a geometric region on the image, is exponentially demagnified compared to the last one according to Lyapunov exponents at each angle on the image (see Johnson et al. 2020; Wong 2021). In physical systems, the increasing optical depth of the source limits how many subrings contribute to the composite image.

Figure 7. Refer to the following caption and surrounding text.

Figure 7. Decomposition of the total intensity into the first three subrings for a snapshot from the high spin a* = 0.94 MAD model with rlow = 1, rhigh = 10, and i = 17°. The intensities in the n = 0, n = 1, and n = 2 subrings are reported in the right three panels, and the location of the critical curve is plotted as a dashed line in green. The total flux density, whose decrease is governed by angle-dependent Lyapunov exponents in the limit of large n for optically thin systems, is given in the top right of each panel. The left panel shows a composite image produced by adding the n = 0 subring in red, n = 1 subring in blue, and n = 2 subring in green. The visual appearance of each subring image is due to the details of the underlying emission, although the geometric extent of each image is due to the spacetime.

Standard image High-resolution image

ipole can be used to track the source of the observed flux density, as seen in Figure 8. All MAD simulations tend to show the same characteristic fragmentary emission structure, which corresponds to the disjoint accretion (see, e.g., the top left panel of Figure 4). Much of the emission in the MAD case is thus produced in the hot, chaotic region of the flow near the horizon. SANE emission is comparatively more structured. Changing the electron temperature model in the SANE simulations can have drastic effects by heating up the jet funnel and shifting emission out of the disk.

Figure 8. Refer to the following caption and surrounding text.

Figure 8. Location of emission around high-spin black holes with a* = 0.94 on a single KS time slice. In this figure, the emissivity is embedded in and represented with respect to a Cartesian coordinate system from simulation coordinates, and all densities are with respect to the Cartesian volume element. Top: 3D rendering of emission source with color and transparency determined by the total emission produced within that region of space. Bottom: the same data as above after summing across the azimuthal dimension ϕ. The left two columns are a representation of emission from a single snapshot from a typical MAD simulation; emission tracks the fragmentary plasma and is relatively insensitive to the electron thermodynamics, parameterized here by rhigh. The right two columns represent emission from the same SANE simulation snapshot but with the different electron temperature prescriptions. Larger values of rhigh shift the emission from the turbulent but steady disk into the funnel region. All simulations have rlow = 1 and are imaged at i = 17°. The total flux density produced by each simulation is the same, so the color scales show the relative concentration of emission in the azimuthal sum.

Standard image High-resolution image

Running igrmonty is significantly more computationally expensive than running ipole, so it is infeasible to generate spectra for every fluid snapshot across every radiation model. The two spectrum constraints considered in the EHT analysis were the overall radiative efficiency of the flow and a boolean determination if the simulation X-ray flux was consistently too high. At minimum, we generate 10 spectra per radiation model, but we checked that producing a denser sampling of spectra in time does not change the result in a statistical sense. Figure 9 shows example spectra produced from one MAD and one SANE snapshot at rlow = 1 and across different values of rhigh.

Figure 9. Refer to the following caption and surrounding text.

Figure 9. Example spectra for MAD (left) and SANE (right) snapshots at low inclination in an rlow = 1 model. The spectra for rhigh = 1, 10, 40, 160 are shown in the top row, and the components of the rhigh = 10 model are shown in the bottom row. All data points are taken from the simultaneous multiwavelength measurement campaign performed coincident with the 2017 EHT observations of M87 and reported in The EHT MWL Science Working Group et al. (2021).

Standard image High-resolution image

5. Future Directions

We now briefly discuss future directions as well as improvements and modifications that can be made to the PATOKA pipeline.

5.1. Radiative Transfer Model

ipole produces images calculated at a single frequency, which neglects the observing bandwidth of the instrument. Extensions to the ray-tracing code could allow for the synthesis of finite bandwidth observations; however, this approximation has been found to be inconsequential in the context of the M87 library.

Our treatment generally assumes that the eDF is thermal, i.e., that it is well described by a Maxwell–Jüttner distribution. This assumption enters through definition of the transfer coefficients, which are calculated from the underlying distribution function (see, e.g., Shcherbakov 2008; Pandya et al. 2016, 2018). The introduction of nonthermal electrons can change both the spectral shape and the image morphology (e.g., Özel et al. 2000; Yuan et al. 2003; Broderick & Loeb 2009; Chael & Narayan 2017; Mao et al. 2017; Davelaar et al. 2018).

GRMHD simulations produce snapshots of the fluid at different, discrete KS times. The ipole and igrmonty data we present were generated under the fast light approximation, where only a single snapshot is used to generate an image or spectrum. This approximation is invalid if the fluid evolves on timescales shorter than the simulation light-crossing time or if one wishes to simulate various time-dependent phenomena, like light echoes from flares or glimmer (e.g., Broderick & Loeb 2005; Moriyama & Mineshige 2015; Hadar et al. 2021; Wong 2021).

The alternative slow light method relies on a high fluid snapshot cadence, has large data-storage requirements, and has a high throughput cost, which is associated with ray tracing through different time slices. Although not presented here, several slow light simulations were generated to confirm that the fast light approximation does not seriously affect the library results (see, e.g., Dexter et al. 2010; Mościbrodzka et al. 2021).

5.2. Radiative Effects

The pipeline we have described uses ideal GRMHD to generate the fluid simulations; we have assumed that M87 can be described by models in which radiative cooling is negligible so that it does not affect the dynamics of the plasma. This assumption was probed in EHTC V, but it is likely that a full radiative treatment of the fluid simulation will be required in future analyses.

The M87 jet funnel may be populated by electron–positron pairs, produced either via cascades (e.g., Beskin et al. 1992; Levinson & Rieger 2011; Broderick & Tchekhovskoy 2015; Hirotani & Pu 2016) or via drizzle (Mościbrodzka et al. 2011; Laurent & Titarchuk 2018; Kimura & Toma 2020; Wong et al. 2021b; Yao et al. 2021). Ideal GRMHD cannot produce unscreened electric fields and therefore cannot track pair cascades. Computing the cross sections for pair drizzle often requires a high-resolution sample of the radiation field, so it is expensive to track in situ in GRMHD simulations and is often evaluated in a postprocessing step. Future study is warranted to investigate the signatures of pair plasma emission and whether or not pairs can populate the jet (see Broderick & Tchekhovskoy 2015; Anantua et al. 2020a; Emami et al. 2021; Yao et al. 2021).

5.3. The Accretion Model

The library presented in Section 4 is not comprehensive. A more dense sampling of spin may enable a better understanding of how spin affects observables, particularly as spin approaches its maximal value . The transitory regime in which the magnetic flux increases from the comfortably SANE state toward the MAD state has been explored, but no dense parameter survey yet exists. The tilted-disk scenario merits further attention and study, even though it increases the dimensionality of the final parameter space by two.

A detailed study of the convergence properties of GRMHD simulations, with respect to both spatial resolution and simulation duration, would be valuable, as would a systematic survey of how the initial conditions affect the statistical properties of the fluid evolution and electromagnetic observables.

5.4. Viscosity

GRMHD simulations typically assume a single ideal fluid and do not model the effects of fluid viscosity or resistivity. While the magnetic Reynolds number in accretion flows is likely so high as to render resistivity irrelevant to the bulk dynamics of the flow, the very long mean free paths of accreting electrons and ions suggest it may be important to include the effect of viscosity, and thus heat conduction and pressure anisotropy. Note that explicit consideration of resistivity may still be important in studies of plasmoids and reconnection, which can contribute to the generation of flares and nonthermal particle injection (e.g., Ripperda et al. 2022).

Simulations treating pressure anisotropies and heat conduction find that, although the pressure anisotropies with respect to the magnetic field can be large in such systems, the overall structure of the accretion state remains similar to ideal simulations (Foucart et al. 2017). Thus, 230 GHz images of so-called extended GRMHD simulations are not expected to be drastically different from images of ideal GRMHD simulations. Nevertheless, including the effects of viscosity may alter the electron thermodynamics and particular dynamics of the flow. Extended GRMHD simulations would have been prohibitively expensive to conduct at the range of parameters required for the current study and are the subject of current development (see Chandra et al. 2015, 2017; Most & Noronha 2021; Most et al. 2021).

5.5. Electron Acceleration

Single-fluid ideal GRMHD simulations assume that the plasma is well described by the bulk variables that represent both the constituent electrons and ions; however, the supermassive black hole accretion flows are likely collisionless, the electrons and ions are unequilibrated, and the dynamical imbalance between electron acceleration and cooling may produce a significantly nonthermal eDF.

The influence of the nonthermal electrons on the observables can be roughly accounted for given a model for the eDF as a function of the local fluid parameters, since the radiative transfer coefficients associated with a given eDF can be computed and applied during postprocessing. Rather than use a local post-hoc model for assigning the eDF, it is also possible to use a model for energy injection rates, e.g., motivated from PIC simulations (e.g., Sironi & Spitkovsky 2014), and track the nonthermal electron component directly in the fluid simulation by evolving the fraction of electrons associated with different energy bins in each simulation zone (e.g., Ball et al. 2016; Chael & Narayan 2017; Petersen & Gammie 2020).

6. Summary

We have described PATOKA, a numerical simulation pipeline that includes the iharm, ipole, and igrmonty codes. PATOKA can be used to generate the spectra and fully polarimetric images of the radiatively inefficient accretion flows around supermassive black holes. We have provided a brief sample of the PATOKA data products that were used in the EHT analysis of the total intensity and linear polarimetric data from the M87 black hole. Simulations and data produced from some or all of the PATOKA pipeline have also been used in other works, including Porth et al. (2019), Johnson et al. (2020), Lin et al. (2020), Palumbo et al. (2020), Gold et al. (2020), Wielgus et al. (2020b), Ricarte et al. (2020), Gelles et al. (2021a), Wong et al. (2021a), Ricarte et al. (2021, 2022), Tiede et al. (2022).

The authors thank Zachary Gelles, Michael Johnson, and Jordy Davelaar for many discussions surrounding the radiative transfer codes, and Christian Fromm and the anonymous referee for their careful reading of and helpful suggestions regarding the manuscript. This work was supported by National Science Foundation grants AST 17-16327, OISE 17-43747, AST 20-07936, AST 20-34306, by a Donald C. and F. Shirley Jones Fellowship to G.N.W., and by the US Department of Energy through Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the US Department of Energy (Contract No. 89233218CNA000001). M.M. acknowledges support by the NWO grant No. OCENW.KLEIN.113. A.R. also acknowledges support from the Black Hole Initiative (BHI), enabled by grants from the Gordon and Betty Moore Foundation and the John Templeton Foundation. G.N.W. also gratefully acknowledges support from the Institute for Advanced Study. This work has been assigned document release number LA-UR-21-28826.

Appendix A: Fishbone–Moncrief Torus Description

Fishbone & Moncrief (1976; hereafter FM) described a two-parameter solution of the relativistic Euler equations for an isentropic, stationary, axisymmetric, and purely azimuthal flow (ur = uθ = 0) in a stationary, axisymmetric background. In Kerr, the Euler equations for the above scenario are

where we write in Boyer–Lindquist (BL) coordinates, h is the specific enthalpy of the fluid, s is its specific entropy, T is temperature,

and the projections of uμ onto the orthonormal vectors of the locally nonrotating frame are (Bardeen et al. 1972)

In our work, we take the solutions with constant l = uϕ ut . For these solutions, the fluid enthalpy can be expressed analytically as a function of l, r, a*, θ, and a parameter defining the location of the inner edge of the disk in the midplane rin,

where the second line is the negation of the first, as evaluated in the midplane at the inner edge of the torus; thus Σin, Δin, and Ain are the expressions defined in Equations (A2)–(A4) with r = rin.

In iharm, we set the angular momentum density of the FM torus in terms of its value at , the radius of maximum pressure in the midplane

The region outside the torus is initialized to zero enthalpy. The remaining fluid variables are computed from the enthalpy using the ideal gas equation of state, the definition of specific enthalpy , and the isentropic condition, which enables us to equate the polytropic index and adiabatic index of the fluid. All grid zones that fall outside the torus are initialized to the density and internal energy floor values.

For u(ϕ), we use Equation (3.3) of FM

so the azimuthal component of the four-velocity in BL coordinates is

The temporal component ut is computed using the four-velocity normalization condition uμ uμ = −1. Figure 10 shows the fluid density and magnetic-field initial conditions for several example configurations.

Figure 10. Refer to the following caption and surrounding text.

Figure 10. Azimuthal slices of initial condition for fluid density (normalized and shown as logarithm) and magnetic field in MAD and SANE flows with different disk sizes (indicated by rin/) and adiabatic indices drawn from Table 1. The black hole corotates with the flow and has dimensionless spin parameter a* = 0.94. The black hole event horizon is plotted as a filled gray circle at x = z = 0.

Standard image High-resolution image

Appendix B: Primitive to Conserved Variables: Conversion and Inversion

Hyperbolic systems of conservation laws can be written in the form

where U is the conserved quantity, Fi is the associated flux in the ith spatial direction, and S is the source term for the quantity. The partial differential equation is homogeneous when there are no source terms, leading to U being conserved with time. Comparing Equation (B1) with Equations (28)-(30), the vector of conserved variables for the governing equations of ideal GRMHD is

Computing the fluxes requires solving a local Riemann problem at the zone faces. For the higher-order reconstruction schemes, computing the fluxes is most convenient using an alternate vector of primitive variables,

as the fluxes can be analytically computed from the primitives as F ( P ), unlike F ( U ). Although some primitives have the added benefit that it is easier to understand what they mean physically, note that has been chosen to improve the numerical stability and is not the plasma three velocity:

where α2 ≡ −1/gtt is the lapse.

The conserved variables are complicated, analytical functions of the primitives. iharm evaluates U ( P ) in two steps:

  • 1.  
    Calculate the fluid four-velocity uμ and magnetic induction four-vector bμ from P ,
    where the Lorentz factor with respect to the normal observer can be computed as and βi gti α2 is the shift.
  • 2.  
    Calculate the stress–energy tensor (Equation (32)) from the primitives and four-vectors, and use it to evaluate U .

The inverse operation for the primitives P ( U ) is performed using the "1DW " scheme (Mignone & McKinney 2007; and see also Noble et al. 2006). The matter conserved variables are nonlinear functions of the corresponding primitives, and there are no known analytic expressions for the inverse. The 1DW scheme treats the conserved variables in the normal observer frame and defines a set of scalars, which can be constructed from U , that reduces the 5D system to a 2D system. Of the two equations to be solved, one of them can be analytically inverted to obtain the Lorentz factor of the fluid. The other, which involves the energy density, is a nonlinear expression, and iharm uses a 1D Newton-Raphson scheme to invert it. The magnetic-field primitives are equal to the corresponding conserved variables up to a factor of .

Appendix C: Interpolation and the σ Cutoff

GRMHD codes are unable to robustly model fluid evolution in regions with high magnetization, so they often rely on limiting "flooring" procedures to ensure the numerical stability of the algorithm. In addition to modifying velocities, the flooring procedures introduce extra, artificial mass and energy, so they must be accounted for when performing ray tracing to generate electromagnetic observables. The standard procedure involves masking each fluid snapshot according to a threshold magnetization value σcut. One would hope that the choice of σcut does not materially affect the simulated observables, and in many cases it does not; however, it has been observed to occasionally alter the morphology of images, e.g., in simulations of hot MAD flows (see Chael et al. 2018).

Depending on the interpolation algorithm used to reconstruct fluid variables at nonzone-centered locations and the step size of the integrator, the use of a σ mask may also introduce ridging in simulated images, which is due to the resolution of the underlying fluid snapshot. Figure 11 shows severe (left) and minor (right) examples of the effect. The left panel shows the former default piecewise constant interpolation scheme for σ that zeroed the zone-centered electron number density and introduces sharp boundaries in the image at low intensities. The right panel shows the same image but when σ is (tri-)linearly interpolated, and the mask is applied directly to the emissivity at each geodesic step. In general, the ridges produce artifacts in the Fourier domain that can completely change the signal at long baselines.

Figure 11. Refer to the following caption and surrounding text.

Figure 11. Comparison of log-scaled total intensity image using zone-based vs. interpolated-σ cutoffs. In both panels, the cutoff is implemented by zeroing the electron number density. The sharp edges are from interpolation artifacts due to how linear interpolation deals with rapidly varying quantities on a nearly-regular grid. Notice that ridges are still present in the interpolated-σ cutoff image; they are particularly visible near the bottom of the image. The ridge features may be accentuated when using an integrator that takes large steps.

Standard image High-resolution image

Even in the right panel, ridges can still be seen near the bottom of the image. This is due to the linear interpolation scheme, which produces values whose first derivatives are discontinuous at zone centers. In multiple dimensions, the interpolation artifacts are particularly clear: Figure 12 shows the interpolated values of a smoothly varying scalar that has been sampled only at zone centers. When the interpolated value is used to mask the transfer coefficients along the geodesics (shown in white in the figure), then neighboring geodesics may have noticeably different path lengths.

Figure 12. Refer to the following caption and surrounding text.

Figure 12. Paths of five neighboring geodesics (white lines) plotted over bilinearly interpolated magnetization σ (color, rapidly varying rainbow color scale). The σ = 1 surface is denoted with a solid black line, and zone boundaries are denoted by the dotted black lines. The magnetization decreases steadily from left to right but is defined at zone centers. The value of σ at all other points is reconstructed using bilinear interpolation. The interpolation scheme produces a jagged transition at zone centers, and the slight deviations in path length produce ridges even when the value of σ is determined from an interpolated value.

Standard image High-resolution image

Appendix D: Varying the Density Scale

The output of the GRMHD simulation must be scaled to a particular choice of absolute, physical units in order to perform the radiative transfer calculation. The mass density scaling is chosen to reproduce a target flux density, which is determined by observation. There is a mild degeneracy between the black hole mass M, the distance between the observer and the black hole dsrc, and the mass density unit . Although keeping M/dsrc constant ensures that the angular size of the black hole on the sky remains fixed, changing M alters the physical scale of the system. Thus at constant , increasing M increases, e.g., the optical depth along a geodesic, which can change the image morphology. Thus, at large deviations of M or dsrc, the degeneracy is broken, since is used to fit the observational flux density constraint. A similar analysis as below but in the context of this degeneracy is presented in Appendix A of Roelofs et al. (2021), where both the black hole mass M and the mass density unit are rescaled at a fixed flux density.

The image flux density Fν is not a simple function of primarily because of the complicated relationship between number density, magnetic-field strength, emission, and absorption: increasing alters the magnetic-field strength, which nonlinearly changes the local emissivity; the local absorptivity along the line of sight increases as well, which can occlude parts of the flow. As a result, the relative brightness in different parts of the image changes.

Figure 13 shows the images produced from a single MAD fluid snapshot with different values of , which are centered around the density scale that produces an M87-like flux density at . Different snapshots produce quantitatively distinct but qualitatively similar results. At low , the flow is optically thin, and most emission is produced near the hole. As is increased, optical depths through the disk and jet surface increase, and the primary observed emission shifts to the jet region. Note that, as is increased, the accretion rate increases, cooling becomes more important, and the nonradiative assumption in the GRMHD is violated. Thus, the images with large are likely not representative of physical accretion states.

Figure 13. Refer to the following caption and surrounding text.

Figure 13. Comparison of images produced from the same fluid snapshot across a range of mass density scales , where is normalized to produce an M87-like flux density Fν = 0.8 Jy. Larger values of have higher optical depths; this effect becomes increasingly obvious at , where the flow begins to occlude the central flux depression. The high density scales chosen for the right two panels light up the funnel wall enough to substantially alter the image morphology, so the field of view for those two panels has been increased. Note that the models with large may not be physical since they are produced from fluid simulations that do not include all dynamically important effects at that mass accretion rate, e.g., cooling.

Standard image High-resolution image

Typically, is determined automatically through a numerical root-finding procedure. Because the locus of emission shifts outward at high , the function may appear to have multiple roots (potentially because the image is too small, the simulation domain is too small, or simply due to the changing structure of the emissivity and absorptivity as magnetic-field strengths increase). It is therefore important to verify that the resultant mass accretion rate is reasonable and that the synthetic observables are consistent with observation, e.g., by checking that the image is compact.

Figure 14 shows the relationship between Fν and for the same snapshot shown in Figure 13 for four different values of rhigh. Although the flux density typically increases with and is relatively well-behaved, it is not possible to analytically determine which corresponds to Fν = 1 in general. Nevertheless, given the correct value of for a particular model, it is often possible to use that value as a seed estimate for the root-finding procedure when it is run on similar models.

Figure 14. Refer to the following caption and surrounding text.

Figure 14. Total flux density Fν at ν = 230 GHz vs. for the same snapshot as Figure 13 but for the four different values of rhigh. As before, note that the values of Fν at large may not be representative of physically meaningful scenarios, since cooling was not included in the fluid calculation.

Standard image High-resolution image

How well do images at a given value compare to nearby images? We quantify the difference between a trial image K and the reference image I according to three different metrics: the L1 "taxicab" metric, the mean-squared error across pixels, and the overall structural dissimilarity between the two images. In order to perform this comparison, we first rescale the trial image intensities so that both images have the same total flux density.

We compute the L1 error according to

and we compute the mean-squared error as

where the sums are taken over all pixels i in the image. We compute the structural similarity index measure between the two images as

where

are measures of the mean pixel intensity, variance, and covariance. We report SSIM in terms of the conventional structural dissimilarity

Figure 15 shows the L1, MSE, and DSSIM measures over images from the rhigh = 40 model as compared to the M87-like reference image at Fν = 0.8 Jy.

Figure 15. Refer to the following caption and surrounding text.

Figure 15. Quantitative comparison between fiducial reference image with chosen for Fν = 0.8 Jy and images with other . All comparison images are drawn from a MAD a* = 0.94 model with rhigh = 40, as in Figure 13. As before, images produced with large may not be representative of physically meaningful scenarios.

Standard image High-resolution image

Appendix E: One-zone Model Comparison

We have described the polarized and approximate "unpolarized" radiative transfer equations in Section 2.2, and in Section 3, we have outlined the observer-to-emitter and emitter-to-observer schemes used by ipole and igrmonty, respectively. We now use a simple one-zone model, like the one used in EHTC V, to compare between the different schemes and codes.

Our model test problem comprises a ball in flat space with uniform electron number density ne = 2 × 105 cm−3, temperature Θe = 10, and magnetic field B = 5 G. The magnetic field has been oriented vertically, such that zero emission is produced along the two polar directions. The ball has radius rout = 6.056 ×1013 cm, corresponding to 100 GM/c2 for M = 4.1 × 106 M.

Both igrmonty and ipole are used to sample the full inclination-dependent synchrotron spectrum of the ball. igrmonty only supports the approximate radiative transfer calculation, but we run ipole in both modes for comparison. Since ipole only produces images, we run it multiple times at different inclinations and frequencies to synthesize the full spectrum one point at a time. The ν Lν spectrum is computed from ipole images by . Here, 10−23 converts from the Jansky of spectral flux density, and dsrc = 2.469 × 1022 cm = 8000 kpc is the distance between the camera and the center of the ball.

In igrmonty, emission is sampled at zone centers, so it is important to adequately resolve the width of the photosphere to allow superphotons that undergo significant extinction to be treated correctly. In ipole, images must both have high enough resolution to accurately reproduce gradients as well as a large enough field of view to see the full extent of the ball.

Figure 16 shows the results of the comparison. The igrmonty and approximate "unpolarized" flux densities produced by ipole agree across both inclination and frequency. The polarized and approximate solutions disagree near ν = 1011 Hz. The two methods disagree because the approximate solution does not account for the αQ,U,V polarized absorptivities, and the effect is most significant in the transition between the regimes where the ball is optically thin (ν ≳ 1011 Hz) and optically thick (ν ≲ 1011 Hz).

Figure 16. Refer to the following caption and surrounding text.

Figure 16. Comparison of three radiative transport schemes for the one-zone model. Top panel shows flux densities at the edge-on 90° inclination. Bottom panel shows flux densities for ν = 230 GHz. The approximate unpolarized transport scheme disagrees with the polarized one in the transition region as the ball becomes optically thick.

Standard image High-resolution image

Appendix F: Coordinates

GRMHD simulations are carried out on logically Cartesian grids, meaning that the internal coordinate representation forms a Cartesian grid, and the non-Euclidean effects are accounted for through explicit treatment of geometric terms. KS coordinates are chosen as the standard base coordinate system for fluid simulations because they are horizon penetrating. The eKS-coordinate system is one of the simplest extensions of KS coordinates and uses an exponential radial coordinate , which increases the number of zones at small radii where both the relevant dynamical timescale is shorter and it is more important to recover the detailed dynamics of the flow.

The simulations in this paper were performed in FMKS coordinates , which are an extension to the MKS coordinates introduced in Gammie et al. (2003), which are themselves an extension of eKS coordinates. Positive integer superscripts in this section should be interpreted as indices, not exponents. Coordinate modifications are generally chosen to both reduce the computational cost and increase effective resolution by concentrating zones in the regions of the domain where more interesting physics occurs—like the midplane and near the horizon at small radii—and through "cylindrification," which expands "unnecessary" small zones—like those at the pole at small radii (see Tchekhovskoy et al. 2011). Each of FMKS, MKS, eKS, and KS is axisymmetric in ϕ.

FMKS makes two modifications to the elevation coordinate x2. The first reproduces MKS and increases the number of zones near the midplane by introducing a sinusoidally varying dependence of Δ(x2) on θ, as

where h is the midplane finification parameter, which we set to h = 0.3. Note that there is no finification when h = 1.

FMKS also introduces a cylindrification in θ whereby zones that are near the poles but at small radii have larger extent in θ. This choice increases the length of the required numerical time step, which is set by the minimum of the signal-crossing time over all zones. The signal-crossing time in the zones near the funnel often approaches the speed of light, and thus this fact combined with the structure of spherical geometry (which keeps the number of azimuthal zones constant regardless of θ) results in many small zones with fast signal-crossing times. Thus, through cylindrification, we increase the size of the smallest zones and similarly gain an increase in time step. The cylindrification is achieved by defining

where α and B are parameters and where

is a normalization term. Finally, the KS colatitudinal coordinate is

where measures the FMKS distance from the inner edge of the simulation. In our simulations, we take s = 0.5, B = 0.82, and α = 14. Figure 17 shows an example of the fluid zone boundaries versus KS coordinates for a grid with the above parameters.

Figure 17. Refer to the following caption and surrounding text.

Figure 17. Mapping between Kerr–Schild (KS) coordinates and funky modified Kerr–Schild (FMKS) coordinates. Left: lines of constant FMKS radial coordinate x1 (vertical) and latitudinal coordinate x2 (left-to-right) plotted vs. KS radius r and elevation θ. Right: same as left but plotted in a Cartesian embedding with and . FMKS coordinates concentrate resolution near the midplane θ = π/2 and away from the poles θ = 0, π at small radii.

Standard image High-resolution image

We are not aware of an analytic inversion for , so a nonlinear root-finding step may be required to find the FMKS coordinates for a particular KS event, which is necessary when, e.g., setting the camera position during ray tracing.

Footnotes

  • 17  

    Assuming a distance of 16.8 Mpc (see EHTC I for more detail).

  • 18  

    The magnetorotational instability (MRI; Balbus & Hawley 1991) may play a crucial role in facilitating the angular momentum transport necessary for accretion in steady disks, whose horizon-scale flows are not choked, cf. MAD flows in Section 2.1.2. Even in MAD flows, the MRI may play a role in transient configurations when weaker magnetic fields thread the bulk accretion flow.

  • 19  

    We let Greek indices run over all four dimensions (0, 1, 2, 3) and use Roman indices for the three spatial ones (1, 2, 3).

  • 20  

    In Lorentz–Heaviside units, the vacuum permittivity and permeability are epsilon0 = μ0 = 1 as in Gaussian units, but the strength of the electric and magnetic fields (and the reciprocal of the fundamental charge) differs from Gaussian units by a factor of . In Gaussian units, ϕc ≈ 50.

  • 21  

    The specific intensity is related to the brightness temperature by Iν = 2kTb ν2/c2 = 2Θb me ν2, where in the last term, we used the electron rest-mass energy to construct a dimensionless measure of brightness temperature . Here we have taken the Rayleigh–Jeans limit with h νkT, which is appropriate for mildly relativistic electrons and frequencies near 230 GHz.

  • 22  

    We use the International Astronomical Union convention (see Hamaker & Bregman 1996) in which positive Q is oriented north–south (vertically), and positive U is oriented along the northeast–southwest direction (top left to bottom right). In this convention, EVPA is measured east of north, i.e., counterclockwise from vertical on the sky.

  • 23  

    Note here that the mass scale is independent of the black hole mass because of scale separation—the accretion flow behaves as a test fluid and is effectively nonself-gravitating.

  • 24  

    Some codes instead set the ion temperature equal to the fluid temperature, which overcounts the energy in the system.

  • 25  

    Notice that y and z are not the conventional Y and Z mass fractions.

  • 26  
  • 27  

    In ideal MHD, the plasma is assumed to have infinite conductivity so that the Lorentz force vanishes, with E + v × B = 0. In covariant form, this condition is uμ Fμ ν = 0.

  • 28  

    Eight variables are saved for basic, ideal GRMHD. Extra variables may be tracked in extensions.

  • 29  

    In 3D, a purely fluid torus is unstable to the Papaloizou & Pringle (1984) instability.

  • 30  

    Earlier camera prescriptions defined the camera in the frame with uμ ∝ (1, 0, 0, 0), which is identical in the limit that the camera radius goes to infinity.

  • 31  

    Because we treat synchrotron radiation, it is convenient to construct the tetrad from the fluid uμ , the geodesic wavevector kμ , and the local magnetic-field orientation bμ .

  • 32  

    The original version can be found on Github at https://github.com/moscibrodzka/ipole. The version described here can be found at https://github.com/afd-illinois/ipole as version 1.4. Both versions are currently maintained.

  • 33  

    We have chosen for orbits to be defined with respect to θ rather than ϕ. This choice leads to a subring index that roughly corresponds to the number of times the geodesic has passed through the midplane, although here n may also increase as the geodesic tracks to infinity.

  • 34  
  • 35  

    Here, negative values of spin imply that the angular momentum of the hole and the disk are antiparallel, i.e., the system is retrograde.

10.3847/1538-4365/ac582e
undefined