Comparison of Polarized Radiative Transfer Codes Used by the EHT Collaboration

In 2019, the Event Horizon Telescope (EHT) Collaboration published images of the central black hole (BH) in the galaxy M87 (hereafter M87*), which measured and interpreted the total intensity of radio emission in two bands near 230 GHz (Event Horizon Telescope Collaboration et al. 2019a, 2019b, 2019c, 2019d, 2019e, 2019f, hereafter EHTC I, EHTC II, EHTC III, EHTC IV, EHTC V, EHTC VI). In 2021, additional results were released measuring the degree and distribution of linear polarization across the image of M87*, measured via the Stokes parameters Q and U (Event Horizon Telescope Collaboration et al. 2021a, 2021b, hereafter EHTC VII, EHTC VIII). Linear polarization results are also expected of the central Milky Way BH Sgr A* accompanying total-intensity results published in 2022 (Event Horizon Telescope Collaboration et al. 2022a, 2022b, 2022c, 2022d, 2022e, 2022f).

In order to interpret polarized observations, the collaboration generated models of the accreting plasma around M87*, usually with general relativistic magnetohydrodynamic (GRMHD) simulations. Simulated images were generated from the models via general relativistic radiative transfer (GRRT) calculations in order to predict the emission visible from Earth from the generated plasma state (EHTC V; see also Wong et al. 2022). The total-intensity images produced by various GRRT codes were validated against analytically defined tests in Gold et al. (2020) and found to be in good agreement. That paper also compared the output of certain pairs of codes, but not all codes, when imaging GRMHD simulation data.

The interpretation of the linear-polarimetric EHT image performed in EHTC VIII also used synthetic images. Polarimetric images are more complicated than total-intensity images, as they predict the linear and circular polarization parameters (Stokes Q, U, V) in addition to the total intensity (Stokes I). Predicting polarized emission involves solving the coupled polarized radiative transfer equations, which can introduce significant additional computational problems, such as the treatment of rapid Faraday rotation and the need to parallel transport the linear polarization direction through curved spacetime. The additional complexity merits a separate comparison of polarized radiative transfer schemes present in several of the codes compared in Gold et al. (2020). That comparison is presented in this paper.

This paper provides brief descriptions of the codes compared, specifications of the tests performed, and measurements of code error (where available) or similarity as a group, compared against parameter changes and estimated detector accuracy.

The paper is structured as follows. In Section 2 we briefly describe all codes participating in the comparison study. In Section 3 we define three test problems used to compare the codes. In Section 4 we define a metric to evaluate light curves and image similarities and present the results of the comparisons. The discussion of the results and limitations of the examined ray-tracing radiative transfer schemes is given in Section 5. We conclude our study in Section 6.

2.1.  BHOSS

2.2.  ipole

2.3.  ipole-IL

2.4.  grtrans

2.5.  Odyssey

2.6.  RAPTOR

Three test problems were used to evaluate the codes. The problems were chosen to reflect tests already present in the literature, highlighting specifically the aspects of code performance related to polarized transport. A previous comparison (Gold et al. 2020) evaluated many of the same codes for similarity and accuracy in producing total-intensity images. An additional goal was to verify data product similarity when imaging the output of GRMHD simulations, a test evaluated only for certain pairs of codes considered in Gold et al. (2020). The tests are described here from least to most complex, with each testing a larger subset of code features.

3.1. Comparison to Analytic Result

The first test problem is a straightforward integration of the nonrelativistic polarized transfer equation using constant coefficients, chosen for the availability of an analytic solution from Landi Degl'Innocenti & Landi Degl'Innocenti (1985), allowing direct evaluation of code accuracy, in addition to code similarity. The test here is taken directly from Dexter (2016), with coefficients as listed in Mościbrodzka & Gammie (2018).

In the Stokes basis I, Q, U, V, the nonrelativistic polarized radiative transfer equation is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d}{{ds}}\left(\begin{array}{l}I\\ Q\\ U\\ V\end{array}\right)=\left(\begin{array}{l}{j}_{I}\\ {j}_{Q}\\ {j}_{U}\\ {j}_{V}\end{array}\right)\\ -\left(\begin{array}{cccc}{\alpha }_{I} & {\alpha }_{Q} & {\alpha }_{U} & {\alpha }_{V}\\ {\alpha }_{Q} & {\alpha }_{I} & {\rho }_{V} & -{\rho }_{U}\\ {\alpha }_{U} & -{\rho }_{V} & {\alpha }_{I} & {\rho }_{Q}\\ {\alpha }_{V} & {\rho }_{U} & -{\rho }_{Q} & {\alpha }_{I}\end{array}\right)\left(\begin{array}{l}I\\ Q\\ U\\ V\end{array}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

In the test, this equation is integrated twice with different subsets of nonzero coefficients. This minimizes the complexity of the analytic comparison functions, isolating any bugs in treating emission and absorption from those in Faraday rotation and conversion. The coefficients for each integration are given in Table 1.

Table 1. Constant Coefficients for the Analytic Comparison Test Integrations

	jI	jQ	jU	jV	αI	αQ	αU	αV	ρQ	ρU	ρV
Emission/Absorption	2.0	1.0	0.0	0.0	1.0	1.2	0.0	0.0	0.0	0.0	0.0
Rotation	0.0	0.1	0.1	0.1	0.0	0.0	0.0	0.0	10.0	0.0	−4.0

Note. These mirror values in Mościbrodzka & Gammie (2018). See Section 3.1.

Download table as:  ASCIITypeset image

3.2. Thin-disk Model

The second test problem consists of imaging emission from a thin opaque disk aligned to the midplane of a near-maximally spinning BH, as described in Novikov & Thorne (1973). This involves solving the geodesic equation in two contexts: first in tracing lines of sight from the camera through the Kerr metric, and second in parallel-transporting the direction of linearly polarized emission from the disk back to the camera. This test closely mirrors a figure from Schnittman & Krolik (2009), which is reproduced as a test in Dexter (2016).

The thin-disk test does not include any diffuse emission, that is, all transport coefficients j_S , α_S , ρ_S = 0 uniformly for all Stokes parameters S. Instead, the initial Stokes parameters are set as a boundary condition at the first midplane crossing of each geodesic when traced backward from the camera. As in Schnittman & Krolik (2009), the total flux F is taken from Page & Thorne (1974), and the intensity at the desired frequency I_ν is obtained by calculating an effective temperature and assuming a blackbody distribution diluted by a hardening factor n = 1.8:

$$\begin{eqnarray}&&{T}_{\mathrm{eff}}\equiv {\left(\displaystyle \frac{F}{\sigma }\right)}^{1/4},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{I}_{\nu }=\displaystyle \frac{1}{{n}^{4}}{B}_{\nu }(n\cdot {T}_{\mathrm{eff}}),\end{eqnarray} \tag{ 3 }$$

where F is the total emitted power per area of the thin disk and B_ν (T) is the blackbody function of the temperature T and the emitted frequency ν in the fluid frame.

The emitted intensity and horizontal polarization fraction in the outgoing direction are determined by assuming scattering from a semi-infinite atmosphere, as in Chandrasekhar (1960), Table 24, with the direction of linear polarization pointing along the plane of the disk. Emission is enabled only between r_ISCO and R_out = 100r_g , where r_g is the system gravitational radius GM_BH/c², with G the gravitational constant and c the speed of light. The fluid orbital angular velocity ${}_{{u}^{t}}^{{u}^{\phi }}$ is assumed to be Keplerian:

$$\begin{eqnarray}&&\frac{{u}^{\phi }}{{u}^{t}}=\frac{1}{{r}^{3/2}+{a}_{* }},\end{eqnarray} \tag{ 4 }$$

where a_* is the dimensionless form of the BH angular momentum J, a_* ≡ Jc/GM² with −1 ≤ a_* ≤ 1. As the test will need to be implemented in many different codes, we simplify the original problem from Schnittman & Krolik (2009) by observing at only a single frequency rather than summing over a range. The full set of parameters used for this image is

$$\begin{eqnarray}&&{a}_{* }=0.99\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}=10\,{M}_{\odot }=1.477\cdot {10}^{6}\mathrm{cm}\displaystyle \frac{{c}^{2}}{G}\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}&&{D}_{\mathrm{source}}=0.05\,\mathrm{pc}\end{eqnarray} \tag{ 5c }$$

