First Sagittarius A* Event Horizon Telescope Results. III. Imaging of the Galactic Center Supermassive Black Hole

The highly sensitive phased-ALMA array participated in three of the five observing days, 2017 April 6, 7, and 11. April 7 is the only day that additionally includes PV observations of Sgr A* and is therefore the day with the longest observation duration, the largest number of detections, and the best overall (u, v)-coverage, as shown in Figure 1. On April 11 an X-ray flare was reported shortly before the start of the EHT observations (Paper II). Strongly enhanced flux density variability is seen in the light curves on that day (Wielgus et al. 2022), possibly posing difficulties for the static imaging of the April 11 data set. These constraints motivate utilizing the less variable 2017 April 7 data as the primary data set for static image reconstruction, with the April 6 observations as a secondary validation data set. Analysis of the remaining 2017 EHT observations of Sgr A* will be presented elsewhere.

Zoom In Zoom Out Reset image size

In Figure 1, the (u, v)-coverage on 2017 April 7 is shown to be asymmetric, with the longest baselines along the north–south direction. The shortest baselines in the EHT are intrasite and sensitive to arcsecond-scale structure (i.e., the SMA and JCMT are separated by 0.16 km; ALMA and APEX are separated by 2.6 km). In contrast, the longest baselines are sensitive to microarcsecond-scale structure (see Table 1). The ∼8.7 Gλ detections on PV–SPT and SMT–SPT baselines are among the longest published projected baseline lengths obtained with ground-based VLBI, alongside the recent EHT 3C 279 results of Kim et al. (2020), slightly longer than the longest baselines in the EHT observations of M87* (8.3 Gλ; M87* Paper IV). Sgr A* was detected on all baselines between stations with mutual visibility, leading to the April 7 (u, v)-coverage approaching the best one theoretically possible with the EHT 2017 array. Table 1 shows the angular resolutions derived from (u, v)-coverage on both April 6 and 7 data.

Table 1. Metrics of EHT Angular Resolution for the 2017 Observations of Sgr A* for the 229.1 GHz Band

	FWHMmaj	${\mathrm{FWHM}}_{\min }$	P.A.
	(μas)	(μas)	(deg)
Minimum Fringe Spacing $1/| {\boldsymbol{u}}{| }_{\max }$ (All Baselines)
Apr. 6	24.2	⋯	⋯
Apr. 7	23.7	⋯	⋯
Minimum Fringe Spacing (ALMA Baselines)
Apr. 6–7	28.6	⋯	⋯
CLEAN Beam (Uniform Weighting)
Apr. 6	24.8	15.3	67.0
Apr. 7	23.0	15.3	66.6
CLEAN Restoring Beam (Used in This Paper)
Apr. 6–7	20	20	⋯

Note. In order to avoid asymmetries introduced by restoring beams, and to homogenize the images among epochs, we adopt a circular Gaussian restoring beam with 20 μas FWHM for all CLEAN reconstructions.

Download table as:  ASCIITypeset image

In the top panel of Figure 2 we show the S/N of the Sgr A* observations as a function of projected baseline length, for the coherent averaging time of 120 s. The split in S/N distributions at various projected baseline lengths is due to the difference in sensitivity for the colocated Chile sites ALMA and APEX, with ALMA baselines yielding detections stronger by about an order of magnitude. In the bottom panel, we show the visibility amplitude (correlated flux density in units of Jy) for Sgr A* as a function of projected baseline length after applying full data calibration.

Zoom In Zoom Out Reset image size

The fully calibrated visibility amplitudes exhibit a prominent secondary peak between two local minima. The first minimum is located at ∼3.0 Gλ and is probed by the Chile–LMT north–south baselines on 2017 April 7. On April 6 recording started about 2 hr later, thus missing the relevant detections, as shown in Figure 3. The second minimum appears at ∼6.5 Gλ, probed by the Chile–Hawai'i baselines on both 2017 April 6 and 7. Overall amplitude structure of the source appears to be consistent across both days, which is particularly well visible in the fully calibrated, light-curve-normalized data sets shown in Figure 3, as the light-curve normalization procedure strongly suppresses the large-scale source intrinsic variability (Broderick et al., 2022). The observed local visibility amplitude minima can be associated with the nulls of the Bessel function J₀, corresponding to the Fourier transform of an infinitely thin ring. For a ring that is 54 μas in diameter we would obtain local amplitude minima at 2.92 and 6.71 Gλ. This is illustrated in the bottom panel of Figure 2, where an analytic Fourier transform of an infinitely thin ring blurred with a 23 μas FWHM Gaussian kernel is shown with dashed lines. ¹⁴⁹ While a blurred ring model roughly captures the dependence of visibility amplitudes on projected baseline length, there is also a clear indication of the source asymmetry, manifesting as amplitude differences between the Chile–LMT and SMT–Hawai'i baselines at the first minimum, probing the same range in projected baseline length (∼2.5–3.5 Gλ) in orthogonal directions. There is also a deficit of flux density with respect to the simple ring model at projected baseline lengths of ∼5–6 Gλ.

Zoom In Zoom Out Reset image size

Finally, we detect very clear and unambiguous nonzero closure phases, which are indicative of source asymmetry. With multiple independent triangles and high S/N, these data sets offer insight into the source phase structure greatly surpassing that of any previous mm-wavelength observations (Fish et al. 2016; Lu et al. 2018). In Figure 4, we show examples of closure phases on several triangles exhibiting various degrees of probed asymmetry and inter/intraday time variability. Closure phases on ALMA–SMT–SMA and ALMA–LMT–SMA triangles immediately show interday variability of the source structure. In the case of Sgr A*, intrinsic source variability is expected also on timescales as short as minutes, adding to the closure phase intraday variability caused by the nontrivial average structure of the source (see Section 3.2). Very long baselines, such as LMT–SPT in the ALMA–LMT–SPT triangle, are additionally affected by the presence of refractive scattering structure (see Section 3.1).

Zoom In Zoom Out Reset image size

