First Sagittarius A* Event Horizon Telescope Results. IV. Variability, Morphology, and Black Hole Mass

As a radio interferometer, the EHT is natively sensitive to the Fourier transform of the sky-plane emission structure. For a source of emission I( x , t), the complex visibility ${ \mathcal V }({\boldsymbol{u}},t)$ is given by

$$\begin{eqnarray}&&{ \mathcal V }({\boldsymbol{u}},t)=\iint {e}^{-2\pi i{\boldsymbol{u}}\cdot {\boldsymbol{x}}}I({\boldsymbol{x}},t){d}^{2}{\boldsymbol{x}},\end{eqnarray} \tag{ 1 }$$

where t is time, x = (x, y) are angular coordinates on the sky, and u = (u, v) are projected baseline coordinates in units of the observing wavelength (see, e.g., Thompson et al. 2017).

The ideal visibilities ${ \mathcal V }$ are not directly observable because they are corrupted by both statistical errors and a variety of systematic effects. For the EHT, the dominant systematics are complex station-based gain corruptions. The relationship between an ideal visibility ${{ \mathcal V }}_{{ij}}$ and the observed visibility V_ij on a baseline connecting stations i and j is given by

$$\begin{eqnarray}&&{V}_{{ij}}={g}_{i}{g}_{j}^{* }{{ \mathcal V }}_{{ij}}+{\sigma }_{\mathrm{th},{ij}}\equiv | {V}_{{ij}}| {e}^{i{\phi }_{{ij}}},\end{eqnarray} \tag{ 2 }$$

where σ_th,ij is the statistical (or "thermal") error on the baseline, g_i and g_j are the station gains, and we have defined the visibility amplitude ∣V_ij ∣ and phase ϕ_ij . The statistical error is well described as a zero-mean circularly symmetric complex Gaussian random variable with a variance determined (per the radiometer equation) by the station sensitivities, integration time, and frequency bandwidth (Thompson et al. 2017). The station gains vary in time at every site and must in general be either calibrated out or determined alongside the source structure.

The presence of station-based systematics motivates the construction and use of "closure quantities" that are invariant to such corruptions. A closure phase ψ_ijk (Jennison 1958) is the sum of visibility phases around a closed triangle of baselines connecting stations i, j, and k,

$$\begin{eqnarray}&&{\psi }_{{ijk}}={\phi }_{{ij}}+{\phi }_{{jk}}+{\phi }_{{ki}}.\end{eqnarray} \tag{ 3 }$$

Closure phases are invariant to station-based phase corruptions, such that the measured closure phase is equal to the ideal closure phase, up to statistical errors. Similarly, a closure amplitude A_ijkℓ (Twiss et al. 1960) is the ratio of pairs of visibility amplitudes on a closed quadrangle of baselines connecting stations i, j, k, and ℓ,

$$\begin{eqnarray}&&{A}_{{ijk}{\ell }}=\displaystyle \frac{| {V}_{{ij}}| | {V}_{k{\ell }}| }{| {V}_{{ik}}| | {V}_{j{\ell }}| }.\end{eqnarray} \tag{ 4 }$$

Analogous with closure phases, closure amplitudes are invariant to station-based amplitude corruptions. Because closure quantities are constructed from nonlinear combinations of complex visibilities, they have correlated and non-Gaussian error statistics; a detailed discussion is provided in Blackburn et al. (2020).

The EHT observed Sgr A* on 2017 April 5, 6, 7, 10, and 11 with the phased Atacama Large Millimeter/submillimeter Array (ALMA) and the Atacama Pathfinder Experiment (APEX) on the Llano de Chajnantor in Chile, the Large Millimeter Telescope Alfonso Serrano (LMT) on Volcán Sierra Negra in Mexico, the James Clerk Maxwell Telescope (JCMT) and phased Submillimeter Array (SMA) on Maunakea in Hawai'i, the IRAM 30 m telescope (PV) on Pico Veleta in Spain, the Submillimeter Telescope (SMT) on Mt. Graham in Arizona, and the South Pole Telescope (SPT) in Antarctica (M87* Paper II). Only the April 6, 7, and 11 observations included the highly sensitive ALMA station, and the April 11 light curve exhibits strong variability (Wielgus et al. 2022) that is presumably associated with an X-ray flare that occurred shortly before the start of the track (Paper II). In this paper, we thus analyze primarily the April 6 and April 7 data sets. We note that while Paper III focuses on the April 7 data set, with the April 6 data set used for secondary validation, most of the analyses carried out in this paper instead focus on a joint data set that combines the April 6 and April 7 data.

At each site the data were recorded in two 1.875 GHz wide frequency bands, centered around sky frequencies of 227.1 GHz (low band; LO) and 229.1 GHz (high band; HI), and in each of two polarization modes. For all telescopes except ALMA and JCMT, the data were recorded in a dual circular polarization mode: right-hand circular polarization (RCP; R) and left-hand circular polarization (LCP; L). ALMA recorded using linear feeds, and the data were later converted to a circular polarization basis during the DiFX (Deller et al. 2011) correlation (Martí-Vidal et al. 2016; Matthews et al. 2018; Goddi et al. 2019). The JCMT observed only a single hand of circular polarization at a time, with the specific handedness (RCP or LCP) changing from day to day. All other stations observed in a standard dual-polarization mode, which allows the construction of RR, RL, LR, and LL correlation products. The analyses in this paper use only the parallel-hand correlations (i.e., RR and LL), which are averaged to form Stokes I data products. Because JCMT records only a single hand at a time, we instead form "pseudo-I" data products for JCMT baselines, using whichever parallel-hand correlation is available as a stand-in for Stokes I. ¹⁵⁰

After correlation, residual phase and bandpass errors are corrected with two independent processing pipelines: EHT-HOPS (Blackburn et al. 2019) producing "HOPS" (Whitney et al. 2004) data and rPICARD (Janssen et al. 2018, 2019) producing "CASA" (McMullin et al. 2007) data. Relative phase gains between RCP and LCP have been corrected based on the assumption of zero circular polarization on baselines between ALMA and other EHT stations. Absolute flux density scales are based on a priori measurements of each station's sensitivity, resulting in a ∼10% typical uncertainty in the amplitude gains (M87* Paper II). The amplitude gains of the colocated ALMA/APEX and SMA/JCMT stations have been further refined via time-variable network calibration (M87* Paper III) using a light curve of the compact Sgr A* flux measured by ALMA and SMA (Wielgus et al. 2022). For the remaining stations, gross amplitude gain errors have been corrected by a transfer of gain solutions from the J1924–2914 and NRAO 530 calibrator sources as described in Paper II.

Following the completion of the above calibration pipelines, additional preprocessing of the data has been carried out as described in Paper III, including calibration of the LMT and JCMT station gains and normalization of the visibility amplitudes by the total light curve. The characterization of residual calibration effects (e.g., polarization leakage) into a systematic error budget, as well as a more comprehensive description of the overall EHT Sgr A* data reduction, is provided in Paper II.

In low-luminosity SMBH systems such as Sgr A*, we expect the emission to originate from the immediate vicinity of the black hole, i.e., on scales comparable to the event horizon size. Here, all characteristic speeds of the hot relativistic gas approach the speed of light. The timescales associated with these processes are therefore set by the gravitational timescale, GM/c³, which is ∼20 s for Sgr A*. This timescale is ∼3 orders of magnitude shorter than the nightly observations carried out by the EHT, so a single observation contains many realizations of the underlying source variability.

GRMHD simulations can model the dynamical processes in Sgr A* and, using ray-tracing and radiative transfer, provide a theoretical expectation for the observed emission. Paper V provides a library of GRMHD simulations and associated movies, which have been scaled to the conditions during the EHT 2017 observations (e.g., the average total 230 GHz flux is set to the EHT measurement). We use the variability characteristics of these simulations as our expectation for the Sgr A* variability seen by the EHT.

GRMHD simulations are universally described by a "red–red" power spectrum, with the largest fluctuations in the emission occurring on the longest timescales and the largest spatial scales (Georgiev et al. 2022). Spatially, the largest scale for variability is limited to the size of the emitting region, which for an observing frequency of 230 GHz is typically several GM/c² and for the EHT Sgr A* data is constrained to be ≲87 μas (Paper II). Temporally, the simulations exhibit a red power spectrum that flattens on timescales ≳ 1000 GM/c³. Observations of the total flux variability in Sgr A* corroborate this expectation, finding a red-noise spectrum extending to timescales of several hours and flattening on longer timescales (Wielgus et al. 2022).

We can, without loss of generality, express the time-variable image structure I in terms of some static mean image I_avg and a zero-mean time-variable component δ I that captures all of the variation,

$$\begin{eqnarray}&&I({\boldsymbol{x}},t)={I}_{\mathrm{avg}}({\boldsymbol{x}})+\delta I({\boldsymbol{x}},t).\end{eqnarray} \tag{ 5 }$$

The linearity of the Fourier transform ensures that an analogous decomposition holds for ${ \mathcal V }$, which is thus simply the sum of an analogous ${{ \mathcal V }}_{0}$ and $\delta { \mathcal V }$. The variation $\delta { \mathcal V }$ represents the component of the data we wish to mitigate.

The EHT stations ALMA and SMA are themselves interferometric arrays capable of separating out extended structure (such as the Galactic center "minispiral"; Lo & Claussen 1983; Goddi et al. 2021; Wielgus et al. 2022) from the Sgr A* light curve,

$$\begin{eqnarray}&&L(t)={ \mathcal V }({\bf{0}},t)=\iint I({\boldsymbol{x}},t){d}^{2}{\boldsymbol{x}},\end{eqnarray} \tag{ 6 }$$

on the largest spatial scales, predicted by GRMHD simulations to be the most variable. Using this motivation, the light-curve-normalized image is defined to be

$$\begin{eqnarray}&&\hat{I}({\boldsymbol{x}},t)\equiv \displaystyle \frac{I({\boldsymbol{x}},t)}{L(t)},\end{eqnarray} \tag{ 7 }$$

with ${\hat{I}}_{\mathrm{avg}}$ and $\delta \hat{I}$ similarly defined; here, the "hat" diacritic denotes light-curve normalization. From GRMHD simulations, the expected noise is well approximated by a broken power law,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{var}}^{2}\equiv \langle \delta {\hat{{ \mathcal V }}}^{2}\rangle \approx \displaystyle \frac{{a}^{2}{\left(| {\boldsymbol{u}}| /{u}_{0}\right)}^{c}}{1+{\left(| {\boldsymbol{u}}| /{u}_{0}\right)}^{b+c}},\end{eqnarray} \tag{ 8 }$$

along any radial direction (Georgiev et al. 2022). This broken power law is described by four parameters: a break at u₀, an amplitude a representing the amount of noise at the break location, and long- and short-baseline power-law indices b and c, respectively. Typically, we expect that c ≳ 2, due to the compact nature of the source.

In Figure 1, red lines show ${\sigma }_{\mathrm{var}}^{2}$ measured for an example GRMHD simulation about average images that have been constructed on observationally relevant timescales. The variability has been averaged in azimuth and across different black hole spin orientations. As the timescale over which the average image is constructed increases, the location of the break u₀ decreases and the amount of power at the break increases. ¹⁵¹ This behavior can intuitively be understood as the GRMHD simulations changing less for short timescales. For comparison, we show the thermal, systematic, and refractive scattering noise. For timescales longer than ∼10 minutes, the variability noise dominates on EHT VLBI baselines.

Zoom In Zoom Out Reset image size