$$\begin{eqnarray}&&\dot{M}=0.1\,{\dot{M}}_{\mathrm{Edd}}\approx 2.218\cdot {10}^{-9}{M}_{\odot }\,{\mathrm{yr}}^{-1}\end{eqnarray} \tag{ 5d }$$

$$\begin{eqnarray}&&h\nu =1\mathrm{keV},\end{eqnarray} \tag{ 5e }$$

where these parameters define observation at a single frequency ν of a BH of mass M_BH at distance D_source, characterized by a Kerr spacetime with BH spin parameter a_*, and accreting at rate $\dot{M}$. Following Dexter and Schnittman, ${\dot{M}}_{\mathrm{Edd}}{\equiv }_{{c}^{2}}^{{L}_{\mathrm{Edd}}}$.

Example output from this test run with ipole-IL is shown in Figure 1. As no circularly polarized emission or Faraday conversion occurs in the problem, the Stokes V flux remains exactly zero.

Zoom In Zoom Out Reset image size

3.2.1. Camera

In this and the following test, the camera tetrad is constructed such that a geodesic at the center of the camera's field of view (FOV) would have zero angular momentum k_ϕ . The polar angle θ_cam is defined relative to the BH angular momentum vector. This is identical to the camera definition from Gold et al. (2020).

The camera is placed at radius 10⁴ r_g in both tests, to reduce the discrepancy between pinhole and planar cameras. Note that R_cam is in this case and in practice much smaller than D_source—so long as R_cam is large enough to eliminate camera effects, the image intensity is invariant with distance.

The FOV of each test is given in two forms: DX, the in-plane distance from one edge of the imaged material to the other in gravitational radii r_g , and FOV, the angular size from Earth in microarcseconds (μs).

Pixels on each image correspond to evenly spaced geodesics, starting from each pixel center. That is, an image with an FOV of 80 μs to a side and a resolution NX of 80 pixels would be calculated using geodesics originating at −39.5 to 39.5 μs away from the FOV center, spaced 1 μs apart in each cardinal direction. All images in these tests are square, with equal FOV from west to east and north to south.

For the thin-disk test, the camera parameters are

$$\begin{eqnarray}&&{\theta }_{\mathrm{cam}}=75^\circ \end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&\mathrm{NX}=80\,\mathrm{pixels}\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}&&\mathrm{DX}=40\,{r}_{g}\end{eqnarray} \tag{ 6c }$$

$$\begin{eqnarray}&&\mathrm{FOV}=78.99\,\mu \mathrm{as}.\end{eqnarray} \tag{ 6d }$$

3.3. GRMHD Snapshot

The last test consists of imaging the relativistic thermal synchrotron emission at 230 GHz from one snapshot from a GRMHD simulation. This test exercises all aspects of the code, as well as code-specific choices, such as interpolation of fluid state recorded at discrete locations and calculation of the transport coefficients via fitting functions. The standard snapshot file used for this test is taken from a SANE simulation with spin a_* = 0.9375, performed using iharm3D (Prather et al. 2021) with a resolution of 288 × 128 × 128 cells in r, θ, and ϕ respectively. Except for the coordinate system, this simulation exactly reflects the simulations performed with iharm3D as a part of the library used in EHTC V and EHTC VIII, as described in Wong et al. (2022). Not all GRRT codes can read all GRMHD output, as coordinate systems and fluid state descriptions can differ from code to code—thus, not every polarized radiative transfer code used in the EHTC can be directly compared with this test. The iharm3D format is chosen since it is readable by a majority of codes used in the EHTC, and in particular those codes relevant to studies in EHTC VIII.

The snapshot is taken at 4500 r_g /c after simulation start, well into the run's quiescent period. The file is available upon request for testing future codes. Since it involves creating an image from just one snapshot of the simulation, the test makes the assumption of "fast light," i.e., that the fluid is static as light propagates from emission to observer.

The parameters of this test are chosen to reflect values for M87*, specifically those used in creating the libraries of simulated images used in EHTC V, EHTC VI and EHTC VIII, hereafter the "EHT image libraries." These are

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}=6.2\times {10}^{9}\,{M}_{\odot }\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&D=16.9\,\mathrm{Mpc}\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&\nu =230\,\mathrm{GHz}.\end{eqnarray} \tag{ 7c }$$

The camera is defined as in Section 3.2.1, with parameters chosen to reflect the angle of the M87* jet and the FOV observed by the EHT:

$$\begin{eqnarray}&&{\theta }_{\mathrm{cam}}=163^\circ \end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&{\phi }_{\mathrm{cam}}=0^\circ \end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&\mathrm{FOV}=160\,\mu \mathrm{as}\end{eqnarray} \tag{ 8c }$$

$$\begin{eqnarray}&&\mathrm{DX}=44.17\,{r}_{g}\end{eqnarray} \tag{ 8d }$$