Typically, the antenna gains and sensitivities as a function of elevation are derived using polynomial fits to opacity-corrected antenna temperature measurements from quasars and solar system objects, tracked over a wide range of elevations. Residual errors in the characterization of these antenna gains lead to corruptions in the flux calibration of the visibilities. Quantifying these effects enables us to disentangle astrophysical variability of Sgr A* from apparent flux variations caused by the imperfect calibration. In this work, we mitigate the systematic gain errors in the Sgr A* data sets based on the analysis of the calibrator sources (J1924–2914 and NRAO 530), which remained stationary in their source structure and flux density on relevant timescales. The detailed procedure to estimate the antenna gains from the two calibrators is described in Section 5.1.3 in Paper II. In particular, it is shown that the a priori gains are 5%–15% for all baselines, except intrasite baselines (∼1%) and those including the LMT (∼35%).

In addition to the thermal noise and the antenna gain uncertainties, we also estimate the nonclosing errors in the data, based on deviations of the trivial closure quantities from zero, as well as from the observed inconsistencies in the distributions of closure quantities between bands and polarizations (M87* Paper III; Paper II). These nonclosing errors are expected to arise from the presence of a small circular polarization component, as well as from uncorrected polarimetric leakages, and other systematic errors, such as residual bandpass effects. For Sgr A*, the nonclosing errors are estimated to be 2° in closure phase and 8% in log closure amplitude (Paper II). Assuming that the errors are baseline independent, these translate to 1° systematic nonclosing uncertainties in visibility phases and 4% systematic nonclosing uncertainties in visibility amplitudes, on the top of the uncertainties related to the amplitude gain calibration. We found that the RR − LL discrepancies in closure quantities are more significant for Sgr A* than in the case of the calibrators. This hints at an intrinsic source property, possibly a contribution from a small circular polarization component. This is consistent with Goddi et al. (2021), reporting ∼ 1% circular polarization in the simultaneous ALMA-only data.

Angular broadening is described by a convolution of an unscattered image with a scattering kernel, or equivalently by a multiplication of unscattering interferometric visibilities by the appropriate Fourier-conjugate kernel. The Fourier-conjugate kernel is given by $\exp \left[-\tfrac{1}{2}{D}_{\phi }({\boldsymbol{b}}/(1+{ \mathcal M }))\right]$, where b is the baseline vector of the interferometer and ${ \mathcal M }$ is the magnification of the scattering screen given by the ratio of the observer-to-screen distance to the screen-to-source distance. Interferometric measurements of Sgr A* with the EHT at the observing wavelength of 1.3 mm are primarily obtained on long baselines of $| {\boldsymbol{b}}| \gtrsim (1+{ \mathcal M }){r}_{\mathrm{in}}$, or equivalently on angular scales of $\theta \lesssim \lambda /(1+{ \mathcal M }){r}_{\mathrm{in}}$, where r_in is the inner scale of the fluctuations. In this regime, the angular broadening is affected by the power-law density fluctuations on scales between the inner and outer scales, giving the phase structure function of D_ϕ ( r ) ∝ λ²∣ r ∣^−α , where α = β − 2.

In Figure 5, we show the scattering kernel in the visibility domains based on the scattering parameters in Johnson et al. (2018). Johnson et al. (2018) imply a near-Kolmogorov power-law spectral index β ∼ 3.38 (or α ∼ 1.38), providing a non-Gaussian kernel more compact than the conventional Gaussian kernel. Consequently, the angular broadening effect, i.e., multiplication of the intrinsic visibilities with the Fourier-conjugate kernel of scattering, causes a slight decrease in visibility amplitudes and therefore also the S/N by a factor of a few at maximum.

Zoom In Zoom Out Reset image size

Angular broadening provides deterministic and multiplicative effects on the observed visibility. Therefore, it is invertible—they can be mitigated by dividing the observed visibility and associated uncertainties by the diffractive kernel visibility (often called deblurring; Fish et al. 2014). However, the actual interferometric measurements of Sgr A* with the EHT have contributions from substructure that arises from the refractive scattering, often referred to as the "refractive noise." The refractive noise is stochastic and additive and therefore not invertible (Johnson & Narayan 2016). In fact, since refractive effects are included in the diffractively blurred image, the refractive noise will be amplified by simply deblurring with the scattering kernel and will likely create artifacts in the reconstructions if not accounted for (Johnson 2016). To avoid effects from refractive noise, we expand the noise budgets of the visibility data prior to deblurring, as described in Section 3.1.2.

The contribution of the refractive noise to the observed visibility is anticipated to be not dominant, except for a small fraction of data beyond ∼6 Gλ (Figure 5). ¹⁵⁰ Signature of the refractive noise, namely, the long, flat tail of the visibility amplitude at long baselines found in recent longer-wavelength observations of Sgr A* (Johnson & Gwinn 2015; Johnson et al. 2018; Issaoun et al. 2019, 2021; Cho et al. 2022), is not clearly seen in the EHT data. In this EHT regime, where the refractive substructure is not unambiguously constrained from data, it is challenging to apply complex strategies that account for the stochastic properties or explicitly recover the refractive screen (e.g., Johnson 2016; Johnson et al. 2018; Issaoun et al. 2019; Broderick et al. 2020a).

We instead mitigate the effect of refractive substructure by introducing error budget models that approximate the anticipated refractive noise. The visibility error budget is increased based on these models prior to the mitigation of angular broadening via deblurring (Section 3.1.1). In this work, we consider four base models that approximate the refractive noise budgets: (i) Const: a constant noise floor (e.g., 10 mJy) for all baselines motivated by the fact that the refractive noise has a mostly flat profile as a function of the baseline length (see Figure 5). (ii) J18model1: (u, v)-dependent noise floor based on the scattering model and parameters described in Psaltis et al. (2018) and Johnson et al. (2018). Since the scattered image is not unique, we have simulated hundreds of scattering realizations and generated corresponding synthetic data that match the (u, v)-coverage of the actual April 7 observation of Sgr A*. The refractive noise values are then computed by taking the standard deviation of the complex visibilities across different realizations. Since the refractive noise is also dependent on the intrinsic source structure, in this case we consider a circular Gaussian model with the second moment that matches the pre-imaging size constraints (see Paper II). (iii) J18model2: Same as J18model1, but using the average refractive noise value of seven geometric models as possible intrinsic source structures (see Section 5). (iv) Considering not only the standard deviation of the refractive effects but also their correlations via a covariance matrix.

Note that using all the information encoded in the covariance matrix (not only in the variance of the refractive noise variables) will provide a better approximation of the refractive noise. However, the short time cadence and the redundant baselines in our data make this covariance matrix noninvertible and thus difficult to use in imaging. For the remaining three refractive noise models, we compute the complex visibility χ² for a suite of synthetic data (corrupted only by thermal noise along with scattering effects) based on seven geometric models of the intrinsic source structure (Section 5),