The intraday variability expected from theoretical considerations can be observed directly in the Sgr A* data. Figure 2 shows the combined baseline coverage for the EHT's 2017 Sgr A* campaign, including the observations on April 5, 6, 7, and 10. The upper limit on the source size of 87 μas (see the second-moment analysis in Paper II) implies that the complex visibilities will be correlated in regions of the (u,v)-plane smaller than ∼ 2 Gλ. In practice, the visibility amplitudes exhibit variations on scales smaller than this and otherwise appear strongly correlated on scales of 1 Gλ (see Paper III, Figure 3). Therefore, among the baseline tracks in Figure 2 there are four regions where the (u,v)-coverage is redundant, i.e., multiple baselines pass within 1 Gλ of the same (u,v)-position. We separate the redundant baseline combinations into "crossing tracks," in which two baseline tracks intersect at a single (u,v)-point, and "following tracks," in which two baselines follow a nearly identical extended track in the (u,v)-plane. Both sets of redundant baselines provide an opportunity to directly probe the degree of intraday variability in the visibilities at specific locations in the (u,v)-plane.

Zoom In Zoom Out Reset image size

Prior to making comparisons, we apply the data preprocessing steps outlined in Section 2.3 to mitigate unphysical sources of variability. To avoid addressing the unknown atmospheric phase delays, we focus exclusively on visibility amplitudes. Because source structure will produce additional variations in the visibility amplitudes that are hard to visualize in projection and obscure the relative degree of variability, we detrend the visibility amplitudes with a linear model. The crossing and following tracks discussed below are shown in the top and bottom subpanels of Figure 2, respectively.

Chile–PV versus Chile–SPT: The first crossing track we consider contains baselines between the Chile stations (ALMA, APEX) and PV and SPT, which both cross near (u,v) = (4 Gλ, 3.5 Gλ) at times separated by 6.2 hr. The concurrent ALMA and APEX baselines are consistent within the reported statistical errors, and thus there is no evidence for unaddressed baseline-specific dominant systematic errors. The normalized visibility amplitudes for the Chile–PV and Chile–SPT baselines individually vary smoothly with time. Nevertheless, they differ significantly at the crossing point, and this difference is consistent in magnitude with the variation found across days (indicated by the gray band in the relevant panel of Figure 2).

Chile–SMT versus Chile–SPT: The second crossing track we consider contains baselines between the Chile stations (ALMA, APEX) and SMT and SPT, which both cross near (u, v) =(3 Gλ, 4.5 Gλ) at times separated by 5.2 hr. Again, we find excellent agreement between ALMA and APEX baselines, individually smooth variations on the Chile–SMT and Chile–SPT baselines, and significant differences in the visibility amplitudes between those baselines.

SMA–SPT versus LMT–SPT: The first following track we consider contains baselines between the SPT, which is located at the South Pole, and SMA and LMT, which have similar latitudes. Because the baseline tracks are coincident across a large range of locations in the (u,v)-plane, this following track permits many direct comparisons at a baseline length of 8 Gλ at times separated by 3.4 hr. As with both crossing tracks, significant differences exist between the two sets of baselines, consistent with the range across multiple days.

SMT–SPT versus PV–SPT: The final following track we consider again involves the SPT, and now the SMT and PV, which also have similar latitudes. This is the longest set of baselines that we consider, with a length of roughly 8.5 Gλ and covering similar regions in the (u,v)-plane at times separated by 6.7 hr. Again, significant variations are exhibited, consistent with those across days.

In summary, intraday variability is observed on multiple baselines with lengths ranging from 5 to 8.5 Gλ and on timescales as short as 3.4 hr. In all cases, this variability is broadly consistent with that observed on interday timescales. Furthermore, the variability behavior is consistent with theoretical expectations from GRMHD simulations and empirical expectations from the Sgr A* light curve, both of which imply that the variable elements of the Sgr A* emission should be uncorrelated beyond a timescale of a few hours (Georgiev et al. 2022; Wielgus et al. 2022). Any average image of Sgr A* reconstructed from data spanning a time range longer than several hours captures the long-timescale asymptotic source structure; the intrinsic image averaged over a single day or multiple days is thus expected to exhibit similar structure.

To quantify the variability observed in the EHT Sgr A* data, we make use of the procedure described in Broderick et al. (2022). This procedure provides an estimate of the excess variability—i.e., the visibility amplitude variance in excess of that caused by known sources, such as average source structure, statistical and systematic uncertainties, and scattering—as a function of baseline length. We apply the same data preparation steps summarized in Section 3.2 and described in Broderick et al. (2022), combining the visibility amplitudes measured on April 5, 6, 7, and 10 in both observing bands. All data points are weighted equally.

The procedure is illustrated in Figure 3. We again make use of the strong correlations induced by the finite source size, and for every location in the (u,v)-plane we consider only those data points falling within a circular region of diameter 1 Gλ centered at that point (red circular region in the top panel of Figure 3). Within each such region containing at least three data points, we linearly detrend the light-curve-normalized visibility amplitudes with respect to u and v to remove variations due to physical structure (bottom panel of Figure 3), and we compute the variance of the residuals. This variance is then debiased to remove the contributions from the reported statistical errors, as described in Broderick et al. (2022). Finally, the variances from all regions having a common baseline length are averaged to produce an azimuthally averaged set of variances. The uncertainty in the variance estimates is obtained via Monte Carlo sampling of the unknown gains, leakage terms, and statistical errors.

Zoom In Zoom Out Reset image size

Figure 4 shows the results of applying this procedure to the Sgr A* data, with the normalized visibility amplitude variance measurements given by the black points. For baselines shorter than 2.5 Gλ, the LMT calibration procedure precludes an accurate estimate of the variance, and thus these baselines have been excluded. For baselines between ∼2.5 and 6 Gλ in length, our empirical estimates of the noise exceed the typical contributions from statistical errors and refractive scattering, indicating the presence of an additional source of structural variability. The degree of inferred variability is consistent with that seen in prior millimeter-VLBI data sets, which is discussed further in Appendix A. For baselines longer than 6 Gλ, our measurements are consistent with the degree of variability expected from the statistical uncertainties in the data; we thus do not directly constrain the source variability on these long baselines.

Zoom In Zoom Out Reset image size

To characterize the variability behavior within the (u,v)-plane, we fit a broken power law of the form in Equation (8) to the normalized variance measurements. As indicated by the filled black circles in Figure 4, significant measurements exist only in the range of baselines with lengths ∼2.5–6 Gλ; on baselines longer than 6 Gλ, we are unable to distinguish the variability from its associated measurement uncertainties. We thus perform the broken power-law fit only to the ∼2.5–6 Gλ range of baselines, where we have significant measurements, and we find no evidence for a break in the power law in this region. As a result, only an upper limit can be placed on u₀, and we are not able to constrain the short-baseline power-law index, c. The range of permitted broken power-law fits is illustrated in Figure 4 by the orange shaded region, with several samples from the posterior distribution explicitly plotted as orange lines.

Because the location of the broken power-law break is poorly constrained, the parameters u₀ and a (describing the location of the break and the amplitude of the power law at the break, respectively) are strongly correlated and highly uncertain. However, it is clear from the orange shaded region in Figure 4 that there is only a narrow range of variances permitted over the 2.5–6 Gλ range of baselines over which the data are constraining. We thus choose to characterize the amplitude of the excess variability noise at ∣ u ∣ = 4 Gλ, which we denote as a₄. Joint posteriors for a₄, the break location u₀, and the long-baseline power-law index b are shown in Figure 5. These constraints are used to inform the prior distributions for the full-track geometric modeling described in Section 7; the associated prior ranges on each parameter are indicated by the purple shaded regions in Figure 5.

Zoom In Zoom Out Reset image size

Having established the existence of structural variability and quantified its magnitude in the Sgr A* data, we now turn to strategies for mitigating its impact on downstream analyses. We employ the light-curve-normalized visibility data, which eliminates large-scale variations and correlations by construction. In principle, there are four methods that we might pursue to address the remaining structural variability:

• 1.  
  Analyze time-averaged data products.
• 2.  
  Employ explicitly time-variable models.
• 3.  
  Analyze short time segments of the data and combine the results afterward to characterize the average source structure.
• 4.  
  Simultaneously reconstruct the average source structure and a statistical characterization of the structural variability.

The first of these options is complicated substantially by the uncertain visibility phases, which limit our ability to coherently average the data on timescales longer than several minutes. The second option can be employed either when a descriptive low-dimensional model for the source structure can be constructed (e.g., Miller-Jones et al. 2019; Kim et al. 2020) or when there is sufficient (u,v)-coverage for nonparametric dynamical imaging algorithms to be successful (e.g., Johnson et al. 2017; Bouman et al. 2018; Arras et al. 2022). The latter approach is explored in the dynamical imaging analyses described in Paper III, ultimately demonstrating that the Sgr A* (u,v)-coverage is insufficient to permit unambiguous reconstructions of the variable source structure.

We dub the third option "snapshot" modeling, whereby a simple geometric model of the source structure is fit to segments of the data that are short enough in duration ( ≲ 3 minutes) for the impact of structural variability to be subdominant to other sources of visibility uncertainty (e.g., refractive scattering; see Figure 1). Though the data sparsity is exacerbated by restricting the reconstructions to only a single snapshot at a time, the model itself is also correspondingly restricted in its parameterization of the source structure. The results of the fits to each individual snapshot are then combined across the entire data set, effectively averaging over the source variability. Details of our snapshot modeling analyses as applied to Sgr A* are presented in Section 6.

The fourth option we refer to as "full-track" modeling, which aims to simultaneously reconstruct both the average source structure and a set of parameters describing the contribution of the structural variability to the visibility data variances (Broderick et al. 2022). In contrast to the snapshot modeling, full-track modeling considers the entire data set at once and uses a parameterized "variability noise" model to appropriately modify the data uncertainties as part of the fitting procedure. In this way, the full-track modeling retains access to sufficient (u,v)-coverage to permit fitting a nonparametric image model to the data (see Paper III), though in Section 7 we also pursue full-track geometric modeling to provide a cross-comparison with the results from the snapshot geometric modeling. Our parameterization of the variability noise follows Equation (8), with the amplitude specified at a baseline length of 4 Gλ as described in Section 3.3. A detailed description of our full-track modeling approach as applied to Sgr A* is presented in Section 7.

Both the snapshot and full-track modeling approaches focus on describing the average source structure and treating the structural variability in a statistical manner. This goal is formally mismatched with what the EHT data measure for a single day, which is instead a collection of complex visibilities that sample different instantaneous realizations of the intrinsic Sgr A* source structure. The nature of this mismatch impacts the full-track analyses significantly.

The variability mitigation scheme employed by full-track modeling presumes that the variability may be modeled as excess uncorrelated fluctuations in the complex visibility data. This assumption is well justified on timescales exceeding a few hours, but significant correlations between visibilities exist on shorter timescales. Within a single day, subhour correlations that are localized in the (u,v)-plane can induce significant biases in the source structure reconstructed from the sparse EHT (u,v)-coverage of Sgr A*. The noise model is thus fundamentally misspecified for EHT data, with the level of misspecification increasing as shorter-in-time segments of data are analyzed; Appendix B describes pathological behavior that can arise when analyzing EHT Sgr A* data from only a single day. While prominent artifacts associated with these subhour correlations are present in the April 6 reconstructions shown in Appendices B and C, we note that the underlying origin of these artifacts is no less present on April 7.

The impact of unmodeled correlations on the reconstructed source structure can be ameliorated by combining multiple days, which provides visibility samples associated with independent realizations of the source structure. This additional sampling rapidly brings the statistical properties of the data into better agreement with the assumptions underpinning the full-track analyses; even the combination of just 2 days is often sufficient to mitigate the subhour correlations in analysis experiments that make use of GRMHD simulation data. For this reason, we combine both the April 6 and April 7 data sets during the analysis of the Sgr A* data. For comparison, Appendix C presents the results of equivalent analyses applied to the April 6 and April 7 data sets individually.

In reconstructing images of Sgr A*, Paper III explores a large space of imaging algorithms and associated assumptions. The resulting "top sets" of images contain primarily "ring-like" image structures, though a small fraction of the images are morphologically ambiguous. These "nonring" images still nominally provide a reasonable fit to the Sgr A* data and so are not ruled out from the Paper III results.

