Hybrid mapping of the Black Hole Shadow in M87

The EHT data are described in detail in Event Horizon Telescope Collaboration et al. (2019c). In brief, observations were made of the nuclear regions of the nearby radio galaxy, M87 (Virgo A), with the goal of imaging the event horizon of the hypothesized supermassive black hole. Observations were made on four days in April (5th to 11th), 2017, at 227.1 GHz and 229.1 GHz, each with a total bandwidth of 1.875 GHz, using an array comprised of seven telescopes spanning the globe, including Europe, South America, continental USA, and Hawaii. The observations each day spanned about 6 hours, in a series of $\sim 7$ min scans, interspersed with observations of other sources.

The EHT collaboration provides public data that have had a priori gain (visibility flux density) calibration applied based on the measured system parameters at each telescope, as well as delay calibration via visibility fringe fitting (Cotton, 1995; Walker, 1999), plus further adjustments based on a few redundant baselines in the array (Event Horizon Telescope Collaboration et al., 2019c). The gain calibration provides reasonable visibility amplitudes (to within $\sim 10\%$). The initial calibration provides enough phase stability to average the data in time to 10 s records, and in frequency to a single 1.875 GHz channel. The EHT collaboration designates these data as the ‘network-calibrated data’, and we follow this designation below.

However, the EHT collaboration emphasize that the initial calibration alone does not allow for phase coherent imaging (Event Horizon Telescope Collaboration et al., 2019a), due to residual phase errors arising from e.g. errors in station clocks, the tropospheric model, or polarization leakage. Subsequent element-based phase self-calibration is required to produce a phase coherent astronomical image. They state: ‘Lack of absolute phase information and a priori calibration uncertainties in EHT measurements require multiple consecutive iterations of CLEAN followed by self-calibration, a routine that solves for station gains to maximize consistency with visibilities of a specified trial image (Pearson & Readhead, 1984)’.

Starting with the network-calibrated data, the EHT collaboration explored three separate imaging and self-calibration processes, including inverse-modeling and forward-modeling approaches. The forward-modeling approach searches through image-plane parameter space for sky surface brightness models which best match the visibility measurements (based on some goodness-of-fit criterion). The model-to-data comparison is made between the visibility amplitudes and closure phases on closed triads of array stations. The antenna complex gains are adjusted to optimize the fit. Closure phases are used because they are insensitive to antenna-based phase errors (Jennison, 1958; Readhead & Wilkinson, 1978; Schwab, 1980; Wilkinson, 1989; Cornwell & Fomalont, 1999; Thompson et al., 2017), although closure phases do not constrain arbitrary sky-plane translations of images made from independent closed triads (Readhead & Wilkinson, 1978; Thyagarajan & Carilli, 2020), thereby requiring some initial alignment procedure using a prior sky-model.

The inverse-modeling approach employed by the EHT team corresponded to a standard hybrid mapping (self-calibration and imaging) approach using DIFMAP²²2https://github.com/eventhorizontelescope/2019-D01-02/tree/master/difmap (Shepherd, 1997). Hybrid imaging uses a simple starting model for initial phase self-calibration of the data. Self-calibration entails making adjustments to antenna-based complex gains (phase, or phase and amplitude), to obtain enough coherence to synthesize initial images of source structure³³3Note that by restricting corrections to antenna-based calibration terms only, these corrections necessarily conserve closure phase measurements on closed triads of antennas (Thompson et al., 2017; Wilkinson, 1989; Readhead & Wilkinson, 1978; Thyagarajan & Carilli, 2020)..

The station-based gain adjustments are derived by minimizing $\chi^{2}$ of the differences between the visibility data and visibilities derived from the Fourier transform of the sky brightness model (or another goodness-of-fit criterion). Subsequent iterations of calibration then use source models derived from the data itself, to derive further corrections of station-based phases and amplitude. The typical process employs the CLEAN imaging and deconvolution algorithm, which effectively breaks down the source structure into point sources, with the self-calibration model then corresponding to the point source CLEAN components (Perley, 1999; Walker, 1999).