$$\begin{eqnarray}&&{\chi }^{2}\left({\sigma }_{\mathrm{ref},i}\right)=\frac{1}{2{N}_{\mathrm{vis}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{vis}}}\frac{| {V}_{i}-{V}_{\mathrm{ea}}{| }^{2}}{{\sigma }_{\mathrm{th},i}^{2}+{\sigma }_{\mathrm{ref},i}^{2}},\end{eqnarray} \tag{ 1 }$$

where V_i are the data visibilities for each scattering realization, V_ea are the ensemble-averaged visibilities (i.e., corresponding to the image experiencing only diffractive scattering), σ_th,i is the thermal noise, and σ_ref,i is the corresponding refractive noise budget for each strategy. The χ² metric provides us with a statistic on how well the ensemble average image represents the synthetic data after taking into account the different modeled refractive noise budgets. Figure 6 shows a comparison of the different χ² distributions for 400 realizations of synthetic data for every scattering mitigation strategy. For J18model1 we have derived a scaling factor to make the median of χ² of all models equal to 1 in order to overcome the dependence of the refractive noise level on the intrinsic source structure (see Appendix B).

Zoom In Zoom Out Reset image size

Figure 6 demonstrates that the (Const, J18model1, and J18model2) refractive noise models result in reasonable χ² values for the simulated scattering realizations we tested, with a majority of the χ² values below 3.0 (ideally χ² = 1). For this reason, in the rest of the paper we focus on the simpler Const and J18model1 strategies.

The EHT Sgr A* data contain unambiguous signatures of an evolving emission structure. On the largest spatial scales, the light curve varies at the ∼10% level on timescales of ∼hours (Wielgus et al. 2022). Variations in excess of those expected from thermal noise are seen on timescales as short as ∼1 minute. Over timescales typical of observation scans, ∼10 minutes, the degree of variation is on the order of 5% (Wielgus et al. 2022).

Direct evidence for short-timescale structural variations may be found in the evolution of closure quantities. Closure phases measured on certain triangles of baselines (e.g., ALMA–SMT–SMA) exhibit significant differences between April 6 and 7 (see Figure 4). The variations seen in the closure phases measured on multiple triangles show significant excesses, relative to thermal noise levels, as captured using the ${ \mathcal Q }$-metric statistic (Roelofs et al. 2017; Paper II).

Nonparametric estimates of the degree of variability as a function of baseline length may be generated by inspecting the visibility amplitudes directly. This is made possible by two fortuitous facts: first, the existence of crossing tracks, and second, that Sgr A* was observed on multiple days. As a result, many presumably independent realizations of the source structure may be compared. Practically, this is obtained by collecting visibility amplitudes in (u, v) bins, linearly detrending to remove the contribution from the static component of the image, and computing the mean and variance of the residuals. We average these estimates azimuthally to improve the significance of variability detection. For details on the procedure and validation examples, we direct the reader to Georgiev et al. (2022) and Paper IV.

In the case of Sgr A*, this nonparametric estimate produces a clear detection of variability that is significantly in excess of the expected thermal noise (Paper IV). Within the range of baseline lengths over which meaningful estimates can be produced, roughly 2 Gλ < ∣u∣ < 6 Gλ, the observed excess variability is broadly consistent with that anticipated by GRMHD simulations, both in magnitude and in dependence on baseline length (Papers IV and V).

Strategies for imaging in the face of source variability can be classified into one of three general categories:

• 1.  
  Variability reconstruction, or "dynamic imaging," in which the evolution of the source emission structure is explicitly recovered during the imaging process. The output of this strategy is a movie of the source emission structure. We refer the reader to Section 4.4 for more discussion on methods for dynamic imaging.
• 2.  
  Variability circumvention, or "snapshot imaging," in which standard image reconstruction is performed on segments of data ("snapshots") that are sufficiently short that the source may be approximated as static across them. The output of this strategy is a time series of static images.
• 3.  
  Variability mitigation, or "variability noise modeling," in which the impact of structural changes in the visibilities is absorbed into an appropriately inflated error budget. The output of this strategy is a single static image of the source, indicative of the time-averaged image over this observation period.

In practice, for segments of data short enough that Sgr A* may be reasonably approximated as static, the (u, v)-coverage of the EHT is insufficient to support reliable snapshot image reconstruction (though more restrictive parameterizations of the source structure, such as permitted using geometric modeling, can still be applied; see Section 9 and Paper IV). However, because dynamic imaging enforces a degree of temporal continuity, it is able to leverage the information provided by densely covered intervals of time to augment the lack of information available during intervals of sparser coverage. Dynamic imaging can thus be thought of as a generalization of both standard (static) imaging and snapshot imaging, with the former being equivalent to dynamic imaging with maximal temporal continuity enforcement and the latter being equivalent to dynamic imaging with no temporal continuity enforcement at all. Because dynamic imaging falls in between these two extremes, it can potentially recover reliable source structure in regions where the data are both too variable for standard imaging and too sparse for snapshot imaging. Our efforts to perform dynamic imaging in the most densely (u, v)-covered regions of data are described in Section 9.

The third strategy listed above—the variability noise modeling approach—permits static images to be reconstructed even from time-variable data. Depending on the specifics of the implementation, the recovered image captures some representation of "typical" source structure. Imaging with variability noise modeling requires that the error budget of the data be first inflated in a way to capture the statistics—or "noise"—of the source variability. The specific form of the noise model we use in this work is a broken power law, for which the variance ${\sigma }_{\mathrm{var}}^{2}$, as a function of baseline length ∣u∣, takes the form

$$\begin{eqnarray}&&{\sigma }_{\mathrm{var}}^{2}(| u| )={a}^{2}{\left(\displaystyle \frac{| u| }{4\,{\rm{G}}\lambda }\right)}^{c}\displaystyle \frac{1+{\left[(4\,{\rm{G}}\lambda )/{u}_{0}\right]}^{b+c}}{1+{\left(| u| /{u}_{0}\right)}^{b+c}}.\end{eqnarray} \tag{ 2 }$$