We can quantify the preference for a ring-like image structure by fitting the data with a set of simple geometrical models. Employing the snapshot geometric modeling technique detailed in Section 6, we compare the Bayesian evidence, ${ \mathcal Z }$, between these different geometric models. The value of ${ \mathcal Z }$ serves as a model comparison metric that naturally balances improvements in fit quality against increases in model complexity, with larger values of ${ \mathcal Z }$ indicating preferred models (see, e.g., Trotta 2008). Figure 6 shows the results of a survey over simple geometric models with varying complexity, captured here by the number of parameters required to specify the model. At all levels of complexity, ring-like models outperform the other tested models. This disparity is most stark for the simplest models but continues to hold as the models increase in complexity.

Zoom In Zoom Out Reset image size

The remainder of this paper proceeds with analyses that presuppose a ring-like emission structure for Sgr A*.

The overall structure of the Sgr A* visibility amplitudes (see the left panel of Figure 7) exhibits at least three distinct regions:

• 1.  
  A "short-baseline" region containing baselines shorter than ∼2 Gλ. The effects of data calibration and preprocessing—particularly the light-curve normalization and LMT calibration procedures (Paper III)—are evident in the unit total flux density and the Gaussian structure of the visibility amplitudes in this region.
• 2.  
  An "intermediate-baseline" region containing baselines between ∼2 and 6 Gλ. The visibility amplitudes in this region exhibit a general rise and then fall with increasing baseline length, peaking at a flux density ∼20% of the total at a baseline length of ∼4 Gλ.
• 3.  
  A "long-baseline" region containing baselines with lengths in excess of ∼6 Gλ. The visibility amplitudes in this region generally rise with increasing baseline length from a deep minimum near ∼6.5 Gλ, approximately flattening out at longer baselines to a level that is ∼3%–10% of the total flux density.

Zoom In Zoom Out Reset image size

The visibility amplitudes exhibit indications of asymmetric source structure, particularly on baselines with lengths of ∼3 Gλ that fall near the first minimum. Here, the baselines between the SMT and Hawai'i stations (oriented approximately in the east–west direction) have systematically higher correlated flux densities than the similar-length baselines between the LMT and Chile stations (oriented approximately in the north–south direction). The implication for the source morphology is that we would expect to see more symmetric structure in the north–south than in the east–west direction. Detailed geometric modeling analyses that are able to capture this asymmetry are described in Sections 6 and 7; here, we consider only a simple azimuthally symmetric toy model that captures some salient features of interest.

We attempt to understand the visibility behavior in light of expectations for a ring-like emitting structure. Specifically, we consider a geometric construction whereby an infinitesimally thin circular ring bordering an inner disk of emission is convolved with a Gaussian blurring kernel. The visibility function V produced by such an emission structure is given by

$$\begin{eqnarray}\begin{array}{rcl}V & = & {F}_{0}{V}_{\mathrm{Gauss}}\left[{f}_{d}{V}_{\mathrm{disk}}+(1-{f}_{d}){V}_{\mathrm{ring}}\right]\\ & = & {F}_{0}\exp \left(-\displaystyle \frac{{w}^{2}{\xi }^{2}}{4\mathrm{ln}(2)}\right)\left[\displaystyle \frac{2{f}_{d}}{\xi }{J}_{1}(\xi )+(1-{f}_{d}){J}_{0}(\xi )\right],\end{array}\end{eqnarray} \tag{ 9 }$$

where F₀ is the total flux in the image, f_d is the fraction of that flux that is contained in the disk component, w = W/d is a fractional ring width, W is the FWHM of the Gaussian convolving kernel, d is the diameter of the ring and disk components, ξ ≡ π∣ u ∣d is a normalized radial visibility-domain coordinate, and J_n (ξ) is a Bessel function of the first kind of order n.

The three regions of Sgr A* data identified above are separated by apparent minima in the visibility amplitudes, and they can be approximately characterized by the baseline locations of those minima and the peak flux density levels achieved at the visibility maxima between them. Figure 7 illustrates how this characterization manifests as constraints on the defining parameters of the geometric toy model. The cyan and purple shaded regions in the left panel indicate the approximate ranges of baseline lengths corresponding to the locations of the first and second visibility minima, respectively. The locations of these minima constrain the diameter of the emitting structure, as shown in the top right panel of Figure 7. To be consistent with both a first minimum falling between ∼2.5 and 3.5 Gλ and a second minimum falling between ∼6 and 7 Gλ, the emitting region must be between ∼50 and 60 μas across. The amplitudes of two visibility maxima—one falling between the first and second visibility minima, and the second following the second minimum—constrain a combination of the fractional disk flux f_d and the fractional ring width w. The bottom right panel of Figure 7 shows the constraints from the first and second visibility maxima in red and orange, respectively, and from the ratio of the two in green.

Taken together, even these few, simple, and only modestly constrained visibility features result in a rather narrow permitted range of model parameter values for d, w, and f_d ; an example of a "best-fit" model from within the permitted range is shown by the gray curve in the left panel of Figure 7. However, we stress that the above constraints only strictly hold within the context of the specific toy model used to derive them. More general and robust constraints on the emission structure require a model that can accommodate more than just the gross features; such models are produced as part of the imaging (Paper III and Section 5) and geometric modeling analyses (see Sections 6 and 7) carried out in this paper series.

The ring-like images reconstructed in Paper III are not azimuthally symmetric, but instead show pronounced azimuthal brightness variations that we would like to capture in our geometric modeling analyses. In this section, we specify the "mG-ring" model that we use in Sections 6 and 7 to quantify the morphological properties of the observed Sgr A* emission.

4.3.1. Image-domain Representation of mG-ring Model

Adopting the construction developed by Johnson et al. (2020), we can model an infinitesimally thin circular ring with azimuthal brightness variations using a sum over angular Fourier modes indexed by integer k,

$$\begin{eqnarray}&&{I}_{\mathrm{ring}}(r,\phi )=\displaystyle \frac{{F}_{\mathrm{ring}}}{\pi d}\,\delta \,\left(r-\displaystyle \frac{d}{2}\right)\sum _{k=-m}^{m}{\beta }_{k}{e}^{{ik}\phi }.\end{eqnarray} \tag{ 10 }$$

Here r is the image radial coordinate, ϕ is the azimuthal coordinate (east of north), d is the ring diameter, {β_k } are the set of (dimensionless) complex azimuthal mode coefficients, and m sets the order of the expansion. Because the image is real, ${\beta }_{-k}={\beta }_{k}^{* };$ we enforce β₀ ≡ 1 so that F_ring sets the total flux density of the ring. Given that the images from Paper III show a ring of radius ∼ 25 μas and the diffraction-limited EHT resolution is ∼20 μas, we expect the data to primarily constrain ring modes with $m\lesssim \pi \left(25/20\right)\approx 4$. We refer to this asymmetric ring as an "m-ring" of order m.

For the purposes of constraining additional image structures, we augment this m-ring in two ways. First, we convolve the m-ring with a circular Gaussian kernel of FWHM W,

$$\begin{eqnarray}&&{I}_{\mathrm{ring}}(r,\phi ;W)={I}_{\mathrm{ring}}(r,\phi )* \left[\displaystyle \frac{4\mathrm{ln}(2)}{\pi {W}^{2}}\exp \,\left(-\displaystyle \frac{4\mathrm{ln}(2){r}^{2}}{{W}^{2}}\right)\right].\end{eqnarray} \tag{ 11 }$$

Second, we add a circular Gaussian component that is concentric with the ring, which serves to provide a nonzero brightness floor interior to the ring. The Gaussian component has a total flux density of F_Gauss and an FWHM of W_Gauss,

$$\begin{eqnarray}&&{I}_{\mathrm{Gauss}}(r,\phi )=\displaystyle \frac{4\mathrm{ln}(2){F}_{\mathrm{Gauss}}}{\pi {W}_{\mathrm{Gauss}}^{2}}\exp \,\left(-\displaystyle \frac{4\mathrm{ln}(2){r}^{2}}{{W}_{\mathrm{Gauss}}^{2}}\right).\end{eqnarray} \tag{ 12 }$$

We refer to the resulting composite model I(r, ϕ), where

$$\begin{eqnarray}&&I(r,\phi )={I}_{\mathrm{ring}}(r,\phi ;W)+{I}_{\mathrm{Gauss}}(r,\phi ),\end{eqnarray} \tag{ 13 }$$

as an "mG-ring." An example mG-ring is shown in Figure 8.

Zoom In Zoom Out Reset image size

An mG-ring of order m has 5 + 2m model parameters: the flux density in the ring (F_ring), the diameter of the ring (d), the flux density in the central Gaussian (F_Gauss), the FWHM of the central Gaussian (W_Gauss), the FWHM of the ring convolving kernel (W), and two parameters for each complex Fourier coefficient β_k with 1 ≤ k ≤ m.

4.3.2. Visibility-domain Representation of mG-ring Model

To aid in efficient parameter space exploration, the mG-ring model is intentionally constructed using components and transformations that permit analytic Fourier transformations. The Fourier transform of the m-ring image (Equation (10)) is given by

$$\begin{eqnarray}&&{V}_{\mathrm{ring}}(| {\boldsymbol{u}}| ,{\phi }_{u})={F}_{\mathrm{ring}}\sum _{k=-m}^{m}{\beta }_{k}{J}_{k}(\pi | {\boldsymbol{u}}| d){e}^{{ik}({\phi }_{u}-\pi /2)},\end{eqnarray} \tag{ 14 }$$

where (∣ u ∣, ϕ_u ) are polar coordinates in the Fourier domain. The convolution with a circular Gaussian in the image plane corresponds to multiplication of this function by the Fourier transform of the convolving kernel,

$$\begin{eqnarray}&&{V}_{\mathrm{ring}}(| {\boldsymbol{u}}| ,{\phi }_{u};W)=\exp \left(-\displaystyle \frac{{\pi }^{2}{W}^{2}| {\boldsymbol{u}}{| }^{2}}{4\mathrm{ln}(2)}\right){V}_{\mathrm{ring}}(| {\boldsymbol{u}}| ,{\phi }_{u}).\end{eqnarray} \tag{ 15 }$$

The Fourier transform of the Gaussian image (Equation (12)) is given by

$$\begin{eqnarray}&&{V}_{\mathrm{Gauss}}(| {\boldsymbol{u}}| ,{\phi }_{u})={F}_{\mathrm{Gauss}}\exp \left(-\displaystyle \frac{{\pi }^{2}{W}_{\mathrm{Gauss}}^{2}| {\boldsymbol{u}}{| }^{2}}{4\mathrm{ln}(2)}\right).\end{eqnarray} \tag{ 16 }$$

By the linearity of the Fourier transform, the visibility-domain representation of the mG-ring model is then simply the sum of these two components,

$$\begin{eqnarray}&&V(| {\boldsymbol{u}}| ,{\phi }_{u})={V}_{\mathrm{ring}}(| {\boldsymbol{u}}| ,{\phi }_{u};W)+{V}_{\mathrm{Gauss}}(| {\boldsymbol{u}}| ,{\phi }_{u}).\end{eqnarray} \tag{ 17 }$$

When interpreting model-fitting results in subsequent sections, we are interested in a number of derivative quantities. We will typically work with the fractional thickness of the ring, w, defined to be

$$\begin{eqnarray}&&w\equiv \displaystyle \frac{W}{d}.\end{eqnarray} \tag{ 18 }$$

Similarly, we are typically interested in fractional representations of flux densities. We define

$$\begin{eqnarray*}&&{F}_{0}\equiv {F}_{\mathrm{ring}}+{F}_{\mathrm{Gauss}}\end{eqnarray*}$$

to be the total flux density, and then

$$\begin{eqnarray}&&{f}_{\mathrm{ring}}\equiv \displaystyle \frac{{F}_{\mathrm{ring}}}{{F}_{0}}\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}&&{f}_{\mathrm{Gauss}}\equiv \displaystyle \frac{{F}_{\mathrm{Gauss}}}{{F}_{0}}\end{eqnarray} \tag{ 20 }$$

