Hybrid Very Long Baseline Interferometry Imaging and Modeling with themis

Generating images from radio wavelength very long baseline interferometric (VLBI) experiments presents a substantial computational and mathematical challenge. Generally, the quantity that is directly observed is related to the spatial Fourier transform of the image, i.e., the complex visibilities, V(u, v), defined in the spatial frequency plane (u, v). In principle, if V were measured at all spatial frequencies (i.e., all points in the u–v plane), and if it were well calibrated (i.e., the full complex V were known), then image reconstruction would reduce to an inverse Fourier transform. In practice neither of the above conditions is met.

First, even for experiments with the most densely sampled u–v coverage possible, its finite extent precludes a unique image reconstruction. That is, the image becomes unconstrained on angular scales smaller than  ∼ λ/u_max—this is simply the interferometric manifestation of the diffraction limit. However, more typical in VLBI is sparse u–v coverage, consisting of densely sampled tracks and significant gaps in which no visibilities are measured. As a result, the lack of uniqueness can become significant, extending even to intermediate- and large-scale features. Second, and particularly within the context of millimeter-wavelength VLBI (see, e.g., Event Horizon Telescope Collaboration et al. 2019a), the full complex visibilities are typically not available. Atmospheric phase delays induce potentially large station-specific phase errors. Imperfect antenna and receiver performance results in station gain amplitude errors. Both of these may be ameliorated by careful treatment of correlations among baselines, e.g., self-calibration, or via the use of closure quantities, e.g., closure phase and/or closure amplitudes, or some combination thereof. Again, these complicate and potentially alter the uniqueness of possible image solutions. Both of these are relevant for observations made by the Event Horizon Telescope (EHT), a global millimeter-VLBI array (Event Horizon Telescope Collaboration et al. 2019b, 2019c, 2019a, hereafter Papers I, II, and III, respectively).

Nevertheless, a number of efficient and high-fidelity methods have been developed to produce image reconstructions from interferometric data. These include the venerable CLEAN algorithm (Högbom 1974; Schwarz 1978; Clark 1980; Schwab 1984), which works directly with the inverse Fourier transform of the calibrated visibility data (the so-called "dirty image"). CLEAN attempts to model the image structure as a sum of point sources, whose locations and intensities are determined by iterative subtraction of appropriately scaled and shifted versions of the point-spread function (the so-called "dirty beam") from the dirty image. This procedure is typically repeated with interleaving "self-calibration" steps that attempt to infer the station-specific gain and phase errors consistent with the modeled image. The resulting collection of point sources is subsequently smoothed to the ostensible diffraction scale, i.e., the radio "beam," for display purposes, though the goodness-of-fit is computed directly from the point sources.

A number of imaging methods have also been developed that attempt to forward model the image directly, traditionally using "maximum entropy" methods (e.g., Frieden 1972; Gull & Daniell 1978; Narayan & Nityananda 1986), and more recently with "regularized maximum likelihood" (RML) methods (e.g., Chael et al. 2016, 2018; Akiyama et al. 2017a, 2017b). These algorithms have two main advantages over the traditional "inverse modeling" approach used by CLEAN: (1) they permit the imposition of general-purpose regularizers on the image structure, such as positivity, smoothness, or sparseness; (2) they can fit directly to closure quantities rather than requiring complex visibilities. As with CLEAN methods, RML modeling may be iterated with self-calibration steps, during which the station-based gain terms may also be subject to prior constraints. Both CLEAN and RML methods have been used with considerable success to reconstruct images of M87 from the 2017 April EHT observations (Event Horizon Telescope Collaboration et al. 2019d, hereafter Paper IV).

In principle, the forward-modeling approach used by RML imaging methods makes them amenable to posterior exploration within a Bayesian framework, generating not just a single optimal fit but the full range of images consistent with the data. Access to the posterior is highly desirable because it enables a quantification of the non-uniqueness of the image reconstruction; i.e., it permits an estimation of the uncertainty in the image in both a formal and practical sense. Posterior inference can also be used to naturally assess choices of resolution and field of view (FOV) of the image, systematically addressing two key parameters often left to the experience of the end user.

Previous efforts to perform Bayesian posterior exploration within the context of radio interferometric imaging have often focused on algorithms that function for connected-element arrays with excellent (u, v)-coverage and large FOVs (e.g., Sutter et al. 2014), motivating the development and use of sparsity priors (e.g., Cai et al. 2018a, 2018b) or resorting to single-point statistics that bypass the computationally taxing sampling (e.g., Greiner et al. 2016; Junklewitz et al. 2016; Naghibzadeh & van der Veen 2018). Nearly all of these efforts construct likelihoods that assume perfectly calibrated complex visibility data subject only to Gaussian thermal noise, and those that do attempt to fit for calibration terms either do not pursue image reconstruction (e.g., Lochner et al. 2015, who fit parameterized source models) or impose simplifying approximations about the structure of the posterior (e.g., Arras et al. 2019, who use variational inference).

In this paper we present an imaging framework implemented within Themis, a flexible and highly parallel model comparison and parameter estimation code developed for the EHT (Broderick et al. 2020). Themis has already been used to compare geometric and general relativistic magnetohydrodynamic (GRMHD) models to EHT data obtained from the 2017 April observing campaign (Event Horizon Telescope Collaboration et al. 2019e, 2019f, hereafter Papers V and VI, respectively). It also provides a powerful set of existing sampling, modeling, and calibration tools, including samplers that can find and explore multi-modal likelihoods, and the ability to trivially generate new models by combining existing models (see Broderick et al. 2020 for details).

In this way image reconstruction is treated in exactly the same standardized form as any other form of model fitting. By leveraging the ability of Themis to explore image posteriors, we gain access to direct estimates of the pixel uncertainties and their spatial correlations throughout the image. Perhaps even more importantly, the implementation of a Themis-based image model permits the self-consistent introduction of additional parameterized model components to the image. Such components can include those on scales smaller than the image pixel size or larger than the image FOV, enabling the search for and characterization of features that may be infeasible to image but which may nevertheless be constrained by the data in the Fourier domain.