A starting model is required for self-calibration to obtain enough coherence to produce a rough image of the source, and launch the hybrid mapping process. As a starting model in the DIFMAP inverse-modeling approach, the EHT collaboration employed the simplest of possible starting models, namely, a point source (Gomez priv. comm.). They also restrict the subsequent CLEAN models to disk-shaped regions in the image plane, for further self-calibration. They state (Event Horizon Telescope Collaboration et al., 2019a): ‘A disk-shaped set of CLEAN windows, or imaging mask, aligned with the previously found geometrical center of the underlying emission structure (e.g.the ring in M87), and with a specified diameter is then supplied; these windows define the area on the image where the CLEAN algorithm searches for point sources. We limited the cleaning windows to image only the compact structure (100 $\mu$as) in order to prevent CLEAN from adding larger-scale emission features that are poorly constrained by the lack of short EHT baselines.’ For the forward-modeling approach starting model, the EHT collaboration denote a Gaussian image as the prior (Event Horizon Telescope Collaboration et al., 2019a): ‘Similar to the EHT-imaging and DIFMAP scripts, the SMILI imaging procedure is iterative, with four stages of imaging and self-calibration. Reconstructions at each stage begin with a circular Gaussian with a FWHM of 20 $\mu$as as the initial image. After the first imaging attempt in each stage, subsequent initializations use the previously obtained image convolved with the initial Gaussian.’

The EHT collaboration explored the robustness of their imaging and self-calibration process by generating source models of very different types (double Gaussian sources, disk, annulus, etc.), adding realistic thermal noise and station-based complex gain errors, and then reconstructing the images via their imaging and self-calibration processes (Event Horizon Telescope Collaboration et al., 2019a). This analysis resulted in choosing imaging and self-calibration procedures that were robust to different morphological features. These procedures were then applied in the imaging of the M87 data, using the simple starting models, as described above.

Here, we present a representative exploration of how changes to the starting models in the self-calibration process itself, can affect the final outcome of the Hybrid imaging process for the EHT observations of the ring in M87 itself. The starting model is an issue because, for relatively sparse VLBI arrays, such as the EHT, and extended source morphologies, it is well known that, depending on the complexity of the source and the density of UV-coverage, it is possible that the iterative self-calibration process can ‘turn the data into the model’ (Perley, 1999; Readhead & Wilkinson, 1978; Walker, 1999), meaning, the optimization problem may be under-, or marginally, constrained. This issue is particularly true for the EHT data, in which the ALMA array, acting as a single element in the array, is close to two orders of magnitude more sensitive than any other element in the array. We return to this issue below.

Figure 1 shows the network-calibrated (ie. a priori flux density calibrated), visibility data from the EHT observations of M87 (Event Horizon Telescope Collaboration et al., 2019c). Plotted in the figure are the visibility amplitudes vs. UV-distance. The EHT collaboration estimate errors on the flux density scale for the visibilities of typically around 10% to 20%, and this error range is supported by subsequent amplitude calibration corrections in the self-calibration process (see Figure 3). Also shown are visibility phases vs. UV-distance, after self-calibration. We discuss this latter plot in Section 3.3.