are the fraction of the total flux density that is contained in the ring and in the Gaussian components, respectively. Note that F₀ is typically close to or fixed to unity as a consequence of normalizing the data by the light curve. We also define a fractional central flux as

$$\begin{eqnarray}&&{f}_{{\rm{c}}}\equiv \displaystyle \frac{{F}_{\mathrm{Gauss}}(r\lt d/2)}{{F}_{\mathrm{ring}}+{F}_{\mathrm{Gauss}}(r\lt d/2)},\end{eqnarray} \tag{ 21 }$$

where F_Gauss(r < d/2) is the integrated flux density of the central Gaussian component interior to the ring radius, given by

$$\begin{eqnarray}&&{F}_{\mathrm{Gauss}}(r\lt d/2)={F}_{\mathrm{Gauss}}\left[1-\exp \left(-\displaystyle \frac{{d}^{2}\mathrm{ln}(2)}{{W}_{\mathrm{Gauss}}^{2}}\right)\right].\end{eqnarray} \tag{ 22 }$$

Following M87* Paper IV, the m-ring position angle η and degree of azimuthal asymmetry A are both determined by the coefficient of the m = 1 mode,

$$\begin{eqnarray}\begin{array}{rcl}\eta & \equiv & \arg \left({\displaystyle \int }_{0}^{2\pi }I(\phi ){e}^{i\phi }d\phi \right)\\ & = & -\arg \left({\beta }_{1}\right),\end{array}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}\begin{array}{rcl}A & \equiv & \displaystyle \frac{\left|{\displaystyle \int }_{0}^{2\pi }I(\phi ){e}^{i\phi }d\phi \right|}{{\displaystyle \int }_{0}^{2\pi }I(\phi ){e}^{i\phi }d\phi }\\ & = & | {\beta }_{1}| .\end{array}\end{eqnarray} \tag{ 24 }$$

A number of these derivative quantities are illustrated in the example mG-ring shown in Figure 8.

The parameters returned by the geometric modeling and feature extraction analyses used in this paper to describe the Sgr A* emission structure do not correspond directly to physical quantities. Instead, the relationship between measured and physical quantities must be calibrated using data for which we know the correct underlying physical system's defining parameters. For ring size measurements, the associated physical quantity of interest is related to the angular size of the gravitational radius,

$$\begin{eqnarray}&&{\theta }_{g}=\displaystyle \frac{{GM}}{{c}^{2}D}.\end{eqnarray} \tag{ 25 }$$

which sets the absolute scale of the system.

Under the assumption that the emission near the black hole originates from some "typical" radius, a measurement of the angular diameter d of the emitting region will be related to θ_g by a scaling factor α,

$$\begin{eqnarray}&&d=\alpha {\theta }_{g}.\end{eqnarray} \tag{ 26 }$$

If the observations were directly sensitive to the critical curve bounding the black hole shadow, then α could be determined analytically and would take on a value ranging from ∼9.6 to 10.4 depending on the black hole spin and inclination (Bardeen 1973; Takahashi 2004). For more realistic emission structures and measurement strategies, the value of α cannot be determined from first principles and must instead be calibrated.

Our α calibration strategy generally follows the procedure developed in M87* Paper VI. Using the library of GRMHD simulations described in Paper V, we generate a suite of 100 synthetic data sets that emulate the cadence and sensitivity of the 2017 EHT observations and that contain a realistic character and magnitude of data corruption; Appendix D describes the generation of these synthetic data sets. In the analyses described in Sections 5, 6, and 7, 90 of these 100 synthetic data sets are used to derive the α calibration for each analysis pathway, while the remaining 10 data sets are used to validate the calibration.

After carrying out ring size measurements on each of the data sets in the suite, we determine α (for each specific combination of data set and measurement technique) by dividing the measured ring diameter by the known value of θ_g (per Equation (26)). For a given measurement technique, the distribution of α values that results from applying this procedure to the entire suite of synthetic data sets then provides a measure of α and its theoretical uncertainty. The α value associated with each measurement technique can then be used to translate Sgr A* ring size measurements into their corresponding θ_g constraints. We note that this calibration strategy assumes that the images contained in the GRMHD library provide a reliable representation of the emission structure in the vicinity of Sgr A*; a separate calibration strategy that relaxes this GRMHD assumption is presented in Paper VI.

Appendix I describes elements of the calibration and validation strategy that are specific to each of the analysis pathways detailed in Sections 5, 6, and 7.

The imaging analyses carried out in Paper III use four different algorithms classified into three categories: one sampling-based posterior exploration algorithm (Themis; Broderick et al. 2020a, 2020b), one CLEAN-based deconvolution algorithm (DIFMAP; Shepherd 1997), and two "regularized maximum likelihood" (RML) algorithms (eht-imaging,Chael et al. 2016, 2018; SMILI, Akiyama et al. 2017a, 2017b). All methods produce image reconstructions using band-combined data (i.e., both low band and high band), and the latter three are run on two versions of the Sgr A* data: a "descattered" version that attempts to deconvolve the effects of the diffractive scattering kernel from the data, and a "scattered" version that applies no such deconvolution. The posterior exploration imaging method Themis instead applies the effects of diffractive scattering as part of its internal forward model, rather than deconvolving the data; the analogous "scattered" and "descattered" versions of the Themis images thus correspond simply to those for which the scattering kernel has been applied or not, respectively. The posterior exploration imaging jointly reconstructs the combined April 6 and April 7 data sets (see Appendix B), while the CLEAN and RML imaging reconstructs each day individually, focusing primarily on the April 7 data and using the April 6 data for cross-validation. Example fits and residuals for each of the imaging pipelines are shown in Figure 9, ¹⁵² and χ² statistics for each image are provided in Appendix E; detailed descriptions of the data preprocessing and imaging procedures for each imaging algorithm are provided in Paper III.

Zoom In Zoom Out Reset image size

For the CLEAN and RML imaging methods, there are a number of tunable hyperparameters associated with each algorithm whose values are determined through extensive "parameter surveys" carried out on synthetic data sets. During a parameter survey, images of each synthetic data set are reconstructed using a broad range of possible values for each hyperparameter. Settings that produce high-fidelity image reconstructions across all synthetic data sets are collected into a "top set" of hyperparameters, and these settings are then applied for imaging the Sgr A* data. The resulting top sets of Sgr A* images capture emission structures that are consistent with the data, and we use these top-set images for the feature extraction analyses in this paper.

The Themis imaging algorithm explores a posterior distribution over the image structure, and there are no hyperparameters that require synthetic data surveys to determine. Rather than producing a top set of images, Themis instead produces a sample of images drawn from the posterior determined from the Sgr A* data. We use these posterior image samples for the feature extraction analyses in this paper.

Given the top-set and posterior images from Paper III, we carry out IDFE analyses using two separate tools: REx and VIDA. An independent cross-validation of both IDFE tools has been carried out in P. Tiede et al. (2022, in preparation). In this section, we provide a brief overview of each method and specify the details relevant for the analyses presented in this paper.

5.2.1.  REx

The Ring Extractor (REx) is an IDFE tool for quantifying the morphological properties of ring-like images. It is available as part of the eht-imaging software library and is described in detail in Chael (2019). REx was the main tool used in M87* Paper IV to extract ring properties from the M87* images, and detailed definitions of the various REx parameters are provided in that paper.

For the majority of the REx-derived ring parameters, we retain the same definitions as used in M87* Paper IV. REx first defines a ring center (x₀, y₀), which is determined to be the point in the image from which radial intensity profiles have a minimum dispersion in their peak intensity radii. The ring radius, r₀, is then taken to be the average of these peak intensity radii over all angles, and the ring thickness w is taken to be the angular average of the FWHM about the peak measured along each radial intensity profile. To avoid biases associated with a nonzero floor to the image brightness outside of the ring, we subtract out the quantity

$$\begin{eqnarray}&&{I}_{\mathrm{floor}}=\frac{1}{2\pi }{\int }_{0}^{2\pi }I({r}_{\max }=60\,\mu \mathrm{as},\phi )d\phi \end{eqnarray} \tag{ 27 }$$

when computing the FWHM, i.e., we compute the average FWHM of I(r, ϕ) − I_floor. For all other ring parameters, the definitions remain the same as those used in M87* Paper IV.

REx defines the ring position angle η and asymmetry A as the argument and amplitude, respectively, of the first circular mode,

$$\begin{eqnarray}&&{\beta }_{1}={\lt\frac{{\int }_{0}^{2\pi }I(\phi )\cos (\phi )d\phi }{{\int }_{0}^{2\pi }I(\phi )d\phi }\gt}_{r},\end{eqnarray} \tag{ 28 }$$

where the angled brackets denote a radial average between r₀ − w/2 and r₀ + w/2. ¹⁵³ These definitions are analogous to those used to define the corresponding position angle and asymmetry of the mG-ring model (Equations (23) and (24), respectively). The fractional central brightness f_c is defined to be the ratio of the mean brightness within 5 μas of the center to the azimuthally averaged brightness along the ring (i.e., along r = r₀).

As in Paper III, we replace the negative pixels in Themis images with zero values before performing REx analyses.

5.2.2.  VIDA

Variational Image-Domain Analysis (VIDA; P. Tiede et al. 2022, in preparation) is an IDFE tool for quantifying the parameters describing a specifiable image morphology; it is written in Julia (Bezanson et al. 2017) and contained in the package VIDA.jl. ¹⁵⁴ VIDA employs a template-matching approach for image analysis, using parameterized templates to approximate an image and adjusting the parameters of the templates until a specified cost function is minimized. Within VIDA, the cost function takes the form of a probability divergence, which provides a distance metric between the image and template; the template parameters that minimize this divergence are taken to provide the best description of the image. The VIDA optimization strategy and additional details are provided in P. Tiede et al. (2022, in preparation).

For the IDFE analyses in this paper, we use VIDA's SymCosineRingwFloor template and the least-squares divergence (for details, see Section 8 of Paper III). This template describes an image structure that is similar to the mG-ring model (Section 4.3), and it is characterized by a ring center (x₀, y₀), a ring diameter d = 2r₀, an FHWM fractional ring thickness w, and a cosine expansion describing the azimuthal brightness distribution S(ϕ),

$$\begin{eqnarray}&&S(\phi )=1-2\sum _{k=1}^{m}{A}_{k}\cos \left[k(\phi -{\eta }_{k})\right].\end{eqnarray} \tag{ 29 }$$

To maintain consistency with the geometric modeling analyses (see Sections 6 and 7), we use m = 4. We also restrict the value of the A₁ parameter to be <0.5 to avoid negative flux in the template. As with the mG-ring model, the orientation η is equal to the first-order phase η₁, and the asymmetry A is equal to the first-order coefficient A₁.

To permit the presence of a central brightness floor, the SymCosineRingwFloor template contains an additional component in the form of a circular disk whose center point is fixed to coincide with that of the ring. The disk radius is fixed to be r₀. A Gaussian falloff is stitched to the outer edge of the disk, such that for radii larger than r₀ the intensity profile becomes a Gaussian with mean r₀ and an FWHM that matches the ring thickness. The flux of this disk component is a free parameter in the template. We then retain the same definition of the fractional central brightness f_c as used by REx.

The output of the IDFE analysis is a set of distributions for the ring parameters from each imaging method; Figure 10 shows an example set of results from applying both IDFE software packages to the descattered Sgr A* posterior and top-set images. However, both REx and VIDA implicitly assume that the images fed into them contain a ring-like emission structure. If the input image does not contain a ring, then the output measurements may not be meaningful. For each input image, we thus wish to determine both whether the image contains a ring-like structure and how sensitive the IDFE results are to the specific manner in which "a ring-like structure" is defined.

Zoom In Zoom Out Reset image size