$$\begin{eqnarray}&&\mathrm{NX}=160\,\mathrm{pixels}.\end{eqnarray} \tag{ 8e }$$

In addition to the system parameters above, imaging a GRMHD simulation necessarily involves setting another scale factor that determines the density of accreting material and the strength of magnetic fields. It is expressed here as a mass unit, ${ \mathcal M }$, which gives units to the unscaled density values from a simulation, ρ_code:

$$\begin{eqnarray}&&{\rho }_{\mathrm{CGS}}=\displaystyle \frac{{ \mathcal M }}{{r}_{g}^{3}}{\rho }_{\mathrm{code}}.\end{eqnarray} \tag{ 9 }$$

${ \mathcal M }$ is not known a priori, but it is highly correlated with the total image brightness. Thus, it is scaled so as to match the total image flux density to the observed compact flux density, usually by employing an iterative solver.

In the EHT image libraries, ${ \mathcal M }$ was fit such that images taken over the course of a full simulation would produce an average of 0.5 Jy of compact flux density (see Wong et al. 2022 for details). For this test, the sample image is fit alone using ipole-IL such that it produces 0.50 Jy when imaged with totally unpolarized transport, or about 0.47 Jy of Stokes I flux density when imaged using polarized transport. The value of ${ \mathcal M }$ and the corresponding accretion rate $\dot{M}$ used for the test are listed below:

$$\begin{eqnarray}&&{ \mathcal M }=1.672\cdot {10}^{26}\,{\rm{g}}\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&\dot{M}=8.644\cdot {10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}=6.285\cdot {10}^{-7}\,{\dot{M}}_{\mathrm{Edd}},\end{eqnarray} \tag{ 10b }$$

where ${\dot{M}}_{\mathrm{Edd}}$ is defined as earlier in this work.

Accretion flows around M87* are strongly suspected to be two-temperature, with little thermal coupling between the ions and electrons (Mahadevan & Quataert 1997; Ryan et al. 2017; Sadowski et al. 2017). Since GRMHD simulations evolve only a single fluid with a single temperature, when simulating images the internal energy must be split between the ions and electrons based on a model. While this process is not well constrained, it is generally documented which electron distribution is being assumed, and the electron energy distribution model is standardized between codes when similar performance is expected.

Thus, for simplicity, in this test we set the electron temperature to a fixed ratio of 1/3 of the ion temperature, derived from the single-fluid GRMHD parameters by holding the total internal energy constant (see EHTC V):

$$\begin{eqnarray}&&{T}_{e}=\displaystyle \frac{2{m}_{p}u}{15k\rho },\end{eqnarray} \tag{ 11 }$$

where u and ρ are the local fluid internal energy and rest-mass density per unit volume, respectively, and m_p and k are the proton mass and Boltzmann constant, respectively. Note that this is equivalent to the so-called "R_high" model of EHTC V with R_low = R_high = 3. In splitting the total internal energy rather than setting the fluid and ion temperatures equal, it differs slightly from the original statement of the R_high model in Mościbrodzka et al. (2016).

The emission, absorption, and rotation coefficients are calculated based on the electron temperature (or more broadly, the electron energy distribution) using fitting functions approximating the full synchrotron emission calculations, which are expensive to compute. Codes in this comparison used a few different sets of fitting functions; further discussion is found in Section 5.3 and Appendix B.

Finally, as in EHTC V and commonly in the literature, emission is tracked only from regions of the simulation with σ ≡ B²/ρ < 1. Regions with higher sigma (largely the polar "jet" regions) can overproduce emission if included, due to hot material in the jet inserted by numerical floors to preserve stability of GRMHD algorithms.

Example output for the GRMHD snapshot test from ipole-IL is provided in Figure 2.

Zoom In Zoom Out Reset image size

4.1. Analytic Comparison Results

Results for the analytic integration tests from RAPTOR, Odyssey, ipole, and ipole-IL are shown in Figure 3. Raw output is plotted in the left panels, and differences from the analytic result are plotted in the right panels. This test verifies that the default accuracy parameters of each code allow them to match an analytic solution to within acceptable errors. Note that this is not a good measure of relative code accuracy or convergence—for convergence tests, see the accompanying code papers cited in Section 2.

Zoom In Zoom Out Reset image size

Note that the results from two integrators are shown for RAPTOR as the "RK4" and "IE" variants. In normal integration, RAPTOR uses the "RK4" integrator, reserving the "IE" integrator for the few zones where Faraday rotation is too strong to take steps of an appropriate size with an explicit scheme (Bronzwaer et al. 2020).

In addition, the ipole scheme (also used in ipole-IL) is semianalytic: it uses the analytic solution for constant coefficients whenever it evolves the nonrelativistic polarized radiative transfer equations. Thus, ipole and ipole-IL will perform this test less accurately when taking more steps: a single step of any size would be exactly accurate, but multiple steps accrue round-off error.

4.2. Metrics

In the following two imaging tests, no exact result is available by which to evaluate code accuracy directly; rather, we evaluate consistency between all codes, both in overall image structure and in several metrics used to compare models with EHT results. In particular, we will use the definitions from EHTC VIII, computed over simulated images and used to compare models to the observed EHT result. These summary statistics include the total flux density F, image-integrated or "zero-baseline" linear and circular polarization fractions ∣m∣_net and v_net, and the average linear polarization fraction over the resolved image 〈∣m∣〉. For an image represented as a vector of emitted flux per pixel in each Stokes parameter I_j , Q_j , U_j , V_j over each pixel j, these values are defined as

$$\begin{eqnarray}&&{\rm{F}}=\displaystyle \sum _{j}{I}_{j}\end{eqnarray} \tag{ 12 }$$

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

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

$$\begin{eqnarray}&&{v}_{\mathrm{net}}=\displaystyle \frac{{\sum }_{j}{V}_{j}}{{\sum }_{j}{I}_{j}}.\end{eqnarray} \tag{ 15 }$$

Additionally, EHTC VIII used a complex coefficient reflecting the degree and angle of azimuthally symmetric linear polarization, β₂ (see Palumbo et al. 2020). It is calculated by first centering the image, e.g., with rex (EHTC IV), and then taking the inner product of the complex linear polarization P = Q + iU with a rotationally symmetric function:

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

$$\begin{eqnarray}&&{I}_{\mathrm{ann}}={\int }_{{\rho }_{\min }}^{{\rho }_{\max }}{\int }_{0}^{2\pi }I(\rho ,\phi )\,\rho \,d\phi d\rho ,\end{eqnarray} \tag{ 17 }$$

