The Role of Adaptive Ray Tracing in Analyzing Black Hole Structure

The simplest interpolation scheme to generate a full image from a discrete set of rays is a nearest-neighbors approach. Given a ray sampled at x ₁, x ₂,..., x _n , the nearest-neighbors intensity distribution simply finds the closest sampled ray at every location: $\widehat{I}({\boldsymbol{x}})\equiv I\left(\arg \ {\min }_{{{\boldsymbol{x}}}_{i}}\left|{{\boldsymbol{x}}}_{i}-{\boldsymbol{x}}\right|\right)$. A second interpolation scheme is a linear/bilinear approach, which is defined so that $\widehat{I}({\boldsymbol{x}})$ is an average of its nearest surrounding rays that is weighted by distance. For a smooth intensity distribution, the errors in these two schemes are given by the first and second image spatial derivatives, respectively.

Specifically, consider a point x lying equidistant from four rays spaced evenly around x . Then, to leading order, the interpolation residuals will be (see Appendices A.1 and A.2 for full derivation)

$$\begin{eqnarray}\left|\widehat{I}({\boldsymbol{x}})-I({\boldsymbol{x}})\right|\sim \left\{\begin{array}{ll}\left|{\rm{\nabla }}\widehat{I}({\boldsymbol{x}})\right|{\rm{\Delta }}x, & \mathrm{Nearest}\\ \tfrac{1}{4}| {{\rm{\nabla }}}^{2}\widehat{I}({\boldsymbol{x}})| {\rm{\Delta }}{x}^{2}, & \mathrm{Linear}.\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where Δx denotes the distance between x and each of its four surrounding rays.

In Figure 1, we plot the gradients and Laplacians for ray-traced images of both a semianalytic model and a model from the EHT GRMHD library (Event Horizon Telescope Collaboration et al. 2019e). In particular, the GRMHD model is a Magnetically Arrested Disk (MAD; Yuan & Narayan 2014) of M87*, with dimensionless spin a_* = +0.94, inclination angle i = 17°, and a field of view of 160 μas (corresponding to ∼ 44.17M/D). We also use the electron heating parameter R_high = 20 (for details, see Mościbrodzka et al. 2016; Event Horizon Telescope Collaboration et al. 2019e).

Zoom In Zoom Out Reset image size

The semianalytic model mirrors Test 1 of Gold et al. (2020) and models a spherically symmetric fluid distribution around a black hole located at a distance D = 7.78 kpc and with mass M = 4 × 10⁶ M_⊙. Additional parameters for this semianalytic model include spin a_* = 0.9, inclination angle i = 60°, and a field of view of ∼ 152.289 μas (30M/D). For both models, we use a camera distance of r_cam = 10⁶ M to mitigate small errors that may arise from the use of a pinhole camera.

In all cases, the gradient and Laplacian of the intensity increase sharply near the photon ring because of strong gravitational lensing. Indeed, Psaltis et al. (2015) proposed using sharp image gradients as a way to localize the photon ring; the Laplacian would also be an effective choice.

The highly localized regions of Figure 1 with large derivatives suggest that judicious sampling can be used to significantly improve image generation speed and quality relative to a uniform grid. We will now describe an efficient recursive sampling procedure that is defined by a pair of error tolerances: ${{ \mathcal R }}_{\mathrm{abs}}$ and ${{ \mathcal R }}_{\mathrm{rel}}$.

To begin, we ray-trace an image on a grid at a coarse initial resolution, n_0x × n_0y , where n_0x and n_0y are any two integers. Next, we selectively populate an image with finer resolution by a factor of two: (2n_0x − 1) × (2n_0y − 1). For each point, we either ray-trace or interpolate based on the expected relative and absolute interpolation errors:

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{rel}}({\boldsymbol{x}})=\left|\displaystyle \frac{I({\boldsymbol{x}})-\widehat{I}({\boldsymbol{x}})}{I({\boldsymbol{x}})}\right|,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{abs}}({\boldsymbol{x}})=\left|\displaystyle \frac{I({\boldsymbol{x}})-\widehat{I}({\boldsymbol{x}})}{\overline{I}}\right|.\end{eqnarray} \tag{ 3 }$$

Both _abs and _rel are dimensionless, but they differ in their normalization: _rel uses the local image intensity I( x ), while _abs uses the image-averaged intensity $\overline{I}$. If ${\epsilon }_{\mathrm{abs}}\gt {{ \mathcal R }}_{\mathrm{abs}}$ and ${\epsilon }_{\mathrm{rel}}\gt {{ \mathcal R }}_{\mathrm{rel}}$, then I( x ) is computed by ray tracing. Otherwise, I( x ) is computed by interpolation. Hence, rays are computed only where interpolation residuals are expected to be large. This procedure can then be repeated arbitrarily, giving an effective final image resolution that is n_x × n_y , with