Here we present a non-parametric image model and example image reconstructions with and without additional geometric model components. In particular, we demonstrate the ability to faithfully reconstruct bright ring-like features within images from GRMHD simulations appropriate for M87, and thus provide both a validation of the general underlying physical picture and a precision measure of the central supermassive black hole mass in M87 that is less sensitive to astrophysical uncertainties than prior methods. We summarize key relevant features of Themis in Section 2, present the image model in Section 3, provide the introduction and validation of additional geometric features in Sections 4 and 5, respectively, and collect conclusions in Section 6.

Themis is a Bayesian model comparison and feature extraction framework developed for the EHT (Broderick et al. 2020). Its modular design enables rapid development in most facets of the model analysis procedure. Two aspects are particularly relevant here: the ability to add any image-based models rapidly and implement new samplers.

All work reported here makes use of the parallel-tempered, differential-evolution Markov chain Monte Carlo (MCMC) sampler. This is an affine invariant, ensemble MCMC scheme, and thus optimizes the proposal distribution. Multiple tempering levels are used to efficiently locate and explore multi-modal likelihood distributions. This sampler has been demonstrated with prescribed/known likelihoods containing thousands of independent peaks (see Section 6.1.2 of Broderick et al. 2020). However, it is well known that this class of samplers can struggle when parameter sets exceed ≳100 (Huijser et al. 2015). While this sampler performance is sufficient for current EHT applications, additional sampler development is warranted prior to applications with a larger separation between the ostensible image resolution and FOV. Currently, the outputs of this sampler are MCMC chains that contain the full posterior probabilities of parameters of a given model.

Themis already has a number of simple geometric models implemented. These include asymmetric Gaussians and slashed rings, i.e., annuli with a linear brightness gradient across the ring (see Section 4.2 of Broderick et al. 2020). Methods for dynamically generating new image-based models from existing models extend these with minimal effort. A key example is the ability to sum images to generate new models; note that this sum is performed after visibility construction, and thus does not impose resolution restrictions from one model on others. Additional image-based models may be simply implemented, requiring only information about the interface and otherwise agnostic regarding the remainder of Themis.

Themis also incorporates the ability to mitigate key systematic uncertainties in EHT observations. Chief among these are the uncertain station gains, associated with atmospheric propagation and absorption. Phase delays are naturally addressed by the use of closure phases. An efficient and accurate method for reconstructing station gain amplitudes, even in large numbers, has been implemented (see Section 5 of Broderick et al. 2020). A simple gain model is effectively marginalized over after the construction of likelihoods, permitting Themis to reconstruct many gains without creating additional demands on the sampler.

Here we present analyses of various simulated data sets produced with u–v coverages identical and systematic error budgets similar to the 2017 observations for M87. The analyses are performed using the visibility amplitudes and closure phases, with the standard signal-to-noise ratio (S/N) restrictions and adopting the minimally covariant closure phase set (Blackburn et al. 2020; Broderick et al. 2020). In all cases we simultaneously reconstruct the individual station gain amplitudes while fitting image models as described in Broderick et al. (2020) and demonstrated in Paper VI.

The line between "imaging" and "modeling" is conceptually blurred, relying on a qualitative distinction between models that prescribe a specific class of features and those that do not. We approach imaging with Themis as a model-fitting endeavor through the use of an effectively non-parametric image model. In this case, we use the term "non-parametric" loosely; the model is formally parametric (see Section 3.1 for specifics), but the form of this parameterization is chosen to be flexible and to avoid steering the model toward any particular image structure.

We begin with a description of the underlying image model, discuss how resolution and FOV are set in a systematic fashion, and then present validation examples.

3.1. Bicubic Raster Model

The image model⁸ is generated in two steps: (1) amplitudes are selected for each of a number of "control points," followed by (2) convolving those amplitudes with a bicubic smoothing kernel.

First, ${N}_{\mathrm{px}}$ control points are placed at the vertices of a uniform square grid contained within a pre-defined FOV, as shown in Figure 1. These control points act as image pixels, the amplitudes of which are the primary model parameters; we enforce positivity at these control points by modeling the amplitudes logarithmically. Given an FOV and a pixel number per side N_{x, y} (such that ${N}_{\mathrm{px}}={N}_{x}\times {N}_{y}$), we set the control point positions (l_i, m_j) and intensities I_ij to be

$$\begin{eqnarray}&&{l}_{i}={s}_{x}i,\,\,{m}_{j}={s}_{y}j,\,\,{I}_{{ij}}={e}^{{p}_{{ij}}},\end{eqnarray} \tag{ 1 }$$

where ${s}_{x,y}=\mathrm{FOV}/({N}_{x,y}-1)$. We restrict the images considered in this paper to having ${N}_{x}={N}_{y}=\sqrt{{N}_{\mathrm{px}}}$, such that the images are square and s_x = s_y.

Zoom In Zoom Out Reset image size

Next, the grid of control points is convolved with a bicubic interpolation kernel to generate intensities that smoothly vary in all directions. From the I_ij values specified in Equation (1), we generate the actual image intensity at every (l, m) using