Here $| u| \equiv \sqrt{{u}^{2}+{v}^{2}}$ is the dimensionless length of the baseline located at (u, v), u₀ is the baseline length corresponding to the break in the power law, a is the variability noise amplitude at a baseline length of 4 Gλ, and b and c are the long- and short-baseline power-law indices, respectively. Equation (2) represents the variance that is associated with structural variability after removing the mean and normalizing by the light curve; see Paper IV for details.

By adding the variability noise given by Equation (2) in quadrature to the uncertainty of every visibility data point, the image becomes constrained to fit each data point to only within the tolerance permitted by the expected source variability. This parameterized variability noise model is generic and can explain well a wide range of source evolution, including complicated physical GRMHD simulations of Sgr A* (Georgiev et al. 2022).

Paper IV presents a nonparametric analysis of Sgr A*'s variability, which is further inspected to provide ranges of broken power-law model parameters that fit Sgr A* data (see Georgiev et al. 2022). As explained in Paper IV, given the baseline coverage of the 2017 Sgr A* campaign, little traction is found on the location of the break, u₀, and the short-baseline power law, c. However, the amplitude a is well constrained with an interquartile range from 1.9% to 2.1%. Similarly, the long-baseline power law, b, is modestly constrained, with interquartile range 2.2–3.2. These interquartile ranges are used to provide approximate priors on the variability noise that should be considered during static imaging reconstruction; values employed are listed in Table 2 for Sgr A*, as well as for a number of synthetic data sets described in Section 5. In the case of the latter, theoretical considerations imply that under very general conditions c ≈ 2, and thus we adopt a general prior of [1.5, 2.5]. For the CLEAN (Section 4.1) and regularized maximum likelihood (RML; Section 4.2) imaging surveys described in Section 6, variability noise models in the identified Sgr A* range are added to the visibility noise budget before static imaging (for both synthetic and real Sgr A* data). For the Themis imaging method (Section 4.3), the parameters of the noise model are fit simultaneously with the image structure, subject to the data-set-specific values in Table 2 being used to define a uniform prior over each data set.

Extending the RML approach from static to dynamic imaging is simple conceptually. Rather than solving for a single image, $\hat{I}$, our new goal is to solve for a series of K images $\hat{\{{I}_{k}\}}$. Each of these images corresponds with small segments of data, which have been divided to have a time duration similar to the expected time variability of the target (typically tens of minutes for Sgr A*). Since the (u, v)-coverage of each data segment is severely limited, we must include an additional term that regularizes the images {I_k } in time rather than just space. A general prescription in terms of the temporal regularizer S_Q can be written mathematically as

$$\begin{eqnarray}\begin{array}{rcl}J(\{{I}_{k}\}) & = & \displaystyle \sum _{\mathrm{data}\ \mathrm{terms}}{\alpha }_{D}{\chi }_{D}^{2}\left(\{{I}_{k}\},V\right)\\ & & -\,\displaystyle \sum _{\mathrm{spatial}\ \mathrm{reg}.}{\beta }_{R}{S}_{R}\left(\{{I}_{k}\}\right)-\displaystyle \sum _{\mathrm{temporal}\ \mathrm{reg}.}{\beta }_{Q}{S}_{Q}\left(\{{I}_{k}\}\right).\end{array}\end{eqnarray} \tag{ 4 }$$

The additional temporal regularization terms, S _Q , encourage smooth evolution of the target over the full observation. Descriptions of temporal regularizers and their application to EHT data are described in Johnson et al. (2017).

In Appendix G the RML dynamic imaging method is used to explore the structure of Sgr A* over the course of a night independent of the variability noise model introduced in Section 3.2.2.

StarWarps is a dynamic imaging method based on a probabilistic graphical model (Bouman et al. 2018). Similar to RML dynamic imaging, StarWarps makes use of temporal regularization to solve for the frames of a movie $\hat{\{{I}_{k}\}}$ over an observation rather than a static image. In contrast to RML, StarWarps independently solves for the marginal posterior of each frame conditioned on all measurements in time; the reconstructed movie $\hat{\{{I}_{k}\}}$ is the mean of each marginal distribution. The advantage of StarWarps with respect to RML is that, when using a linearized measurement model, StarWarps can solve for the frames of a video $\hat{\{{I}_{k}\}}$ with exact inference—resulting in a better-behaved optimization problem that is less likely to get stuck in local minima when compared to RML dynamic imaging.

StarWarps defines a dynamic imaging model for observed data using the following potential functions:

$$\begin{eqnarray}&&{\psi }_{{V}_{k}| {I}_{k}}={{ \mathcal N }}_{{V}_{k}}({f}_{k}({I}_{k}),{\sigma }_{k}^{2})\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\psi }_{{I}_{k}}={{ \mathcal N }}_{{I}_{k}}(\mu ,{\rm{\Lambda }})\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\psi }_{{I}_{k}| {I}_{k-1}}={{ \mathcal N }}_{{I}_{k}}({I}_{k-1},{\beta }_{Q}^{-1}{\mathbb{1}}),\end{eqnarray} \tag{ 7 }$$

for a normal distribution ${ \mathcal N }(m,C)$ with mean m and matrix covariance or scalar standard deviation C. Each set of observed data V_k taken at time k is related to the underlying image, I_k , through the measurement model, f_k (I_k ) (e.g., visibility model, closure phase model). Spatial regularization is imposed through the second potential; I_k is encouraged to be a sample from a multivariate Gaussian distribution with mean μ and covariance Λ. In this work, we define Λ to encourage spatial smoothness with a spectral distribution profile ${({u}^{2}+{v}^{2})}^{-a/2}$ controlled by hyperparameter a, as described in detail in Bouman et al. (2018). The third potential describes how images evolve over time; as β_Q increases, the temporal regularization increases, and vice versa. Although more complex evolution models are described in Bouman et al. (2018), in this paper we restrict ourselves to a simplified evolution model that encourages only small changes between adjacent frames I_k−1 and I_k .

The joint probability distribution of this dynamic model can be written as

$$\begin{eqnarray}&&p(\{{I}_{k}\},\{{V}_{k}\})=p({I}_{1})\prod _{k=1}^{K}p({V}_{k}| {I}_{k})\prod _{k=2}^{K}p({I}_{k}| {I}_{k-1}),\end{eqnarray} \tag{ 8 }$$