$$\begin{eqnarray}\begin{array}{rcl}{n}_{x} & = & {2}^{N}({n}_{0x}-1)+1\\ {n}_{y} & = & {2}^{N}({n}_{0y}-1)+1.\end{array}\end{eqnarray} \tag{ 4 }$$

Our recursive approach relies on computing estimates ${\widehat{\epsilon }}_{\mathrm{rel}}({\boldsymbol{x}})$ and ${\widehat{\epsilon }}_{\mathrm{rel}}({\boldsymbol{x}})$ for the interpolation residuals, and these estimates change for different configurations of pixels. As shown in Figure 2, each point x lying on the (n + 1)th grid will fall into one of four categories:

• 1.  
  Category 1 (Black): x had its intensity computed at the nth refinement level; there is no interpolation residual.
• 2.  
  Category 2 (Blue): x lies in between two points located above and below, each of whose intensities were computed at the nth refinement level.
• 3.  
  Category 3 (Red): x lies in between two points located to left and right, each of whose intensities were computed at the nth refinement level.
• 4.  
  Category 4 (Green): x lies in between four equidistant corner points, each of whose intensities were computed at the nth refinement level.

Zoom In Zoom Out Reset image size

For cases 2–4, we use finite differences to estimate derivatives and then take the appropriate leading term in the Taylor series approximation to evaluate interpolation uncertainties (see, e.g., Appendix A of Pedrola 2015). For example, consider estimating the intensity at the central point of Figure 2, whose location we will call x . Defining ${\widehat{I}}_{i}\equiv \widehat{I}({{\boldsymbol{x}}}_{i})$ and labeling the points as in Figure 3, we find