$$\begin{eqnarray}&&I(l,m)=\displaystyle \sum _{i,j}w(l-{l}_{i})w(m-{m}_{j}){I}_{{ij}},\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}w(x)=\left\{\begin{array}{ll}0 & -2\geqslant x\\ b\left[-{x}^{3}-5{x}^{2}-8x-4\right] & -1\geqslant x\gt -2\\ -(b+2){x}^{3}-(b+3){x}^{2}+1 & 0\geqslant x\gt -1\\ (b+2){x}^{3}-(b+3){x}^{2}+1 & 1\geqslant x\gt 0\\ b\left[{x}^{3}-5{x}^{2}+8x-4\right] & 2\geqslant x\gt 1\\ 0 & x\geqslant 2\end{array}\right.\end{eqnarray} \tag{ 3 }$$

is the 1D bicubic interpolation kernel and b is a control parameter that affects the monotonicity of the interpolation. For all of the images presented in this paper, we have set b = −0.5. We note that this bicubic interpolation can in principle result in negative values for I at some (l, m), despite the strictly positive nature of I_ij. These departures from positivity are controlled by the value of b, and in practice we have found that setting b = −0.5 ensures that any negative excursions are typically very small.

In practice, we perform the convolution described by Equation (2) in the visibility domain. That is, we set

$$\begin{eqnarray}&&V(u,v)=\displaystyle \sum _{i,j}W(2\pi {s}_{x}u)W(2\pi {s}_{y}v){e}^{2\pi i({l}_{i}u+{m}_{j}v)}{I}_{{ij}},\end{eqnarray} \tag{ 4 }$$

in which

$$\begin{eqnarray}\begin{array}{rcl}W(k) & = & -\displaystyle \frac{4}{{k}^{3}}\sin (k)\left[2b\cos (k)+(4b+3)\right]\\ & & +\,\displaystyle \frac{12}{{k}^{4}}\left\{b\left[1-\cos (2k)\right]+2\left[1-\cos (k)\right]\right\},\end{array}\end{eqnarray} \tag{ 5 }$$

is the Fourier transform of w(x) (see Appendix A). For modest values of N_x, such as those presented here and relevant for the EHT, V(u, v) is most rapidly computed via the direct sum in Equation (4); fast Fourier transforms become advantageous only for N_x ≳ 10².

In the absence of absolute phase information, as occurs here, the image is defined only up to an arbitrary translation. However, the sparse nature of the available control points coupled with the assumed lack of flux outside of the prescribed FOV typically forces the image toward a single location, alleviating this concern. Thus, we make no attempt to fix the location of the image, e.g., by setting the center of light or imposing a prior on the spatial distribution of flux in the image.

3.2. Pixel Size and FOV Selection

Within the context of VLBI, the smallest and largest recoverable image structures are approximately determined by the longest and shortest baselines in the array, respectively. Information about the source structure on spatial scales that are larger than those probed by the shortest baselines, on scales that are smaller than those probed by the longest baselines, or indeed on any scales that are not directly measured by the array, must therefore be guided by the image model. For our purposes, the relevant model hyperparameters are the control point separation (i.e., the pixel size) and the FOV, neither of which has a unique a priori specification. Other imaging methods have to contend with this same issue, but it is common practice to simply employ some rule of thumb when selecting these hyperparameters.

In our image modeling procedure, we would ideally pursue a strategy for selecting unique and data-driven values for both N_x and FOV via Bayesian evidence maximization. In practice, we approximate the Bayes factors using the Bayesian information criterion (BIC), which explicitly penalizes additional parameters:

$$\begin{eqnarray}&&\mathrm{BIC}=-2{{ \mathcal L }}_{\max }+k\mathrm{log}({N}_{\mathrm{data}}).\end{eqnarray} \tag{ 6 }$$

Here, ${{ \mathcal L }}_{\max }$ is the logarithm of the maximum likelihood, k is the number of model parameters, and ${N}_{\mathrm{data}}$ is the number of data points. We also make use of the Akaike information criterion (AIC), defined by

$$\begin{eqnarray}&&\mathrm{AIC}=-2{{ \mathcal L }}_{\max }+2k+\displaystyle \frac{2k(k+1)}{{N}_{\mathrm{data}}-k-1}.\end{eqnarray} \tag{ 7 }$$

Both the BIC and the AIC are listed in Table 1 for a simulated data set constructed to approximate the 2017 April EHT observations of M87. Both information criteria converge on the same optimal values for FOV and N_x, yielding a natural choice for both hyperparameters. We note that the selected FOV and N_x values remain subject to our assumption of a raster geometry. For example, we do not explore asymmetric or disconnected grids here. The inferred resolution (10 μas) and extent of the emission region are very similar to those found by other forward-modeling approaches to image reconstruction applied to the same simulated data set Paper IV. This implies a modest degree of superresolution, by roughly a factor of two, relative to the ostensible resolution of the EHT, approximately 20 μas.

Table 1.  Dependence of the BIC (Left Entry) and AIC (Right Entry) on N_x and FOV for a Simulated Data Set Designed to Approximate the 2017 April 11 EHT Observations of M87

	FOV (μas)
Nx	40	60	80	100	120
4	222/244	(2.02/2.04) × 103	(3.24/3.24) × 104	(5.38/5.38) × 105	(2.12/2.12) × 106
6	31.0/31.0	0/0	106/106	862/862	(5.81/5.81) × 103
8	217/207	178/168	175/165	169/159	246/236
10	446/471	408/433	402/427	394/419	390/415
12	731/902	689/859	689/860	670/840	672/843

Note. Both information criteria are quoted relative to their minima, which occur for N_x = 6, $\mathrm{FOV}=60\,\mu$as (in bold).

Download table as:  ASCIITypeset image

3.3. Image Reconstructions from Synthetic Data

In Figure 2 we present image reconstructions from a variety of simulated data sets with realistic u–v coverage and noise for the 2017 April EHT observations of M87. To facilitate comparison with other EHT imaging methods, we use the same simulated data sets for which reconstructions were presented in Figure 10 of Paper IV, which include realistic sources of systematic error, including large gain amplitude variations at a single station (Large Millimeter Telescope). These images were selected to both characterize the ability of the imaging procedures to reconstruct and probe ringlike structures (e.g., Ring, Crescent, GRMHD) and to assess their ability to distinguish ringlike from non-ringlike structures (e.g., Disk and Double).

Zoom In Zoom Out Reset image size

Prior to fitting, the data were pre-processed in the manner described in Paper VI; that is, we scan-averaged the visibility data and added a 1% systematic error contribution (consistent with the magnitude of the non-closing errors measured in Paper III) in quadrature to the thermal errors prior to constructing closure quantities. We jointly fit to the closure phase and debiased visibility amplitude data products from both the high and low frequency bands, and we excluded all data points having a S/N below 2. To accommodate potential structures that are resolved on all but the intrasite baselines (e.g., those confined to Hawaii; Papers IV; VI), all fits included a large-scale Gaussian component (see Section 4.1).

In all cases we initialize the image model with a broad Gaussian emission structure, though subsequent exploration renders the details of the initial parameter values moot. All tests were run multiple times, and their MCMC chains have reached the well-mixed regime.

Good fits are found for all reconstructions, with the maximum likelihood fit achieving a χ² per degree of freedom of 387/360, 385/360, 472/360, 370/360, and 330/360 for the Ring, Crescent, Disk, Double, and GRMHD tests, respectively. An example set of residual plots is shown in Figure 3 for the GRMHD test, and we see no obvious structure or poorly fit sections of data. To facilitate comparisons with other methods used for EHT image reconstruction, for each test we computed the equivalent blurring kernel defined and presented in Table 4 of Paper IV. We find similar values to those reported there: 14.8 μas, 14.6 μas, 28.0 μas, and 12.5 μas, for the Ring, Crescent, Disk, and Double tests, respectively.

Zoom In Zoom Out Reset image size

Despite the rectilinear nature of the control point gridding, a surprising variety of images is possible, as can be seen in Figure 2 (e.g., the reproduction of nearly circular Gaussian components in the Double test). On the smallest angular scales accessible to the model—i.e., the 10 μas inter-pixel spacing—we see imperfections in the maximum a posteriori reconstruction (e.g., knots on the Crescent, asymmetry in the components of the Double) associated with the limited (u, v)-coverage. Smoothing the models with a 15 μas FWHM Gaussian kernel mutes these small-scale structures and highlights the larger scales, for which the agreement between model and truth is more visually apparent (see the bottom two rows of Figure 2). However, we emphasize that even the apparent discrepancies in the unsmoothed model are statistically consistent with the truth, within the confidence levels claimed by the posterior. That is, the locations and magnitudes of the various small "artifacts" are no more than expected given the data-driven uncertainties in each pixel's amplitude, and the myriad realizations of all possible small-scale discrepancies are fully captured by the posterior distribution (see Section 3.4).

All of the image reconstructions shown in Figure 2 recover the scale and morphology of the truth images well. This includes the sizes and separations between features, a brightness gradient when present and its absence when not present, and the overall image orientations. Of particular relevance for detecting black hole shadows with the EHT is the ability to detect ring-like structures with high fidelity, that is, to reproduce rings when present in the underlying truth image but convincingly rule them out when they are not.

The overall flux normalization is a function of the station gain amplitude reconstruction, and thus subject to the associated priors. This limits the flux reconstructions to a precision of roughly 20% and is responsible for differences in the total brightness of some image reconstructions, e.g., the Disk test.

3.4. Image Posteriors

In practice, image parameter exploration as described in the previous sections is far more computationally expensive than the methods that have been used to date to generate images from EHT observations, which typically perform a constrained optimization of some comparison metric between image and data (see Paper IV). This increased expense is a consequence of the general parameter estimation paradigm, which considers the best fit to be less important than the full set of compatible fits. Thus, a key feature of image reconstruction within our modeling framework is the generation of a posterior distribution for all image parameters, which can be usefully conceptualized as an "image posterior." The image posterior can be used to quantitatively address a variety of questions about a given image reconstruction. Note that, despite the absence of absolute phase information, we do find that the image centroids are well fixed by the FOV combined with the modest number of pixels employed. As such, we neglect the potential contributions to the statistics we present here from translations of the reconstructed image.

A common question in radio interferometric imaging pertains to the reliability (or "believability") of specific image features. In Figure 4 we show a selection of 15 reconstructed images taken from the MCMC chain for the GRMHD test image (fifth column of Figure 2). Ten of these images (shown in the middle two rows of the figure) are sampled directly from the posterior distribution, such that the magnitude and frequency of any variations seen across the images mirror those present in the image posterior (e.g., a feature seen in two of the 10 images indicates that ∼20% of the image posterior space is consistent with that feature). Additionally, the bottom row of Figure 4 shows the five most extreme⁹ outlier images contained within 10⁴ elements randomly drawn from the MCMC chain; these images highlight image morphologies that live on the fringe of what the data will permit.

Zoom In Zoom Out Reset image size

Single-point statistics computed from the image posterior can be used to produce "typical" images and to provide pixel-specific error budgets. The top row of Figure 4 shows examples of such statistics in the form of the mean and standard deviation of the image posterior. We find that the standard deviation map is highly inhomogeneous, varying substantially across the image domain; such behavior is both expected and frequently observed in radio interferometric imaging (see, e.g., Figure 17 of Paper IV), despite often being presented as a single global rms value that is assumed to be shared by all pixels, but it is potentially critical for accurate image assessment.

Perhaps most striking is the clear ringing present within the standard deviation map, associated with the control point structure. The origin of this ringing is elucidated by the pairwise marginal posterior distributions for pixels across the image, shown in Figure 5 for a single row and column. The nodes in the standard deviation map are associated with a set of sinusoidal oscillations with a wavelength of  ∼ 20 μas, corresponding to the maximum spatial frequencies accessible to the 2017 EHT (u, v)-coverage. For M87 in particular, limited north–south (u, v)-coverage (see Figure 12 in Paper III) increases the uncertainty in the mode amplitude in that direction. In contrast, we find little to no correlation between pixels in the east–west direction, where the 2017 EHT (u, v)-coverage was more complete for M87. That is, the joint pixel flux distributions indicate both the presence of the expected resolution limit and the magnitude of the uncertainties this limit induces in the reconstructed image (roughly 10% in this case).

Zoom In Zoom Out Reset image size

With access to n-point statistics provided by the image posterior, arbitrarily sophisticated image interrogation procedures could be developed, such as metrics to assess how likely it is that the image contains a ring-like feature. The ability to generate statistically rigorous posteriors is also critical for applications such as presented in Section 4, where rigidly parameterized modeling is performed simultaneously with image reconstructions.

Visibility-domain modeling of VLBI data can be a powerful tool for robust parameter estimation, but it assumes that the underlying source structure is well-described by the model under consideration. This assumption is rarely satisfied for simple models when the data have high S/N, at which point it becomes necessary to include additional and flexible model complexity to account for deviations between the desired model structure and the true source structure. For a primary EHT science target, M87, the additional complexities have been associated with some combination of stochastic fluctuations in the emission region (as caused by, e.g., turbulence) and uncalibrated systematic uncertainties in the data. The flexibility provided by imaging has played a critical role in constraining the nature of these deviations and in guiding the introduction of ad hoc model features (e.g., the "nuisance Gaussian components" used alongside the geometric crescent models in Paper VI) to address them.

The implementation in the previous section of imaging within a modeling framework enables a novel path forward toward realizing the best elements of both approaches: (a) modeling known features and thus accessing substantial improvements in their parameter estimation, and (b) imaging the remainder of the source structure, yielding an otherwise agnostic view of the stochastic aspects of the image. We refer to this simultaneous modeling and imaging approach as "hybrid imaging."¹⁰

There are two key elements to successful hybrid imaging. First, because a primary goal is the accurate reconstruction of the parameters of the model components, the ability to explore the posterior distribution for the full image/model combination is required. Where simultaneous exploration is not possible—e.g., when modeling separately from or iteratively with finding best-fit images—the sampler will only have restricted access to the parameter space and the reconstructed posteriors will not necessarily be reliable. Thus, the ability to construct full image priors as described in the previous section is critical.

Second, some ordering parameter in which the model and image components can be distinguished is highly desirable if not strictly necessary. This ordering may be accomplished by enforcing a separation in angular scale, variability, polarization properties, wavelength, etc. Where features in the image component can be reconstructed independently by either the model or image components, a degeneracy is induced between the two. There are many problems for which these components are well separated in a natural way, and in this section we describe two such cases that are of particular relevance to EHT applications.

4.1. Large-scale Flux

The primary EHT targets exhibit structure over a wide range of angular scales, including on those much larger than the ∼200 μas FOV of the EHT images. This large-scale image structure impacts only the shortest (intrasite, <10 km) baselines while having little effect on intermediate and long baselines (see, e.g., Papers IV; VI). Nevertheless, these short baselines provide important constraints on station gains and removing them would thus be undesirable.

To address the intrasite baseline flux excess, it has been found sufficient to include a large-scale (FWHM ≳ 1 mas) Gaussian model component in the image generation and model comparison exercises. A Gaussian is chosen simply to have a structure-agnostic model component with a specifiable size scale. This component is naturally separated in spatial scale from the much more compact (FWHM < 100 μas) image features produced in the reconstructions (see Figure 6). We include such a Gaussian in all of the examples presented in this and the following section for this purpose.

Zoom In Zoom Out Reset image size

4.2. Photon Rings

In this section we describe a number of imaging experiments designed to assess and demonstrate the ability of existing and future EHT observations to detect and constrain photon rings. We focus on several example simulations from the GRMHD library presented in Paper V, all of which contain a bright, narrow ringlike feature. These simulations differ substantially in their other properties, however, including the brightness and structure of the extended emission, presence of additional ringlike and spiral features, asymmetric components, and dynamical structures. For GRMHD simulations appropriate for M87 we are able to faithfully reconstruct the properties of the bright ringlike feature, and recover its radius.

5.1. Ring Reconstruction from GRMHD Simulations

Five example GRMHD movie snapshots that are broadly consistent with the EHT 2017 observations of M87 were selected, listed in Table 2. All of these come from simulations that were deemed acceptable in Paper V and were resized to best match the EHT data and rescaled to a total flux of 0.6 Jy. Data were prepared with the u–v coverage from the 2017 EHT observations of M87 from Paper III. The fitting procedure is identical to that described in Section 3.3, but now with the inclusion of a ring model.

Table 2.  Properties of GRMHD Simulations in Figure 7 (from Left to Right) and Associated Reconstructions

a	i	Rhi	Flow	rph	rpeak	rring	Ftotal	Fring/Ftotal	w/rring
			Type	(μas)	(μas)	(μas)	(Jy)
0.5	158°	10	MAD	21.2	21.4	${21.8}_{-0.4}^{+0.3}$	0.50 ± 0.03	${0.57}_{-0.03}^{+0.04}$	0.02 ± 0.01
0.94	158°	40	MAD	18.7	19.3	19.3 ± 1.0	0.42 ± 0.02	0.44 ± 0.05	${0.02}_{-0.01}^{+0.02}$
−0.94	17°	20	SANE	19.6	20.5	20.0 ± 0.9	0.37 ± 0.02	0.42 ± 0.04	0.03 ± 0.02
0.94	163°	160	SANE	16.9	17.0	${17.9}_{-0.4}^{+0.3}$	0.47 ± 0.02	${0.52}_{-0.05}^{+0.04}$	0.03 ± 0.02
0.94	163°	160	SANE	16.9	16.6	${16.6}_{-0.6}^{+0.5}$	${0.47}_{-0.01}^{+0.02}$	0.40 ± 0.03	0.03 ± 0.02

Note. Listed are the simulation parameters (spin, inclination, R_hi, flow type) and the associated average asymptotic photon ring radius (r_ph), radius of the peak of the azimuthally averaged flux (r_peak), and reconstructed ring radius (median and 68% confidence interval), total compact flux (F_total), fraction of the flux contained in the ring (F_ring/F_total), and the fractional width of the ring (w/r_ring).

Download table as:  ASCIITypeset image

Figure 7 shows the reconstructions in comparison to the truths both before and after smoothing. In each instance high-quality fits exist with χ²/DoF of 130/165, 447/453, 276/253, 618/592, and 539/596, for the models as shown from left to right.¹¹

Zoom In Zoom Out Reset image size

The slashed ring model component accurately reproduces the bright ringlike features in each image. The remaining image component reconstructs the low-brightness extended emission. There is little traction on the ring thickness, for which the posterior distribution remains dominated by the imposed uniform prior restricting the fractional width to below 5%. The weakness of this constraint is consistent with the fundamental resolution limits imposed by the u–v coverage. In contrast, the total flux and ring diameter are well constrained. The combined ring component and image component together do an excellent job of reproducing the image morphology. In all cases the ring component is modestly subdominant, contributing between 40% and 60% of the total flux in the image.

The quantitative success of the ring reconstruction is presented in Table 2 and Figure 8, where the azimuthally averaged intensity profiles from the ground truth image are compared against the posteriors on the ring radii. In each case, the posterior covers in part or in its entirety the bright excess associated with the n = 1 photon ring. In the case of the MAD, a = 0.94 model the posterior is double peaked, associated with the two apparent ringlike structures in the underlying truth image.

Zoom In Zoom Out Reset image size

The fourth simulation standard and normal evolution (SANE; a = 0.94) has multiple ringlike features appearing at significantly different radii. Of these, only one is associated with the gravitational, n = 1 photon ring feature, the others with orbiting, turbulent features. This results in multiple peaks in the azimuthally averaged brightness profile, limiting the formal applicability of the single-ring model. Despite this, the ring radius posterior has a significant overlap with the brightness peak associated with the n = 1 photon ring.

However, these additional ringlike features are highly dynamic and expected to evolve significantly on timescales of weeks in M87. To determine how the ring reconstruction depends on the dynamical state of the emitting gas, we repeated the reconstruction exercise with a frame taken from the same GRMHD simulation taken at a sufficiently different time so as to be effectively uncorrelated. This time, the resulting radius posterior matches the location of the bright peak with a similar fidelity to the other GRMHD simulation experiments. More importantly, the two ring radii reconstructions are inconsistent with each other at the 2.7σ level. This suggests that repeated measurements separated by many dynamical times presents a natural strategy for distinguishing transient turbulent structures that arise due to astrophysical processes from the anticipated bright lensing features caused by gravity.

5.2. Gravitational Inferences from Ring Reconstructions

Encoded within the size and shape of the asymptotic photon ring is the structure of the underlying spacetime (Johannsen & Psaltis 2010). However, the ability to infer the gravitational properties of the black hole, e.g., the mass or spin, depends on the relationship between the location of the bright ringlike feature, presumably the n = 1 photon ring, and the asymptotic photon ring. In all of the GRMHD simulations listed in Table 2 the peak of the bright ring is biased away from the asymptotic photon ring. Three factors are relevant to the interpretation and magnitude of this bias.

First, note that the underlying resolution of the images from which the simulated data were constructed is 1 μas. This dominates the widths of the flux profiles in Figure 8. Thus, even in principle, the precision with which the bright ring radius can be constrained is ≳0.1 μas, limited fundamentally by the number of pixels covered by the ring. In all instances in Table 2 the disparity is larger than this ostensible precision limit, though less than the underlying pixel resolution. The role of GRMHD simulation resolution will be explored elsewhere.

Second, a dominant contribution to the ring flux is from the n = 1 photon ring, which is produced by photons that have executed half orbits from emission to detection (see, e.g., Johnson et al. 2020). These are biased away from the asymptotic photon ring ($n\to \infty$). Typical values for this bias from GRMHD simulations are ≈0.5 μas (0.13 GM/c²D) (Johnson et al. 2020) and are consistent with the biases observed in the first three GRMHD simulations in Table 2. The direction and precise value of this bias does depend on the distribution of emission near the black hole. However a strong upper limit is imposed by the condition that polar photon trajectories deflect around the black hole sufficiently to return to the equatorial plane (see Figure 9 and Appendix B), implying that the bias is less than 4.2 μas (1.1 GM/c²D) generally and typically considerably less than 1 μas (0.26 GM/c²D) for emission regions similar to those inferred for M87 from the geometric crescent fits described in Paper VI and found in GRMHD simulations of M87.

Zoom In Zoom Out Reset image size

Third, dynamical features will present transient flux excesses inside and outside the asymptotic ring that bias the radial location of the flux peak. This bias is clearly exemplified in the last two entries of Table 2, where the radius of the flux peak within a single simulations varies by 0.4 μas over long times (>20 days).

However, despite the potential sources of bias in the location of flux peak, the posteriors of the reconstructed radii for the ring model component are consistent with the expected asymptotic photon ring radius. In the first three entries of the Table 2, for which a ringlike structure is present, the ring radius estimates are consistent with the asymptotic photon ring size at the 1σ level. For the final GRMHD simulation, the reconstructed ring radius estimates are consistent at the 2σ level despite the presence of confounding features. An exhaustive study of the accuracy of the ring reconstructions for different underlying emission region morphologies will be presented in future publications.

Fully Bayesian image reconstruction from VLBI observations is now practical. This method has been implemented in a fashion appropriate for EHT applications within the Themis analysis framework by treating the image reconstruction process identically with existing model-fitting tasks. Employing Themis for this has exploited two key advantages: first, it permits the immediate and efficient application of high-performance computing resources to the reconstruction of radio images with realistic observational parameters and complications; the Themis image reconstruction tools can efficiently make use of thousands of cores simultaneously. Second, it ensures that additional development within Themis is immediately available, including the implementation of additional samplers and mitigation of new sources of systematic error. These image reconstruction tools have been successfully applied to a number of simulated EHT data sets. Thus, Bayesian image reconstruction is now a practical tool for radio image analysis, though additional work on samplers may be required to fully realize its potential.

Performing image reconstruction in a Bayesian fashion presents two unique opportunities that arise due to treating the image reconstruction problem as a parameter estimation exercise in a statistically rigorous manner. First, a consequence of this approach is the construction of a full image posterior, i.e., not a single best-fit image but rather a statistically meaningful collection of all of the images that are consistent with the data. In parameter estimation applications, this posterior is the primary product. In this case, it permits a quantitative assessment of the sources and magnitudes of errors. For example, two-point functions on the image posterior elucidate the origin and size of the error induced by known holes in the u–v coverage of the 2017 April EHT observations of M87. More complex n-point functions may be developed in the future to address specific questions about image uniqueness and fidelity.

Second, treating image reconstruction and model parameter estimation identically enables seamless integration of the two methods of data analysis, a procedure we term hybrid imaging. Importantly, the reconstructed posteriors on model parameter estimates are meaningful in the normal way.

We also take this opportunity to demonstrate the ability to make a data-informed choice of the FOV and resolution, removing a frequent uncertainty in the image reconstruction process and source of subjective operator input.

We demonstrate hybrid imaging with two explicit applications. The first is the inclusion of large-scale Gaussian components to model structure on scales much larger than the ostensible FOV, a procedure that has been utilized in previous EHT analyses (Papers IV; VI). The second and more exciting application is the inclusion of narrow ring model components to extract and characterize the n = 1 photon rings from EHT observations of M87.

Hybrid imaging of GRMHD simulations with realistic simulated EHT observations of M87 are able to faithfully reconstruct the properties of the prominent photon ring feature in the underlying truth images. These produce ring size posteriors that typically cover the true value within accuracies of 2%–5%. Where turbulent features are mistaken for photon rings, the ring size estimates are variable on timescales of weeks in M87. Thus, the benefits of non-parametric reconstructions combined with the strong restrictions on potential additional image features are already poised to significantly improve EHT constraints on the properties of the supermassive black hole in M87.

Photon ring characterization in Sagittarius A*, the supermassive black hole in the Galactic center, brings the additional challenges of addressing the interstellar scattering screen and short-timescale source variability. Both of these will be addressed in a future publication. One considerable benefit of doing so is that high-precision mass and distance estimates exist for Sgr A* (Gravity Collaboration et al. 2019), enabling precision tests of general relativity using the photon ring.

There are many potential future applications of hybrid imaging. All of these exploit its ability to constrain the freedom of image components via model specification, where such constraints are required due to experimental considerations (e.g., baseline coverage) or desirable because of high sensitivity. Such applications within the context of the EHT include the connection of small-scale features with large-scale model structures, e.g., the modeling of the millimeter-wavelength milliarcsecond jet in M87 while imaging the horizon-scale millimeter-wavelength image structure. As the EHT adds stations and increases in sensitivity, hybrid imaging may be able to constrain the shape of the photon ring, producing an independent measure of black hole spin, and to detect higher-order rings, verifying a key predicted consequence of gravitational lensing. At the same time, hybrid imaging effectively separates the otherwise unresolved ring component from the larger-scale diffuse emission, and thus will enable dynamical studies on timescales as short as a week in M87.

Finally, hybrid imaging permits an extension of the class of sources for which shadows may be detected and characterized beyond the two primary targets of the EHT, for which the shadow is formally resolvable. By assuming that a crescent-like structure exists and contributes substantially to the emission in the bright radio cores, hybrid imaging extends the ability to estimate black hole masses to thousands of systems.

We thank Mareki Honma, Michael Johnson, Maciek Wielgus, and Katherine Boumann for useful comments. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET:www.sharcnet.ca) and Compute/Calcul Canada (www.computecanada.ca). Computations were made on the supercomputer Mammouth Parallèle 2 from University of Sherbrooke, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the ministère de l'Économie, de la science et de l'innovation du Québec (MESI) and the Fonds de recherche du Québec—Nature et technologies (FRQ-NT). This work was supported in part by 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 Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade. A.E.B. thanks the Delaney Family for their generous financial support via the Delaney Family John A. Wheeler Chair at Perimeter Institute. A.E.B. and P.T. receive additional financial support from the Natural Sciences and Engineering Research Council of Canada through a Discovery Grant. R.G. receives additional support from the ERC synergy grant "BlackHoleCam: Imaging the Event Horizon of Black Holes" (grant No. 610058). D.W.P. was supported in part by the Black Hole Initiative at Harvard University, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation to Harvard University.