To determine whether the images we are analyzing with REx and VIDA contain ring-like structures, we use metronization, ¹⁵⁵ a software that preprocesses the images into a form suitable for topological analysis and extracts topologically relevant features with the help of the open-source computational topology code Dionysus 2 ¹⁵⁶ (Morozov 2017). A detailed description of metronization can be found in Christian et al. (2022).

The metronization preprocessing procedure consists of the following steps:

• 1.  
  First, the image undergoes a "robust" thresholding step, in which the pixels are sorted by brightness in a cumulative sequence, and all pixels below a certain threshold in this sequence have their values set to zero and the rest are set to a value of one.
• 2.  
  Next, in a process called "skeletonization," the Boolean image produced in the first step is reduced to its topological skeleton that preserves the topological characteristics of the original shape. This step thins large contiguous areas of flagged pixels and enlarges the "holes."
• 3.  
  The topological skeleton is rebinned and downsampled. Holes smaller than the rebinning resolution are preserved by the skeletonization in the previous step.
• 4.  
  The downsampled image undergoes skeletonization once more.

The resulting output is a low-resolution image that preserves the topologically relevant information from the original image, speeding up the application of computationally expensive topological algorithms that follow. A technique known as persistent homology is then used to convert this low-resolution image into a topological space that preserves features that are topologically invariant. It computes a quantity known as the first Betti number that provides a metric for measuring the number of holes present in the image.

The metronization software contains a number of tunable parameters that determine how closely the emission structure in the input image must resemble that of a topological ring, and for how many cumulative threshold levels it must persist, for it to be classified as a ring. We identify three modes for these parameters—a "permissive" mode, a "moderate" mode, and a "strict" mode—and explore the impact on the REx and VIDA measured parameter distributions when the input top sets and posterior images are restricted only to those that are classified as containing rings. We compare these results with those of a fourth, default setting in which the top sets and posterior images are not filtered by the classification prescribed by metronization.

We note that metronization differs from the ring identification methods presented in Paper III in that it searches for the presence of a topological ring in the input image. Figure 11 compares the mean ring and nonring descattered images for each imaging pipeline as classified by metronization in the "permissive" mode and the clustering analysis from Paper III. Both methods classify all the posterior imaging samples as rings, while the top-set imaging samples contain both ring and nonring images. We find that the mean ring and nonring images for each imaging pipeline are broadly consistent between the two classification methods.

Zoom In Zoom Out Reset image size

The definition of what constitutes a ring is subjective, and there will always be images that are ambiguous to the human eye. Different automated methods will classify these images differently. Hence, it is important to verify that the ring parameters measured by REx and VIDA are robust against the specifics of the ring identification scheme used. Figure 12 shows the resulting diameter distributions from ring fitting to the descattered Sgr A* images from all imaging pipelines, split out by metronization setting. As we move from the most to the least permissive classification scheme, the tails in the distributions are diminished while the primary peaks are sharpened, but the mean and general shape of the distribution remain largely unchanged. This trend indicates that while metronization penalizes images with emission structures deviating from a topological ring, the distributions of the REx and VIDA measurements are robust against the choice of metronization mode employed.

Zoom In Zoom Out Reset image size

Prior to fitting the mG-ring model to real or synthetic data, we process the data using the pre-imaging pipeline described in Paper III. This preprocessing procedure entails light-curve normalization and an inflation of the error budget to account for residual calibration uncertainties and the effects of refractive scattering in the interstellar medium toward Sgr A*. Specifically, the total error budget σ_sb for a visibility measured on the baseline b during snapshot s is given by

$$\begin{eqnarray}&&{\sigma }_{{sb}}^{2}={\sigma }_{\mathrm{th},{sb}}^{2}+{f}^{2}| V{| }_{{sb}}^{2}+{\sigma }_{\mathrm{ref},{sb}}^{2}.\end{eqnarray} \tag{ 30 }$$

Here the first term corresponds to the baseline-specific thermal noise (see Equation (2)), the second term is a component that is multiplicative in the visibility amplitude and is intended to capture residual (nongain) calibration errors (e.g., residual polarization leakage), and the third term is the J18model1 refractive scattering noise from Paper III. For the snapshot modeling, we fix f = 0.02 per the analyses carried out in Paper II. The pre-imaging pipeline also mitigates the impact of diffractive scattering by "deblurring" the data using the Johnson et al. (2018) model.

Following the application of the pre-imaging pipeline, we split the data into 120 s segments, or "snapshots," and coherently average the visibilities in each snapshot over the 120 s window. Finally, we flag snapshots that contain fewer than four unique stations, so as to retain snapshots during which closure amplitudes can be formed.

The first step of our snapshot modeling procedure is to determine the posterior distribution for the mG-ring model parameters on each snapshot of data. The observation is divided up into N_s independent snapshots, which we label using a snapshot index s. Within each snapshot we fit the mG-ring model described in Section 4.3, whose parameter vector we denote as θ _s . For a single snapshot, the posterior is given by Bayes's theorem,

$$\begin{eqnarray}&&{P}_{s}({{\boldsymbol{\theta }}}_{s}| {{\boldsymbol{D}}}_{s})\displaystyle \frac{{{ \mathcal L }}_{s}({{\boldsymbol{D}}}_{s}| {{\boldsymbol{\theta }}}_{s}){\pi }_{s}({{\boldsymbol{\theta }}}_{s})}{{{ \mathcal Z }}_{s}({{\boldsymbol{D}}}_{s})},\end{eqnarray} \tag{ 31 }$$

where D _s denotes the data available on snapshot s, ${{ \mathcal L }}_{s}$ is the likelihood, π_s is the prior distribution, and ${ \mathcal Z }$ is the Bayesian evidence.

In our snapshot modeling analyses we make use of three different classes of interferometric data products: visibility amplitudes ∣V∣, log closure amplitudes $\mathrm{ln}A$, and closure phases ψ. Each analysis uses only a single amplitude data product (either visibility amplitudes or log closure amplitudes) along with the closure phases. For analyses that use visibility amplitudes and closure phases, the likelihood is given by

$$\begin{eqnarray}&&{{ \mathcal L }}_{s}={{ \mathcal L }}_{| V| ,s}{{ \mathcal L }}_{\psi ,s},\end{eqnarray} \tag{ 32 }$$

while for those that use log closure amplitudes, we instead have

$$\begin{eqnarray}&&{{ \mathcal L }}_{s}={{ \mathcal L }}_{A,s}{{ \mathcal L }}_{\psi ,s}.\end{eqnarray} \tag{ 33 }$$

Here ${{ \mathcal L }}_{| V| ,s}$, ${{ \mathcal L }}_{\psi ,s}$, and ${{ \mathcal L }}_{A,s}$ are components of the likelihood on snapshot s associated with the visibility amplitudes, closure phases, and log closure amplitudes, respectively. We assume Gaussian likelihood functions for the amplitude data components and a von Mises likelihood function for the closure phases; the detailed expressions for each likelihood function are provided in Appendix F.

The output of a snapshot fitting analysis is a set of posterior samples for the model parameters from each individual snapshot; Figure 13 shows an example set of posterior distributions for the mG-ring diameter parameter on each snapshot in the April 6 and April 7 data sets. To arrive at a single posterior on these parameters that combines the information from all snapshots across both days, we use a Bayesian hierarchical model similar to the one used in Baronchelli et al. (2020). This approach treats the model fit to each snapshot as a realization from some average model or "hypermodel."

Zoom In Zoom Out Reset image size

6.3.1. Averaging Procedure

We denote the parameters of the average model as $\displaystyle \bar{{\boldsymbol{\theta }}\,}$ and the distribution of the snapshot model conditioned on the average model by ${\pi }_{s}({{\boldsymbol{\theta }}}_{s}| \displaystyle \bar{{\boldsymbol{\theta }}})$. Given this conditional probability, the joint snapshot and average parameter posterior is given by

$$\begin{eqnarray}&&P({\boldsymbol{\Theta }},\displaystyle \bar{{\boldsymbol{\theta }}\,}\,| {\boldsymbol{ \mathcal D }})=\displaystyle \frac{\pi (\displaystyle \displaystyle \bar{{\boldsymbol{\theta }}\,}\,)}{{ \mathcal Z }({\boldsymbol{ \mathcal D }})}\prod _{s}{{ \mathcal L }}_{s}({{\boldsymbol{D}}}_{s}| {{\boldsymbol{\theta }}}_{s}){\pi }_{s}({{\boldsymbol{\theta }}}_{s}| \displaystyle \displaystyle \bar{{\boldsymbol{\theta }}\,}\,),\end{eqnarray} \tag{ 34 }$$

where $\pi (\displaystyle \bar{{\boldsymbol{\theta }}\,}\,)$ is the prior distribution for the hypermodel parameters (the "hyperprior"), ${\boldsymbol{\Theta }}=({{\boldsymbol{\theta }}}_{1},{{\boldsymbol{\theta }}}_{2},...,{{\boldsymbol{\theta }}}_{{N}_{s}})$ is the parameter vector across all snapshots, and ${\boldsymbol{ \mathcal D }}=({{\boldsymbol{D}}}_{1},{{\boldsymbol{D}}}_{2},...,{{\boldsymbol{D}}}_{{N}_{s}})$ is the data vector across all snapshots. To find the marginal average parameter posterior, we integrate this expression over all the snapshot parameters,

$$\begin{eqnarray}&&P(\displaystyle \bar{{\boldsymbol{\theta }}\,}\,| {\boldsymbol{ \mathcal D }})=\displaystyle \frac{\pi (\displaystyle \bar{{\boldsymbol{\theta }}\,}\,)}{{ \mathcal Z }({\boldsymbol{ \mathcal D }})}\int \left(\prod _{s}{{ \mathcal L }}_{s}({{\boldsymbol{D}}}_{s}| {{\boldsymbol{\theta }}}_{s}){\pi }_{s}({{\boldsymbol{\theta }}}_{s}| \displaystyle \bar{{\boldsymbol{\theta }}\,}\,)\right){\rm{d}}{\boldsymbol{\Theta }}.\end{eqnarray} \tag{ 35 }$$

In general, this integral is analytically intractable. However, a bit of manipulation permits us to use the posterior samples from the individual snapshot fits to make headway. Because the snapshots are independent, we can swap the order of the integral and product in Equation (35) and use Bayes's theorem to substitute in for the snapshot likelihood (Equation (31)), giving

$$\begin{eqnarray}&&P(\displaystyle \bar{{\boldsymbol{\theta }}\,}\,| {\boldsymbol{ \mathcal D }})=\pi (\displaystyle \bar{{\boldsymbol{\theta }}\,}\,)\prod _{s}\left(\int P({{\boldsymbol{\theta }}}_{s}| {{\boldsymbol{D}}}_{s})\displaystyle \frac{{\pi }_{s}({{\boldsymbol{\theta }}}_{s}| \displaystyle \bar{{\boldsymbol{\theta }}\,}\,)}{{\pi }_{s}({{\boldsymbol{\theta }}}_{s})}{\rm{d}}{{\boldsymbol{\theta }}}_{s}\right).\end{eqnarray} \tag{ 36 }$$

Note that the evidence term from the prefactor denominator in Equation (35) has now been subsumed into the posterior term P( θ _s ∣ D _s ) inside of the integral. To evaluate Equation (36), we make use of the fact that the snapshot posterior samples ${{\boldsymbol{\theta }}}_{s}^{(i)}$ permit us to approximate the integral by a sum,

$$\begin{eqnarray}&&P(\displaystyle \bar{{\boldsymbol{\theta }}\,}\,| {\boldsymbol{ \mathcal D }})\approx {\pi }_{s}(\displaystyle \bar{{\boldsymbol{\theta }}\,}\,)\prod _{s}\sum _{i}\displaystyle \frac{{\pi }_{s}({{\boldsymbol{\theta }}}_{s}^{(i)}| \displaystyle \bar{{\boldsymbol{\theta }}\,}\,)}{\pi ({{\boldsymbol{\theta }}}_{s}^{(i)})}.\end{eqnarray} \tag{ 37 }$$