$$\begin{eqnarray}{\widehat{\epsilon }}_{\mathrm{abs}}({\boldsymbol{x}})=\left\{\begin{array}{ll}\left|\displaystyle \frac{1}{{\overline{I}}_{\mathrm{int}}}\sqrt{{\left(\displaystyle \frac{{\widehat{I}}_{4}-{\widehat{I}}_{1}}{2}\right)}^{2}+{\left(\displaystyle \frac{{\widehat{I}}_{3}-{\widehat{I}}_{2}}{2}\right)}^{2}}\right|, & \mathrm{Nearest}\\ \left|\displaystyle \frac{-({\widehat{I}}_{1}+{\widehat{I}}_{2}+{\widehat{I}}_{3}+{\widehat{I}}_{4})+({\widehat{I}}_{5}+{\widehat{I}}_{6}+{\widehat{I}}_{7}+{\widehat{I}}_{8})}{16{\overline{I}}_{\mathrm{int}}}\right|, & \mathrm{Linear},\end{array}\right.\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{\widehat{\epsilon }}_{\mathrm{rel}}({\boldsymbol{x}})=\left\{\begin{array}{ll}\left|\displaystyle \frac{1}{{\widehat{I}}_{1}}\sqrt{{\left(\displaystyle \frac{{\widehat{I}}_{4}-{\widehat{I}}_{1}}{2}\right)}^{2}+{\left(\displaystyle \frac{{\widehat{I}}_{3}-{\widehat{I}}_{2}}{2}\right)}^{2}}\right|, & \mathrm{Nearest}\\ \left|\displaystyle \frac{-({\widehat{I}}_{1}+{\widehat{I}}_{2}+{\widehat{I}}_{3}+{\widehat{I}}_{4})+({\widehat{I}}_{5}+{\widehat{I}}_{6}+{\widehat{I}}_{7}+{\widehat{I}}_{8})}{4({\widehat{I}}_{1}+{\widehat{I}}_{2}+{\widehat{I}}_{3}+{\widehat{I}}_{4})}\right|, & \mathrm{Linear}.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

Here, the nearest neighbor is chosen (arbitrarily) to be ${\widehat{I}}_{1}$. Additionally, ${\overline{I}}_{\mathrm{int}}$ is defined as the interpolated average intensity after the first pass of ray tracing (with resolution n_0x × n_0y ). For a detailed derivation of Equations (5) and (6), see Appendices A.1 and A.2.

Zoom In Zoom Out Reset image size

While we have derived estimates for ${\widehat{\epsilon }}_{\mathrm{abs}}({\boldsymbol{x}})$ and ${\widehat{\epsilon }}_{\mathrm{rel}}({\boldsymbol{x}})$ using the simplified approximation of a smooth image with interpolation residuals dominated by low-order derivatives, these estimates are also useful for the more general case of images with small-scale structure and sharp gradients. For example, if the small-scale image structure is stochastic with a power-law spectrum, then we can compare the ensemble-average properties of the exact (Equation (2)) and estimated (Equation (5)) interpolation residuals. For 1D interpolation, we find

$$\begin{eqnarray}\displaystyle \frac{{\left\langle {\widehat{\epsilon }}_{\mathrm{abs}}^{2}\right\rangle }^{\tfrac{1}{2}}}{{\left\langle {\epsilon }_{\mathrm{abs}}^{2}\right\rangle }^{\tfrac{1}{2}}}=\left\{\begin{array}{ll}{2}^{\beta /2-2}, & \mathrm{Nearest}\\ {2}^{\tfrac{\beta }{2}-4}\sqrt{\displaystyle \frac{1+{2}^{\beta -1}-{3}^{\beta -2}}{1-{2}^{\beta -4}}}, & \mathrm{Linear},\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where β is the power-law exponent. GRMHD simulations often have β ∼ 2.5, giving ratios ${\left\langle {\widehat{\epsilon }}_{\mathrm{abs}}^{2}\right\rangle }^{\tfrac{1}{2}}/{\left\langle {\epsilon }_{\mathrm{abs}}^{2}\right\rangle }^{\tfrac{1}{2}}$ of 0.59 and 0.27 for the 1D nearest and linear interpolation, respectively. Thus, even in this generalized case of stochastic image fluctuations with a power-law spectrum (for which a local Taylor series expansion is poor), our approximate estimates for interpolation error will still provide useful refinement criteria. For additional discussion of the interpolation errors for images with power-law spectra, see Appendix B.

Previous uses of adaptive ray tracing have used more complicated refinement criteria than the ones we have selected. Parkin (2011) also implemented a second-order criterion (the Löhner 1987 criterion), which uses multiple cross derivatives. Similarly, VRT ² employs a bicubic interpolator. However, we expect the benefit of these higher-order schemes to be marginal in the case of an image with a turbulent power spectrum (see Appendix B). For the remainder of this paper, we use the linear interpolation scheme.

We implement the recursive ray-tracing scheme described above into the polarized general relativistic radiative transfer (GRRT) code ipole ⁶ (Mościbrodzka & Gammie 2017). This implementation allows us to assess the performance of our algorithm and subsequently generate extremely high-resolution images of black holes.

In addition to the MAD model and spherical semianalytic model used in Figure 1, we also generate images using a Standard and Normal Evolution (SANE) model from the EHT GRMHD library (Yuan & Narayan 2014; Event Horizon Telescope Collaboration et al. 2019e), as well as a semianalytic model of a geometrically thin, rotating disk (Test 5 of Gold et al. 2020).

Before generating images, however, we must first determine suitable error tolerances ${{ \mathcal R }}_{\mathrm{abs}}$ and ${{ \mathcal R }}_{\mathrm{rel}}$. To evaluate the effect of different tolerances on the resultant image, we use the following three metrics (the second and third of which are also used in Gold et al. 2020):

$$\begin{eqnarray}&&\mathrm{Interpolation}\,\mathrm{Fraction}\equiv \frac{\#\mathrm{pixels}\,\mathrm{interpolated}}{{n}_{x}\times {n}_{y}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\rm{Flux}}\,{\rm{Error}}\equiv \left|\displaystyle \frac{F-\hat{F}}{F}\right|,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\rm{Mean}}\,{\rm{Squared}}\,{\rm{Error}}\equiv \displaystyle \frac{{\sum }_{{\rm{pixels}}i}| I({{\boldsymbol{x}}}_{i})-\hat{I}({{\boldsymbol{x}}}_{i}){| }^{2}}{{\sum }_{i}| I({\boldsymbol{x}}){| }^{2}}.\end{eqnarray} \tag{ 10 }$$

Here, F ≡ ∫d² x I( x ) is the image flux associated with the model intensity distribution, and $\widehat{F}\equiv \int {d}^{2}{\boldsymbol{x}}\,\widehat{I}({\boldsymbol{x}})$ is the image flux associated with the estimated intensity distribution. We note that while a high flux error necessarily implies a high MSE, we include the former metric because it has a clear physical interpretation and bears more relevance in potential applications to light curves.

Our goal is to maximize the interpolation fraction (IF) while minimizing the flux error (FE) and MSE, as this achieves both efficiency and accuracy. To this end, we run our adaptive scheme on a single GRMHD snapshot (the same snapshot/parameters used in Figure 1) with a wide variety of error tolerances, and we evaluate the IF/FE/MSE for each. The contours in Figure 4 show how the error metrics each depend on ${{ \mathcal R }}_{\mathrm{rel}}$ and ${{ \mathcal R }}_{\mathrm{abs}}$. For this comparison, we take n_0x = n_0y = 65 pixels and n_x = n_y = 1025 pixels, and we compute the model intensity distribution I( x ) by ray-tracing each pixel on a 1025 × 1025 grid.

Zoom In Zoom Out Reset image size

As expected, as ${{ \mathcal R }}_{\mathrm{rel}}$ and ${{ \mathcal R }}_{\mathrm{abs}}$ increase, the adaptively sampled image $\widehat{I}({\boldsymbol{x}})$ increasingly deviates from the fully sampled I( x ). As ${{ \mathcal R }}_{\mathrm{rel}}\to \infty$ and ${{ \mathcal R }}_{\mathrm{abs}}\to \infty$, all adaptive refinement will cease and the image will be given by the initial sampling. For the initial 65 × 65 sampling of this snapshot, MSE ≈ 0.22 ≈ 10^−0.65 and FE ≈ 0.029 ≈ 10^−1.54. These serve as reference values to compare against iterations of the scheme with more stringent tolerances.

The contours in Figure 4 can be used for practical image generation purposes—for the desired accuracy, one can select ${{ \mathcal R }}_{\mathrm{abs}}$ and ${{ \mathcal R }}_{\mathrm{rel}}$ such that the interpolation fraction and hence the efficiency are maximized. For our subsequent analysis of GRMHD images, we choose ${{ \mathcal R }}_{\mathrm{abs}}=0.025$ and ${{ \mathcal R }}_{\mathrm{rel}}=0.001$ (the red star in Figure 4), as this set of tolerances allows for >90% interpolation while constraining the MSE and FE to <0.1%.

For subsequent analysis of the semianalytic models, we choose ${{ \mathcal R }}_{\mathrm{abs}}={{ \mathcal R }}_{\mathrm{rel}}=0.001$, as we find that the semianalytic models are slightly more sensitive to the absolute tolerance.

Using the tolerances selected described above, we adaptively ray-trace in ipole at resolutions of 1025 × 1025 and 4097 × 4097 pixels. Images of the four models for the latter resolution are shown in Figure 5, along with their respective "ray-density plots," illustrating where the program is concentrating rays in each image. The ray-density plots are generally consistent with Figure 1—the majority of ray tracing takes place in an annulus about the photon ring, while the background and shadow are predominantly interpolated. As with Figure 1, the annuli for the semianalytic models are significantly thinner and more pronounced.

Zoom In Zoom Out Reset image size

The corresponding residual images (fully ray-traced images subtracted from adaptively ray-traced images and normalized by the image-averaged intensity: $\left|\tfrac{\widehat{I}({\boldsymbol{x}})-I({\boldsymbol{x}})}{\overline{I}}\right|$) are shown in Figure 6. In all images, the residual amplitudes are relatively constant, which reflects the refinement goal of _abs. In portions of the image that correspond to rays that have been traced, the residuals are zero. In particular, all images contain a base grid of 65 × 65 evenly spaced rays, giving evenly spaced nulls in the residuals that align with this grid. This feature is most evident in the semianalytic models, where the smooth intensity distribution allows significant interpolation.

Zoom In Zoom Out Reset image size

The resultant error statistics (Equations (8)–(10)) for these adaptively ray-traced images are presented in Table 1. The 1025 × 1025 images match the predictions from the contours in Figure 4: both the MSE and FE are lower than 0.1%, while the interpolation fraction exceeds 90%. In Table 1, we also list the total number of rays traced in each image: (1 − IF) × n_x × n_y , which is roughly proportional to computational expense.

While these statistics were generated for a specific choice of image parameters, we obtain similar results for other models. We ray-trace with the same tolerances on all combinations of MAD/SANE, a_* = −0.94, 0, + 0.94, and R_high = 20, 40, 80. We find that among these images, IF ranges from 0.877 to 0.929, FE ranges from 3.02 × 10⁻⁵ to 8.59 × 10⁻⁴, and MSE ranges from 2.82 × 10⁻⁴ to 5.55 × 10⁻⁴.

We may further vary the magnetization (σ) cut (see, e.g., Chael et al. 2019), which is a quantity designed to restrict emission to regions where the fluid density has not been physically invalidated. Regenerating the MAD image from Figure 1 with σ = 1, 5, 10, and 50, we find that IF ranges from 0.917 to 0.931, FE ranges from 2.50 × 10⁻⁴ to 7.50 × 10⁻⁴, and MSE ranges from 2.79 × 10⁻⁴ to 4.69 × 10⁻⁴. Thus, for these GRMHD simulations, the accuracy and efficacy of the adaptive refinement criteria are relatively insensitive to the choice of accretion and radiation model parameters.

The semianalytic models have lower errors than the GRMHD models and higher interpolation fractions. This is likely because the smooth underlying structure in semianalytic models is better suited to interpolation.

While IF is a useful metric to quantify sampling efficiency, it does not directly correspond to the reduction in image generation time. Namely, rays near the photon ring are more expensive to trace, as their geodesics have longer paths through the emitting material and thus require more calculations to perform the radiative transfer. Nevertheless, even though our adaptive scheme predominantly samples in this computationally expensive region, we find that the reduction in image generation time is similar to the interpolation fraction (see Figure 7). In general, the relationship between interpolation fraction and computational expense will depend on details of the underlying model and of the GRRT implementation.

Zoom In Zoom Out Reset image size

We note additionally from Table 1 that as resolution increases, IF increases as well. This behavior is expected because the error estimates decrease at smaller separations, allowing more pixels to satisfy the error tolerances required for interpolation. Indeed, Figure 5 shows that the density of rays in the diffuse parts of the image quickly saturates—upon increasing the spatial resolution from 1025 × 1025 to 4097 × 4097, the sampled ray density only increases in the photon ring region.

Unlike the IF, however, the MSE and FE do not depend strongly on resolution. These quantities instead depend on the tolerances ${{ \mathcal R }}_{\mathrm{rel}}$ and ${{ \mathcal R }}_{\mathrm{abs}}$. Thus, even though more pixels are interpolated at a resolution of 4097 × 4097 compared to a resolution of 1025 × 1025, the output images have comparable accuracy by these two metrics.

VLBI directly measures interferometric visibilities, which correspond to complex Fourier components of the sky image (Thompson et al. 2017). Hence, visibility-domain tests are appropriate to assess suitability for direct comparisons with observables. In this section, we analyze the visibility spectra of our adaptively ray-traced images and show that they exhibit the universal properties expected for black hole images.

4.3.1. Expected Visibility Signatures

4.3.2. Visibility Spectra of High-resolution Images

We now explore the expected long-baseline visibility signatures using our high-resolution images computed with adaptive ray tracing. While short interferometric baselines will have a complex visibility structure that depends on the overall image morphology, the visibility signatures from the photon ring and from turbulence will emerge on baselines that heavily resolve the image. Specifically, to accurately estimate the visibility on a baseline with dimensionless length u requires angular resolution ${\rm{\Delta }}\theta \sim 1/u$ (Thompson et al. 2017). For our images of both the GRMHD and semianalytic models, the visibility spectra can thus be computed to baselines of ${u}_{\max }\gtrsim 5000\,{\rm{G}}\lambda$. This value is approximately 600 times larger than the longest current EHT baselines (Event Horizon Telescope Collaboration et al. 2019b). However, in the following analysis, to minimize errors from finite pixel size, we only analyze baselines shorter than 1000 Gλ.

Figure 8 shows visibility amplitudes for the adaptively ray-traced MAD model, the spherically symmetric semianalytic model, and the disk semianalytic model. Because the semianalytic images have significant emission extending beyond our specified FOV, the sharp artificial cutoffs at the edges of the FOV produce spurious high-frequency visibility power. To suppress this artificial power, we double the FOV and pixel number (keeping the image resolution fixed) and then apply a Gaussian taper with FWHM of FOV/4 before computing visibilities. Figure 8 also shows the visibility errors resulting from the adaptive ray tracing, demonstrating that these errors are a small fraction of the visibility amplitudes on all baselines. Specifically, these errors correspond to the Fourier amplitudes of the residual images defined in Section 4.2.

Zoom In Zoom Out Reset image size

The MAD visibility spectrum (top row) possesses many kinks and does not display a clear-cut pattern. Although a characteristic periodicity may be evident, small-scale image power from turbulence exceeds that of the lensed emission. This turbulent power gives a visibility "noise" that decays on long baselines approximately as V(u) ∼ u^−p with p ≈ 1.2.

The spherical semianalytic model (middle row), on the other hand, displays a smooth, turbulence-free spectrum. Just as distinct subrings are not visible in the images of this model due to the spherical symmetry of the fluid distribution (e.g., Narayan et al. 2019), distinct subrings are not visible in the visibilities of this model.

In contrast, the disk semianalytic model (bottom row) does show clear signs of distinct photon subrings in both the images and visibilities. The spectrum falls steeply around ∼ 500 Gλ before flattening again shortly thereafter, corresponding to the transition between the n = 1 and n = 2 subrings.

Resolving the photon subrings in the image domain is difficult in practice due to the exponential radial demagnification of each successive subimage (Darwin 1959; Luminet 1979; Ohanian 1987). In particular, per Johnson et al. (2020), the width of each subring on the screen scales as w_n ∼ w₀ e^{−γ n} , where the Lyapunov exponent γ is defined as

$$\begin{eqnarray}&&\gamma \equiv \sqrt{\displaystyle \frac{2{ \mathcal R }^{\prime\prime} ({r}_{{\rm{c}}})}{-{u}_{-}{a}^{2}{E}^{2}}}\,K\left(\displaystyle \frac{{u}_{+}}{{u}_{-}}\right).\end{eqnarray} \tag{ 11 }$$

Here, K is the complete elliptic integral of the first kind (with squared modulus u₊/u₋), ${ \mathcal R }(r)$ is the radial effective potential for null geodesics in Kerr, E is the conserved energy, r_c is the corresponding radius of the spherical orbit, and u_± denote roots of the angular effective potential. Furthermore, the radial curve ρ_n of each successive subring exponentially approaches the critical curve ρ_c : δ ρ_n ≈ ρ_c e^{−γ n} . Hence, the image resolution required to resolve each successive subring thus grows exponentially for a fixed FOV.

To explore the properties of photon subrings explicitly, we combine our adaptive scheme with a subring decomposition code. The decomposition code can generate images corresponding to the nth subring by only including emission from the appropriate segment of the full geodesic. Notice that this definition of subrings differentiates between photons that may follow the same geodesic: the geodesic for a photon that makes n half-orbits will overlap with the geodesic for some other photon that makes n − 1 half-orbits. Thus, a pixel that is illuminated in the n = 2 image may also be illuminated in the n = 1 one, and the total pixel brightness will have contributions from the second and first subrings, respectively (see, e.g., Figure 3 of Johnson et al. 2020).

Figure 9 shows the visibility spectra of the disk semianalytic model decomposed into contributions from individual subrings. These reveal a new feature in the subring visibilities: the different subring widths between the top and bottom of the image lead to intermittent beating along the v axis between the n and n + 1 subrings. When the former is resolved, the beating vanishes, and the visibilities smoothly decay until reaching the level of power from the next subring. This phenomenon does not appear on the u axis, as the ring does not have significant thickness asymmetry on the horizontal axis.

Zoom In Zoom Out Reset image size

To illustrate the presence of subrings at high resolutions, we again use the adaptive subring decomposition code to generate 32769 × 32769 images of the n = 1, 2 and 3 subrings of the MAD model, now with a FOV of 80 μas. We show these images along with a lower-resolution image of n = 0 in Figure 10. Figure 11 shows stacked cross sections of the brightness profiles and zoom in on four points of the image (bottom, top, left, right). Individual subrings become thinner, with approximately constant peak brightness. Thus, the peak brightness of the sum of the first n subimages is approximately proportional to (n + 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can use adaptive ray tracing to generate high-resolution movies, which can be used to study how various averaging techniques reduce the turbulent noise present in GRMHD visibility spectra. Reducing turbulent noise is necessary to reveal universal signatures of the lensed emission surrounding the black hole, which reflect the purely geometrical properties of the spacetime. Additionally, by quantifying the amount of turbulence present in these simulations, we can obtain a better understanding of the accretion dynamics surrounding black holes and the timescales over which they vary.

In Figure 12, we show the visibility spectra for the MAD model with t_avg = 0M (i.e., a single snapshot), as well as t_avg = 100M and t_avg = 500M. The images were generated with a resolution of 1025 × 1025 and are thus capable of resolving the n = 0 and n = 1 subrings, whose individual visibility spectra are shown in green and blue, respectively.

Zoom In Zoom Out Reset image size

We see that by an averaging scale of t_avg = 100M, the visibility spectrum of the full image has mostly converged to that of the n = 1 subring for u ≲ 100 Gλ. However, the n = 0 spectrum still neighbors that of n = 1, indicating that turbulent noise is still significant. By 500M, the n = 0 spectrum falls below the n = 1 spectrum in the region where we expect n = 1 to dominate the time-averaged spectrum.

Suppose we wish to interpolate the intensity at x directly from its nearest neighbor at x ₁. Because I( x ) and $\widehat{I}({\boldsymbol{x}})$ are smooth in between x and x ₁, then we may apply Taylor's theorem (or in this case, just the mean value theorem) to see that

$$\begin{eqnarray}&&\widehat{I}-I={\rm{\nabla }}\widehat{I}({\boldsymbol{p}})\cdot ({\boldsymbol{x}}-{{\boldsymbol{x}}}_{1}),\end{eqnarray} \tag{ A1 }$$

for a point p lying on the line segment connecting x to x ₁. To leading order in ∣ x − x ₁∣, we may replace p with x , giving

$$\begin{eqnarray}&&\widehat{I}-I\approx {\rm{\nabla }}\widehat{I}({\boldsymbol{x}})\cdot ({\boldsymbol{x}}-{{\boldsymbol{x}}}_{1}).\end{eqnarray} \tag{ A2 }$$

Let us now restrict our attention to the rectangular gridding scheme presented in Figure 2. On this grid, Equation (A2) is ambiguous, as x will be adjacent to multiple equidistant pixels, rendering the location of x ₁ ill defined. Regardless of where we choose to set x ₁, however, the interpolation residual will be extremized by the quantity $| {\rm{\nabla }}\widehat{I}({\boldsymbol{x}})| | {\boldsymbol{x}}-{{\boldsymbol{x}}}_{1}|$. So, defining Δx ≡ ∣ x − x ₁∣, we take

$$\begin{eqnarray}&&| I({\boldsymbol{x}})-\widehat{I}({\boldsymbol{x}})| \sim | {\rm{\nabla }}\widehat{I}({\boldsymbol{x}})| {\rm{\Delta }}x,\end{eqnarray} \tag{ A3 }$$

which is the first half of Equation (1). This expression is now symmetric—it does not depend on which of the equidistant pixels we define to be x 's nearest neighbor.

Because we only have access to I( x ) at discretely sampled rays, we approximate the gradient using finite differences, which requires an examination of each of the four categories of pixel locations listed in Section 3.2. Let us suppose that x falls into Category 4, as Categories 1–3 are just simplifications thereof. The arrangement of pixels for Category 4 is shown in Figure 3.

Because ∣∇I∣ is rotationally invariant, we may evaluate the partial derivatives along rotated axes $x^{\prime}$ and $y^{\prime}$ that are aligned with the corner pixels. Then adopting the same notation as Figure 3, the central difference approximations to the derivatives are (e.g., Appendix A of Pedrola 2015)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \widehat{I}}{\partial x^{\prime} }{\left|\right.}_{{\boldsymbol{x}}} & \approx & \displaystyle \frac{{\widehat{I}}_{4}-{\widehat{I}}_{1}}{2{\rm{\Delta }}x},\\ \displaystyle \frac{\partial \widehat{I}}{\partial y^{\prime} }{\left|\right.}_{{\boldsymbol{x}}} & \approx & \displaystyle \frac{{\widehat{I}}_{3}-{\widehat{I}}_{2}}{2{\rm{\Delta }}x}.\end{array}\end{eqnarray} \tag{ A4 }$$

Because $\widehat{I}({\boldsymbol{x}})\equiv {\widehat{I}}_{1}$, the nearest-neighbor error estimates become

$$\begin{eqnarray}\begin{array}{rcl}{\epsilon }_{\mathrm{abs}}^{\mathrm{NN}}({\boldsymbol{x}}) & \approx & \left|\displaystyle \frac{1}{{\overline{I}}_{\mathrm{int}}}\sqrt{{\left(\displaystyle \frac{{\widehat{I}}_{4}-{\widehat{I}}_{1}}{2}\right)}^{2}+{\left(\displaystyle \frac{{\widehat{I}}_{3}-{\widehat{I}}_{2}}{2}\right)}^{2}}\right|,\\ {\epsilon }_{\mathrm{rel}}^{\mathrm{NN}}({\boldsymbol{x}}) & \approx & \left|\displaystyle \frac{1}{{\widehat{I}}_{1}}\sqrt{{\left(\displaystyle \frac{{\widehat{I}}_{4}-{\widehat{I}}_{1}}{2}\right)}^{2}+{\left(\displaystyle \frac{{\widehat{I}}_{3}-{\widehat{I}}_{2}}{2}\right)}^{2}}\right|,\end{array}\end{eqnarray} \tag{ A5 }$$

which are the first half of Equations (5) and (6).

For linear interpolation on an arbitrary 2D grid, the error estimate does not reduce to an equation as simple as A2. However, by restricting our attention again to the specific gridding scheme in Figure 2, we are able to derive a straightforward error estimate as follows.

Suppose that x falls into Category 4, as once again, Categories 1–3 will just be simplifications thereof. Adopting the same labels as Figure 3, the linear interpolation residual is explicitly given by

$$\begin{eqnarray}&&\begin{array}{l}\widehat{I}-I=\displaystyle \frac{1}{4}\left(\displaystyle \sum _{i=1}^{4}\widehat{I}({{\boldsymbol{x}}}_{i})\right)-I({\boldsymbol{x}})\\ =\displaystyle \frac{\left(\tfrac{1}{2}[\widehat{I}({{\boldsymbol{x}}}_{1})+\widehat{I}({{\boldsymbol{x}}}_{4})]-I({\boldsymbol{x}})\right)+\left(\tfrac{1}{2}[\widehat{I}({{\boldsymbol{x}}}_{2})+\widehat{I}({{\boldsymbol{x}}}_{3})]-I({\boldsymbol{x}})\right)}{2}.\end{array}\,\end{eqnarray} \tag{ A6 }$$

We recognize this quantity as the average of two 1D linear interpolation residuals along the $x^{\prime}$ and $y^{\prime}$ axes. In 1D, linear interpolation residuals scale with the second spatial derivative (see, e.g., Theorem 4.3 of Epperson 2013), with the leading order expansion $| {\widehat{I}}_{1{\rm{D}}}-I| \approx \tfrac{1}{2}| \widehat{I}^{\prime\prime} ({\boldsymbol{x}})| {\rm{\Delta }}{x}^{2}$. Plugging this into the numerator of (A6), the 2D residual becomes

$$\begin{eqnarray}\begin{array}{rcl}\widehat{I}-I & = & \displaystyle \frac{\tfrac{1}{2}{\widehat{I}}_{x^{\prime} x^{\prime} }({\boldsymbol{x}}){\rm{\Delta }}{x}^{2}+\tfrac{1}{2}{\widehat{I}}_{y^{\prime} y^{\prime} }({\boldsymbol{x}}){\rm{\Delta }}{x}^{2}}{2}\\ & = & \displaystyle \frac{1}{4}{{\rm{\nabla }}}^{2}\widehat{I}({\boldsymbol{x}}){\rm{\Delta }}{x}^{2},\end{array}\end{eqnarray} \tag{ A7 }$$

where in the last step, we used the rotational invariance of the Laplacian. This is the second half of Equation (1).

The final task is now to approximate the Laplacian with finite differences. In doing so, we must be careful not to break the symmetry of the error approximation, and we must rely only on the pixels in Figure 3. To this end, we approximate the second derivatives by averaging the second-order central differences on either side of x . This gives

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\partial }^{2}\widehat{I}}{\partial x{{\prime} }^{2}} & \approx & \displaystyle \frac{{\widehat{I}}_{8}-{\widehat{I}}_{4}-{\widehat{I}}_{1}+{\widehat{I}}_{5}}{4{\rm{\Delta }}{x}^{2}},\\ \displaystyle \frac{{\partial }^{2}\widehat{I}}{\partial y{{\prime} }^{2}} & \approx & \displaystyle \frac{{\widehat{I}}_{7}-{\widehat{I}}_{3}-{\widehat{I}}_{2}+{\widehat{I}}_{6}}{4{\rm{\Delta }}{x}^{2}}.\end{array}\end{eqnarray} \tag{ A8 }$$

And because $\widehat{I}({\boldsymbol{x}})\equiv \tfrac{1}{4}({\widehat{I}}_{1}+{\widehat{I}}_{2}+{\widehat{I}}_{3}+{\widehat{I}}_{4})$, the error estimates become

$$\begin{eqnarray}\begin{array}{rcl}{\epsilon }_{\mathrm{abs}}^{\mathrm{lin}}({\boldsymbol{x}}) & \approx & \left|\displaystyle \frac{-({\widehat{I}}_{1}+{\widehat{I}}_{2}+{\widehat{I}}_{3}+{\widehat{I}}_{4})+({\widehat{I}}_{5}+{\widehat{I}}_{6}+{\widehat{I}}_{7}+{\widehat{I}}_{8})}{16{\overline{I}}_{\mathrm{int}}}\right|,\\ {\epsilon }_{\mathrm{rel}}^{\mathrm{lin}}({\boldsymbol{x}}) & \approx & \left|\displaystyle \frac{-({\widehat{I}}_{1}+{\widehat{I}}_{2}+{\widehat{I}}_{3}+{\widehat{I}}_{4})+({\widehat{I}}_{5}+{\widehat{I}}_{6}+{\widehat{I}}_{7}+{\widehat{I}}_{8})}{4({\widehat{I}}_{1}+{\widehat{I}}_{2}+{\widehat{I}}_{3}+{\widehat{I}}_{4})}\right|,\end{array}\end{eqnarray} \tag{ A9 }$$

which is the second half of Equations (5) and (6).

Appropriate modifications to the scheme are made when the pixels are close to the edge of the image. In this region, we may not have access to ${\widehat{I}}_{1}-{\widehat{I}}_{8}$.