Here we derive the Fourier space representation of the bicubic convolution kernel given in Equation (5). The weight function in Equation (3) is C₂. However it does have discontinuities in the second and third derivatives, being identically zero beyond that. Thus it may be represented as a triple and quadruple integral over a small number of Dirac δ functions.

The second derivative of the cubic interpolation kernel is

$$\begin{eqnarray}\begin{array}{c}{w}^{{\rm{{\prime} }}{\rm{{\prime} }}}(x)=\left\{\begin{array}{ll}0 & -2\geqslant x\\ b\left[-6x-10\right] & -1\geqslant x\gt -2\\ -6(b+2)x-2(b+3) & 0\geqslant x\gt 1\\ 6(b+2)x-2(b+3) & 1\geqslant x\gt 0\\ b\left[6x-10\right] & 2\geqslant x\gt 1\\ 0 & x\geqslant 2.\end{array}\right.\end{array}\end{eqnarray} \tag{ A1 }$$

Note that this is not continuous, with the discontinuities of ±(−8b+6) and ∓2b at x = ± 1 and x = ± 2, respectively.

The third derivative has two components, the first describing the smooth pieces of w''(x) and the second consisting of δ functions associated with the discontinuities. That is,

$$\begin{eqnarray}&&w\prime\prime\prime (x)=g(x)+{F}_{3}(x)\end{eqnarray} \tag{ A2 }$$

where

$$\begin{eqnarray}\begin{array}{c}g(x)=\left\{\begin{array}{ll}0 & -2\geqslant x\\ -6b & -1\geqslant x\gt -2\\ -6(b+2) & 0\geqslant x\gt -1\\ 6(b+2) & 1\geqslant x\gt 0\\ 6b & 2\gt \geqslant x\gt 1\\ 0 & x\geqslant 2,\end{array}\right.\end{array}\end{eqnarray} \tag{ A3 }$$

and

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{F}_{3}(x) & = & 2b\delta (x+2)+(8b+6)\delta (x+1)\\ & & -\,(8b+6)\delta (x-1)-2b\delta (x-2).\end{array}\end{array}\end{eqnarray} \tag{ A4 }$$

Again, there are discontinuities in the non-singular part of w‴, i.e., g(x), though now at x = 0, x = ± 1, and x = ± 2; g(x) is constant otherwise. Thus, the fourth derivative consists solely of δ functions, w^(iv)(x) = F₄(x) where

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{F}_{4}(x) & = & -6b\delta (x+2)-12\delta (x+1)\\ & & +\,12(b+2)\delta (x)-12\delta (x-1)-6b\delta (x-2).\end{array}\end{array}\end{eqnarray} \tag{ A5 }$$

Then, noting that all constants of integration vanish (because $f(x)=0$ for $| x| \gt 2$), we can write

$$\begin{eqnarray}&&w(x)=\iiint {F}_{3}(x)+\iiint \int {F}_{4}(x).\end{eqnarray} \tag{ A6 }$$

The frequency-space representation of F₃ and F₄ are

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{F}_{3}(k) & = & \displaystyle \int {{dxe}}^{-{ikx}}{F}_{3}(x)\\ & = & 4i\sin (k)\left[2b\cos (k)+(4b+3)\right]\\ {F}_{4}(k) & = & \displaystyle \int {{dxe}}^{-{ikx}}{F}_{4}(x)\\ & = & 12a\left[1-\cos (2k)\right]+24\left[1-\cos (k)\right],\end{array}\end{array}\end{eqnarray} \tag{ A7 }$$

where k = 2πΔu, where Δ is the spacing between control points. Therefore, applying the integrals,

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}W(k) & = & \displaystyle \frac{{F}_{3}(k)}{{\left({ik}\right)}^{3}}+\displaystyle \frac{{F}_{4}\left(k\right)}{{\left({ik}\right)}^{4}}=-\displaystyle \frac{4}{{k}^{3}}\sin (k)\\ & & \times \,\left[2b\cos (k)+(4b+3)\right]\\ & & +\,\displaystyle \frac{12}{{k}^{4}}\left\{b\left[1-\cos (2k)\right]+2\left[1-\cos (k)\right]\right\}.\end{array}\end{array}\end{eqnarray} \tag{ A8 }$$

In Figure 9 we show the mapping between points on the equatorial and image planes for observers on the polar axis. Here we collect information about how this mapping is generated practically and what it implies generally.

The mapping is generated by direct integration of the null geodesic equations for photons leaving a distant screen. We simplify the problem by considering an observer on the polar axis, for which axisymmetry reduces the function to a one-dimensional ${\rm{\Delta }}{r}_{n}(R)$, which is a mapping of the Boyer–Lindquist radius R in the equatorial plane (at the (n + 1)th photon crossing) to the radial location in the image plane. While this does not generalize to arbitrary inclinations, the insensitivity of the asymptotic photon ring radius to inclinations below 30° suggests that this is sufficient to estimate the magnitudes of the shifts of the low-order photon ring radii from the asymptotic value (Johannsen & Psaltis 2010).

This mapping is monotonic, a fact that has two important consequences. First, the entirety of the equatorial plane is remapped at every order, and thus the photon rings present a sequence of copies of the direct emission. Second, and more relevant here, is that the peak of the intensity within the nth photon ring is fully determined by the location peak of the emission, R_pk, and the map. This arises immediately from the conservation of I_ν/ν³ along the photon trajectories. Thus, for optically thin emission within the equatorial plane, the total intensity at a radial position of r in the image plane is

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{\nu }}{{\nu }^{3}}(r)=\displaystyle \sum _{m=0}^{n}\displaystyle \frac{{I}_{g\nu }}{{\left(g\nu \right)}^{3}}({R}_{m}),\end{eqnarray} \tag{ B1 }$$