We can use this expression to sample from the posterior distribution over just the hypermodel parameters $\displaystyle \bar{{\boldsymbol{\theta }}}$, having fully marginalized over the parameters from each individual snapshot.

We note that the averaging procedure described here is simply a generalization of standard inverse-variance weighting. If we consider a delta-function hypermodel that contains only a single parameter (i.e., the to-be-determined mean value) for each snapshot model parameter, then in the limit where the individual snapshot posteriors P( θ _s ∣ D _s ) are Gaussian and the priors on $\displaystyle \bar{{\boldsymbol{\theta }}\,}$ and θ _s are uninformative, the posterior maximum for $P(\displaystyle \bar{{\boldsymbol{\theta }}\,}\,| {\boldsymbol{ \mathcal D }})$ is equal to the mean of the snapshot posterior means weighted by their inverse posterior variances. However, because the model we employ (described in the following section) does not conform to the necessary conditions (see, e.g., the non-Gaussian snapshot posteriors shown in Figure 13), we proceed with the more general averaging procedure.

6.3.2. Hypermodel Specification

We now need to specify the hypermodel $\pi ({{\boldsymbol{\theta }}}_{s}| \displaystyle \bar{{\boldsymbol{\theta }}\,}\,)$ that determines the distribution from which the individual snapshot models are drawn; for simplicity, we choose a hypermodel that is approximately Gaussian. Let $\displaystyle \bar{{\boldsymbol{\theta }}\,}\,=({\boldsymbol{\mu }},{\boldsymbol{\sigma }})$, where μ is a vector of the mean parameter values and σ is a vector containing their standard deviations across scans. We assign most hypermodel parameters to be distributed according to a truncated normal distribution,

$$\begin{eqnarray}&&{\pi }_{s,\mathrm{tN}}({\theta }_{s,i}| {\mu }_{i},{\sigma }_{i})={ \mathcal T }({\theta }_{s,i}| {\mu }_{i},{\sigma }_{i};{a}_{i},{b}_{i}).\end{eqnarray} \tag{ 38 }$$

Here ${ \mathcal T }(\theta | \mu ,\sigma ;a,b)$ denotes the density for a truncated normal distribution with mean μ and standard deviation σ, and whose lower and upper bounds are given by a and b, respectively; we index the separate parameters by i. This truncation is necessary to ensure that the support of the hypermodel parameters matches that of the individual snapshot model parameters. However, for angular parameters—i.e., those with values that are periodic in [0, 2π)—we instead use a von Mises distribution,

$$\begin{eqnarray}&&{\pi }_{s,\mathrm{vM}}({\theta }_{s,j}| {\mu }_{j},{\sigma }_{j})=\displaystyle \frac{1}{2\pi {I}_{0}({\sigma }_{j}^{-2})}\exp \left[\displaystyle \frac{\cos ({\theta }_{s,j}-{\mu }_{j})}{{\sigma }_{j}^{2}}\right].\end{eqnarray} \tag{ 39 }$$

The subscripted "tN" in Equation (38) indicates that the corresponding parameters use a truncated normal prior, while the subscripted "vM" in Equation (39) similarly indicates that the corresponding parameters use a von Mises prior. The total hypermodel is then given by

$$\begin{eqnarray}&&\pi ({{\boldsymbol{\theta }}}_{s}| {\boldsymbol{\mu }},{\boldsymbol{\sigma }})=\prod _{i}{\pi }_{s,\mathrm{tN}}({\theta }_{s,i}| {\mu }_{i},{\sigma }_{i})\prod _{j}{\pi }_{s,\mathrm{vM}}({\theta }_{s,j}| {\mu }_{j},{\sigma }_{j}),\end{eqnarray} \tag{ 40 }$$

where the first product runs over nonangular parameters and the second runs over angular parameters.

We set the hyperpriors for μ to be equal to the corresponding snapshot priors, which are specified in Table 1. For σ we instead use a half-normal hyperprior,

$$\begin{eqnarray}&&\pi ({\boldsymbol{\sigma }})=\prod _{i}{ \mathcal T }({\sigma }_{i}| 0,{L}_{i}/4;\ 0,{L}_{i}),\end{eqnarray} \tag{ 41 }$$

where L_i = (b_i − a_i ) is the breadth of support for parameter θ_s,i . Appendix G describes the level of consistency between these selected hyperpriors and the priors for the individual snapshot model parameters.

Table 1. Snapshot Modeling mG-ring Priors

Parameter	Comrade	DPI
F0	${ \mathcal U }(0.8,1.2)$	δ(1)
d (μas)	${ \mathcal U }(25,85)$	...
${d}^{{\prime} }$ (μas)	...	${ \mathcal U }(25,85)$
W (μas)	${ \mathcal U }(1,40)$	...
${W}^{{\prime} }$ (μas)	...	${ \mathcal U }(1,40)$
∣βm ∣	${ \mathcal U }(0,0.5)$	${ \mathcal U }(0,0.5)$
$\arg ({\beta }_{m})$ (deg)	${ \mathcal U }(-180,180)$	${ \mathcal U }(-180,180)$
fGauss	${ \mathcal U }(0,1)$	${ \mathcal U }(0,1)$
WGauss (μas)	${ \mathcal U }(40,200)$	${ \mathcal U }(40,200)$

Note. Prior distributions for Comrade and DPI snapshot geometric modeling analyses. ${ \mathcal U }(a,b)$ denotes a uniform prior on the interval [a, b], and δ(a) denotes a delta-function (i.e., fixed-value) prior, with the parameter value fixed at a. For the definitions of the parameters see Section 4.3.

Download table as:  ASCIITypeset image

We use three different software packages to carry out snapshot geometric modeling on the Sgr A* data and a fourth software to perform the hypermodel sampling. In this section, we specify the relevant implementation specifics for these different tools. Cross-validation tests are detailed in Appendix H.

6.4.1.  Comrade

6.4.2.  eht-imaging

6.4.3.  DPI

The third software we use for snapshot geometric modeling is the Python code Deep Probabilistic Imaging/Inference (DPI/α -DPI; Sun & Bouman 2021; Sun et al. 2022). DPI approximates the posterior over all model parameters by fitting a normalizing flow neural network (Rezende & Mohamed 2015) to the data using a Rényi α-divergence variational inference technique (Li & Turner 2016). DPI is an optimization-based posterior estimation framework, and it uses the auto-differentiation package PyTorch (Paszke et al. 2017) to optimize the neural network weights. The posterior estimation accuracy is further improved post-optimization through importance reweighting of the samples generated by the normalizing flow neural network.

DPI supports fitting to multiple data products, including visibility amplitudes and closure quantities, but it does not currently support the inclusion of station gain amplitudes as model parameters. We thus use DPI to fit to closure phases and log closure amplitudes; the snapshot likelihood is given by Equation (33). Prior to fitting, all time stamps that are unable to form at least one closure phase and at least one closure amplitude are flagged.

DPI differs from both Comrade and eht-imaging in that it defines geometric models in the image domain rather than in the visibility domain, and it uses a nonuniform fast Fourier transform (NFFT) to compute the necessary data products. For the analyses carried out in this paper, we discretize the model as an image containing 32 × 32 pixels spanning a 160 μas field of view.

Because the pixel size is finite, DPI cannot support a model containing infinitesimally thin rings such as that in Equation (10); furthermore, convolutions in the image domain are computationally expensive. The DPI fits in this paper thus employ a modified version of the mG-ring model specification,

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\mathrm{ring}}^{{\prime} }(r,\phi ) & = & {F}_{\mathrm{ring}}\displaystyle \frac{\exp \left[\tfrac{4\mathrm{ln}2}{{{W}^{{\prime} }}^{2}}{\left(r-\tfrac{{d}^{{\prime} }}{2}\right)}^{2}\right]\sum _{k=-m}^{m}{\beta }_{k}{e}^{{ik}\phi }}{2\pi {\displaystyle \int }_{0}^{\infty }r\exp \left[\tfrac{4\mathrm{ln}2}{{{W}^{{\prime} }}^{2}}{\left(r-\tfrac{{d}^{{\prime} }}{2}\right)}^{2}\right]\,{dr}},\\ {I}^{{\prime} }(r,\phi ) & = & {I}_{\mathrm{ring}}^{{\prime} }(r,\phi )+{I}_{\mathrm{Gauss}}(r,\phi ),\\ {F}_{0} & = & {F}_{\mathrm{ring}}+{F}_{\mathrm{Gauss}}=1,\end{array}\end{eqnarray} \tag{ 42 }$$

where we note that ${d}^{{\prime} }$ and ${W}^{{\prime} }$ are conceptually distinct from d and W. The quantities d and W in the mG-ring model from Section 4.3 determine the diameter of the infinitesimally thin ring and the FWHM of its convolving kernel, respectively. In contrast, the quantities ${d}^{{\prime} }$ and ${W}^{{\prime} }$ determine the intensity peak and FWHM, respectively, of a radial Gaussian function. These two specifications converge only in the limit of large d and small W. In addition, the total flux of the DPI model implementation is fixed to be 1 Jy because DPI fits only to closure quantities and closure amplitudes are not sensitive to the absolute flux scale.

6.4.4. Sampling the Hypermodel Posterior

The mG-ring model described in Section 4.3 is not actually a single model but rather a class of models, delineated by the order m (Equation (10)). To determine the m-order that is preferred by the Sgr A* data, we carry out a series of snapshot mG-ring fits to the Sgr A* using different values of m and compare Bayesian evidence estimates. Given a set of log-evidences $\mathrm{ln}{{ \mathcal Z }}_{s}$ computed for every snapshot s in a single observation, the total evidence for the entire observation is simply given by their sum,

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal Z }=\sum _{s}\mathrm{ln}{{ \mathcal Z }}_{s}.\end{eqnarray} \tag{ 43 }$$

The Comrade snapshot fitting analyses directly estimate the Bayesian evidence on every snapshot, and so the total evidence across an entire observation can be computed directly using Equation (43). The results of a Comrade m-order survey covering m = {1, 2, 3, 4, 5} are shown in Figure 14. We find that the m = 4 order is preferred in both bands and across both calibration pipelines. We thus use the m = 4 mG-ring as our fiducial model for all Comrade Sgr A* analyses in this paper.

Zoom In Zoom Out Reset image size

Unlike Comrade, DPI does not directly estimate the Bayesian evidence during each fit. Instead, we use the evidence lower bound (ELBO) to determine the m-order preference. The ELBO is a combination of the true evidence modified by a relative entropy term that encodes the performance of the variational approximation,

$$\begin{eqnarray}&&\mathrm{ELBO}(m)=\mathrm{log}p({{\boldsymbol{D}}}_{s}\parallel m)-{D}_{{KL}}[q({{\boldsymbol{\theta }}}_{s})\parallel p({{\boldsymbol{\theta }}}_{s}| {{\boldsymbol{D}}}_{s},m)],\end{eqnarray} \tag{ 44 }$$

where D_KL[A∥B] is the Kullback–Leibler divergence of A from B, and q( θ _s ) is the optimized DPI normalizing flow distribution. The relative entropy term is zero when the DPI distribution q( θ _s ) and the true posterior p( θ _s ∣ D _s ) are identical, so the ELBO provides a rough estimate of the log-evidence. The results of a DPI m-order survey covering m = {1, 2, 3, 4} indicate that either m = 1 or m = 2 is preferred, depending on the day and band. We choose to err on the side of increased model flexibility and use the m = 2 mG-ring as our fiducial model for all DPI Sgr A* analyses in this paper.

Figure 15 shows representative mG-ring fits to the Sgr A* HOPS low-band data for both the Comrade and DPI pipelines. In all cases, we find that the normalized residuals are distributed around a value of zero with a subunity variance, and there is no evidence of systematic structure. The χ² statistics for each of these fits are provided in Appendix E.