where $p({I}_{1})={\psi }_{{I}_{1}}$, $p({V}_{k}| {I}_{k})={\psi }_{{V}_{k}| {I}_{k}}$, and $p({V}_{k}| {I}_{k})\propto {\psi }_{{I}_{k}}{\psi }_{{I}_{k}| {I}_{k-1}}$. In the case of a linear measurement model, f(I) (e.g., complex visibility model), the expected value of every I_k conditioned on all data V = {V_k } can be solved in closed form efficiently using the elimination algorithm. However, in the case of complex gain errors the measurement model is no longer linear. By linearizing the model, we can solve in closed form for a linearized solution, $\hat{\{{I}_{k}\}}$. We then iterate between linearizing the measurement model around our current solution and solving the linearized solution in closed form until convergence.

The StarWarps method is used in Section 9 alongside snapshot geometric modeling methods to help analyze the short-timescale variations of Sgr A* over a ∼100-minute region of time on April 6 and 7.

We use the following three ring models motivated by the morphology identified in many "first images" presented in Appendix C: symmetric and asymmetric ring models (henceforth Ring and Crescent, respectively), and a symmetric ring model with a bright hot spot that rotates in the counterclockwise direction with a period of 30 minutes (henceforth Ring+hs). The first two models are designed to test whether our imaging methods can identify a symmetric versus asymmetric ring, while the latter hot spot model tests the effects of a fast-moving localized emission on the reconstructed images. Besides the ring models, we use four nonring images. To assess the robustness of the central depression seen in ring reconstructions, we adopt a uniform circular and an elliptical disk model (henceforth Simple Disk and Elliptical Disk, respectively). Finally, motivated by the nonring images recovered in Appendix C, we adopt a point-source and double-point-source model (henceforth Point and Double, respectively).

The parameters (e.g., diameter, width) of each geometric model are selected to be broadly consistent with representative properties of Sgr A*'s deblurred visibility amplitudes. We use the following four criteria: (1) the first null traced by Chile–LMT baselines is located at the baseline length of 3.25–3.65 Gλ and position angle (PA) of ∼50°, with an amplitude of ∼0.1 Jy; (2) the peak of the visibility amplitudes between the first two nulls has ∼0.3 Jy; (3) the second null traced by Chile–Hawai'i and/or Chile–PV baselines is located at the (u, v)-distance of 6 Gλ, with an amplitude less than 0.1 Jy; and (4) for asymmetric models, visibility amplitudes on Chile–LMT baselines are ∼1.5 times larger than on SMT–Hawai'i baselines. Figure 7 shows the comparison of visibility amplitudes between Sgr A* data on April 7 and corresponding synthetic data (after adding temporal variability—see Section 5.1.2), demonstrating qualitative agreement between the synthetic data and Sgr A* visibility amplitudes.

To mimic the temporal variability of Sgr A*, the geometric models, denoted by I_geo(x, y), are modulated by a temporal evolution sampled from a statistical model: inoisy (Lee & Gammie 2021). This model enables sampling random spatiotemporal fields, ρ(t, x, y), according to specified local correlations. Lee & Gammie (2021) and Levis et al. (2021) showed that inoisy is able to generate random fields that capture the statistical properties of accretion disks (see Figure 8). Using this model, we modulate the static geometric models according to

$$\begin{eqnarray}&&I(t,x,y)={I}_{\mathrm{geo}}(x,y)\exp \left\{\rho (t,x,y)\right\}.\end{eqnarray} \tag{ 9 }$$

Figure 7 shows both the time average and a single snapshot of each model highlighting the effect of temporal evolution.

Zoom In Zoom Out Reset image size

The inoisy model parameters were selected to generate a similar degree of time variability to that of Sgr A* measurements. We impose the following conditions to match metrics of temporal variations: (1) the mean of each movie's total flux is 2.3 Jy, consistent with the ALMA light curve on April 7 (Section 2; Wielgus et al. 2022); (2) the standard deviation of the total flux is in the range of 0.09–0.28 Jy, as seen in the ALMA and SMA light curves and intrasite baselines; (3) the ${ \mathcal Q }$-metric (Roelofs et al. 2017) of the intrinsic closure phase variability is comparable to Sgr A* data on all triangles.

Figure 9 shows temporal variations in the total flux density and closure phases with comparison to Sgr A* data. The synthetic data variability is able to capture the real data light-curve variablity. Moreover, the power spectrum density distributions of the light curves from synthetic models are broadly consistent with Sgr A* data. For closure phases, inoisy produces data with visible time variability seen in high-S/N triangles, such as the ALMA–SMT–LMT triangle. These synthetic movies also roughly match the Sgr A* variability amplitudes, averaged over all triangles, as evaluated using the ${ \mathcal Q }$-metric. We note that while these synthetic data are in good agreement with the above aspects of the EHT data, their variability amplitudes in Fourier domain are slightly less than Sgr A* data (see Table 2 in Section 3.2).

Zoom In Zoom Out Reset image size

Similar to previous EHT work (M87* Papers IV and VII), we designed three scripted imaging pipelines utilizing the DIFMAP, eht-imaging, and SMILI software packages. After completing the common pre-imaging processing of data (Section 6.1), each pipeline reconstructs images using a broad parameter space (weights for the regularization functions, mask sizes, station gain constraints, variability noise budget parameters, etc.). We describe each pipeline in detail in Appendices D.1, D.2, and D.3.

Each pipeline explored on the order of 10³–10⁵ parameter combinations, as summarized in Tables 3, 4, and 5 for DIFMAP, eht-imaging, and SMILI, respectively. Each pipeline has some unique choices that are fixed (e.g., the pixel size, or the convergence criterion) and surveyed (e.g., the regularizer weights), while some parameters are commonly explored (e.g., parameters for the scattering and intraday variations in Section 6.1).

While all imaging pipelines adopt the common preprocessing of data described in Section 6.1, there are some differences in data processing. For instance, the noise budgets for refractive scattering and intraday variability are updated during self-calibration rounds in SMILI. RML imaging pipelines (eht-imaging and SMILI) adopt the same prior and initial images across all synthetic and real data sets. The DIFMAP pipeline uniformly explores a library of initial models for a first-phase self-calibration, selecting the one that provides the best fit to the closure phases after a first run of cleaning (see Appendix D.1). All three pipelines use combined low- and high-band data for imaging without any data flagging (including the intrasite baselines).