where ρ/ϕ are polar coordinates in the image plane, measured from/about the image center. This metric is expected to be useful only for images with a relatively low observer angle, as they will be more symmetric; thus, it is computed and compared only for the GRMHD snapshot test, which uses the low observer angle i = 17° expected for M87* and used for libraries of simulated images of that object.

The quantities 〈∣m∣〉 and β₂ are sensitive to image resolution. In order to mirror EHT measurements and reflect how simulated images were used in EHTC VIII, a circular Gaussian blur with an FWHM of 20 μas was applied to all images before computing either resolution-dependent quantity.

In addition to the quantities used for direct comparison in EHTC VIII, we measure a point-source linear polarization direction or electric vector position angle (EVPA) east of due north on the sky, and thus in our Stokes convention defined as

$$\begin{eqnarray}&&\mathrm{EVPA}=\displaystyle \frac{1}{2}\arg \left(\displaystyle \frac{\displaystyle \sum _{j}{U}_{j}}{\displaystyle \sum _{j}{Q}_{j}}\right).\end{eqnarray} \tag{ 18 }$$

As this metric is potentially volatile for images with low net linear polarization, and a corresponding measurement has not been made to which we might compare, we follow EHTC VIII in omitting this as a comparison metric—rather, we use it only in image summaries.

When evaluating image similarity, we use the normalized mean squared error (NMSE):

$$\begin{eqnarray}&&\mathrm{NMSE}(A,B)=\displaystyle \frac{\displaystyle \sum _{j}| {A}_{j}-{B}_{j}{| }^{2}}{\displaystyle \sum _{j}| {A}_{j}{| }^{2}},\end{eqnarray} \tag{ 19 }$$

where A_j and B_j are the intensities of a particular Stokes parameter in two images at pixel j. Regardless of exact nomenclature, all NMSE values listed in this work and in Gold et al. (2020) use this normalization and are thus comparable.

Note that this definition of the NMSE is not symmetric under the ordering of A and B, and in particular as images get dimmer, NMSE(I, 0) = 1, whereas NMSE(0, I) = ∞ . As it is normalized against the sum of squared pixel intensities, the NMSE becomes more volatile when evaluating dimmer images.

The NMSE was one of two metrics used to gauge similarity in Gold et al. (2020). We omit the other, the structural dissimilarity (DSSIM), since for the case of very similar images, values of the DSSIM are highly correlated with the NMSE (see Appendix C).

4.3. Thin-disk Test Results

Each code's output for the thin-disk test is plotted in Figure 4. The images are visually indistinguishable except in a few particular pixels, and this similarity is borne out in the comparison metrics. Stokes V is omitted from plots and comparisons for this test, as all codes produce exactly zero Stokes V across the entire image.

Zoom In Zoom Out Reset image size

Recall that this test involves accurate evaluation of the geodesic equation and accurate parallel transport of the linear polarization vector from emission to camera. As the codes' similarity in tracing geodesics was evaluated extensively in Gold et al. (2020), we focus here on demonstrating the latter through comparison of the resulting Stokes Q, U images and the relevant image-integrated metrics. Table 2 lists total fluxes and net polarization parameters for each image in the test. Figure 5 presents comparisons of each metric between each pair of images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Image-integrated Values for the Thin-disk Test

Code	Flux (Jy)	∣m∣net (%)	〈∣m∣〉 (%)	EVPA (deg)
ipole	6.841e+06	2.3341	2.378	88.123
ipole-IL	6.8699e+06	2.3224	2.3707	87.974
grtrans	6.8227e+06	2.3246	2.3709	88.01
RAPTOR	6.8689e+06	2.3265	2.3726	88.025
Odyssey	6.7338e+06	2.3527	2.3971	88.146
BHOSS	6.6949e+06	2.3439	2.3845	88.131

Note. v_net is omitted since it is uniformly zero, and β₂ magnitude and angle are omitted since the image is not symmetric and thus the magnitude is very small.

Download table as:  ASCIITypeset image

4.4. GRMHD Snapshot Test Results

Results for the GRMHD snapshot test for ipole, ipole-IL, grtrans, Odyssey, and RAPTOR are listed in Table 3 and presented in Figures 6 and 7. Results are formatted similarly to the results of the thin-disk test, with the addition of a Stokes V component and circular polarization fraction in images and tables and the addition of the rotationally symmetric linear polarization coefficient β₂ in tables of integrated values.

Zoom In Zoom Out Reset image size

Table 3. Image-integrated Values for the GRMHD Snapshot Test

Code	Flux (Jy)	∣m∣net (%)	vnet (%)	〈∣m∣〉 (%)	EVPA (deg)	$\left|{\beta }_{2}\right|$	∠β2 (deg)
ipole-IL	0.47976	1.5232	0.64637	31.264	−78.241	0.28165	−20.307
ipole	0.47335	1.5904	0.62402	31.27	−77.243	0.28171	−20.622
RAPTOR	0.49396	1.6057	0.67381	31.893	−77.797	0.28346	−21.31
grtrans	0.45346	1.5135	0.69909	31.197	−77.647	0.28063	−17.903
Odyssey	0.46671	1.3928	0.70068	31.54	−71.432	0.28325	−20.713

Note. Definitions for all values are given in Section 4.2. Note that the EVPA is not used as a comparison metric in Figure 7.

Download table as:  ASCIITypeset image

Except for the total flux density, the color bars in Figure 7 reflect the 1σ values used to make cuts when evaluating models in EHTC VIII.

5.1. Comparison to Observational Constraints

The values for each comparison metric used as cuts in EHTC VIII are listed in Table 4. The table values are based on 1σ ranges for measurements of the same quantities in EHT data, described in EHTC VII. These ranges provide a comparison to evaluate code interchangeability—if the differences between codes are substantially less than the range of measurement uncertainties, the analysis is agnostic to the choice of code employed. The 1σ ranges listed are also used as the color bar ranges in the colored table listings in Figure 7.

Table 4. Detector Uncertainty Compared Against Code Uncertainty in Observed Parameters

Parameter	1σ	Code Uncertainty	Uncertainty/σ
∣m∣net	1.35 %	0.21 %	0.16
vnet	0.4 %	0.08 %	0.20
〈∣m∣〉	2.5 %	0.70 %	0.28
$\left|{\beta }_{2}\right|$	0.015	0.0026	0.17
∠β2	17°	3.4°	0.20

Note. The "1σ" column lists the 1σ detector uncertainty of each parameter, as estimated in EHTC VII and used as an allowable range when womparing models in EHTC VIII. The measurement of 〈∣m∣〉 was an upper bound; this upper bound was doubled to select the cut value. The "Code Uncertainty" column lists the greatest observed difference between codes when computing each parameter, and the final column lists this value as a proportion of the 1σ value.