Zoom In Zoom Out Reset image size

The data preparation for the full-track geometric modeling analyses is similar to that used for snapshot geometric modeling analyses (see Section 6.1). The data are first processed through the pre-imaging pipeline described in Paper III, which applies light-curve normalization and performs some a priori gain calibration. However, unlike in Paper III and Section 6.1, we do not modify the data uncertainties at all beyond their thermal noise values; neither a systematic error term nor a refractive scattering term is added to the error budget. Additionally, no "deblurring" is applied to the data; instead, the blurring is applied directly to the model as described in the next section.

Following the application of the pre-imaging pipeline, we coherently average the visibilities from each baseline on a per-scan basis. A scan length (∼10 minutes) is approximately the amount of time over which we expect structural variability to be subdominant to other sources of uncertainty (see Section 3.1, in particular Figure 1). Furthermore, the station gains are expected to be constant in time across a single scan but not from one scan to the next (M87* Paper III; Paper II), meaning that a scan length is also the longest coherent integration time that the a priori calibration can support.

While the full-track modeling is necessarily focused on reconstructing a time-averaged image structure, the underlying data remain a collection of complex visibilities that sample different instantaneous realizations of the intrinsic Sgr A* source structure. As a consequence, the Sgr A* data exhibit subhour correlations that over a single day are localized in the (u,v)-plane. As previously noted in Section 3.4 and detailed in Appendix B, these unmodeled correlations can result in significant biases in the reconstructed properties of Sgr A*. However, by fitting to multiple days of Sgr A* data, and thus combining multiple samplings of the variable source structure at each location in the (u,v)-plane, we better match the statistical properties of the data to those assumed by the full-track analysis. An additional benefit of combining days is that the multiday analyses more clearly emphasize the static signatures of gravitational lensing from the spurious astrophysical variability. For these reasons, all full-track analyses presented in Section 8 make use of the combined April 6 and April 7 Sgr A* data. For comparison we provide single-day analysis results in Appendix C.

The goal of the full-track geometric modeling procedure is to determine the posterior distribution for the parameters of the static mG-ring model and parameterized noise model that best describe an entire Sgr A* data set. Our specification for the mG-ring model is described in Section 4.3, and we retain the same notation and terminology in this section. Additionally, we incorporate the blurring effects of scattering in the same manner described in Section 4 of Paper III, through multiplication of the mG-ring visibilities by the Fourier transform of the scattering kernel. We note, however, that all images shown in the figures in this paper correspond to the underlying (i.e., nonscattered) image.

7.2.1. Parameterized Noise Model

Our parameterized noise model for a complex visibility V_i measured on a baseline u _i is given by

$$\begin{eqnarray}&&{\sigma }_{i}^{2}={\sigma }_{\mathrm{th},i}^{2}+{f}^{2}| V{| }_{i}^{2}+{\sigma }_{\mathrm{ref}}^{2}+{\sigma }_{\mathrm{var}}^{2}(| {{\boldsymbol{u}}}_{i}| ).\end{eqnarray} \tag{ 45 }$$

Here the first term is the thermal noise in the measurement (see Equation (2)), the second term is a component that is multiplicative in the visibility amplitude ∣V∣_i and is intended to capture residual (nongain) calibration errors (e.g., residual polarization leakage), the third term is a component that is additive and is intended to account for refractive scattering noise, and the fourth term is a component that is a function of the baseline length ∣ u _i ∣ and is intended to capture the effects of source variability. With the exception of the variability term, Equation (45) is similar to the noise budget used in the snapshot modeling (see Equation (30)); the only difference is that now f and σ_ref enter into the model as free parameters.

The variability noise σ_var is described in Section 3 (see Equation (8)) and consists of a broken power law in ∣ u ∣ specified by four parameters: an overall amplitude a₄ specified at a baseline length of 4 Gλ, a falling long-baseline power-law index b, a rising short-baseline power-law index c, and a baseline length u₀ at which the power-law breaks. Informative prior bounds for each of these parameters are determined from the model-agnostic variability quantification described in Section 3.3, and these bounds are listed in Table 2.

Table 2. Full-track Modeling mG-ring Priors

Parameter	Prior
fring	${ \mathcal U }(0.05,4.0)$
d (μas)	${ \mathcal U }(20,85)$
W (μas)	${ \mathcal U }(1,40)$
∣βm ∣	${ \mathcal U }(0.0,0.5)$
$\arg ({\beta }_{m})$ (deg)	${ \mathcal U }(-180,180)$
fGauss	${ \mathcal U }(0.05,4.0)$
WGauss (μas)	${ \mathcal U }(20,200)$
σref	${{ \mathcal N }}_{L}(\mathrm{log}(0.004),1)$
f	${{ \mathcal N }}_{L}(\mathrm{log}(0.01),1)$
a4	${{ \mathcal U }}_{L}({10}^{-4.39784},{10}^{-4.27339})$
b	${ \mathcal U }(2.35213,3.37849)$
c	${ \mathcal U }(1.5,2.5)$
u0 (Gλ)	${{ \mathcal U }}_{L}({10}^{-1.09771},{10}^{0.236534})$

Note. Prior distributions for Themis full-track geometric modeling analyses. The top section lists priors for the mG-ring model parameters, and the bottom section lists priors for the parameterized noise model. ${ \mathcal U }(a,b)$ denotes a uniform prior on the interval [a, b], ${{ \mathcal U }}_{L}(a,b)$ denotes a log-uniform prior on the interval [a, b], and ${{ \mathcal N }}_{L}(\mu ,{\sigma }^{2})$ denotes a lognormal prior with mean μ and variance σ². Priors for the variability noise parameters a₄, b, c, and u₀ are informed by the model-agnostic variability quantification analysis described in Section 3.3.

Download table as:  ASCIITypeset image

7.2.2.  Themis Implementation

We have implemented the combined mG-ring plus noise parameterization as a model within the sampling-based parameter estimation framework Themis developed for the EHT (Broderick et al. 2020a, 2020b). Given a model specification and a data set, Themis works within a Bayesian formalism to produce a set of samples from the posterior distribution of the model parameters. Themis uses a MCMC sampling scheme to explore the posterior space, employing a parallel tempering scheme (Syed et al. 2022) to ensure traversal over the entire prior volume and the Hamiltonian Monte Carlo sampling kernel from the Stan package (Carpenter et al. 2017) to efficiently sample within each tempering level. A detailed description of the Themis sampling framework can be found in Tiede (2021).

The full-track geometric modeling analyses carried out in this paper fit to complex visibility data. Given a vector of geometric model parameters p and a vector of noise model parameters n , the Themis likelihood function for complex visibilities is Gaussian,

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}\sum _{i}\left(\displaystyle \frac{| {V}_{i}-\displaystyle \hat{V}({{\boldsymbol{u}}}_{i};{\boldsymbol{p}}){| }^{2}}{{\sigma }_{i}^{2}({\boldsymbol{n}})}+\mathrm{ln}\left[2\pi {\sigma }_{i}^{2}({\boldsymbol{n}})\right]\right).\end{eqnarray} \tag{ 46 }$$

Here V_i is a measured visibility, u _i is its baseline vector, $\displaystyle \hat{V}$ is the corresponding modeled visibility, and the sum is taken over all data points i. The noise σ_i ( n ) in each visibility is specified as in Equation (45), with the noise model parameters n = {f, σ_ref, a, b, c, u₀}. Themis internally solves for and marginalizes over the full set of complex gain parameters (i.e., one complex gain per station per time stamp) at every sampling step using a Laplace approximation (see Broderick et al. 2020a). It also applies the Johnson et al. (2018) diffractive scattering kernel directly to the model prior to computing visibilities. A validation test of the Themis mG-ring plus noise model implementation is described along with other tests in a dedicated paper on the noise modeling approach (Broderick et al. 2022).

We assess MCMC convergence through both visual inspection of the traces and a number of quantitative chain statistics, including the integrated autocorrelation time, split-$\displaystyle \hat{R}$, and parameter rank distributions (Vehtari et al. 2021). The number of tempering levels is selected to ensure efficient communication between the highest and lowest levels (per Syed et al. 2022), which typically requires about 20 levels. We run the sampler for between 5 × 10⁴ and 10⁵ steps per tempering level.

To compute the Bayesian evidence, Themis uses thermodynamic integration (e.g., Lartillot & Philippe 2006), which computes the log-evidence through

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal Z }={\int }_{0}^{1}{\rm{d}}\beta {\left\langle \mathrm{ln}{ \mathcal L }\right\rangle }_{\beta },\end{eqnarray} \tag{ 47 }$$

where ${\left\langle \mathrm{ln}{ \mathcal L }\right\rangle }_{\beta }$ is the expectation of the log-likelihood taken over the distribution

$$\begin{eqnarray}&&{P}^{(\beta )}({\boldsymbol{p}},{\boldsymbol{n}})={\left(\displaystyle \frac{P({\boldsymbol{p}},{\boldsymbol{n}})}{{\pi }_{\mathrm{ref}}({\boldsymbol{p}},{\boldsymbol{n}})}\right)}^{\beta }{\pi }_{\mathrm{ref}}({\boldsymbol{p}},{\boldsymbol{n}}).\end{eqnarray} \tag{ 48 }$$

Note that Themis does not take π_ref to be the prior distribution. Instead, Themis uses a uniform distribution whose support matches the support of the priors given in Table 2. To compute Equation (47), we compute the average log-likelihood for each tempering level and then use trapezoidal integration to numerically compute the integral.

Priors for all mG-ring model and noise model parameters are listed in Table 2. We impose mean-zero lognormal priors for all station gain amplitudes, with a log standard deviation of 0.01 for all network calibrated stations (ALMA, APEX, JCMT, SMA), 0.2 for the LMT, and 0.1 for the remaining stations (PV, SMT, SPT); these gain priors are motivated by the expected performance of each station after the post-processing described in Section 7.1 (see also Paper II). All station gain phase priors are uniform on the unit circle.

As with the snapshot geometric analyses (see Section 6.5), the mG-ring model used for the full-track analyses is really a class of models that increases in complexity with m. The results of a Themis m-order survey covering m = {2, 3, 4, 5, 6} are shown in Figure 16. In contrast to the snapshot geometric modeling (see Section 6.5), we find that the full-track analysis is able to support more complex model specifications, exhibiting a strong preference for m > 4 over m = 4. However, we find that increasing the m-order does not significantly impact the values of the primary morphological parameters of choice (see the similar diameter measurements and uncertainties for m ≥ 4 in Figure 16). Thus, to maintain consistency among model specifications and to facilitate comparison with the snapshot geometric modeling analyses, we proceed with the m = 4 mG-ring as our fiducial model for all full-track Sgr A* analyses in this paper.

Zoom In Zoom Out Reset image size

A representative m = 4 mG-ring fit to the Sgr A* HOPS low-band data is shown in Figure 17. We find that the normalized residuals are distributed around a value of zero with near-unity variance and that there is no evidence of systematic structure. The χ² statistics for this fit are discussed in Appendix E.

Zoom In Zoom Out Reset image size

The model-agnostic variability quantification analysis carried out in Section 3 demonstrates that the Sgr A* data exhibit variability—quantified here in terms of a normalized visibility amplitude variance—that is significantly in excess of that expected from thermal noise, station gains, and refractive scattering. As illustrated in Figure 4, the measured variability can be broken down into three regions with qualitatively distinct behavior:

• 1.  
  On short baselines with lengths ∣ u ∣ ≲ 2.5 Gλ, corresponding to spatial scales ≳100 μas, limitations in our calibration and the subsequent choices made in our preprocessing procedure preclude meaningful constraints on the variability. The light-curve normalization procedure removes all variability on intrasite baselines and suppresses variability on short intersite baselines that are highly correlated with the light curve; the variability of the light curve itself is thoroughly characterized in Wielgus et al. (2022). The source size constraint used to perform gain calibration of the LMT-SMT baseline (see Paper II; Paper III) imposes a further, more artificial suppression of the variability on this baseline. We thus do not obtain any variability measurements for baselines shorter than 2.5 Gλ.