Large imaging surveys on synthetic data facilitate the evaluation of different potential parameter combinations. Following M87* Paper IV, the principal output from each parameter survey is a Top Set: a set of parameter combinations that produce acceptable images on the suite of synthetic data presented in Section 5.

The fidelity of synthetic image reproduction is measured using the normalized cross-correlation between the reconstructed images and ground-truth images. We define the normalized cross-correlation of two images X and Y made of N pixels as

$$\begin{eqnarray}&&{\rho }_{\mathrm{NX}}=\frac{1}{N}\displaystyle \sum _{i}^{N}\frac{({X}_{i}- \langle X \rangle )({Y}_{i}- \langle Y \rangle )}{{\sigma }_{X}{\sigma }_{Y}}.\end{eqnarray} \tag{ 10 }$$

Here X_i and Y_i denote the image intensity at the ith pixel, 〈X〉 and 〈Y〉 denote the mean pixel values of the images, and σ_X and σ_Y are the standard deviations of pixel values. The position offset between the frames is corrected by shifting one frame relative to the other along R.A. and decl. and identifying the shift coordinate that corresponds with the largest ρ_NX.

In order to recover a Top Set, a threshold for ρ_NX is defined for each synthetic data set. Imaging parameter combinations that recover images that score above that threshold for all synthetic data sets are selected as Top Sets. The threshold values of ρ_NX are determined in a manner similar to M87* Paper IV. For each "ground-truth" image (obtained by time-averaging the synthetic movie), we evaluate ρ_NX between the ground-truth image and the same image blurred with a Gaussian beam of FWHM equivalent to the maximum fringe spacing of the Sgr A* observations, 24 μas. This value of ρ_NX quantifies the potential loss of the image fidelity due to the limited angular resolution. Figure 10 shows examples of the ρ_NX curves between unblurred and blurred ground-truth images as a function of the blurring size; the critical value corresponds to those at α = 24 μas. Note that this value depends on the true source structure. Unlike in M87* Paper IV, we find that a relaxation of the ρ_NX threshold is required to account for the fact that we reconstruct static images from time-variable data sets. Hence, the critical ρ_NX values for all training data sets are multiplied by the relaxation factor of 0.95. In other words, the threshold for each synthetic data set is set to 0.95 × ρ_NX for ρ_NX evaluated at α = 24 μas. This relaxed threshold allows for a large enough number of Top Set parameters to be identified for all epochs and imaging pipelines. We ensure that the relaxed threshold still reconstructs all representative ground-truth morphologies; in Appendix E, the worst ρ_NX images are shown for Top Set reconstructions of each model with each imaging pipeline to demonstrate that the representative features are recovered even in the worst-fidelity Top Set images.

Zoom In Zoom Out Reset image size

In Tables 3, 4, and 5, we summarize the parameters and surveyed values in each pipeline's survey. These tables indicate the fraction of images corresponding to that parameter in each pipeline's Top Set for April 7 data. The results on April 6 data are summarized in Appendix E. The tables also provide the total number of surveyed parameter combinations, as well as the number of combinations selected for each Top Set. As seen in each table, there are more than 1000 parameter combinations in each pipeline's Top Set for April 7, which we find is sizable enough for downstream analysis. In contrast, we find that the April 6 Top Set sizes are much smaller, likely due to the poorer (u, v)-coverage.

In Figure 11, we show the time-averaged ground-truth images and corresponding image reconstructions for each of the synthetic data sets with April 7 (u, v)-coverage. Images from all pipelines are obtained with and without scattering mitigation (henceforth descattered and on-sky reconstructions). These descattered and on-sky reconstructions are compared to the time-averaged ground-truth images of the intrinsic and scattered structure, respectively.

Zoom In Zoom Out Reset image size