where g is the standard Doppler factor, ${\rm{\Delta }}{r}_{m}({R}_{m})\equiv r-{r}_{\mathrm{ph}}$ defines R_m, and n is the highest-order image present at r. This is peaked when

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d({I}_{\nu }^{\mathrm{obs}}/{\nu }^{3})}{d{\rm{\Delta }}r} & = & \displaystyle \sum _{m=0}^{n}\displaystyle \frac{d\left[{I}_{g\nu }^{\mathrm{em}}/{\left(g\nu \right)}^{3}\right]}{{{dR}}_{m}}\displaystyle \frac{{{dR}}_{m}}{d{\rm{\Delta }}r}\\ & \approx & \displaystyle \frac{d\left[{I}_{g\nu }^{\mathrm{em}}/{\left(g\nu \right)}^{3}\right]}{{{dR}}_{n}}\displaystyle \frac{{{dR}}_{n}}{d{\rm{\Delta }}r}=0,\end{array}\end{eqnarray} \tag{ B2 }$$

where we have used the fact that the ring widths decrease exponentially with n, and thus for smoothly varying emission regions dR_n/dΔ r increases exponentially with n (see, e.g., Johnson et al. 2020). Therefore, if R_pk is the location of the peak of ${I}_{g\nu }^{\mathrm{em}}/{(g\nu )}^{3}$, i.e., the peak of the emission after inclusion of the Doppler shift and Doppler beaming, then the peak of the nth-order photon ring emission is ${\rm{\Delta }}{r}_{n,\mathrm{pk}}\approx {\rm{\Delta }}{r}_{n}({R}_{\mathrm{pk}})$.