• 2.  
  On intermediate baselines with lengths between 2.5 Gλ ≲ ∣ u ∣ ≲ 6 Gλ, corresponding to spatial scales between ∼30 and ∼100 μas, we measure significant variability that exhibits an approximately power-law decline with increasing baseline length. The power-law index is between ∼2 and 3, and the magnitude of the variability ranges from a peak rms of ∼5% of the total flux density (∼120 mJy) near 2.5 Gλ down to ∼1% of the total flux density (∼25 mJy) near 6 Gλ.
• 3.  
  On long baselines with lengths ∣ u ∣ ≳ 6 Gλ, corresponding to spatial scales ≲30 μas, the measured variability is comparable in magnitude to that expected from statistical errors and refractive scattering. These measurements thus do not contain statistically significant detections of structural variability.

These measurements describe the level of excess variance that the data exhibit about an underlying average source model. The parameters describing a broken power-law noise model fit to these measurements are thus used as a variability noise budget during image reconstruction (Paper III) and to define priors on the corresponding parameters in the full-track modeling analyses (Section 7).

Determining the intrinsic (i.e., infinite-time) source variability from these measurements requires an additional debiasing step to remove the impact of correlations between data points that are closely spaced in time. The analysis carried out in Section 3 involves binning the visibility data in the (u,v)-plane for the purpose of computing variances. However, many data points within a single bin are from measurements taken close in time, which can introduce correlations that bias the computed variance. A procedure for removing this bias is detailed in Broderick et al. (2022), whereby the factor relating the measured variability to the intrinsic variability at every baseline length is calibrated using synthetic measurements of GRMHD simulations. In practice, Broderick et al. (2022) derive the debiasing function using the same 90 GRMHD simulations we use in this paper for θ_g calibration (see Section 4.4 and Appendix D), and the resulting debiasing factor is close to unity everywhere except between ∼6 and 7.5 Gλ (see Figure 5 in Broderick et al. 2022). Applying this debiasing function to the variability measurements from the Sgr A* data yields the results shown in Figure 18 and reported in Table 3.

Zoom In Zoom Out Reset image size

Due to the near-unit debiasing factor, the variability measurements shown in Figure 18 are similar to those from Figure 4. Quantitative constraints on the parameters of the noise model that are well constrained are presented in Figure 19 and Table 3. Where strong constraints on the excess noise exist (i.e., on baseline lengths between ∼2 and 6 Gλ), it continues to be well described by a single power law with index $b={2.4}_{-0.8-1.6}^{+0.8+2.1}$ and amplitude (evaluated at a baseline length of 4 Gλ) of ${a}_{4}={1.9}_{-0.2-0.3}^{+0.2+0.5} \%$, where in each value the 1σ and 2σ ranges are indicated. The measured variability magnitude is between ∼2 and 10 times higher than that expected from refractive scattering alone. The lack of an observable break places an upper limit on its location of u₀ < 1.3 Gλ at 1σ and u₀ < 3.1 Gλ at 2σ.

Zoom In Zoom Out Reset image size

The excess variability is broadly consistent with that due to structural fluctuations anticipated by the GRMHD simulations discussed in Paper V and Georgiev et al. (2022). The magnitude of the excess variability lies within the range of that predicted by GRMHD simulations for ∣ u ∣ > 2 Gλ, though it does appear to marginally favor less variable models. The long-baseline power-law index is consistent with all GRMHD simulations. A detailed discussion of the implications for GRMHD models is contained within Paper V.

Both the snapshot (Section 6) and full-track (Section 7) geometric modeling analyses produce reconstructions of the Sgr A* emission structure, and the posterior distributions determined for the parameters describing these geometric model reconstructions provide a quantification of the morphological properties directly from the EHT interferometric data. Similarly, the IDFE analyses carried out in Section 5 quantify the morphological properties of the top-set and posterior images reconstructed in Paper III.

Figure 20 compares the geometric models and image reconstructions determined for Sgr A*; for the geometric modeling analyses we show posterior mean images (i.e., the mean of many images sampled from the posterior distribution), while for the image reconstructions we show averages over the top sets (for eht-imaging, SMILI, and DIFMAP) or posterior means (for Themis). We see that both the snapshot and full-track geometric modeling analyses each recover a grossly similar overall structure across frequency bands and that this structure is also similar between the snapshot and full-track analyses. The image reconstructions permit much more flexibility in the permitted image structure, and so we see correspondingly more variation both within the imaging methods and between the imaging and geometric modeling. The primary point of consistency between the ring-like structures recovered from imaging and the rings fit via geometric modeling seems to be their sizes.

Zoom In Zoom Out Reset image size

We use the geometric modeling and IDFE analyses to quantify a number of morphological parameters of interest, which are shown in Figure 21 and listed in Table 4.

Zoom In Zoom Out Reset image size

Table 4. Parameters Describing Sgr A* Image Morphology

				CASA						HOPS
Analysis class	Software	Day	Band	$\displaystyle \hat{d}$ (μas)	W (μas)	A	η (deg)	fc		$\displaystyle \hat{d}$ (μas)	W (μas)	A	η (deg)	fc
Imaging	DIFMAP+ REx	April 7	HI+LO	${50.6}_{3.0}^{+1.1}$	${31.9}_{3.1}^{+4.2}$	${0.09}_{0.04}^{+0.07}$	48.4±57.1	${0.43}_{0.10}^{+0.12}$		${50.1}_{2.2}^{+1.2}$	${31.8}_{2.7}^{+3.3}$	${0.05}_{0.02}^{+0.09}$	−22.8±69.6	${0.42}_{0.09}^{+0.10}$
	DIFMAP+ VIDA	April 7	HI+LO	${51.2}_{2.8}^{+1.2}$	${33.4}_{1.8}^{+4.1}$	${0.10}_{0.04}^{+0.11}$	49.8±61.1	${0.40}_{0.09}^{+0.14}$		${51.3}_{1.4}^{+1.3}$	${33.1}_{1.1}^{+2.2}$	${0.07}_{0.03}^{+0.18}$	−28.1±57.1	${0.39}_{0.08}^{+0.15}$
	eht-imaging+ REx	April 7	HI+LO	${53.0}_{2.5}^{+1.9}$	${25.7}_{3.2}^{+3.6}$	${0.13}_{0.06}^{+0.08}$	−33.9±56.3	${0.33}_{0.18}^{+0.12}$		${53.6}_{1.1}^{+1.8}$	${27.2}_{2.6}^{+1.6}$	${0.11}_{0.03}^{+0.05}$	−81.0±62.6	${0.21}_{0.10}^{+0.11}$
	eht-imaging+ VIDA	April 7	HI+LO	${52.7}_{2.1}^{+5.6}$	${28.0}_{3.2}^{+16.3}$	${0.19}_{0.11}^{+0.26}$	−40.5±65.2	${0.30}_{0.16}^{+0.22}$		${53.5}_{1.5}^{+3.6}$	${27.4}_{2.1}^{+2.5}$	${0.11}_{0.05}^{+0.23}$	−94.4±60.9	${0.19}_{0.05}^{+0.13}$
	SMILI+ REx	April 7	HI+LO	${52.0}_{3.0}^{+5.5}$	${25.1}_{4.4}^{+1.1}$	${0.16}_{0.06}^{+0.04}$	130.9±70.3	${0.29}_{0.20}^{+0.09}$		${51.4}_{4.8}^{+4.0}$	${26.7}_{3.0}^{+1.2}$	${0.12}_{0.03}^{+0.04}$	135.3±90.0	${0.23}_{0.13}^{+0.08}$
	SMILI+ VIDA	April 7	HI+LO	${49.3}_{1.7}^{+13.6}$	${28.5}_{3.7}^{+34.0}$	${0.15}_{0.04}^{+0.21}$	124.4±77.0	${0.26}_{0.07}^{+0.48}$		${50.6}_{1.1}^{+6.1}$	${27.2}_{1.3}^{+0.8}$	${0.07}_{0.01}^{+0.11}$	163.0±90.4	${0.19}_{0.05}^{+0.09}$
	Themis + REx	April 6+7	HI+LO	...	...	...	...	...		${51.7}_{0.3}^{+0.4}$	${23.6}_{0.3}^{+0.3}$	${0.07}_{0.01}^{+0.08}$	−152.7±20.9	${0.02}_{0.01}^{+0.01}$
	Themis + VIDA	April 6+7	HI+LO	...	...	...	...	...		${53.7}_{1.7}^{+0.4}$	${24.7}_{0.4}^{+0.4}$	${0.08}_{0.01}^{+0.05}$	−137.1±15.7	${0.08}_{0.01}^{+0.01}$
Snapshot	Comrade	April 6+7	HI	${53.9}_{0.5}^{+0.6}$	${16.4}_{0.2}^{+0.3}$	${0.23}_{0.02}^{+0.01}$	3.8±5.3	${0.23}_{0.01}^{+0.01}$		${52.3}_{0.6}^{+0.7}$	${17.7}_{0.2}^{+0.2}$	${0.23}_{0.04}^{+0.03}$	−8.1±6.7	${0.15}_{0.01}^{+0.01}$
		April 6+7	LO	${51.4}_{0.6}^{+0.5}$	${17.0}_{0.3}^{+0.2}$	${0.22}_{0.03}^{+0.02}$	−11.4±9.0	${0.20}_{0.01}^{+0.01}$		${53.9}_{0.7}^{+0.8}$	${17.7}_{0.2}^{+0.2}$	${0.22}_{0.02}^{+0.02}$	−0.2±6.6	${0.18}_{0.01}^{+0.01}$
	DPI	April 6+7	HI	${50.9}_{1.0}^{+1.0}$	${16.1}_{0.7}^{+0.7}$	${0.21}_{0.14}^{+0.18}$	−75.6±43.2	${0.18}_{0.02}^{+0.02}$		${47.7}_{1.1}^{+1.1}$	${19.1}_{0.7}^{+0.7}$	${0.31}_{0.20}^{+0.13}$	−44.9±11.7	${0.12}_{0.01}^{+0.02}$
		April 6+7	LO	${51.6}_{0.9}^{+0.9}$	${16.2}_{0.7}^{+0.7}$	${0.14}_{0.10}^{+0.07}$	−97.5±9.5	${0.14}_{0.01}^{+0.02}$		${45.3}_{1.7}^{+1.6}$	${21.6}_{0.9}^{+0.9}$	${0.31}_{0.19}^{+0.13}$	−73.2±8.2	${0.10}_{0.01}^{+0.01}$
Full-track	Themis	April 6+7	HI	${52.1}_{0.4}^{+0.4}$	${19.7}_{0.4}^{+0.4}$	${0.16}_{0.02}^{+0.02}$	−4.0±11.8	${0.09}_{0.02}^{+0.02}$		${52.0}_{0.5}^{+0.4}$	${21.5}_{0.5}^{+0.5}$	${0.14}_{0.02}^{+0.02}$	8.9±8.0	${0.05}_{0.02}^{+0.02}$
		April 6+7	LO	${52.7}_{0.5}^{+0.5}$	${20.3}_{0.4}^{+0.4}$	${0.15}_{0.02}^{+0.02}$	−74.0±10.5	${0.16}_{0.03}^{+0.03}$		${51.4}_{0.5}^{+0.4}$	${22.5}_{0.5}^{+0.4}$	${0.14}_{0.02}^{+0.02}$	1.6±7.5	${0.03}_{0.01}^{+0.02}$

Note. Median values and 68% credible intervals for the morphological quantities of interest, measured from the EHT Sgr A* data. Because medians and quantiles are not well defined for angular variables, for the position angle η we instead quote the circular mean and standard deviation.

Download table as:  ASCIITypeset image