DIFMAP, eht-imaging, and SMILI images for each geometric model in Figure 11 are obtained using cross-validation: the parameter combination that provides the best mean ρ_NX across other geometric models (i.e., except for the model being tested) is selected. The cross-validation images in Figure 11 (which are contained in each pipeline's Top Set) successfully recover the representative morphology of each geometric model, demonstrating the capability of a single imaging parameter combination to identify various source structures. The manifestation of structure significantly different from the ground-truth morphology is only seen in a small fraction of the cross-validated Top Set parameters. For instance, using the ring classification method described in Appendix F, only 2%, <1%, and <1% of the cross-validated Top Set reconstructions for the Ring model are identified as not having a ring feature for DIFMAP, eht-imaging, and SMILI, respectively; the reconstruction of a nonring morphology is also found to be limited to small fractions of 3%, 5%, and 1% for the Crescent model. Similarly, only a small fraction of the reconstructions from nonring models are reconstructed with a ring morphology—in particular, for the DIFMAP, eht-imaging, and SMILI pipelines, 7%, 11%, and 3% and <1%, <1%, and <1% of the cross-validated Top Sets reconstructed a ring feature from the Point and Double source models, respectively.

For Themis reconstructions, the mean posterior image is shown for each model. Figure 11 shows that the posterior images from Themis reconstructions identify the general morphology of each geometric model. Most of the Themis posterior images satisfy the criteria based on a minimum threshold of ρ_NX used in Top Set selections for DIFMAP, eht-imaging, and SMILI pipelines (see Section 6.2.2), with the exception of a few models discussed below. For April 7 images, all posterior images show higher ρ_NX than the threshold except for 1% and 7% of descattered reconstructions for Double source and Ring+hs models, respectively. For April 6, the images below the threshold are limited to 1% and 22% of descattered reconstructions for the Point and Double source models, respectively. However, 95% of on-sky reconstructions of the Point source model are above the ρ_NX threshold on April 6 for Themis. The high fractions of images beyond the ρ_NX threshold for all models on April 7 demonstrate the capability of the Themis pipeline to recover various representative morphologies at an acceptable fidelity.

Figure 11 also shows the resiliency of the reconstructed morphology to the scattering prescriptions. While the on-sky reconstructions without scattering mitigation tend to be slightly blurrier than those with descattering, there are not many other notable differences in their appearance. In particular, the refractive substructure, which adds spatial distortion of images on scales finer than the angular resolution of the EHT, is not well constrained in our EHT data and therefore does not strongly appear in any reconstructions. We further discuss the effects of scattering prescriptions for Sgr A* images in Section 7.5.2.

We show example GRMHD reconstructions on April 7 in the rightmost panels in Figure 11. This GRMHD simulation contains a ring with a diameter of ∼51 μas. GRMHD images from DIFMAP, eht-imaging, and SMILI pipelines are reconstructed with the Top Set parameter combinations that correspond with the largest average ρ_NX value across all seven geometric models. For Themis reconstructions, on the other hand, the expected (i.e., mean) posterior image is shown. ¹⁵²

Unlike with the simple geometric synthetic data sets, the distribution of GRMHD reconstructions shows wide variations in the image appearance. ¹⁵³ Although GRMHD reconstructions in Figure 11 commonly show a ring-like morphology with a diameter of ∼50 μas, the azimuthal intensity distribution is not consistent across the Top Sets. In the top panels of Figure 12, we show images from GRMHD April 7 data from each pipeline, which are randomly selected from DIFMAP, eht-imaging, and SMILI Top Sets and posterior images from Themis. Most of the images have clear asymmetric ring features, but a few images show nonring structures. The diameters of ring features are broadly consistent across different reconstructions, comparable to the ground-truth image. On the other hand, the azimuthal distributions are not uniquely constrained by the Top Set images—different PAs are seen in the randomly sampled images from each Top Set.

Zoom In Zoom Out Reset image size

To visualize the distribution of images with different morphology, we categorize all images into three major groups: ring images peaked at PAs within (a) the range of −180° ≤ PA < 0°, where the ground-truth value of −124° is located; (b) the range of 0° ≤ PA < 180°; and (c) the remaining images comprising nonring or other ring-like images with much less consistency. The definition of a ring used in this paper to classify ring versus nonring images is described in Appendix F. Note that the particular definition of a ring will influence the quoted percentages of ring and nonring images recovered. We find that the particular definition chosen in this paper results in classification that largely aligns with human perception. However, the classification of images that are borderline between ring and nonring classification is sensitive to the exact criteria used. Therefore, ring definitions that make use of different criteria can lead to classification that still largely aligns with human perception but varies somewhat in the ring classification percentages quoted in this paper.

In the middle and bottom panels of Figure 12, we summarize the distribution of images from each pipeline with and without scattering mitigation, respectively. The GRMHD images within each pipeline's Top Set are clustered into image modes. The upper subpanel shows the mean image of each cluster, indicating a common or representative morphology recovered. The middle subpanel shows the distributions of the azimuthally averaged radial intensity profiles, where the vertical order of images is sorted by the peak radii of the profiles. The lower subpanel shows the azimuthal distribution of the radial peak intensity within a radius of 10–40 μas.

Figure 12 demonstrates that most of the Top Set or posterior images reconstruct ring images from the GRMHD data set. In particular, the radial profiles of the first two ring modes show a broad consistency of the peak radius around ∼25 μas, which implies a diameter of ∼50 μas consistent with that of ∼50 μas for the ground-truth model. The capability of identifying a ring with the consistent diameter does not appear to depend on the scattering mitigation. Similarly to geometric synthetic data (see Section 6.4.1), the scattering mitigation does not significantly affect the resulting source morphology in the reconstructed images, except that the images without scattering mitigation tend to be slightly broader. Refractive substructure does not appear in the reconstructed on-sky images, again indicating that the appearance of the refractive substructure is not strongly constrained with EHT 2017 data.

As seen in the azimuthal profiles, each pipeline provides at least one asymmetric ring mode with the peak PA roughly consistent with that of the mean ground-truth image of −124°. For this particular GRMHD example, this mode is found to be the most popular mode in most pipelines. However, Figure 12 indicates that the ring mode with the correct orientation is not always identified as the most popular mode among image samples—for instance, the correct orientation is not identified as the most popular ring mode for the eht-imaging on-sky pipeline. Therefore, caution should be taken, as the popularity of a mode in an RML or the CLEAN pipeline's Top Set does not necessarily always correspond with the true underlying structure.

We find that key takeaways from the GRMHD example are consistent with the "best-bet" GRMHD models presented in Paper V, identified based on various criteria using both EHT and non-EHT data. In Appendix H, we show example reconstructions of a "best-bet" GRMHD model. We find the identification of a ring morphology for the vast majority of the Top Set reconstructions; however, multiple ring modes and nonring images are still reconstructed. Note also that the peak PA of the ground-truth "best-bet" GRMHD model is not necessarily identified as the most popular mode reconstructed in each Top Set. These results indicate that the same trends are seen across multiple GRMHD models that are broadly consistent with EHT data.

The first possible explanation for the small percentage of nonring reconstructions of Sgr A* is simply that Sgr A* does not possess a ring morphology. Although this possibility cannot be ruled out completely, we explore the possibility of a potential bias in our imaging approach toward recovering ring images from underlying nonring sources.

There is always a possibility that the parameters initially explored by the RML and CLEAN parameter surveys were inadvertently biased to produce mostly ring images. This would result in an artificially high percentage of ring reconstructions that could give overconfidence in an incorrect ring solution. To explore this hypothesis, we inspect the percentage of ring and nonring images of Sgr A* that are in the initial parameter survey versus the restricted Top Set on April 7. We find that although ring images make up a large percentage of the initial parameter surveys, 91%, 84%, and 60% for DIFMAP, eht-imaging, and SMILI, respectively, for the descattered reconstructions, this number becomes significantly larger after Top Set selection. In particular, the percentage of ring images rises to values of 95%, 97%, and 98%, respectively. Therefore, when we restrict to those parameters that can best disambiguate between the different ring and nonring source morphologies contained in our synthetic data sets (Section 5.1), the number of ring images increases. We also note that Themis posterior samples for April 7 only contain ring images.

To further test the possibility that we are biased to recover primarily rings from underlying nonring sources, we have explored the performance of our imaging methods on nonring synthetic data sets to see how they compare to Sgr A* results. In particular, we performed a Top Set analysis on the Point and Double source models. Two "cross-validation" Top Sets were identified by excluding the performance of either the Point or Double source model and only using the remaining six geometric models in selection of the imaging parameters. These two cross-validation Top Set parameter sets were then applied to the Point and Double synthetic data sets. Although each cross-validation Top Set does incorrectly reconstruct some ring images, each image set primarily contains reconstructions that look similar to the ground-truth Point and Double source structure. In particular, for descattered reconstructions, only 9% and <1% of the cross-validation Top Set reconstructions possess a ring morphology, compared to 97% in Sgr A* reconstructions of April 7. The Themis pipeline produces no ring images in the posterior samples for both the Point and Double data sets. Thus, in the nonring synthetic data sets we have explored we find that the number of incorrect ring reconstructions is much less than the number that we recover for Sgr A*. In summary, our imaging pipelines do not appear prone to artificially create a majority of ring structures in sources that do not possess an intrinsic ring morphology.

The second possible explanation for the nonring reconstructions of Sgr A* is that interstellar scattering causes a nonring morphology. This could take one of two forms: (a) Sgr A*'s ring morphology has been distorted to a nonring morphology by an interstellar scattering screen, or (b) our imaging algorithms have incorporated an incorrect scattering model that reconstructs a corrupted nonring morphology.

To address the possibility that interstellar scattering is causing Sgr A*'s intrinsic structure to be distorted to a nonring, we inspect differences in the descattered versus on-sky reconstructions. As outlined in Section 3.1, descattered reconstructions attempt to mitigate two primary effects of interstellar scattering: diffractive scattering that causes angular broadening, and refractive scattering that introduces small-scale structure to the on-sky image. Sgr A* reconstructions in Figure 14 show that on-sky images are systematically blurrier than descattered images, as expected owing to deblurring that is performed before recovering a descattered image (refer to Section 3.1.1). This systematic difference between on-sky and descattered images is also seen in synthetic data reconstructions presented in Section 6.4. This angular broadening causes the central dip of a ring to be less prominent. Nonetheless, we note that the vast majority of the on-sky (i.e., without any scattering mitigation prescription) images are still rings, with percentages of 93%, 98%, 90%, and 98% for CLEAN, eht-imaging, SMILI, and Themis, respectively. We note that it is highly unlikely that a nonring morphology was distorted into an on-sky ring morphology by interstellar scattering.

To address the possibility that our assumed model of the interstellar scattering is reconstructing a corrupted image of the descattered source, we explore the use of multiple refractive noise models to mitigate effects of refractive substructure before deblurring (Section 6.2). The contribution of the included refractive noise budget used to mitigate scattering (for descattered reconstructions) is nonnegligible for long-baseline data (see Section 3.1.2 and Figure 5). Nontheless, the RML and CLEAN imaging surveys show that the choice of the refractive noise model does not significantly affect the resulting distributions of image structures. Figure 20 shows the average descattered Sgr A* images (clustered into the same four morphologies as are presented in Figure 13) for each refractive noise prescription explored. The comparable fractions of images in each cluster across the different descattering strategies indicate the resiliency of the recovered image structures to different refractive noise models. Another piece of evidence suggesting that the scattering prescription does not strongly affect results can be seen in the Themis results; in Themis a constant refractive noise floor is able to vary in posterior sampling of descattered images. Although an essentially unlimited refractive noise component is permitted by the Themis model, the posterior estimation instead typically chooses a noise level of <25 mJy, only 1% of the total flux of the source.

Zoom In Zoom Out Reset image size

In summary, results indicate that for EHT measurements of Sgr A*, interstellar scattering does not significantly affect the recovered morphology. Both on-sky and descattered images contain a majority of ring morphologies. In addition, we find that the particular choice of scattering mitigation strategy only minimally affects the recovered image structures.

The third possible origin for nonring Sgr A* images is poor reconstruction quality of the imaging methods. Imaging is solving an ill-posed inverse problem due to sparse (u, v)-coverage, which always will have the possibility of recovering an incorrect image. This is further exasperated in imaging Sgr A* by the challenges that come with recovering an evolving source. In this section we explore this possibility and find that our imaging methods often reconstruct nonring sources for variable ring sources with a comparable small percentage as found for Sgr A*.

The first natural question is how our imaging methods perform on reconstructing ring sources with similar data properties to Sgr A*. We find that our methods produce nonring images, even when the underlying source structure is ring-like. As an example, recall that for the variable GRMHD ring-like sources analyzed in Section 6 and Appendix H our imaging methods produced mostly rings, but also a very small fraction of nonring images. Therefore, we expect that for a variable ring source we would recover some nonrings that fit the data well. However, note that the number of nonrings was still fairly small in the case of the GRMHDs: 5% for the GRMHD presented in Figure 12, and 4% for the "best-bet" GRMHD presented in Figure 37. These values are comparable to the 3% of nonrings reconstructed for Sgr A* descattered images across all pipelines. We also find that cross-validated Top Set images reconstructed of the Crescent and Ring geometric sources produce a small fraction of nonring images. Note that these nonring percentages for the GRMHD, Crescent, and Ring sources are much less than the percentage of nonrings reconstructed for the variable Point and Double sources, as discussed in Section 7.5.1.

Due to Sgr A*'s intraday evolution, mitigation of temporal variability in the data is important for reconstructing a single static image of Sgr A*. We next explore how the variability mitigation approach affects the proportion of ring versus nonring images reconstructed. Section 3.2 introduced an approach to model the temporal variability as an additional noise budget that could be added in quadrature to the visibility thermal noise budget. The dependence of the temporal variability noise model parameters on resulting Sgr A* reconstructions for the RML and CLEAN pipelines was investigated. Although a variability noise budget generally helps with imaging (as evidenced by a higher percentage of Top Set parameters selected—89%–95% of the Top Set parameter combinations include a variability noise model on April 7), similar to scattering, we find that there are no significant differences between the images recovered under different variability noise parameters. We suspect that this is partly caused by the complex interplay between regularizers and different noise model parameters (e.g., variability, scattering, and systematic), although no significant trends were identified.

In Appendix G, we demonstrate how our results are insensitive to a different method of variability mitigation. In particular, we show time-averaged morphologies identified by a large survey over full-track RML dynamical imaging parameters (refer to Section 4.4.1) that do not rely on the variability noise model presented in Section 3.2.2. We again find the ring modes to be dominantly reconstructed, while a small fraction of nonring structures are also identified. The broad consistency indicates that our results are resilient to at least two different methods to recover time-averaged morphology.

In summary, we find that our imaging methods do reconstruct a small percentage of nonring images from ring sources with comparable variability to that seen in Sgr A* data. Similar behavior is observed across a variety of different imaging approaches to mitigate variability. Therefore, we conclude that we would expect our methods to produce a small fraction of nonring images from an underlying variable ring morphology.