Download table as:  ASCIITypeset image

This comparison provides evidence that model evaluations as in EHTC VIII remain similar regardless of which of the included codes is employed. As recorded in Table 4, maximum code variation is universally less than 30% of the detector 1σ range: in ∣m∣_net (0.21% vs. 1.5% ), v_net (0.08% vs. 0.4%), 〈∣m∣〉 (0.70% vs. 2.5%), $\left|{\beta }_{2}\right|$ (0.0026 vs. 0.015), and ∠β₂ (3.4° vs. 17°). In any analysis based on cuts, code differences can shift a few particular images into or out of the final consideration. However, at these uncertainties no image from outside the 1σ detector uncertainty would be consistent with the central observed value.

Broadening the comparison to different images and models shows promising similarities. Appendix D presents distributions of the image differences between ipole-IL and grtrans when run over thousands of snapshots of a very different model from the example: they show a wide variance but a smaller difference on average than in the example image, suggesting that image differences, at least between these codes, are mostly stochastic, further suppressing any potential effect on a cuts-based analysis as in EHTC VIII.

5.2. Potential Measurement of System Parameters

To translate the NMSE into a measure of code accuracy in testing model parameters, we define an "error budget" for each Stokes parameter, consisting of the largest NMSE between code results: 0.02 in I, 0.04 in Q and U, and 0.13 in V. Assuming perfect detector accuracy and modeling, this error characterizes which images are too similar to be effectively distinguished above code-to-code variations.

We then translate this error budget into constraints on the input parameters M_BH, ${ \mathcal M }$, R_high, and the viewing angle by varying these parameters around the nominal values and calculating the resulting NMSE versus the nominal image, using ipole-IL. As illustrated in Figure 8, the required parameter changes are very modest; that is, the possible constraints on system parameters are very precise. In imaging the example model, codes agree well enough to constrain the mass of M87* to within 0.4% (2.5 × 10⁷ M_⊙), ${ \mathcal M }$ to 9% (1.5 × 10²⁵), the observer angle to 0.8°, and the R_high parameter to within 0.18. These values are dependent on the base image—in particular, constraints on R_high will also depend on R_low, which was set differently in this case than for images in EHTC VIII.

Zoom In Zoom Out Reset image size

5.3. Caveats and Limitations

In each of the tests conducted for this comparison, the several GRRT codes used within the EHT Collaboration have produced sufficiently similar results that they are functionally interchangeable for the collaboration's uses. This is true both when measured in terms of image similarity (NMSE) and when measured directly in terms of the image metrics used to compare simulated polarized images to the EHT result in EHTC VIII.

Using their default accuracy parameters, codes match the analytic result for the case of constant transport coefficients to better than 1 part in 10⁻⁵. They agree to within NMSE of 0.015 when imaging an analytically defined problem requiring parallel transport of the polarization vector.

In the more complex task of interpolating, translating, and imaging GRMHD output, codes agree to within an NMSE of 0.13 at worst, when measuring specifically the circular polarization map (NMSE of 0.045 in linear polarization, 0.02 in total intensity). Based on image similarity, the choice of imaging code will matter in model comparisons only when trying to determine the BH mass to within 0.4%, accretion rate within 9%, or observer angle to within 0.8°. These values significantly outclass both the detection and modeling uncertainties available in the near future.

When measured with the image metrics used for model comparison in EHTC VIII, all comparison images agree to much better precision than the detector uncertainty. Further, much of the difference that does appear is shown to be stochastic in nature. Thus, the choice of code is verified directly to have little effect on the analysis performed in that paper.