This uv-plot already tells us much about the source. First, the source is clearly well resolved, with rapidly decreasing flux density versus baseline length. Second, the ‘depression’, or possible null, in the uv-plot suggests something non-Gaussian, with a fairly circularly symmetric surface brightness distribution and a hard-edge, since, for instance, a Gaussian source does not produce any nulls, while a double Gaussian or double point source distribution would have a series of nulls that depend on the position angles of the visibility fringes relative to the source separation. More complex assemblies of point sources will generally have UV-distributions that are broad in amplitude across UV-distance (see Pearson (1999) and Purcell & Wilkinson (2019)⁴⁴4http://www.jb.man.ac.uk/ira4/IRA4%20SuppMat_Chapter9.htm). The UV distribution has the characteristics of a $J_{1}$ Bessel function, corresponding to the Fourier transform of a disk or annulus (an annulus corresponds to the difference of two $J_{1}$ Bessel functions of different diameters). The position of the depression at about $4\times 10^{9}$ wavelengths implies a characteristic size scale $\sim 50~{}\mu$as. We will use this information to guide our starting models, as was also true for the processing employed by EHT team.

We performed our hybrid mapping process using the Astronomical Image Processing System (AIPS; Greisen, 2003). Again, hybrid mapping is a standard ‘inverse modeling’ method to derive element based gains, and the sky brightness, via iterative image generation and self-calibration using the resulting images themselves in each iteration (Pearson & Readhead, 1984; Purcell & Wilkinson, 2019).

There are a number of challenges with the EHT data. First, the coverage of the Fourier plane (uv-coverage) is relatively sparse (see Fig. 12 in Event Horizon Telescope Collaboration et al., 2019c), and the source has significant structure. Second, the initial sky surface brightness model derived from a Fourier transform of the network-calibrated data shows no coherent structure, and hence does not provide a useful starting model.

And third, the ALMA array, as one element in the VLBI array, has roughly 60 times the sensitivity of any other element (Event Horizon Telescope Collaboration et al., 2019c). Hence, if the weights in the visibility data are used in the self-calibration process, ALMA acts as an unmovable anchor, around which the other station gains are adjusted, hence leading to a situation where supposed antenna-based gain corrections approach baseline-based gain corrections. The latter leads to the clear danger of ‘turning the data into the model’ (Perley, 1999). Moreover, when required gain corrections are dominated by systematic errors, and not by thermal noise, then weighting by the thermal sensitivity of a given station may not be optimal (Carilli et al., 1991).⁵⁵5Further, while the signal-to-noise of baselines that contain ALMA are much higher than for non-ALMA baselines, in procedures that employ model fitting to closure phases, every closure triad that includes ALMA must also include one baseline that does not involve ALMA. Hence, in terms of the closure phase, triads that include ALMA are at most root(3) lower noise than those without ALMA.

We address these difficulties as follows. First, we re-weight the visibility data, to lower the dominance of ALMA on the calibration process, as was also done by the EHT collaboration (Event Horizon Telescope Collaboration et al., 2019a), by taking the square root of the data weights for all baselines⁶⁶6We explored other weighting schemes, such as uniform weights across the EHT array, but found a square root filter provided the best results.. We have flagged the very short baselines (ALMA to APEX and JCMT to SMA) as was done in Event Horizon Telescope Collaboration et al. (2019c) before any processing.

The first iteration of self-calibration employs the model image, and phase-only self-calibration. In every self-calibration iteration, we derive scan-averaged solutions, impose a minimum number of array stations of four, and set a signal-to-noise threshold of 2.5 $\sigma$ for acceptable station gain solutions. We use a linear (L1) minimization criterion, and not $\chi^{2}$, to reduce the impact of outlier measurements⁷⁷7L1 minimization implies a weighted sum of the moduli of the residuals is minimized (Schwab, 1981). After the first model-based self-calibration, we then derive source models from the data using the CLEAN imaging algorithm. For the imaging, we employed IMAGR in AIPS, with robust weighting ($R=0$), a cell size of 2 $\mu$as, and a loop gain of 0.03. The resulting CLEAN Gaussian restoring beam has a FWHM of about $22\,\mu$as$\times 16\,\mu$as, major axis position angle of $25^{\circ}$, in all cases, and we restore the final CLEAN images with this beam.

We perform one iteration of phase self-calibration using the a priori model, and then proceed through a second iteration of phase self-calibration using the CLEAN component sky model derived from these data. We then employ two iterations of amplitude and phase self-calibration, with the improving CLEAN models. There is an inevitable subjectivity in the processing due to the fact that the EHT collaboration have already presented robust images of the ring in M87. However, it remains true that the characteristic angular scale of the source is already fairly well determined by Figure 1, as was also employed as a guide in the EHT imaging paper. To systemitize the model generation process, in each iteration of CLEAN we set CLEAN boxes to cover regions with $\geq 5\sigma$ surface brightness, and we stop the CLEAN process at 1 mJy component flux density.

A total of four iterations of self-calibration was adequate to achieve convergence (two phase, two amplitude and phase), where convergence was measured by a change of less than 10% in the rms and peak negative sidelobe, in the final image.

We begin by showing two examples of the calibration solutions that are derived at key steps in the hybrid mapping process. We look at the solutions for two different models: the annulus and the point source. Comparative results for the other models are similar.

Figure 2 shows the phase solutions after the first phase self-calibration using the a priori sky model, for baselines involving ALMA. There is clear coherence in time for the phase solutions: smooth phase gradients over time are seen, and the gradients from the point source and the annulus models parallel each other with time. There are also clear phase offsets between the point source and the annular models. These systematic offsets between models demonstrate that self-calibration does not preserve absolute astrometry. In particular, phase calibration with a point source model at the phase center will tend to shift the brightest feature in the sky brightness to the phase center. We shall see that trend in the final images below.

Figure 3 shows the amplitude solutions after the first iteration of phase and amplitude self-calibration. The amplitude corrections are generally $\pm$10% to 20%, with the LMT showing the largest corrections, as was also found in the EHT collaboration processing (Event Horizon Telescope Collaboration et al., 2019a). Again, the amplitude corrections for the point-source starting model and the annular starting model, track each other reasonably. More uv-data are flagged due to the signal-to-noise criteria for solutions with the point source starting model, but only at the start of the observation on some baselines to the Pico Veleta antenna.

Figure 4 shows the final images at 229.1 GHz for the different starting models, as well as for the network-calibrated data. The network-calibrated data clearly shows little/no coherence, as expected. We re-emphasize that, in all cases, the hybrid mapping process followed a uniform sequence of (i) setting CLEAN boxes to cover regions with $\geq 5\sigma$ surface brightness, (ii) cleaning to a uniform residual level, (iii) stopping after 4 iterations (two phase only, two amplitude and phase self-calibration).

The self-calibrated images all show an extended source with a scale $\sim 40~{}\mu$as, consistent with the UV plot shown in Figure 1. Table 1 shows the results for the maximum surface brightness, rms noise, maximum negative artifact, and total flux density in the final images.

The hybrid mapping processes starting with the annulus and disk models converge on a ring-like structure, with a clear central depression. These two images have the best image-quality metrics, meaning off-source noise and artifacts.

In the images starting with a point source, Gaussian, and an asymmetric double source, the source does appear extended on a similar scale of $\sim 50~{}\mu$as. However, the images themselves are substantially noisier, by 30% to 60% (see Table 1), and with larger peak residuals by similar factors. The asymmetric double starting model produces the poorest of the final images, in terms of noise and artifacts. In these cases, the self-calibration process shifts the peak of the final image to the position of the peak on the input model, since phase self-calibration does not preserve absolute astrometry (Readhead & Wilkinson, 1978). And while depressions do exist within the source perimeter, any possible ring is not as clearly defined in the final images.

Figure 5 shows the closure phases on six representative triads for one day of data, compared to the model derived from the CLEAN components in the final imaging stage. In all cases, we CLEAN to the same number of iterations. The five different starting models in the self-calibration process are shown in different colors. Column 6 in Table 1 lists the rms deviation (represented by square-root of the normalized $\chi^{2}$) of closure phases between the final CLEAN models and the data for all days, and all triads in which the time on-source was more than 50 min, and the rms deviations of the data itself were $<50^{\circ}$. The rms deviations between data and model are normalized for each triad by the rms of the data on each triad, and then summed to give a global, normalized rms value for all the data.

The normalized rms closure phase deviations between model and data are very similar for the annulus, disk, and Gaussian starting models, while higher for the point source starting model, and highest for the asymmetric double model, although only by about 50%. Closure phase is a useful diagnostic in that the values are robust of antenna-based calibration terms. However, the images themselves contain added information from the visibility amplitudes, and hence provide potentially deeper insight into the quality of the final product from the self-calibration process.

Overall, hybrid imaging of EHT interferometric visibility data, following a uniform set of steps in the process, leads to an extended source on a scale of $\sim 40~{}\mu$as. Starting the hybrid mapping process with a disk or annulus model converges to the best images, in terms of image noise and artifacts, and closure phase residuals. These images show well defined rings. Starting with point, Gaussian, and asymmetric double models show sources with internal depressions, but any ring itself is less clearly defined.

We emphasize that in no case does the deepest depression within the source go to zero surface brightness. The deepest depression still has a surface brightness greater than $10\sigma$, and only a factor 2.5 below the minimum surface brightness along the ring itself in e.g. the annulus starting model image. Indeed, whether there is a true annulus, with zero intensity in the center, as opposed to, e.g. an edge-brightened disk, is impossible to tell from these data for two reasons. First is that the FWHM of the synthesized beam is only a factor two smaller than the ring diameter. And second is, again, the sparse UV-coverage.

Readhead & Wilkinson (1978) showed early in the history of VLBI that imaging an edge-brightened, and disk-like (their figure 8), sources with internal structure using a sparse array is challenging: simply going from a four to five element array makes substantial improvements in recovered structure. The EHT has 7 elements, although the mutual visibilities are not continuous, meaning the UV-coverage changes substantially over the course of the observation. Sometimes only four antennas are on source, in particular at the start of the observation. Moreover, the Readhead & Wilkinson (1978) study employed large antennas and bright sources at low frequency, and hence were of uniformly high signal-to-noise, as shown in Wilkinson et al. (1979). The EHT has large differences in signal-to-noise for different baselines, and many low signal-to-noise baselines that do not include ALMA. Hence, the challenges cited by Readhead & Wilkinson (1978) may be even greater for the EHT, and it is not surprising that the EHT hybrid imaging results, while suggestive, are not conclusive as to the ring structure in M87.

One simple argument for a ring or disk remains the UV-plot shown in Figure 1. For the visibility amplitudes vs. UV-distance, we plot the amplitudes with only network-calibration (a priori gains) applied, i.e., before any self-calibration. Again, the depression, or possible null, in the amplitude vs. distance distribution, argues for an annular or disk model (meaning, a UV-distance vs. amplitude plot more like a $J_{1}$ Bessel function), as opposed to e.g. a double source, which will have nulls that vary with the relative projection of the fringe spacing along the double source position angle. Generally, on a UV-plot involving multiple antennas and just the UV-distance (i.e. projected length of baseline, but no position angle information), such as Figure 1, collections of point sources will typically not show a distinct null (Pearson, 1999).

In Figure 1 we show two model visibility distributions. The models are for a disk and an annulus. The sizes for the models were guided by the results from the hybrid mapping process, such that the position of the null corresponds roughly to the lowest visibility amplitudes (see Section 3.4. The model amplitudes were then adjusted to best fit the visibilities (lowest unnormalized $\chi^{2}$). The data are noisy, and the position of the null is gross, but the models do show one significant difference: the second peak of the disk model is distinctly lower than that for the annular model (Pearson, 1999; Purcell & Wilkinson, 2019). Specifically, the flat disk model falls well below the data itself. This difference in the visibility amplitude vs. UV-distance for a disk vs. annulus vs. Gaussian was also noted in Figure 1 in Event Horizon Telescope Collaboration et al. (2019b).

For the phases, the lack of phase coherence in the network-calibrated data necessitates the use of the self-calibrated data. The visibility phase vs. UV-distance plot for the models shows either zero, or $\pm 180^{o}$ phases, since the models are circularly symmetric, and centered on the phase center. The simple point we want to make from Figure 1 is that the disk and the annulus models show very similar behaviour in phase vs. UV-distance. The implication is that the visibility phase distribution is not a strong diagnostic for differentiating between an annulus or disk.

We then continue the processing of data using the annular starting model. First, we process the 227.1 GHz data in an identical manner as was employed at 229.1 GHz, specifically, limiting the CLEAN boxes for regions around the ring/disk. At both frequencies, a modest extension appears to the Southwest (see Appendix A). As a check, we also processed the combined frequency data, but for the four days separately, and again see the extension in all days, although the detailed structure varies from day to day and frequency to frequency, indicating the limitations of the data and process, and/or real time variability (see Appendix A).

We then combined the two frequency datasets, and continued the self-calibration process, but now extending a single square clean box to 76 $\mu$as in full width, and centered on the disk/ring. This will include the disk/ring area, and just beyond the disk/ring in all directions, including the modest extension to the southwest. Using a larger, symmetric CLEAN box, centered on the disk/ring, avoids bias toward generating structure in any particular direction around the disk/ring. We performed one more iteration of amplitude and phase self-calibration. The result is shown in the left panel of Figure 6. Note that the color scale has been adjusted to emphasize that the surface brightness at the center of the disk/ring is still well above the noise ($>10\sigma$).

The ring itself becomes better defined, with a diameter between maxima around the annulus of $\sim 40~{}\mu$as to $\sim 44~{}\mu$as, consistent with the imaging by the EHT collaboration (Event Horizon Telescope Collaboration et al., 2019a). In all data sets (days and frequencies; Section A), the southern half of the ring shows about a factor two higher surface brightness, although, again, the details are not exactly identical across days and frequencies. The total flux density in this image is 243 mJy, with a peak surface brightness of 50.6 mJy beam^-1, an rms noise of 0.4 mJy beam^-1, and a peak negative sidelobe of $-1.5$ mJy beam^-1. The rms noise and peak negative sidelobe in this image is about a factor 1.5 to 2 lower than those shown in Figure 4. Of course, this image includes two times more data, and further self-calibration. The modest extension to the southwest becomes more prominent, extending about 60 $\mu$as from the center of the ring at a surface brightness level of $\sim 4$ mJy beam^-1.

Initially, we considered the southwest emission to be a likely artifact of the processing, since the main jet on arcsecond (ie. kpc) scales in M87 is seen at a position of about $+18^{\circ}$ north of west ($=288^{\circ}$ east of north; Walker et al. (2018)). However, we reinspected the latest 90 GHz VLBI images of Kim et al. (2018a) (see also EHT MWL Science Working Group et al. (2021)). In these images, the jet on scales below a few milliarcseconds is strongly edge-brightened, with an opening angle $\geq 40^{\circ}$. In particular, on scales at the limit of their resolution (beam FWHM $=120~{}\mu$as$\,\times\,50~{}\mu$as, major axis position axis north-south Kim et al. (2018a)), they find in the inner $\sim 200~{}\mu$as, the brightest limb of the jet extends in the southwest direction from the core, at a position angle of $\approx 15^{\circ}$ south of west ($255^{\circ}$ east of north), similar to the observed direction of the extension in Figure 6 (left panel). For reference, we produce an overlay of the publicly available 90 GHz VLBI image (Kim et al., 2018a), in Figure 7.

Given the very limited Fourier coverage of the EHT data, we emphasize caution concerning any jet. The right panel of Figure 8 shows an overlay of the EHT image with a contour plot of the synthesized beam (PSF) for these data, where the PSF has been shifted to the peak surface brightness location around the ring/disk. The PSF has peak sidelobes of about 59%, located in two directions: one roughly North-South, and a second along a position angle that corresponds roughly to the modest extension to the southwest. Hence, it is possible the extension is an artefact of the PSF shape for these observations.

To explore this latter possibility further, we start with a uniform brightness annulus model, with characteristics (size and uniform surface brightness) guided by the observational results. We then Fourier transform and sample this model with the exact UV-sampling of the EHT data itself. We add thermal noise appropriate for the array, and we add systematic phase errors at levels consistent with the large phase corrections we found in the first iteration of phase self-calibration, with phase offsets of up to $\pm 100^{o}$ (Figure 2). We then process the corrupted simulated data through the hybrid mapping procedure in an identical manner as for the real data. The resulting image is also shown in Figure 8 (right panel).

The image resulting from this simulated data has a peak surface brightness, rms noise, and off-source artifacts, similar to the real data. However, the image does not show any extension in any direction comparable to what is seen with the actual data itself. At the lowest contour level (about 4$\sigma$), there are weak ‘ears’ noticeable in the direction of the southwest-northeast PSF sidelobes, but again, at a level and size much reduced relative to the southwest extension seen in the real data.

Overall, the southwest extension is not reproduced in the simulated data, but it does repeat from day to day, and frequency to frequency (although not identically) in the real data. On the other hand, the PSF at both frequencies and all days is comparable (meaning, all the data have comparable UV-coverage), and one of the main sidelobes of the PSF is along this southwest direction. Hence, we cannot be certain of reality of this extension as a core jet in M87, and encourage further observations with improved UV-coverage, to verify or discover any possible core-jet in this important target.

We note recent work by Arras et al. (2020), who present evidence for a similar southwest extension, independent of (and mutually unknown to) our analysis. This is not surprising, given we employ the same data. However, they follow a forward-modeling inference approach, as opposed to our inverse-modeling approach. They converge on a ring, with a possible extension to the southwest, and possibly also to the northwest. They reach the interesting conclusion that the detailed structure around the ring, and of the possible extensions beyond the ring, may be time variable on timescales of days, although the essential ring-geometry itself is robust from day to day (see also Satapathy et al. (2021)). Such variability is physically reasonable, given the size of the regions, and likely relativistic characteristic velocities (Wielgus et al., 2020). Strong variability could affect any final image made when combining all the data in the hybrid mapping sequence, and could be a limiting factor, depending on the level of variability. We discuss the time variability in Appendix A.