The Event Horizon Telescope Collaboration thanks the following organizations and programs: the Academia Sinica; the Academy of Finland (projects 274477, 284495, 312496, 315721); the Agencia Nacional de Investigación y Desarrollo (ANID), Chile via NCN19₀₅₈ (TITANs) and Fondecyt 1221421, the Alexander von Humboldt Stiftung; an Alfred P. Sloan Research Fellowship; Allegro, the European ALMA Regional Centre node in the Netherlands, the NL astronomy research network NOVA and the astronomy institutes of the University of Amsterdam, Leiden University, and Radboud University; the ALMA North America Development Fund; the Black Hole Initiative, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation (although the opinions expressed in this work are those of the author(s) and do not necessarily reflect the views of these Foundations); the Brinson Foundation; Chandra DD7-18089X and TM6-17006X; the China Scholarship Council; the China Postdoctoral Science Foundation fellowships (2020M671266, 2022M712084); Consejo Nacional de Ciencia y Tecnología (CONACYT, Mexico, projects U0004-246083, U0004-259839, F0003-272050, M0037-279006, F0003-281692, 104497, 275201, 263356); the Consejería de Economía, Conocimiento, Empresas y Universidad of the Junta de Andalucía (grant P18-FR-1769), the Consejo Superior de Investigaciones Científicas (grant 2019AEP112); the Delaney Family via the Delaney Family John A. Wheeler Chair at Perimeter Institute; Dirección General de Asuntos del Personal Académico-Universidad Nacional Autónoma de México (DGAPA-UNAM, projects IN112417 and IN112820); the Dutch Organization for Scientific Research (NWO) for VICI award (grant 639.043.513), grant OCENW.KLEIN.113 and the Dutch Black Hole Consortium (with project No. NWA 1292.19.202) of the research program the National Science Agenda; the Dutch National Supercomputers, Cartesius and Snellius (NWO grant 2021.013); the EACOA Fellowship awarded by the East Asia Core Observatories Association, which consists of the Academia Sinica Institute of Astronomy and Astrophysics, the National Astronomical Observatory of Japan, Center for Astronomical Mega-Science, Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute; the European Research Council (ERC) Synergy Grant "BlackHoleCam: Imaging the Event Horizon of Black Holes" (grant 610058); the European Union Horizon 2020 research and innovation program under grant agreements RadioNet (No. 730562) and M2FINDERS (No. 101018682); the Horizon ERC Grants 2021 program under grant agreement No. 101040021; the Generalitat Valenciana postdoctoral grant APOSTD/2018/177 and GenT Program (project CIDEGENT/2018/021); MICINN Research Project PID2019-108995GB-C22; the European Research Council for advanced grant "JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales" (grant No. 884631); the Institute for Advanced Study; the Istituto Nazionale di Fisica Nucleare (INFN) sezione di Napoli, iniziative specifiche TEONGRAV; the International Max Planck Research School for Astronomy and Astrophysics at the Universities of Bonn and Cologne; DFG research grant "Jet physics on horizon scales and beyond" (grant No. FR 4069/2-1); Joint Columbia/Flatiron Postdoctoral Fellowship (research at the Flatiron Institute is supported by the Simons Foundation); the Japan Ministry of Education, Culture, Sports, Science and Technology (MEXT; grant JPMXP1020200109); the Japanese Government (Monbukagakusho:MEXT) Scholarship; the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for JSPS Research Fellowship (JP17J08829); the Joint Institute for Computational Fundamental Science, Japan; the Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS, grants QYZDJ-SSW-SLH057, QYZDJSSW-SYS008, ZDBS-LY-SLH011); the Leverhulme Trust Early Career Research Fellowship; the Max-Planck-Gesellschaft (MPG); the Max Planck Partner Group of the MPG and the CAS; the MEXT/JSPS KAKENHI (grants 18KK0090, JP21H01137, JP18H03721, JP18K13594, 18K03709, JP19K14761, 18H01245, 25120007); the Malaysian Fundamental Research Grant Scheme (FRGS) FRGS/1/2019/STG02/UM/02/6; the MIT International Science and Technology Initiatives (MISTI) Funds; the Ministry of Science and Technology (MOST) of Taiwan (103-2119-M-001-010-MY2, 105-2112-M-001-025-MY3, 105-2119-M-001-042, 106-2112-M-001-011, 106-2119-M-001-013, 106-2119-M-001-027, 106-2923-M-001-005, 107-2119-M-001-017, 107-2119-M-001-020, 107-2119-M-001-041, 107-2119-M-110-005, 107-2923-M-001-009, 108-2112-M-001-048, 108-2112-M-001-051, 108-2923-M-001-002, 109-2112-M-001-025, 109-2124-M-001-005, 109-2923-M-001-001, 110-2112-M-003-007-MY2, 110-2112-M-001-033, 110-2124-M-001-007, and 110-2923-M-001-001); the Ministry of Education (MoE) of Taiwan Yushan Young Scholar Program; the Physics Division, National Center for Theoretical Sciences of Taiwan; the National Aeronautics and Space Administration (NASA, Fermi Guest Investigator grant 80NSSC20K1567, NASA Astrophysics Theory Program grant 80NSSC20K0527, NASA NuSTAR award 80NSSC20K0645); NASA Hubble Fellowship grants HST-HF2-51431.001-A, HST-HF2-51482.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555; the National Institute of Natural Sciences (NINS) of Japan; the National Key Research and Development Program of China (grant 2016YFA0400704, 2017YFA0402703, 2016YFA0400702); the National Science Foundation (NSF, grants AST-0096454, AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-1126433, AST-1140030, DGE-1144085, AST-1207704, AST-1207730, AST-1207752, MRI-1228509, OPP-1248097, AST-1310896, AST-1440254, AST-1555365, AST-1614868, AST-1615796, AST-1715061, AST-1716327, AST-1716536, OISE-1743747, AST-1816420, AST-1935980, AST-2034306); NSF Astronomy and Astrophysics Postdoctoral Fellowship (AST-1903847); the Natural Science Foundation of China (grants 11650110427, 10625314, 11721303, 11725312, 11873028, 11933007, 11991052, 11991053, 12192220, 12192223); the Natural Sciences and Engineering Research Council of Canada (NSERC, including a Discovery Grant and the NSERC Alexander Graham Bell Canada Graduate Scholarships-Doctoral Program); the National Youth Thousand Talents Program of China; the National Research Foundation of Korea (the Global PhD Fellowship Grant: grants NRF-2015H1A2A1033752, the Korea Research Fellowship Program: NRF-2015H1D3A1066561, Brain Pool Program: 2019H1D3A1A01102564, Basic Research Support Grant 2019R1F1A1059721, 2021R1A6A3A01086420, 2022R1C1C1005255); Netherlands Research School for Astronomy (NOVA) Virtual Institute of Accretion (VIA) postdoctoral fellowships; Onsala Space Observatory (OSO) national infrastructure, for the provisioning of its facilities/observational support (OSO receives funding through the Swedish Research Council under grant 2017-00648); the Perimeter Institute for Theoretical Physics (research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science); the Princeton Gravity Initiative; the Spanish Ministerio de Ciencia e Innovación (grants PGC2018-098915-B-C21, AYA2016-80889-P, PID2019-108995GB-C21, PID2020-117404GB-C21); the University of Pretoria for financial aid in the provision of the new Cluster Server nodes and SuperMicro (USA) for a SEEDING GRANT approved toward these nodes in 2020; the Shanghai Pilot Program for Basic Research, Chinese Academy of Science, Shanghai Branch (JCYJ-SHFY-2021-013); the State Agency for Research of the Spanish MCIU through the "Center of Excellence Severo Ochoa" award for the Instituto de Astrofísica de Andalucía (SEV-2017- 0709); the Spinoza Prize SPI 78-409; the South African Research Chairs Initiative, through the South African Radio Astronomy Observatory (SARAO, grant ID 77948), which is a facility of the National Research Foundation (NRF), an agency of the Department of Science and Innovation (DSI) of South Africa; the Toray Science Foundation; the Swedish Research Council (VR); the US Department of Energy (USDOE) through the Los Alamos National Laboratory (operated by Triad National Security, LLC, for the National Nuclear Security Administration of the USDOE, contract 89233218CNA000001); the YCAA Prize Postdoctoral Fellowship; and the Astrophysics and High Energy Physics program by MCIN (with funding from European Union NextGenerationEU, PRTR-C17I1).

We thank the staff at the participating observatories, correlation centers, and institutions for their enthusiastic support. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.01154.V. ALMA is a partnership of the European Southern Observatory (ESO; Europe, representing its member states), NSF, and National Institutes of Natural Sciences of Japan, together with National Research Council (Canada), Ministry of Science and Technology (MOST; Taiwan), Academia Sinica Institute of Astronomy and Astrophysics (ASIAA; Taiwan), and Korea Astronomy and Space Science Institute (KASI; Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, Associated Universities, Inc. (AUI)/NRAO, and the National Astronomical Observatory of Japan (NAOJ). The NRAO is a facility of the NSF operated under cooperative agreement by AUI. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC05-00OR22725. We also thank the Center for Computational Astrophysics, National Astronomical Observatory of Japan. The computing cluster of Shanghai VLBI correlator supported by the Special Fund for Astronomy from the Ministry of Finance in China is acknowledged. This work was supported by FAPESP (Fundacao de Amparo a Pesquisa do Estado de Sao Paulo) under grant 2021/01183-8.

APEX is a collaboration between the Max-Planck-Institut für Radioastronomie (Germany), ESO, and the Onsala Space Observatory (Sweden). The SMA is a joint project between the SAO and ASIAA and is funded by the Smithsonian Institution and the Academia Sinica. The JCMT is operated by the East Asian Observatory on behalf of the NAOJ, ASIAA, and KASI, as well as the Ministry of Finance of China, Chinese Academy of Sciences, and the National Key Research and Development Program (No. 2017YFA0402700) of China and Natural Science Foundation of China grant 11873028. Additional funding support for the JCMT is provided by the Science and Technologies Facility Council (UK) and participating universities in the UK and Canada. The LMT is a project operated by the Instituto Nacional de Astrófisica, Óptica, y Electrónica (Mexico) and the University of Massachusetts at Amherst (USA). The IRAM 30 m telescope on Pico Veleta, Spain is operated by IRAM and supported by CNRS (Centre National de la Recherche Scientifique, France), MPG (Max-Planck-Gesellschaft, Germany), and IGN (Instituto Geográfico Nacional, Spain). The SMT is operated by the Arizona Radio Observatory, a part of the Steward Observatory of the University of Arizona, with financial support of operations from the State of Arizona and financial support for instrumentation development from the NSF. Support for SPT participation in the EHT is provided by the National Science Foundation through award OPP-1852617 to the University of Chicago. Partial support is also provided by the Kavli Institute of Cosmological Physics at the University of Chicago. The SPT hydrogen maser was provided on loan from the GLT, courtesy of ASIAA.

This work used the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562, and CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442. XSEDE Stampede2 resource at TACC was allocated through TG-AST170024 and TG-AST080026N. XSEDE JetStream resource at PTI and TACC was allocated through AST170028. This research is part of the Frontera computing project at the Texas Advanced Computing Center through the Frontera Large-Scale Community Partnerships allocation AST20023. Frontera is made possible by National Science Foundation award OAC-1818253. This research was done using services provided by the OSG Consortium (Pordes et al. 2007; Sfiligoi et al. 2009), which is supported by the National Science Foundation award Nos. 2030508 and 1836650. This research was carried out using resources provided by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy Office of Science. Additional work used ABACUS2.0, which is part of the eScience center at Southern Denmark University. Simulations were also performed on the SuperMUC cluster at the LRZ in Garching, on the LOEWE cluster in CSC in Frankfurt, on the HazelHen cluster at the HLRS in Stuttgart, and on the Pi2.0 and Siyuan Mark-I at Shanghai Jiao Tong University. The computer resources of the Finnish IT Center for Science (CSC) and the Finnish Computing Competence Infrastructure (FCCI) project are acknowledged.This research was enabled in part by support provided by Compute Ontario (http://computeontario.ca), Calcul Quebec (http://www.calculquebec.ca), and Compute Canada (http://www.computecanada.ca).

The EHTC has received generous donations of FPGA chips from Xilinx Inc., under the Xilinx University Program. The EHTC has benefited from technology shared under open-source license by the Collaboration for Astronomy Signal Processing and Electronics Research (CASPER). The EHT project is grateful to T4Science and Microsemi for their assistance with hydrogen masers. This research has made use of NASA's Astrophysics Data System. We gratefully acknowledge the support provided by the extended staff of the ALMA, from the inception of the ALMA Phasing Project through the observational campaigns of 2017 and 2018. We would like to thank A. Deller and W. Brisken for EHT-specific support with the use of DiFX. We thank Martin Shepherd for the addition of extra features in the Difmap software that were used for the CLEAN imaging results presented in this paper. We acknowledge the significance that Maunakea, where the SMA and JCMT EHT stations are located, has for the indigenous Hawaiian people.

In all imaging ray-tracing codes radiative transfer equations are solved along null geodesics that terminate at a "camera" at some large distance from the supermassive BHs and where a polarization map at a chosen observing frequency is constructed. By integrating Stokes parameters over entire images made for different frequencies, one can also construct a polarized synchrotron spectral energy distribution of any model. However, instead of comparing model spectra from the discussed imaging codes, here we carry out an alternative comparison. Namely, we compare spectra produced by the ipole code to polarized spectra generated via Monte Carlo scheme radpol. In the Monte Carlo code the polarized radiative transfer integration scheme is conceptually distinct from all discussed imaging codes (for detailed description see Mościbrodzka 2020). Showing a convergence of two different approaches is an independent validation of emission produced by imaging codes. Hence, we compare spectra produced by ipole and radpol codes using the plasma model setup described in Section 3.3 (the test with ${{ \mathcal M }}_{\mathrm{low}}$). In Figure 9, we show radio−millimeter spectra of Stokes I luminosity, fractional linear polarization, and circular polarizations. The relative difference between luminosities is less than 10%, except for high-frequency emission. Both codes show consistent amplitude of fractional linear and circular polarizations and agree on handedness of circular polarization.

Zoom In Zoom Out Reset image size

The emission, absorption, and rotation coefficients of a fluid are well determined for a particular distribution of electron energies. However, the integrations involved are numerically expensive; since the coefficients must be calculated at every step when integrating the radiative transfer equations, fitting functions have been developed to approximate the coefficients quickly. A few sets of such fitting functions exist applicable to our regime: one outlined in Dexter (2016), the other in Pandya et al. (2016) and Marszewski et al. (2021).

The differences between these functions at various points within a representative set of input parameters are given in their respective papers, but we wish to provide some intuition concerning the differences these functions make in practice, and consequently whether fitting accuracy might be a driving factor in code differences.

All of the images in this comparison were created using the coefficient fits from Dexter (2016). Substituting the coefficient fits from Pandya et al. (2016) produces results more dissimilar than the disagreement between codes on three metrics: the net circular polarization at 0.13 points rather than 0.08 points, and average linear polarization fraction at 1.8 points rather than 0.7 points, and the $\left|{\beta }_{2}\right|$ coefficient at 0.016 rather than 0.0026. This is due to significant differences between the fits in computing the emission coefficient for circularly polarized light, j_V .

Unlike most emission coefficients, Faraday rotation coefficient ρ_V does not go to zero with low temperature—therefore, the low-temperature behavior of fitting functions is important. In particular, the expression from Dexter (2016) for ρ_V should not be used at low temperature Θ_e < 1, since it can produce a catastrophic cancellation not matching the desired limiting behavior of one of its quotients. Due to this instability, most codes either switch to the Shcherbakov fit at low temperature (e.g., grtrans) or use the Shcherbakov fit exclusively (e.g., ipole-IL). As illustrated in Table 5, these approaches produce nearly identical results, different by an NMSE less than 10⁻⁵ in the worst case.

Table 5. Comparison of Changes to the ipole-IL Result under Different Emission Coefficient Fitting Functions

ipole-IL Image Run...	NMSE I	NMSE Q	NMSE U	NMSE V	Δ∣m∣net (%)	Δ〈∣m∣〉 (%)	Δvnet (%)	${\rm{\Delta }}\left|{\beta }_{2}\right|$	Δ∠β2 (deg)
...with Pandya jS	4.7e-05	0.0055	0.0059	0.0087	0.0246	1.78	0.127	0.0168	0.353
...with Dexter ρV	1.4e-09	2.6e-06	2.1e-06	6.8e-06	−0.00746	0.0186	−0.0031	0.000131	0.00225
...with approximate Kn (x)	3.9e-08	5.6e-05	2.5e-05	0.00014	−0.000404	0.0239	−0.00524	0.000123	−0.00336

Note. Each row lists a change made to the default emissivity values and the resulting differences between the new image and the one used in the comparison. Substituting emissivities from Pandya et al. (2016) changes the result by more than the overall level of code agreement, whereas substituting the rho_V fit from Dexter (2016) changes almost nothing about the image. Approximating the Bessel functions K_n does not badly affect this image but can be a substantial source of error in images with higher Faraday rotation.

Download table as:  ASCIITypeset image

In addition, both the expressions from Shcherbakov (2008) and Dexter (2016) involve Bessel functions, which are tempting to approximate by assuming that emission is exclusive to the regime Θ_e > 1 in order to avoid unnecessary computation. However, this approximation produces clearly incorrect limiting behavior for ρ _V at low temperature. Thus, the Faraday rotation is misapplied, producing too little rotation in the EVPA, or in some cases, rotation in the wrong direction. In the sample image this is a minor effect due to an overall small Faraday rotation, but it can severely affect SANE disks seen from larger observer angles.

In addition to the NMSE, several other metrics could be used to gauge image dissimilarity between codes. Three additional metrics were evaluated in the context of this comparison: the normalized mean linear error (NMLE), structural dissimilarity (DSSIM), and inverse zero-normalized cross-correlation (DZNCC). These are defined as follows:

$$\begin{eqnarray}&&\mathrm{NMLE}(A,B)=\displaystyle \frac{{\sum }_{j}| {A}_{j}-{B}_{j}| }{{\sum }_{j}| {A}_{j}| }\end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{SSIM}(A,B) & = & \left(\displaystyle \frac{2{\mu }_{A}{\mu }_{B}}{{\mu }_{A}^{2}+{\mu }_{B}^{2}}\right)\\ & & \times \left(\displaystyle \frac{\tfrac{2}{N}{\sum }_{j}\left({A}_{j}-{\mu }_{A}\right)\left({B}_{j}-{\mu }_{B}\right)}{{\sigma }_{A}^{2}+{\sigma }_{B}^{2}}\right)\end{array}\end{eqnarray} \tag{ C2 }$$

$$\begin{eqnarray}&&\mathrm{DSSIM}(A,B)=\displaystyle \frac{1}{\left|\mathrm{SSIM}(A,B)\right|}-1\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&\mathrm{ZNCC}(A,B)=\displaystyle \frac{1}{N}\displaystyle \sum _{j}\displaystyle \frac{1}{{\sigma }_{A}{\sigma }_{B}}\left({A}_{j}-{\mu }_{A}\right)\left({B}_{j}-{\mu }_{B}\right)\end{eqnarray} \tag{ C4 }$$

$$\begin{eqnarray}&&\mathrm{DZNCC}(A,B)=\displaystyle \frac{1}{\left|\mathrm{ZNCC}(A,B)\right|}-1,\end{eqnarray} \tag{ C5 }$$

where μ_X is the average pixel value of an image and σ_X is the standard deviation of the pixel values.

In this comparison, most images were very similar—for this limited case, the various similarity metrics were found to correlate strongly in ordering of similarity between all images, and usually even in relative magnitude, as shown in Figure 10. The absolute values of the constraint metrics mean little without the context of Section 5.2, so any particular metric could fill the role of a similarity gauge to compare to image variation from other sources.

Zoom In Zoom Out Reset image size

While the GRMHD snapshot file used for the test in Section 3.3 reflects the simulations and imaging parameters used in practice by the EHTC in studies of M87*, there is always the chance that the snapshot itself is a particularly simple case, not reflective of average code differences in practice. Furthermore, in characterizing the impact of code differences on metric-based model comparisons, it would be useful to have an idea of what portion of code differences in metrics are due to systematic errors, versus stochastic products of limited accuracy parameters or sampling differences.

To measure the variation in results of this test over a typical variety of GRMHD states, two codes with substantially different algorithms, ipole-IL and grtrans, were compared across 2000 snapshots of a GRMHD simulation used in generating the EHT image libraries. This particular simulation represented a magnetically arrested disk (MAD) state about a BH of spin a_* = 0.9375, and the 2000 snapshots shown represent the entire quiescent portion of the simulation from 5000 r_g /c to 15,000 r_g /c. Details of the initial conditions, resolution, etc., are available in Wong et al. (2022).

Figure 11 compares the total unpolarized flux computed by ipole-IL and grtrans over the entire window. Figure 12 provides histograms of the differences in all image-integrated values over the window, along with Gaussian functions following their means and standard deviations. Table 6 lists the mean (i.e., average difference) and standard deviation (i.e., span of differences) between codes in each metric. The parameters 〈∣m∣〉 and ∠β₂ appear to be almost entirely stochastic. Flux, ∣m∣_net, and v_net are approximately a quarter systematic, and $\left|{\beta }_{2}\right|$ is half systematic (though the error itself is minuscule).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 6. Mean Values and Standard Deviations of the Distributions in Figure 12

Variable	μ	σ	μ/σ
$\tfrac{{\rm{\Delta }}\mathrm{Flux}}{\mathrm{Flux}}$ (%)	−0.197	0.841	0.234
Δ∣m∣net (%)	−0.0627	0.27	0.232
Δ〈∣m∣〉 (%)	−0.0182	0.206	0.0884
Δ∣v∣net (%)	0.00168	0.0116	0.145
${\rm{\Delta }}\left|{\beta }_{2}\right|$	−0.00178	0.00361	0.493
Δ∠β2 (deg)	−0.0757	1.31	0.0576

Note. The last column lists the proportion μ/σ, which can be taken as an estimate of the relevance of systematic versus stochastic errors in describing differences between these codes.

Download table as:  ASCIITypeset image