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  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 and 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 .

We note that  differs from the ring identification methods presented in  in that it searches for the presence of a topological ring in the input image. <ref> compares the mean ring and nonring descattered images for each imaging pipeline as classified by  in the “permissive” mode and the clustering analysis from . 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.

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 and are robust against the specifics of the ring identification scheme used. <ref> shows the resulting diameter distributions from ring fitting to the descattered images from all imaging pipelines, split out by  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  penalizes images with emission structures deviating from a topological ring, the distributions of the and measurements are robust against the choice of  mode employed.

Diameter distributions determined by and for all descattered images from the HOPS pipeline, organized by  mode. Each panel shows the fraction of images found to possess topological ring structure. For posterior imaging we use the combined April 6+7 results, and for top-set imaging we use the April 7 results.

§ SNAPSHOT GEOMETRIC MODELING

Because the data are observed to be time-variable (see <ref>), a static model cannot reproduce the observed data.  As described in <ref>, one method for mitigating the effects of this variability on the reconstructed source structure is through the use of an inflated variability noise budget, as pursued during the imaging (; <ref>) and full-track geometric modeling (<ref>) analyses.  In this section we instead pursue “snapshot” geometric modeling, whereby we fit a geometric model – the model described in <ref> – for which the parameters are allowed to vary as a piecewise constant function of time. To this end, we divide the data into many independent and short “snapshots” over which the source is assumed to be static. In this section, we detail our formalism for fitting the model to snapshots of data and for combining the fits from across snapshots into a global posterior distribution.

§.§ Data preparation

Prior to fitting the model to real or synthetic data, we process the data using the pre-imaging pipeline described in . 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 .  Specifically, the total error budget $\sigma_{sb}$ for a visibility measured on the baseline $b$ during snapshot $s$ is given by

\begin{equation}
\sigma_{sb}^2 = \sigma_{\text{th},sb}^2 + f^2 |V|_{sb}^2 + \sigma_{\text{ref}, sb}^2 . \label{eqn:SnapshotNoise}
\end{equation}

Here the first term corresponds to the baseline-specific thermal noise (see <ref>), 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  refractive scattering noise from .
For the snapshot modeling, we fix $f=0.02$ per the analyses carried out in . The pre-imaging pipeline also mitigates the impact of diffractive scattering by “deblurring” the data using the <cit.> 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.

§.§ Snapshot fitting procedure

The first step of our snapshot modeling procedure is to determine the posterior distribution for the 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 model described in <ref>, whose parameter vector we denote as $\boldsymbol{\theta}_s$.  For a single snapshot, the posterior is given by Bayes's theorem,

\begin{equation}\label{eqn:SnapshotBayes}
\begin{aligned}
P_s(\boldsymbol{\theta}_s | \boldsymbol{D}_s)
&= \frac{\mathcal{L}_s(\boldsymbol{D}_s | \boldsymbol{\theta}_s)\pi_s(\boldsymbol{\theta}_s)}{\mathcal{Z}_s(\boldsymbol{D}_s)} ,
\end{aligned}
\end{equation}

where $\boldsymbol{D}_s$ denotes the data available on snapshot $s$, $\mathcal{L}_s$ is the likelihood, $\pi_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 $\ln A$, and closure phases $\psi$.  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{equation}
\mathcal{L}_s = \mathcal{L}_{|V|,s} \mathcal{L}_{\psi, s} , \label{eqn:ShapshotLikelihoodVisamp}
\end{equation}

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

\begin{equation}
\mathcal{L}_s = \mathcal{L}_{A, s} \mathcal{L}_{\psi, s} . \label{eqn:ShapshotLikelihoodClosures}
\end{equation}

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 <ref>.

§.§ Averaging the snapshot results

The output of a snapshot fitting analysis is a set of posterior samples for the model parameters from each individual snapshot; <ref> shows an example set of posterior distributions for the 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 <cit.>. This approach treats the model fit to each snapshot as a realization from some average model or “hypermodel.”

Example snapshot modeling results and averaging scheme applied to the April 6 and 7 low-band HOPS data sets. The blue filled regions show the posterior distribution of the diameter parameter for an $m=4$ model fit to each 120 s snapshot. We find the tightest diameter posteriors from 12.5 to 14 UTC, which corresponds to the best time region from . For visual clarity we only show the distributions for every second snapshot. The black curves at the top show the diameter distribution corresponding to 100 random draws from the hypermodel posterior (<ref>).

§.§.§ Averaging procedure

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

\begin{equation}\label{eq:snapshot:topset}
\begin{aligned}
P(\boldsymbol{\Theta},\bm{\bar{\theta}}| \bm{\mathcal{D}})
&= \frac{\pi(\bm{\bar{\theta}})}{\mathcal{Z}(\bm{\mathcal{D})}}\prod_s \mathcal{L}_s(\bm{D}_s|\bm{\theta}_s) \pi_s(\bm{\theta}_s|\bm{\bar{\theta}}),
\end{aligned}
\end{equation}

where $\pi(\bm{\bar{\theta}})$ is the prior distribution for the hypermodel parameters (the “hyperprior”), $\boldsymbol{\Theta} = (\boldsymbol{\theta}_1, \boldsymbol{\theta}_2, \dots, \boldsymbol{\theta}_{N_s})$ is the parameter vector across all snapshots, and $\boldsymbol{\mathcal{D}} = (\boldsymbol{D}_1, \boldsymbol{D}_2, \dots, \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{equation}\label{eq:snapshot:topset_marginal}
\begin{aligned}
P(\bm{\bar{\theta}}| \bm{\mathcal{D}})
&= \frac{\pi(\bm{\bar{\theta}})}{\mathcal{Z}(\bm{\mathcal{D}})}\int \left( \prod_s \mathcal{L}_s(\bm{D}_s|\bm{\theta}_s) \pi_s(\bm{\theta}_s|\bm{\bar{\theta}}) \right) \mathrm{d}\boldsymbol{\Theta} .
\end{aligned}
\end{equation}

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 <ref> and use Bayes's theorem to substitute in for the snapshot likelihood (<ref>), giving

\begin{equation}\label{eq:snapshot:avgposterior}
P(\bm{\bar{\theta}}| \bm{\mathcal{D}})
= \pi(\bm{\bar{\theta}}) \prod_s \left( \int  P(\bm{\theta}_s|\bm{D}_s)\frac{\pi_s(\bm{\theta}_s|\bm{\bar{\theta}})}{\pi_s(\bm{\theta}_s)} \mathrm{d}\bm{\theta}_s \right) .
\end{equation}

Note that the evidence term from the prefactor denominator in <ref> has now been subsumed into the posterior term $P(\bm{\theta}_s|\bm{D}_s)$ inside of the integral. To evaluate <ref>, we make use of the fact that the snapshot posterior samples $\bm{\theta}_s^{(i)}$ permit us to approximate the integral by a sum,

\begin{equation}\label{eq:snapshot:avgposterior_approximate}
P(\bm{\bar{\theta}}| \bm{\mathcal{D}})
\approx \pi_s(\bm{\bar{\theta}})\prod_s \sum_{i} \frac{\pi_s(\bm{\theta}^{(i)}_s|\bm{\bar{\theta}})}{\pi(\bm{\theta}^{(i)}_s)} .
\end{equation}

We can use this expression to sample from the posterior distribution over just the hypermodel parameters $\bm{\bar{\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(\bm{\theta}_s|\bm{D}_s)$ are Gaussian and the priors on $\bar{\bm{\theta}}$ and $\bm{\theta}_s$ are uninformative, the posterior maximum for $P(\bm{\bar{\theta}}| \bm{\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 <ref>), we proceed with the more general averaging procedure.

§.§.§ Hypermodel specification

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

\begin{equation}
\pi_{s, \text{tN}}(\theta_{s,i} | \mu_i, \sigma_i) =  \mathcal{T}(\theta_{s,i} | \mu_i, \sigma_i;\, a_i, b_i) . \label{eqn:hypermodel:nonangular}
\end{equation}

Here $\mathcal{T}(\theta | \mu, \sigma;\, a, b)$ denotes the density for a truncated normal distribution with mean $\mu$ and standard deviation $\sigma$, 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\pi)$ – we instead use a von Mises distribution,

\begin{equation}
\pi_{s, \text{vM}}(\theta_{s,j} | \mu_j, \sigma_j) = \frac{1}{2\pi I_0(\sigma_j^{-2})}\exp\left[\frac{\cos(\theta_{s,j} - \mu_j)}{\sigma_j^2}\right] . \label{eqn:hypermodel:angular}
\end{equation}

The subscripted “tN” in <ref> indicates that the corresponding parameters use a truncated normal prior, while the subscripted “vM” in <ref> similarly indicates that the corresponding parameters use a von Mises prior.
The total hypermodel is then given by

\begin{equation}
\pi(\bm{\theta}_s | \bm{\mu}, \bm{\sigma}) = \prod_{i} \pi_{s, \text{tN}}(\theta_{s,i} | \mu_i, \sigma_i) \prod_{j}\pi_{s, \text{vM}}(\theta_{s,j} | \mu_j, \sigma_j) ,
\end{equation}

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

We set the hyperpriors for $\bm{\mu}$ to be equal to the corresponding snapshot priors, which are specified in <ref>. For $\bm{\sigma}$ we instead use a half-normal hyperprior,

\begin{equation}
\pi(\bm{\sigma}) =  \prod_i \mathcal{T}(\sigma_i| 0, L_i/4;\; 0, L_i),
\end{equation}

where $L_i = (b_i - a_i)$ is the breadth of support for parameter $\theta_{s,i}$.
<ref> describes the level of consistency between these selected hyperpriors and the priors for the individual snapshot model parameters.

Snapshot modeling priors
$F_0$     $\mathcal{U}(0.8, 1.2)$     $\delta(1)$

$d$ ()     $\mathcal{U}(25, 85)$     …

$d^{\prime}$ ()     …    $\mathcal{U}(25, 85)$

$W$ ()     $\mathcal{U}(1, 40)$      …

$W^{\prime}$ ()     …    $\mathcal{U}(1, 40)$

$|\beta_m|$     $\mathcal{U}(0, 0.5)$     $\mathcal{U}(0, 0.5)$

$\arg(\beta_m)$ (deg)     $\mathcal{U}(-180, 180)$     $\mathcal{U}(-180, 180)$

$f_{\rm Gauss}$     $\mathcal{U}(0, 1)$     $\mathcal{U}(0, 1)$

$W_{\rm Gauss}$ ()     $\mathcal{U}(40, 200)$     $\mathcal{U}(40, 200)$

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

§.§ Software implementations

We use three different software packages to carry out snapshot geometric modeling on the 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 <ref>.

Our primary snapshot fitting software is the modeling framework <cit.>, which is written in the dynamic programming language Julia <cit.>. does not natively include functionality for constructing a joint probability describing both observations and model. Instead, it interfaces with existing probabilistic programming languages present in Julia. For the analyses presented in this paper, we use the probabilistic programming package Soss[<https://github.com/cscherrer/Soss.jl>] to construct the joint probability.  This interface is specified in the package ComradeSoss.jl[<https://github.com/ptiede/ComradeSoss.jl>]. To sample from the posterior, we use the nested sampling package , which also produces estimates of the Bayesian evidence <cit.>.

Given a model specification, can fit a variety of interferometric data products, including visibility amplitudes, closure phases, and log closure amplitudes. Unless otherwise specified, for the snapshot modeling analyses performed in this paper, we use to fit to visibility amplitudes and closure phases to the mG-ring model; the snapshot likelihood is thus given by <ref>.  Prior to fitting, all time stamps that contain fewer than four baselines are flagged.

When fitting to visibility amplitudes, we include the station gain amplitudes as model parameters alongside the geometric parameters that describe the model. For the gain amplitudes we use a lognormal prior with a log-mean of zero (i.e., corresponding to unit gain amplitude) and a log standard deviation of 0.1 on all stations except for LMT, for which we use a log standard deviation of 0.2 to accommodate its larger variations .

We also utilize the geometric model-fitting tools developed within the Python library <cit.>. This library enables visibility-domain fitting to arbitrary combinations of simple analytic models, including the model, and it can do so using a variety of interferometric data products, including visibility amplitudes, closure phases, and log closure amplitudes.  is also able to interface with a variety of external packages to perform parameter optimization or posterior exploration.

For the snapshot modeling analyses performed in this paper, we match the operation of with that of .  Unless otherwise specified, we use to fit to visibility amplitudes and closure phases, so that the snapshot likelihood is given by <ref>, and we use  <cit.> for posterior exploration and evidence estimation. We also specify the same priors for the station gain parameters as used in the fits.

Given that both and use identical model specifications, priors, and samplers, we expect all results produced by these softwares to be identical up to sampling precision.  We thus use only fits for all snapshot geometric analyses in this paper.

The third software we use for snapshot geometric modeling is the Python code Deep Probabilistic Imaging/Inference <cit.>.
approximates the posterior over all model parameters by fitting a normalizing flow neural network <cit.> to the data using a Rényi $\alpha$-divergence variational inference technique <cit.>.
is an optimization-based posterior estimation framework, and it uses the auto-differentiation package PyTorch <cit.> 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.

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  to fit to closure phases and log closure amplitudes; the snapshot likelihood is given by <ref>. Prior to fitting, all time stamps that are unable to form at least one closure phase and at least one closure amplitude are flagged.

differs from both and 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 \times 32$ pixels spanning a 160 field of view.

Because the pixel size is finite,  cannot support a model containing infinitesimally thin rings such as that in <ref>; furthermore, convolutions in the image domain are computationally expensive.  The  fits in this paper thus employ a modified version of the model specification,

\begin{equation} \label{eq:dpi-mGring}
\begin{split}
I_{\text{ring}}^{\prime}(r, \phi) &= F_{\text{ring}} \frac{\exp\left[ \frac{4\ln{2}}{{W^{\prime}}^2} \left( r-\frac{d^{\prime}}{2} \right)^2 \right] \sum_{k=-m}^{m} \beta_k e^{ik\phi}}{2 \pi \int_0^{\infty} r \exp\left[ \frac{4\ln{2}}{{W^{\prime}}^2} \left( r-\frac{d^{\prime}}{2} \right)^2 \right] \,dr}, \\
I^{\prime}(r, \phi) &= I_{\text{ring}}^{\prime}(r, \phi) + I_{\text{Gauss}}(r, \phi) , \\
F_0 &= F_{\text{ring}} + F_{\text{Gauss}} = 1 ,
\end{split}
\end{equation}

where we note that $d^{\prime}$ and $W^{\prime}$ are conceptually distinct from $d$ and $W$.  The quantities $d$ and $W$ in the model from <ref> 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  model implementation is fixed to be 1 Jy because  fits only to closure quantities and closure amplitudes are not sensitive to the absolute flux scale.

§.§.§ Sampling the hypermodel posterior

To sample from the hypermodel posterior $P(\bm{\bar{\theta}}| \bm{\mathcal{D}})$ (<ref>), we use the adaptive Metropolis sampler from <cit.> via its implementation in the Julia package RobustAdaptiveMetropolisSampler.jl[<https://github.com/anthofflab/RobustAdaptiveMetropolisSampler.jl>]. The sampler is initialized by first running an adaptive genetic algorithm from the Julia package BlackBoxOptim.jl[<https://github.com/robertfeldt/BlackBoxOptim.jl>], which provides a starting point near the maximum posterior density. We run the sampler for a minimum of 2 million Markov Chain Monte Carlo (MCMC) steps or until we have effective sample sizes of 500 for all parameters.

§.§ Model selection

Relative Bayesian evidence and best-fit model diameter vs. $m$-order, from the snapshot+averaging geometric modeling from the pipeline applied to the data.  Because the absolute values of $\ln\mathcal{Z}$ can be substantially different for each data set (i.e., each combination of observing day, frequency band, and calibration pipeline), and because only the relative values carry information about model specification preferences, we reference all $\ln\mathcal{Z}$ values to the maximum value attained at any $m$ for that data set.

Representative examples of snapshot modeling results for both the (top row) and  (bottom row) pipelines.  The top row of panels shows results from fitting an $m=4$ with to the HOPS low-band data on both April 6 and 7, while the bottom row of panels shows results from fitting an $m=2$ with  to the same dataset.  Each panel is arranged analogously to the individual panels of <ref>, though the plotted data products are different.  For the results, the visibility amplitudes and closure phases are used during fitting, so these are the data products shown in the top left and top right panels, respectively.  For the  results, the log closure amplitudes and closure phases are used during fitting, so these are the data products shown in the bottom left and bottom right panels, respectively.  All closure quantities are plotted as a function of the perimeter of the relevant polygon (i.e., triangles for closure phases, quadrangles for log closure amplitudes); the evident zero-valued closure phases and log closure amplitudes primarily correspond to so-called “trivial” polygons, i.e., those with near-zero area .  We note that the data-flagging procedures differ slightly between the two pipelines (see Sections <ref> and <ref>), resulting in small differences in the fitted data sets.

The model described in <ref> is not actually a single model but rather a class of models, delineated by the order $m$ (<ref>). To determine the $m$-order that is preferred by the data, we carry out a series of snapshot fits to the using different values of $m$ and compare Bayesian evidence estimates.  Given a set of log-evidences $\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{equation}
\ln \mathcal{Z} = \sum_s \ln \mathcal{Z}_s . \label{eqn:EvidenceCombination}
\end{equation}

The 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 <ref>.  The results of a $m$-order survey covering $m=\{1,2,3,4,5\}$ are shown in <ref>.  We find that the $m=4$ order is preferred in both bands and across both calibration pipelines.  We thus use the $m=4$ as our fiducial model for all analyses in this paper.

Unlike ,  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{equation}
\label{eq:elbo}
\text{ELBO}(m) = \log p(\bm{D}_s \| m ) - D_{KL}[q(\bm{\theta}_s)\|p(\bm{\theta}_s|\bm{D}_s, m)],
\end{equation}

where $D_{\rm{KL}}[A\|B]$ is the Kullback-Leibler divergence of $A$ from $B$, and $q(\bm{\theta}_s)$ is the optimized  normalizing flow distribution.
The relative entropy term is zero when the  distribution $q(\bm{\theta}_s)$ and the true posterior $p(\bm{\theta}_s|\bm{D}_s)$ are identical, so the ELBO provides a rough estimate of the log-evidence. The results of a  $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$ as our fiducial model for all  analyses in this paper.

<ref> shows representative fits to the HOPS low-band data for both the and  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 $\chi^2$ statistics for each of these fits are provided in <ref>.

§ FULL-TRACK GEOMETRIC MODELING

The snapshot modeling analysis presented in the previous section addresses the variability of the data by explicitly permitting the source structure to vary in time.  As described in <ref>, an alternative approach to fitting variable data is to statistically capture the impact of variability, treating it as an additional source of uncertainty modifying data that otherwise describe a static (or average) source structure.  We pursue such an approach here in the form of “full-track” geometric modeling, whereby we fit the model (see <ref>) to an entire data set at once and account for the variability by simultaneously fitting a parameterized noise model.  In this section we detail our formalism for fitting the geometric model alongside a model that captures the noise budget inflation associated with source variability.

§.§ Data preparation

The data preparation for the full-track geometric modeling analyses is similar to that used for snapshot geometric modeling analyses (see <ref>).  The data are first processed through the pre-imaging pipeline described in , which applies light-curve normalization and performs some a priori gain calibration.  However, unlike in  and <ref>, 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 ($\sim$10 minutes) is approximately the amount of time over which we expect structural variability to be subdominant to other sources of uncertainty (see <ref>, in particular <ref>).  Furthermore, the station gains are expected to be constant in time across a single scan but not from one scan to the next , 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 source structure.  As a consequence, the data exhibit subhour correlations that over a single day are localized in the -plane.  As previously noted in <ref> and detailed in <ref>, these unmodeled correlations can result in significant biases in the reconstructed properties of .  However, by fitting to multiple days of data, and thus combining multiple samplings of the variable source structure at each location in the -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 <ref> make use of the combined April 6 and April 7 data.  For comparison we provide single-day analysis results in <ref>.

§.§ Model specification and implementation

The goal of the full-track geometric modeling procedure is to determine the posterior distribution for the parameters of the static model and parameterized noise model that best describe an entire data set.  Our specification for the model is described in <ref>, 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 , through multiplication of the 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.

§.§.§ Parameterized noise model

Our parameterized noise model for a complex visibility $V_{i}$ measured on a baseline $\boldsymbol{u}_i$ is given by

\begin{equation}
\sigma_{i}^2 = \sigma_{\text{th},i}^2 + f^2 |V|_{i}^2 + \sigma_{\text{ref}}^2 + \sigma_{\text{var}}^2(|\boldsymbol{u}_i|) . \label{eqn:FullTrackNoise}
\end{equation}

Here the first term is the thermal noise in the measurement (see <ref>), the second term is a component that is multiplicative in the visibility amplitude $|V|_{i}$ and which is intended to capture residual (nongain) calibration errors (e.g. residual polarization leakage), the third term is a component that is additive and which is intended to account for refractive scattering noise, and the fourth term is a component that is a function of the baseline length $|\boldsymbol{u}_i|$ and which is intended to capture the effects of source variability.  With the exception of the variability term, <ref> is similar to the noise budget used in the snapshot modeling (see <ref>); the only difference is that now $f$ and $\sigma_{\text{ref}}$ enter into the model as free parameters.

The variability noise $\sigma_{\text{var}}$ is described in <ref> (see <ref>) and consists of a broken power law in $|\boldsymbol{u}|$ specified by four parameters: an overall amplitude $a_4$ specified at a baseline length of 4 G$\lambda$, a falling long-baseline power-law index $b$, a rising short-baseline power-law index $c$, and a baseline length $u_0$ at which the power-law breaks.  Informative prior bounds for each of these parameters are determined from the model-agnostic variability quantification described in <ref>, and these bounds are listed in <ref>.

Full-track modeling priors
Parameter     Prior
$f_{\text{ring}}$     $\mathcal{U}(0.05, 4.0)$

$d$ ()     $\mathcal{U}(20, 85)$

$W$ ()     $\mathcal{U}(1, 40)$

$|\beta_m|$     $\mathcal{U}(0.0, 0.5)$

$\arg(\beta_m)$ (deg)     $\mathcal{U}(-180, 180)$

$f_{\rm Gauss}$     $\mathcal{U}(0.05, 4.0)$

$W_{\rm Gauss}$ ()     $\mathcal{U}(20, 200)$

$\sigma_{\text{ref}}$     $\mathcal{N}_L(\log(0.004), 1)$

$f$     $\mathcal{N}_L(\log(0.01), 1)$

$a_4$     $\mathcal{U}_L(10^{-4.39784},10^{-4.27339})$

$b$     $\mathcal{U}(2.35213,3.37849)$

$c$     $\mathcal{U}(1.5, 2.5)$

$u_0$ (G$\lambda$)     $\mathcal{U}_L(10^{-1.09771},10^{0.236534})$

Prior distributions for full-track geometric modeling analyses.  The top section lists priors for the 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 $\mu$ and variance $\sigma^2$. Priors for the variability noise parameters $a_4$, $b$, $c$, and $u_0$ are informed by the model-agnostic variability quantification analysis described in <ref>.

Relative Bayesian evidence and median posterior model diameter vs. $m$-order, from the full-track geometric modeling applied to the data; error bars on the diameter measurements show $1\sigma$ credible intervals.  Because the absolute values of $\ln\mathcal{Z}$ can be substantially different for each data set (i.e., each combination of frequency band and calibration pipeline), and because only the relative values carry information about model specification preferences, we reference all $\ln\mathcal{Z}$ values to the maximum value attained at any $m$ for that data set.

§.§.§ implementation

We have implemented the combined plus noise parameterization as a model within the sampling-based parameter estimation framework developed for the EHT <cit.>.  Given a model specification and a data set, works within a Bayesian formalism to produce a set of samples from the posterior distribution of the model parameters.  uses a MCMC sampling scheme to explore the posterior space, employing a parallel tempering scheme <cit.> to ensure traversal over the entire prior volume and the Hamiltonian Monte Carlo sampling kernel from the
Stan package <cit.> to efficiently sample within each tempering level.  A detailed description of the sampling framework can be found in <cit.>.

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

\begin{equation}
\ln \mathcal{L} = - \frac{1}{2} \sum_i \left( \frac{| V_i - \hat{V}(\boldsymbol{u}_i ; \boldsymbol{p}) |^2}{\sigma_i^2(\boldsymbol{n})} + \ln\left[ 2 \pi \sigma_i^2(\boldsymbol{n}) \right] \right) .
\end{equation}

Here $V_i$  is a measured visibility, $\boldsymbol{u}_i$ is its baseline vector, $\hat{V}$ is the corresponding modeled visibility,
and the sum is taken over all data points $i$.  The noise $\sigma_i(\boldsymbol{n})$ in each visibility is specified as in <ref>, with the noise model parameters $\boldsymbol{n} = \{ f, \sigma_{\text{ref}}, a, b, c, u_0 \}$.  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 <cit.>.  It also applies the <cit.> diffractive scattering kernel directly to the model prior to computing visibilities.  A validation test of the plus noise model implementation is described along with other tests in a dedicated paper on the noise modeling approach  <cit.>.

We assess MCMC convergence through both visual inspection of the traces and a number of quantitative chain statistics, including the integrated autocorrelation time, split-$\hat{R}$, and parameter rank distributions <cit.>.  The number of tempering levels is selected to ensure efficient communication between the highest and lowest levels <cit.>, which typically requires about 20 levels.  We run the sampler for between $5 \times 10^4$ and $10^5$ steps per tempering level.

To compute the Bayesian evidence, uses thermodynamic integration <cit.>, which computes the log-evidence through

\begin{equation}\label{eq:TI}
\ln \mathcal{Z} = \int_{0}^1{\rm d}\beta \left<\ln \mathcal{L}\right>_{\beta},
\end{equation}

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

\begin{equation}
P^{(\beta)}(\bm{p},\bm{n}) = \left(\frac{P(\bm{p},\bm{n})}{\pi_{\rm ref}(\bm{p},\bm{n})}\right)^\beta \pi_{\rm ref}(\bm{p}, \bm{n}).
\end{equation}

Note that does not take $\pi_{\rm ref}$ to be the prior distribution. Instead, uses a uniform distribution whose support matches the support of the priors given in <ref>. To compute <ref>, we compute the average log-likelihood for each tempering level and then use trapezoidal integration to numerically compute the integral.

Priors for all model and noise model parameters are listed in <ref>.  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 <ref> .  All station gain phase priors are uniform on the unit circle.

§.§ Model selection

As with the snapshot geometric analyses (see <ref>), the model used for the full-track analyses is really a class of models that increases in complexity with $m$.  The results of a $m$-order survey covering $m=\{2,3,4,5,6\}$ are shown in <ref>.  In contrast to the snapshot geometric modeling (see <ref>), 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 \geq 4$ in <ref>).  Thus, to maintain consistency among model specifications and to facilitate comparison with the snapshot geometric modeling analyses, we proceed with the $m=4$ as our fiducial model for all full-track analyses in this paper.

A representative $m=4$ fit to the HOPS low-band data is shown in <ref>.  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 $\chi^2$ statistics for this fit are discussed in <ref>.

Results of full-track modeling using an $m=4$ fit to the HOPS low-band data on April 6 and 7, arranged analogously to the individual panels of <ref>.  As in <ref>, $\hat{V}$ denotes light-curve-normalized complex visibilities.

§ RESULTS

In this section we aggregate and present the results from the analyses described in Sections <ref> through <ref> as applied to the 2017 EHT data.

Similar to <ref>, but after the direct visibility estimates have been debiased to account for the short-time temporal correlations as described in <cit.>.  These estimates are directly comparable to the power spectra implied by GRMHD simulations <cit.>. A single example and the range associated with the library presented in  are shown by the red line and band, respectively.

Joint posteriors of the constrained parameters after fitting a broken power law to the debiased model-agnostic normalized variance estimates, specifically the normalization at $|\boldsymbol{u}|=4\,{\rm G}\lambda$, denoted by $a_4$, the break location $u_0$, and the long-baseline power-law index $b$.
Contours enclose 50%, 90%, and 99% probability.

Debiased variability parameters

Baseline length ()     Normalized variance     Fit

1c2.75     $(17.2 \pm 14.3)\times10^{-4}$     Y

1c3.25     $(4.9 \pm 2.3)\times10^{-4}$     Y

1c3.75     $(3.7 \pm 1.8)\times10^{-4}$     Y

1c4.25     $(2.6 \pm 1.6)\times10^{-4}$     Y

1c4.75     $(2.0 \pm 0.9)\times10^{-4}$     Y

1c5.25     $(3.2 \pm 1.7)\times10^{-4}$     Y

1c5.75     $(1.5 \pm 0.8)\times10^{-4}$     Y

1c6.25     $(0.5 \pm 0.4)\times10^{-4}$     N

1c6.75     $(0.5 \pm 0.4)\times10^{-4}$     N

1c7.25     $(0.4 \pm 0.4)\times10^{-4}$     N

1c7.75     $(0.3 \pm 0.2)\times10^{-4}$     N

1c8.25     $(0.3 \pm 0.2)\times10^{-4}$     N

1c8.75     $(0.3 \pm 0.3)\times10^{-4}$     N

Quantity     Symbol     Estimatea

Excess noise at $|\boldsymbol{u}|=4$     $a_4$       $1.90\pm0.2~\%$

Long-baseline power-law index     $b$         $2.4\pm0.8$

Break baseline length             $u_0$       $<1.3~{\rm G}\lambda$

Measurements of intrinsic variability from the EHT data.  The top portion of the table lists the visibility amplitude variance measurements from <ref> after applying the debiasing procedure described in <cit.>; these values correspond to the black points plotted in <ref>.  The bottom portion of the table lists the best-fit broken power-law model parameters for the variability measurements.
aQuoted are the median values and $1\sigma$ ranges.  Upper limits are $1\sigma$ limits.

§.§ Structural variability measurements

The model-agnostic variability quantification analysis carried out in <ref> demonstrates that the 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 <ref>, the measured variability can be broken down into three regions with qualitatively distinct behavior:

* On short baselines with lengths $|\boldsymbol{u}| \lesssim 2.5~{\rm G}\lambda$, corresponding to spatial scales $\gtrsim$100 , 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 <cit.>.  The source size constraint used to perform gain calibration of the LMT-SMT baseline  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~{\rm G}\lambda$.

* On intermediate baselines with lengths between $2.5~{\rm G}\lambda \lesssim |\boldsymbol{u}| \lesssim 6~{\rm G}\lambda$, corresponding to spatial scales between $\sim$30 and $\sim$100 , we measure significant variability that exhibits an approximately power-law decline with increasing baseline length.  The power-law index is between $\sim$2 and 3, and the magnitude of the variability ranges from a peak RMS of $\sim$5% of the total flux density ($\sim$120 mJy) near $2.5~{\rm G}\lambda$ down to $\sim$1% of the total flux density ($\sim$25 mJy) near $6~{\rm G}\lambda$.

* On long baselines with lengths $|\boldsymbol{u}| \gtrsim 6~{\rm G}\lambda$, corresponding to spatial scales $\lesssim$30 , 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  and to define priors on the corresponding parameters in the full-track modeling analyses (<ref>).

Mean images for each geometric modeling and image reconstruction pipelines applied to the data, showing both low band and high band separately for the geometric modeling and the combined bands for image reconstruction .  The geometric modeling and imaging pipelines have been applied to the combined April 6 and 7 data, while the , , and imaging pipelines have been applied to the April 7 descattered data; <ref> and <ref> show single-day results for all pipelines on April 6 and 7, respectively.  The upper group of images have been produced from the HOPS calibration pipeline data, while the bottom group of images correspond to the CASA calibration pipeline. All of the images share a common brightness color scale; the absolute brightness scale is arbitrary because each image has been normalized to have unit total flux density.

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 <ref> involves binning the visibility data in the -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 <cit.>, 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, <cit.> derive the debiasing function using the same 90 GRMHD simulations we use in this paper for $\theta_g$ calibration (see <ref> and <ref>), and the resulting debiasing factor is close to unity everywhere except between $\sim$6 and 7.5 ${\rm G}\lambda$ (see Fig. 5 in ).  Applying this debiasing function to the variability measurements from the data yields the results shown in <ref> and reported in <ref>.

Due to the near-unit debiasing factor, the variability measurements shown in <ref> are similar to those from <ref>. Quantitative constraints on the parameters of the noise model that are well constrained are presented in <ref> and <ref>. Where strong constraints on the the excess noise exist (i.e., on baseline lengths between $\sim$2 and 6 G$\lambda$), it continues to be well described by a single power law with index
and amplitude (evaluated at a baseline length of $4\,{\rm G}\lambda$) of
$a_4 = 1.9^{+0.2+0.5}_{-0.2-0.3}\%$,
where in each value the $1\sigma$ and $2\sigma$ ranges are indicated.  The measured variability magnitude is between $\sim$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_0<1.3\,{\rm G}\lambda$ at $1\sigma$ and $u_0<3.1\,{\rm G}\lambda$ at $2\sigma$.

The excess variability is broadly consistent with that due to structural fluctuations anticipated by the GRMHD simulations discussed in  and <cit.>.
The magnitude of the excess variability lies within the range of that predicted by GRMHD simulations for $|\boldsymbol{u}|>2~{\rm G}\lambda$, 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 .

Parameters describing image morphology
5cCASA         5cHOPS
5-9 11-15
Analysis class     Software     Day     Band     $\hat{d}$ ()     $W$ ()     $A$     $\eta$ (deg)     $f_c$         $\hat{d}$ ()     $W$ ()     $A$     $\eta$ (deg)     $f_c$
8*Imaging     +     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

+     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
+     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

+     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
+     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

+     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
+     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

+     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
2-15 4*Snapshot     2*     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
2*     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
2-15 2*Full-track     2*     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
Median values and 68% credible intervals for the morphological quantities of interest, measured from the EHT data. Because medians and quantiles are not well defined for angular variables, for the position angle $\eta$ we instead quote the circular mean and standard deviation.

§.§ Image morphology measurements

Both the snapshot (<ref>) and full-track (<ref>) geometric modeling analyses produce reconstructions of the 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 <ref> quantify the morphological properties of the top-set and posterior images reconstructed in .

<ref> compares the geometric models and image reconstructions determined for ; 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 , , and ) or posterior means (for ).  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.

We use the geometric modeling and IDFE analyses to quantify a number of morphological parameters of interest, which are shown in <ref> and listed in <ref>.

Morphological parameters measured from data using the geometric modeling (left columns), posterior imaging IDFE (middle column), and top-set imaging IDFE (right columns) analyses; each marker denotes a median value, and the error bars indicate 68% credible intervals. For geometric and posterior imaging we show the combined April 6 and 7 results, while for top-set imaging we show just descattered April 7. Each row of panels shows the results for a single morphological parameter, with the markers colored according to the method used to make the measurement (per the legend on the right).  Measurements made using data from both calibration pipelines are shown using circular and square markers for HOPS and CASA, respectively.  IDFE measurements made using and are indicated by filled and open markers, respectively; no metronization cuts have been applied to either top-set or posterior images.  The IDFE analyses have been applied to images reconstructed using both frequency bands simultaneously (LO+HI), while the geometric modeling analyses have been applied to each frequency band separately (LO and HI).

§.§.§ Ring size

The parameter of most interest for gravitational studies  is the ring size, which we quantify using its diameter.  For the modeling results from both snapshot and full-track analyses, we report a debiased diameter $\hat{d}$, given by

\begin{equation}\label{eq:moran_bias}
\hat{d} = d - \frac{1}{4\ln(2)}\frac{W^2}{d} ,
\end{equation}

where $d$ and $W$ are the ring diameter and thickness, respectively (see <ref>).  This debiasing corrects for the lowest-order impact of the Gaussian blurring kernel on the radial location of the peak intensity – which is shifted inward with respect to the radius of the pre-convolved ring  – and thus aids a more direct comparison of the geometric modeling diameter values with those obtained from IDFE.

The diameter measurements from the geometric modeling and IDFE analyses are compared in the top row of <ref> across frequency bands and calibration pipelines.  We find that the diameter is the most well constrained of the geometric parameters we attempt to quantify, with both a typical measurement uncertainty and a scatter between measurement types that is substantially smaller than the magnitude of the value itself. An average of the geometric modeling results across both frequency bands and calibration pipelines yields a debiased diameter of $51.9$ , with a corresponding symmetrized uncertainty of $2.0$ . The quoted error is the  68% (i.e., approximately $1\sigma$) probability and corresponds to the samples from each measurement weighted equally. The IDFE measurements are broadly consistent with the results from geometric modeling, yielding $\hat{d} = 51.8 \pm 2.6$ . The corresponding joint constraint from both geometric modeling and IDFE analyses yields $\hat{d} = 51.8 \pm 2.3$ .

We note that even after debiasing the diameter measurements from different analysis pathways remain interpretationally distinct quantities.  To ensure mutual consistency between different measurement methods, in <ref> we calibrate the diameter measurements to a common physical scale using the GRMHD synthetic data sets generated for this purpose (see <ref>).

§.§.§ Ring thickness

The thickness of the ring is of interest for its ability to constrain the location and size of the emitting region near the black hole <cit.>.  The ring thickness measurements from the geometric modeling and IDFE analyses are compared in the second row of <ref> across frequency bands and calibration pipelines.  We find that the geometric modeling methods recover similar ring thicknesses, with the snapshot analyses obtaining $\sim$16–22 and the full-track analyses obtaining $\sim$19–23 .  The IDFE analyses obtain consistently thicker rings, with $\sim$30 being a more typical value and a somewhat larger scatter (from $\sim$25 to 35 ) seen both across and within pipelines.  The increased ring thickness recovered from the IDFE analyses likely arises in part from image smoothing introduced by some reconstruction algorithms (e.g., the CLEAN algorithm used by ).

The thickness of the 87 ring was a parameter that the geometric modeling analyses carried out in  had difficulty constraining, and it showed substantial variation across days and between modeling methods; only a relatively weak upper limit for the fractional thickness of $W/\hat{d} \lesssim 0.5$ was obtained.  In contrast, we find for that the thickness parameter is relatively well constrained by geometric modeling approaches. The fractional thickness is $W/\hat{d} = 0.35 \pm 0.05$, where the uncertainties quoted are symmetrized $1\sigma$. IDFE analyses obtain systematically larger fractional thicknesses, finding $W/\hat{d} = 0.53 \pm 0.1$.

Unlike with the diameter measurements, we do not debias the ring thicknesses obtained from different analysis pathways or attempt to calibrate them to a common scale.  The ring thicknesses from geometric modeling and IDFE thus represent two interpretationally distinct quantities, and we do not produce an analysis-agnostic measurement of the ring thickness.

§.§.§ Position angle and asymmetry

The magnitude and orientation of any asymmetry in the azimuthal brightness distribution around the ring is of interest because it can be related to the spin and inclination of the black hole .  As described in <ref>, the position angle ($\eta$) and degree of azimuthal asymmetry ($A$) are both determined by the coefficient of the $m=1$ mode as specified in <ref> and <ref>, respectively.  These definitions match closely the corresponding IDFE quantities defined in <ref>.  The asymmetry and position angle measurements we obtain for are shown in the third and fourth rows of <ref>, respectively.

Unlike for 87, where the image morphology exhibits a clearly defined asymmetry axis whose magnitude and orientation can be consistently quantified using either geometric modeling or IDFE analyses , the image structure for is less amenable to such a description.  The asymmetry magnitude measurements show a large scatter across methods, spanning $\sim$0.15–0.3 for the geometric modeling methods and $\sim$0.04–0.20 for the IDFE methods.  The IDFE methods recover systematically smaller median levels of asymmetry than the geometric modeling methods, but the uncertainties are large; several of the measurement methods have statistical uncertainties that cover nearly the entire 0–0.5 range of the prior distribution for $A$.

The position angle measurements show similarly little consistency between different analysis methods, spanning essentially the full $(-180^{\circ},180^{\circ})$ range when compared across all data sets and measurement techniques.  The geometric modeling analyses find position angles that are loosely confined to a region between ${\sim}-100^{\circ}$ and $0^{\circ}$ across, but the IDFE analyses show a large ($>$100$^{\circ}$) scatter between methods and a similar magnitude of uncertainty for individual measurements.

§.§.§ Brightness depression

The depth of the brightness depression interior to the ring is a key signature of the presence of a black hole.  Additionally, it can be used to constrain the presence of an emitting or reflecting surface, as a potential alternative to a horizon (e.g., , ).  For the model, the fractional central flux $f_c$ is given by <ref>.  For the IDFE analyses, we retain the definitions from <ref> for $f_c$.

The $f_c$ measurements for are shown in the bottom row of <ref>.  We find a large spread in values across analysis methods, ranging from $\sim$0.1 to 0.25 for the geometric modeling analyses and from $\sim$0.0 to 0.5 for the IDFE analyses.  Compared against the constraints from geometric modeling of the 87 ring structure, which consistently found $f_c \lesssim 0.1$ , the results obtained here for allow for the possibility of substantially more emission interior to the ring.

§.§ Gravitational radius and mass

The ring size measurements presented in <ref> have been made using a variety of different analysis techniques, with different inherent assumptions and biases.  To bring these otherwise disparate measurement techniques to a common scale, we follow a strategy similar to that developed in 87 in  and calibrate the diameter measurements using simulations from the  GRMHD library.  As described in <ref>, our calibration suite consists of synthetic data sets constructed from 90 GRMHD simulations spanning a range of accretion flow and black hole parameters, and for which an absolute reference size scale (i.e., the angular size of the gravitational radius $\theta_g$) is known.  We apply the same data processing and ring diameter measurement strategies as used for the data to each of these synthetic data sets, and we use the resulting distribution of diameter measurements to derive the value and uncertainty in the scaling factor ($\alpha$) between $\theta_g$ and $d$ for every method (<ref>).  We note that a conceptually similar calibration is carried out in the companion , in which the GRMHD assumption is relaxed and a more diverse set of spacetimes and accretion flow models is used to calibrate ring size measurements.

§.§.§ Calibrated scaling factors

When applied to the calibration suite datasets, each diameter measurement technique produces a discrete distribution of $\alpha$ scaling factors.  We use a kernel density estimator (KDE) from the  package <cit.> to produce a nonparametric estimate of the continuous distribution corresponding to these discrete samples, and we use this KDE to construct the $\theta_g$ distribution for described below (<ref>).

<ref> lists the derived $\alpha$ value and its uncertainty – as computed from the KDE distribution – for each of the ring size measurement methods used in this paper.  The uncertainty in $\alpha$ contains two main components: a statistical uncertainty $\sigma_{\alpha}^{\text{(stat)}}$ associated with the fidelity of the ring measurement from each data set, and a theoretical uncertainty $\sigma_{\alpha}^{\text{(theory)}}$ associated with the intrinsic scatter in $\alpha$ as measured across different GRMHD calibration data sets.  The total uncertainty $\sigma_{\alpha}^{\text{(tot)}}$ is a combination of both the statistical and theoretical uncertainties, meaning that in practice we do not have access to $\sigma_{\alpha}^{\text{(theory)}}$ in isolation, and thus we do not report it in <ref>.  Nevertheless, the relative values of $\sigma_{\alpha}^{\text{(tot)}}$ and $\sigma_{\alpha}^{\text{(stat)}}$ indicate that the theoretical uncertainty is typically the dominant component.

The calibrated $\alpha$ values span a range of $\sim$10–12, depending on the specific measurement technique.  Statistical uncertainties in the $\alpha$ values determined using the geometric modeling techniques are a few percent, while for the IDFE techniques the statistical uncertainty is typically larger and reaches as high as $\sim$20% in the worst cases.  Folding in both the statistical and the theoretical components, the total uncertainties are more comparable across methods, though they still exhibit a large range spanning $\sim$15–35%.  Overall, the calibrated $\alpha$ values show similar magnitudes to what  derived from fits to 87, but the calibration uncertainty in the case of is substantially larger.  This increased uncertainty reflects the increased flexibility that has been built into the ring size measurement techniques to capture structural variability in the source, as well as the increased morphological diversity of the GRMHD calibration suite that is necessary to accommodate the a priori unknown inclination of .

$\alpha$ calibration parameters
Analysis Class     Software     Day     $\alpha$     $\sigma^{\rm (stat)}_\alpha$     $\sigma^{\rm (tot)}_\alpha$
April 6+7     12.0     (+0.2, -0.2)    (+1.6, -1.4)

April 6+7     11.0     (+0.8, -0.8)    (+2.2, -4.3)

April 6+7     11.7     (+0.1, -0.1)    (+1.3, -1.3)

April 7     10.5     (+0.9, -1.4)    (+2.0, -2.3)

April 7     10.6     (+1.0, -1.3)    (+1.7, -3.1)

April 7     11.0     (+1.4, -1.3)    (+2.1, -2.5)

April 7     11.0     (+1.2, -1.3)    (+1.7, -3.2)

April 7     10.3     (+2.4, -2.1)    (+2.8, -4.4)

April 7     10.4     (+1.4, -1.4)    (+1.8, -3.7)

April 6+7     10.3     (+0.5, -0.4)    (+1.5, -2.7)

April 6+7     10.6     (+0.4, -0.4)    (+1.2, -3.9)

Median values and 68% credible intervals for the calibrated $\alpha$ values, averaged over frequency bands and calibration pipelines.

$\theta_g$ distributions for all of the analyses in this paper. The solid color curves show the fitted generalized lambda distribution. The transparent curves show the kernel density estimate for the distribution. The plotted results have been obtained using the HOPS calibration pipeline. For reference, we overlay the $\theta_g$ constraints measured from stellar orbits with the Keck <cit.> and Very Large Telescope <cit.> facilities.

§.§.§ Sgr A* angular gravitational radius

We apply the calibrated $\alpha$ values to the measured ring diameters for each measurement technique, which produces a distribution of $\theta_g$ values that captures the uncertainties in both the ring measurements and in the GRMHD calibration. The resulting $\theta_g$ distribution exhibits sampling noise from KDE of the $\alpha$ distribution (due to the finite size of the calibration suite), and it can exhibit secondary low-probability modes at large values of $\theta_g$ (see <ref>). To provide a smooth unimodal estimate of the $\theta_g$ distribution, we fit a generalized lambda distribution (G$\lambda$D) to the KDE distribution for each measurement technique. The G$\lambda$D is a unimodal distribution representing a diverse family of probability density functions.  We use the G$\lambda$D parameterization from <cit.>, and we use the GLDEX package in R <cit.> to carry out the fitting.

The resulting $\theta_g$ distributions are shown in <ref> and listed in <ref>.  An average across all methods and datasets yields a joint constraint of $\theta_g = $ , where the uncertainties are quoted at the 68% (i.e., $1\sigma$) probability level and the systematic uncertainty is taken to be the standard deviation across all measurement methods.  This value is consistent with the considerably more precise constraints obtained from measurements of stellar orbits <cit.>, and the gravitational implications of this consistency are explored in .

$\theta_g$ values and uncertainties across analyses
2cHOPS        2cCASA

Analysis     Software     Day     LO     HI         LO     HI
April 6+7    $ 4.6_{-0.5}^{+0.7} $    $ 4.5_{-0.5}^{+0.7} $         $ 4.4_{-0.4}^{+0.7} $    $ 4.6_{-0.5}^{+0.7} $
April 6+7    $ 4.5_{-0.8}^{+1.7} $    …        $ 4.9_{-0.8}^{+1.8} $    …
April 6+7    $ 4.7_{-0.4}^{+0.8} $    $ 4.7_{-0.4}^{+0.8} $         $ 4.7_{-0.4}^{+0.8} $    $ 4.7_{-0.4}^{+0.8} $
April 7    2c$4.8_{-0.9}^{+1.2}$          2c$4.8_{-0.9}^{+1.2}$
April 7    2c$4.8_{-0.7}^{+1.4}$          2c$4.8_{-0.8}^{+1.5}$
April 7    2c$4.9_{-0.8}^{+1.2}$          2c$4.9_{-0.9}^{+1.3}$
April 7    2c$5.0_{-0.8}^{+1.6}$          2c$5.0_{-0.9}^{+1.7}$
April 7    2c$5.0_{-1.1}^{+2.8}$          2c$5.2_{-1.2}^{+2.8}$
April 7    2c$5.1_{-0.9}^{+1.9}$          2c$5.2_{-1.3}^{+2.6}$
April 6+7    2c$5.1_{-0.7}^{+1.1}$          2c…
April 6+7    2c$5.0_{-0.5}^{+1.3}$          2c…

The angular gravitational radius, $\theta_g$, measured across analysis schemes and data sets; all values are quoted as medians and 68% credible intervals in units of . Entries that straddle both low-band and high-band columns have been computed from fits to band-combined data sets.

§.§.§ Sgr A* mass

The angular size of the gravitational radius is proportional to the ratio of the mass and distance to , per <ref>.  Our constraints on $\theta_g$ can thus be mapped directly to constraints on the black hole mass $M$ by incorporating an independent distance measurement to .  <cit.> report a distance of $D = 8.15 \pm 0.15$ (at $1\sigma$ probability) to the Galactic center, measured using trigonometric VLBI parallaxes of a large number ($\sim$200) of masers.  Using this distance measurement along with our measurement of $\theta_g$ from <ref> yields a constraint on the mass of of $M = $ , where we again quote $1\sigma$ uncertainties and the systematic component is taken to be the weighted standard deviation across methods.  This measurement is once again consistent with the more precise constraints obtained from stellar orbits <cit.>, and the uncertainty in the mass remains dominated by our uncertainty in $\theta_g$.

§ SUMMARY AND CONCLUSIONS

In this paper we quantify the temporal variability and morphological properties of the horizon-scale emission from , using data taken by the EHT in 2017 April.  Our primary morphological quantity of interest is the diameter of the observed ring of emission, which we quantify using multiple independent analysis pathways.  We then use the ring diameter to place constraints on the angular size of the gravitational radius ($\theta_g$) and on the mass ($M$) of .  The analyses presented here have been carried out using data taken on April 6 and April 7, across two frequency bands and two data calibration pipelines.

Motivated by theoretical expectations that the dynamical timescales in should be much shorter than the duration of EHT observing tracks, we employ a new method developed in <cit.>
for quantifying the time variability observed in the visibilities in a manner that is agnostic to the specifics of the average underlying source structure.
We find that the visibility amplitudes
exhibit a light-curve-normalized variance that is in excess of that expected from thermal noise, station gains, or refractive scattering effects, and we attribute this excess variance to intrinsic structural changes in the source.  The detected variability is most statistically significant on baselines with lengths between 2.5 and 6 G$\lambda$, where it exhibits an approximately power-law decline with increasing baseline length, with a power-law index of $\sim$2–3.  The magnitude of this variability on baselines near 3 G$\lambda$ in length exceeds 0.1 Jy, which is roughly equal to the value of the correlated flux density on these same baselines.

Through an exploration of potential simple geometric source structures, we demonstrate that the EHT data statistically prefer ring-like morphologies over other morphologies with comparable complexity.
We develop and deploy two new methods for fitting the time-variable data using static geometric ring models with azimuthally modulated brightness structures.  In the first method, called “snapshot” geometric modeling, we first fit the models to short segments of data over which the source variability is subdominant to other sources of uncertainty.  The fits from individual segments are then combined via a hierarchical model averaging scheme to provide parameter constraints across the entire observation.  In the second method, called “full-track” geometric modeling, we fit a static geometric model to the entire data set alongside parameters that describe the statistical fluctuations that time variability induces in the data.  Our parameterization for the variability “noise” is motivated by the work of <cit.> and takes the form of a broken power law in baseline length that contributes to the data uncertainties.

We compare the results from snapshot and full-track geometric modeling, both with one another and with the results of IDFE from the images reconstructed in , to constrain the horizon-scale morphology of .  The ring diameter is well constrained and stable across both frequency bands and calibration pipelines, with geometric modeling and IDFE techniques jointly determining a value of (68% credible intervals).  We find that the magnitude and orientation of the ring asymmetry, as well as the depth of its central brightness depression, are poorly constrained and have values that can depend sensitively on the measurement method employed.  The thickness of the ring is well measured by individual analysis methods but takes on a value that depends on the specifics of each method; geometric modeling methods find a FWHM ring thickness of $35 \pm 5\%$ of the ring diameter, while image-domain feature extraction methods find a FWHM of $53 \pm 10\%$.

Using a suite of synthetic data sets generated from the  GRMHD simulation library, we calibrate the diameter measurements from both geometric modeling and image-domain feature extraction methods to a common physical scale.  The resulting constraint on the angular size of the gravitational radius, combined across all methods and data sets, is $\theta_g = $ .  This large uncertainty arises from both the model flexibility necessary to capture structural variability in the source and the broad morphological diversity of the GRMHD calibration suite that reflects the a priori unknown inclination of .  Combining our $\theta_g$ constraint with an independent distance measurement from <cit.>, we determine the mass of to be $M = $ .  Though the uncertainties are large compared to those derived using other techniques (e.g., stellar orbit modeling), our measurement represents the first time that the mass of has been constrained by observations of light bending near the horizon.

We thank an anonymous referee for insightful and constructive comments that helped improve the quality of this paper.

The Event Horizon Telescope Collaboration thanks the following
organizations and programs: the Academia Sinica; the Academy
of Finland (projects 274477, 284495, 312496, 315721); the Agencia Nacional de Investigación
y Desarrollo (ANID), Chile via NCN$19\_058$ (TITANs) and Fondecyt 1221421, the Alexander
von Humboldt Stiftung; an Alfred P. Sloan Research Fellowship;
Allegro, the European ALMA Regional Centre node in the Netherlands, the NL astronomy
research network NOVA and the astronomy institutes of the University of Amsterdam, Leiden University and Radboud University;
the ALMA North America Development Fund; the Black Hole Initiative, which is funded by grants from the John Templeton Foundation and the Gordon
and Betty Moore Foundation (although the opinions expressed in this work are those of the author(s)
and do not necessarily reflect the views of these Foundations);
Chandra DD7-18089X and TM6-17006X; the China Scholarship
Council; China Postdoctoral Science Foundation fellowship (2020M671266); Consejo Nacional de Ciencia y Tecnología (CONACYT,
Mexico, projects  U0004-246083, U0004-259839, F0003-272050, M0037-279006, F0003-281692,
104497, 275201, 263356);
the Consejería de Economía, Conocimiento,
Empresas y Universidad
of the Junta de Andalucía (grant P18-FR-1769), the Consejo Superior de Investigaciones
Científicas (grant 2019AEP112);
the Delaney Family via the Delaney Family John A.
Wheeler Chair at Perimeter Institute; Dirección General
de Asuntos del Personal Académico-Universidad
Nacional Autónoma de México (DGAPA-UNAM,
projects IN112417 and IN112820);
the Dutch Organization for Scientific Research (NWO) VICI award
(grant 639.043.513) and grant OCENW.KLEIN.113; the Dutch National Supercomputers, Cartesius and Snellius
(NWO Grant 2021.013);
the EACOA Fellowship awarded by the East Asia Core
Observatories Association, which consists of the Academia Sinica Institute of Astronomy and
Astrophysics, the National Astronomical Observatory of Japan, Center for Astronomical Mega-Science,
Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute;
the European Research Council (ERC) Synergy
Grant “BlackHoleCam: Imaging the Event Horizon
of Black Holes" (grant 610058);
the European Union Horizon 2020
research and innovation programme under grant agreements
RadioNet (No 730562) and
M2FINDERS (No 101018682);
the Generalitat
Valenciana postdoctoral grant APOSTD/2018/177 and
GenT Program (project CIDEGENT/2018/021); MICINN Research Project PID2019-108995GB-C22;
the European Research Council for advanced grant `JETSET: Launching, propagation and
emission of relativistic jets from binary mergers and across mass scales' (Grant No. 884631);
the Institute for Advanced Study; the Istituto Nazionale di Fisica
Nucleare (INFN) sezione di Napoli, iniziative specifiche
the International Max Planck Research
School for Astronomy and Astrophysics at the
Universities of Bonn and Cologne;
DFG research grant “Jet physics on horizon scales and beyond” (Grant No. FR 4069/2-1);
Joint Princeton/Flatiron and Joint Columbia/Flatiron Postdoctoral Fellowships,
research at the Flatiron Institute is supported by the Simons Foundation;
the Japan Ministry of Education, Culture, Sports, Science and Technology (MEXT; grant JPMXP1020200109); the Japanese Government (Monbukagakusho:
MEXT) Scholarship;
the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for JSPS
Research Fellowship (JP17J08829); the Joint Institute for Computational Fundamental Science, Japan; the Key Research
Program of Frontier Sciences, Chinese Academy of
Sciences (CAS, grants QYZDJ-SSW-SLH057, QYZDJSSW-SYS008, ZDBS-LY-SLH011);
the Leverhulme Trust Early Career Research
Fellowship; the Max-Planck-Gesellschaft (MPG);
the Max Planck Partner Group of the MPG and the
CAS; the MEXT/JSPS KAKENHI (grants 18KK0090, JP21H01137,
JP18H03721, JP18K13594, 18K03709, JP19K14761, 18H01245, 25120007); the Malaysian Fundamental Research Grant Scheme (FRGS) FRGS/1/2019/STG02/UM/02/6; the MIT International Science
and Technology Initiatives (MISTI) Funds;
the Ministry of Science and Technology (MOST) of Taiwan (103-2119-M-001-010-MY2, 105-2112-M-001-025-MY3, 105-2119-M-001-042, 106-2112-M-001-011, 106-2119-M-001-013, 106-2119-M-001-027, 106-2923-M-001-005, 107-2119-M-001-017, 107-2119-M-001-020, 107-2119-M-001-041, 107-2119-M-110-005, 107-2923-M-001-009, 108-2112-M-001-048, 108-2112-M-001-051, 108-2923-M-001-002, 109-2112-M-001-025, 109-2124-M-001-005, 109-2923-M-001-001, 110-2112-M-003-007-MY2, 110-2112-M-001-033, 110-2124-M-001-007, and 110-2923-M-001-001);
the Ministry of Education (MoE) of Taiwan Yushan Young Scholar Program;
the Physics Division, National Center for Theoretical Sciences of Taiwan;
the National Aeronautics and
Space Administration (NASA, Fermi Guest Investigator
grant 80NSSC20K1567, NASA Astrophysics Theory Program grant 80NSSC20K0527, NASA NuSTAR award
NASA Hubble Fellowship
grant HST-HF2-51431.001-A awarded
by the Space Telescope Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Inc., for NASA, under contract NAS5-26555;
the National Institute of Natural Sciences (NINS) of Japan; the National
Key Research and Development Program of China
(grant 2016YFA0400704, 2017YFA0402703, 2016YFA0400702); the National
Science Foundation (NSF, grants AST-0096454,
AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-1126433, AST-1140030,
DGE-1144085, AST-1207704, AST-1207730, AST-1207752, MRI-1228509, OPP-1248097, AST-1310896, AST-1440254,
AST-1555365, AST-1614868, AST-1615796, AST-1715061, AST-1716327,  AST-1716536, OISE-1743747, AST-1816420, AST-1935980, AST-2034306);
NSF Astronomy and Astrophysics Postdoctoral Fellowship (AST-1903847);
the Natural Science Foundation of China (grants 11650110427, 10625314, 11721303, 11725312, 11873028, 11933007, 11991052, 11991053, 12192220, 12192223);
the Natural Sciences and Engineering Research Council of
Canada (NSERC, including a Discovery Grant and
the NSERC Alexander Graham Bell Canada Graduate
Scholarships-Doctoral Program); the National Youth
Thousand Talents Program of China; the National Research
Foundation of Korea (the Global PhD Fellowship
Grant: grants NRF-2015H1A2A1033752, the Korea Research Fellowship Program:
NRF-2015H1D3A1066561, Brain Pool Program: 2019H1D3A1A01102564,
Basic Research Support Grant 2019R1F1A1059721, 2021R1A6A3A01086420, 2022R1C1C1005255);
Netherlands Research School for Astronomy (NOVA) Virtual Institute of Accretion (VIA) postdoctoral fellowships;
Onsala Space Observatory (OSO) national infrastructure, for the provisioning
of its facilities/observational support (OSO receives
funding through the Swedish Research Council under
grant 2017-00648);  the Perimeter Institute for Theoretical
Physics (research at Perimeter Institute is supported
by the Government of Canada through the Department
of Innovation, Science and Economic Development
and by the Province of Ontario through the
Ministry of Research, Innovation and Science); the Spanish Ministerio de Ciencia e Innovación (grants PGC2018-098915-B-C21, AYA2016-80889-P,
PID2019-108995GB-C21, PID2020-117404GB-C21);
the University of Pretoria for financial aid in the provision of the new
Cluster Server nodes and SuperMicro (USA) for a SEEDING GRANT approved towards these
nodes in 2020;
the Shanghai Pilot Program for Basic Research, Chinese Academy of Science,
Shanghai Branch (JCYJ-SHFY-2021-013);
the State Agency for Research of the Spanish MCIU through
the “Center of Excellence Severo Ochoa” award for
the Instituto de Astrofísica de Andalucía (SEV-2017-
0709); the Spinoza Prize SPI 78-409; the South African Research Chairs Initiative, through the
South African Radio Astronomy Observatory (SARAO, grant ID 77948),  which is a facility of the National
Research Foundation (NRF), an agency of the Department of Science and Innovation (DSI) of South Africa;
the Toray Science Foundation; Swedish Research Council (VR);
the US Department
of Energy (USDOE) through the Los Alamos National
Laboratory (operated by Triad National Security,
LLC, for the National Nuclear Security Administration
of the USDOE (Contract 89233218CNA000001); and the YCAA Prize Postdoctoral Fellowship.

We thank
the staff at the participating observatories, correlation
centers, and institutions for their enthusiastic support.
This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2016.1.01154.V. ALMA is a partnership
of the European Southern Observatory (ESO;
Europe, representing its member states), NSF, and
National Institutes of Natural Sciences of Japan, together
with National Research Council (Canada), Ministry
of Science and Technology (MOST; Taiwan),
Academia Sinica Institute of Astronomy and Astrophysics
(ASIAA; Taiwan), and Korea Astronomy and
Space Science Institute (KASI; Republic of Korea), in
cooperation with the Republic of Chile. The Joint
ALMA Observatory is operated by ESO, Associated
Universities, Inc. (AUI)/NRAO, and the National Astronomical
Observatory of Japan (NAOJ). The NRAO
is a facility of the NSF operated under cooperative agreement
by AUI.
This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National
Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract
No. DE-AC05-00OR22725. We also thank the Center for Computational Astrophysics, National Astronomical Observatory of Japan.
The computing cluster of Shanghai VLBI correlator supported by the Special Fund
for Astronomy from the Ministry of Finance in China is acknowledged.

APEX is a collaboration between the
Max-Planck-Institut für Radioastronomie (Germany),
ESO, and the Onsala Space Observatory (Sweden). The
SMA is a joint project between the SAO and ASIAA
and is funded by the Smithsonian Institution and the
Academia Sinica. The JCMT is operated by the East
Asian Observatory on behalf of the NAOJ, ASIAA, and
KASI, as well as the Ministry of Finance of China, Chinese
Academy of Sciences, and the National Key Research and Development
Program (No. 2017YFA0402700) of China
and Natural Science Foundation of China grant 11873028.
Additional funding support for the JCMT is provided by the Science
and Technologies Facility Council (UK) and participating
universities in the UK and Canada.
The LMT is a project operated by the Instituto Nacional
de Astrófisica, Óptica, y Electrónica (Mexico) and the
University of Massachusetts at Amherst (USA). The
IRAM 30-m telescope on Pico Veleta, Spain is operated
by IRAM and supported by CNRS (Centre National de
la Recherche Scientifique, France), MPG (Max-Planck-Gesellschaft, Germany)
and IGN (Instituto Geográfico
Nacional, Spain). The SMT is operated by the Arizona
Radio Observatory, a part of the Steward Observatory
of the University of Arizona, with financial support of
operations from the State of Arizona and financial support
for instrumentation development from the NSF.
Support for SPT participation in the EHT is provided by the National Science Foundation through award OPP-1852617
to the University of Chicago. Partial support is also
provided by the Kavli Institute of Cosmological Physics at the University of Chicago. The SPT hydrogen maser was
provided on loan from the GLT, courtesy of ASIAA.

This work used the
Extreme Science and Engineering Discovery Environment
(XSEDE), supported by NSF grant ACI-1548562,
and CyVerse, supported by NSF grants DBI-0735191,
DBI-1265383, and DBI-1743442. XSEDE Stampede2 resource
at TACC was allocated through TG-AST170024
and TG-AST080026N. XSEDE JetStream resource at
PTI and TACC was allocated through AST170028.
This research is part of the Frontera computing project at the Texas Advanced
Computing Center through the Frontera Large-Scale Community Partnerships allocation
AST20023. Frontera is made possible by National Science Foundation award OAC-1818253.
This research was carried out using resources provided by the Open Science Grid,
which is supported by the National Science Foundation and the U.S. Department of
Energy Office of Science.
Additional work used ABACUS2.0, which is part of the eScience center at Southern Denmark University.
Simulations were also performed on the SuperMUC cluster at the LRZ in Garching,
on the LOEWE cluster in CSC in Frankfurt, on the HazelHen cluster at the HLRS in Stuttgart,
and on the Pi2.0 and Siyuan Mark-I at Shanghai Jiao Tong University.
The computer resources of the Finnish IT Center for Science (CSC) and the Finnish Computing
Competence Infrastructure (FCCI) project are acknowledged. This
research was enabled in part by support provided
by Compute Ontario (http://computeontario.ca), Calcul
Quebec (http://www.calculquebec.ca) and Compute
Canada (http://www.computecanada.ca).

The EHTC has
received generous donations of FPGA chips from Xilinx
Inc., under the Xilinx University Program. The EHTC
has benefited from technology shared under open-source
license by the Collaboration for Astronomy Signal Processing
and Electronics Research (CASPER). The EHT
project is grateful to T4Science and Microsemi for their
assistance with Hydrogen Masers. This research has
made use of NASA's Astrophysics Data System. We
gratefully acknowledge the support provided by the extended
staff of the ALMA, both from the inception of
the ALMA Phasing Project through the observational
campaigns of 2017 and 2018. We would like to thank
A. Deller and W. Brisken for EHT-specific support with
the use of DiFX. We thank Martin Shepherd for the addition of extra features in the Difmap software
that were used for the CLEAN imaging results presented in this paper.
We acknowledge the significance that
Maunakea, where the SMA and JCMT EHT stations
are located, has for the indigenous Hawaiian people.

, , , , , , , Stan, ,

§ COMPARISON OF SGR A* VARIABILITY IN 2013 AND 2017

Estimate of the variance of the visibility amplitudes from the 2013 1.3 mm VLBI observations reported in <cit.>, indicated by the red open circles.  The associated statistical errors are shown by the red triangles.  For comparison, all of the elements of <ref> are reproduced as shown there, including the 2017 EHT variance estimates (filled and open black circles).

Additional prior epochs of millimeter-VLBI observations of provide a means to explore a wider range of baseline lengths and assess the consistency of the variability across many years.  <cit.> reported 1.3 mm VLBI observations of from 2013 March 21, 22, 23, 26, and 27 with an array that included JCMT, SMA, SMT, and the Combined Array for Research in Millimeter-wave Astronomy (CARMA, which has since been decommissioned).  By virtue of the small number of participating stations, these observations were much more limited in their -coverage than the 2017 EHT observations.  Nevertheless, the 2013 observations provide a second multiday data set for which normalized variance estimates may be produced and compared to those reported in this paper.

<cit.> reported only the visibility amplitudes, which we average in time on a per-scan basis.  This averaging presumes that there is no intrinsic phase evolution in the visibilities over the approximately 10-minute scan lengths; this approximation is well justified for the source sizes inferred from the 2017 EHT observations.  The visibility data are normalized by the intrasite SMA-JCMT, SMA-SMA, and CARMA-CARMA baselines.  Note that this normalization eliminates the need to perform a phase calibration like that applied to JCMT and LMT for the 2017 EHT data.

<ref> shows the reconstructed normalized variances from the 2013 observations in comparison to those from the 2017 EHT observations.
Baselines between 2.5 and 4 provide variability estimates that are broadly consistent between the two sets of observations.  This agreement suggests that the degree of structural variability exhibited by during the 2017 EHT campaign is not anomalous.

Similarly, because the station gains were well characterized for all stations during the 2013 observations, it is not necessary to make assumptions about the source size on short baselines.  Therefore, the 2013 observations provide estimates of the normalized variance on baselines shorter than 2 . The large statistical errors of these measurements preclude strong constraints on the variability below 1 , but there are nevertheless hints of a turnover in the variability power between 1 and 2 .

§ ORIGIN AND MITIGATION OF BIASES IN SINGLE-DAY ANALYSES

Visibility amplitudes from the HOPS low-band data set, averaged coherently over 120 s segments, on Apr 6 (red), April 7 (blue), and April 5 and 10 (gray) on the Chile-LMT baselines as functions of baseline length (left) and observing time (right).  Error bars indicate the error implied by the mean noise model and are intended to account for fluctuations due to variability in addition to statistical and known systematic error components.

Visibility amplitudes from the HOPS low-band GRMHD validation synthetic data set (data set 092 in <ref>), averaged coherently over 120 s segments, on April 10 (red), April 7 (blue), and April 5 and 10 (gray) on the Chile-LMT baselines as functions of baseline length (left) and observing time (right).  Error bars indicate the error implied by the mean noise model and are intended to account for fluctuations due to variability in addition to statistical and known systematic error components.  Vertical dotted green lines indicate the positions at which frames are shown in <ref>.

Average (left) and snapshot (right) images from the GRMHD simulation used to generate the validation synthetic data set shown in <ref>.  Three frames from the portion of the simulation used to generate the April 10 (top, outlined in red) and April 7 (bottom, outlined in blue) synthetic data sets are shown at the GMST times specified (corresponding to the times indicated by the vertical dotted green lines in <ref>).  All of the images share a common brightness color scale; the absolute brightness scale is arbitrary because each image has been normalized to have unit total flux density, and a modest amount of saturation has been permitted in the brightest regions to enhance the visibility of low-brightness features.

As discussed in <ref>, there is compelling evidence that exhibits structural variability on timescales ranging from minutes to a full observation night.  The degree of variability is estimated in <ref> and <ref>, and found to be the dominant contribution to the difference between the observed visibilities and those associated with a mean image for baselines with lengths between $\sim$3 and 6 .  Underlying the noise modeling mitigation method (see <ref>) is the assumption that the added “variability noise” modifies the data in a stochastic manner, i.e., coherent deviations do not persist throughout large patches of the -plane.  However, Earth-rotation aperture synthesis naturally results in correlated variability between visibilities that are nearby in the -plane, because visibilities that have small separations tend to also have small temporal separations.  Furthermore, the sparsity of the EHT array prevents most locations in the -plane from being sampled more than once in a single observation (though see <ref> for several exceptions), meaning that multiple observing days must be combined to access more than a single instantiation of the source variability.

The impact of structural variability can be seen most prominently on the Chile-LMT baselines, which exhibit coherent deviations on $\sim$1 hr timescales that are evident in the visibility amplitudes presented in <ref>.  These deviations are most pronounced near 1 h GMST ($\sim$12 UTC), where they are evident on both April 6 and April 7.  The residual gain uncertainties of $\sim$10%-20%  are insufficient to explain the dramatic drop near 4 around 1 h GMST on April 6.  At the same GMST on April 7, the visibility amplitudes fluctuate upward by a similar amount. On the remaining days the amplitudes at this GMST lie between the April 6 and April 7 values, indicating that the variations are associated with a process that is uncorrelated on interday timescales.

Similar coherent deviations are also observed in the synthetic data sets produced from GRMHD simulations for the purposes of calibrating ring size measurements to a common physical scale (see <ref> and <ref>). <ref> shows a similar set of visibility amplitudes on the Chile-LMT baselines for one of these synthetic data sets (data set 092 in <ref>).  While the date and time of the spurious feature differ from those seen in the real data – arising in the simulated data set near 3 h GMST ($\sim$14 UTC) – a dramatic deviation is present and persists for $\sim$0.5 hr.  In both the data and GRMHD simulations, these observed deviations are ${\gtrsim}3\sigma$ outliers (per the variance expected from the variability quantification scheme detailed in <ref>), they are nearly exclusively confined to the Chile-LMT baseline, and they are rare.  However, due to their coherent nature – i.e., many data points on a given day are similarly displaced – these fluctuations violate the assumption of statistical independence made by the noise modeling mitigation scheme on a single day.

The origin of the visibility amplitude deviations for the GRMHD simulations can be identified with coherent variable structures moving about the ring.  As shown in <ref>, instantaneous images from the GRMHD simulations can differ qualitatively from the average image, with the former sometimes dominated by small, bright patches of emission.  When aligned in the NW-SE direction and separated by the $\sim$50 ring diameter, these bright emission regions significantly impact the visibilities on the $|\boldsymbol{u}| \approx 4$ Chile-LMT baselines.  Days that do not exhibit large variations correspond to periods less impacted by such patchy emission structures.

Comparison of the results from imaging (top row), full-track geometric modeling (middle row), and snapshot modeling (bottom row) across different combinations of April 6 and 7 data sets.  The first column shows results from fitting to the April 6 data, the second column shows results from fitting to the April 6 data after flagging baselines between Chile and LMT, the third column shows results from fitting to the April 7 data, and the fourth column shows results from fitting to the combined April 6 and 7 data.  In all panels, we show an image corresponding to the posterior mean; for the imaging results, each sample image has been shifted during averaging so as to maximize the normalized cross-correlation computed with respect to a reference sample.  The full-track and snapshot modeling results show fits to the HOPS pipeline low-band data, while the imaging results show fits to the HOPS pipeline combined low- and high-band data.  All of the images share a common brightness color scale; the absolute brightness scale is arbitrary because each image has been normalized to have unit total flux density, and a modest amount of saturation has been permitted in the brightest regions to enhance the visibility of low-brightness features.

The coherent deviations seen on the Chile-LMT baselines can manifest as pathological behaviors in the single-day analyses carried out in this paper and in , particularly when applied to the sparser April 6 data set (see left column of <ref>).  The specific manner in which the visibility fluctuations impact the reconstructed source structure depends on the details of the analysis scheme and the freedom each underlying image model has to accommodate a subset of visibilities that are discrepant with the time-averaged structure.

Images reconstructed from the April 6 data exhibit parallel NE-SW streaks, typical of baseline artifacts associated with miscalibration of a single baseline, and a natural consequence of visibility amplitudes on the Chile-LMT baseline that are discrepant with the time-averaged image.  The top left panel of <ref> shows such artifacts in a image reconstruction; similar image artifacts are observed in the RML (and ) and, to a lesser extent, the CLEAN () reconstructions presented in .

By virtue of their specification, the source models are not capable of introducing streak-like features into the image structure.  When applied to , full-track fits to the April 6 data instead exhibit smaller ring sizes than those applied to the April 7 data (see the first and third panels in the second row of <ref>), with the April 6 fits preferring rings with a $\sim$40 diameter and the April 7 fits preferring rings with a $\sim$55 diameter.  This discrepancy in ring size may be associated with a shift in the location of the visibility minimum from $\sim$3 to $\sim$4 between April 7 and April 6, respectively.  However, we note that the behavior of the full-track model when applied to the April 6 data is a sensitive function of the $m$-order; changing the $m$-order can cause the model to prefer a $\sim$55 diameter.  Snapshot fits exhibit qualitatively similar behavior to the full-track fits, as shown in the third row of <ref>.  We again find that the April 6 data prefer a somewhat smaller ring diameter than the April 7 fits, though for snapshot fits the posterior distributions for the diameter parameter are consistent between the 2 days.

Despite their disparate forms, the various artifacts observed in April 6 reconstructions are effectively ameliorated by flagging the Chile-LMT baselines prior to carrying out the analyses (see the second column in <ref>), providing further evidence that the origins of the artifacts are confined to (or at least dominated by) the Chile-LMT baselines.  However, though flagging of these baselines is successful in preventing the specific analysis pathologies discussed above, such flagging is not otherwise motivated; there is no evidence for atypical data calibration issues on these baselines, and thus no reason to believe that the observed variability excess is anything other than intrinsic to the source.  Rather than flagging data, we proceed instead with the noise modeling scheme described in <ref>, which is itself intended to mitigate the effects of intrinsic variability on the source reconstructions.

As described above, the apparent failure of the noise modeling method to produce consistent results on some individual single-day analyses can be attributed to the coherent (rather than stochastic) nature of the variability sampled at any single location in the -plane, which violates the key assumption underlying the noise modeling approach that each data point represent an independent sample of the source variability.  This assumption would be more faithfully adhered to by a data set containing a larger number of independent variability realizations; within the context of the EHT observations, additional variability realizations are most naturally acquired by combining data sets across observing days.  Multiday analyses are also more consistent in spirit with the model-agnostic variability estimates in <ref> and <ref>.  Combining the independent realizations of the structural variability from multiple days improves the estimate of the mean visibilities, as evident in <ref>, in which the multiday mean of the visibility amplitudes is both smoother than and intermediate between those of April 6 and 7 individually.

The fourth column in <ref> shows reconstructions made using the combined April 6 and 7 data sets.  The improved behavior of the visibility means is reflected in better consistency across methods and a reduction in image artifacts.

§ SINGLE-DAY FITS

In this section we present the results from analyses carried out on the April 6 and April 7 data sets individually.  <ref> and <ref> show single-day images from each of the analysis pipelines (analogous to those shown in <ref>), and <ref> shows the corresponding measurements of morphological properties (analogous to those shown in <ref>).

Same as <ref>, but showing results obtained from using only the April 6 data.

Same as <ref>, but showing results obtained from using only the April 7 data.

Similar to <ref>, but showing results obtained from fitting to the April 6 and 7 data sets separately.  Each column from <ref> has now been split into a pair of columns showing April 6 results on the left and April 7 results on the right.

§ GRMHD SYNTHETIC DATA SET GENERATION

To calibrate measurements of the angular gravitational radius $\theta_g$ (see <ref>), we rely on the library of GRMHD simulations and associated GRRT synthetic movies produced and described in .
In this appendix we provide an overview of our model selection and data generation procedures, which are conceptually similar to the calibration analysis carried out in .

We select 90 simulations from the GRMHD library to be used for $\theta_g$ calibration and another 10 simulations to be used to validate this calibration.  The 90 calibration data sets uniformly grid a range of GRMHD parameters: every combination of the two MAD and SANE accretion states, five black hole spin values of $[-0.94, -0.5, 0, 0.5, 0.94]$, three inclinations of $[10^{\circ}, 50^{\circ}, 90^{\circ}]$, and three $R_{\text{high}}$ values of $[10, 40, 160]$ are represented.[$R_{\text{high}}$ is a parameter that sets the ratio of ion to electron temperatures in the simulated images; see , , and <cit.> for details.  For the GRMHD calibration suite, we have not included models with $R_{\text{high}}=1$ because these models tend to produce images with substantial extended emission that is large compared to the black hole shadow size; we note that the $R_{\text{high}}=1$ models are also rejected by the model selection constraints applied in .]  The 10 validation data sets are split evenly between MAD and SANE, but the black hole spins are randomly selected from $[-0.94, -0.5, 0, 0.5, 0.94]$, the inclinations are randomly selected from $[30^{\circ}, 70^{\circ}, 110^{\circ}, 130^{\circ}, 150^{\circ}]$, and the $R_{\text{high}}$ values are randomly selected from $[10, 40, 160]$.  The resulting model images contain a large variety of possible image morphologies, along with self-consistent dynamics as governed by the equations of GRMHD and GRRT.  The GRMHD parameters corresponding to each selected calibration and validation model are listed in <ref>, and some example average images are shown in <ref>.

GRMHD synthetic data set parameters
Data set index     Accretion state     Spin     Inclination     $R_{\text{high}}$     Position angle     Input $\theta_g$
000     MAD      -0.94                             10      10      155       4.477

001     MAD      -0.94                             10      40      -114      4.140

002     MAD      -0.94                             10      160     24        5.317

003     MAD      -0.94                             50      10      168       5.493

004     MAD      -0.94                             50      40      90        5.493

005     MAD      -0.94                             50      160     89        6.423

006     MAD      -0.94                             90      10      -177      3.843

007     MAD      -0.94                             90      40      -73       5.503

008     MAD      -0.94                             90      160     112       6.152

009     MAD      -0.5                  10      10      168       5.704

010     MAD      -0.5                  10      40      51        5.684

011     MAD      -0.5                  10      160     -12       4.502

012     MAD      -0.5                  50      10      -22       5.724

013     MAD      -0.5                  50      40      178       5.563

014     MAD      -0.5                  50      160     105       4.039

015     MAD      -0.5                  90      10      -171      5.935

016     MAD      -0.5                  90      40      146       3.596

017     MAD      -0.5                  90      160     -3        5.110

018     MAD      0.0       10      10      35        3.677

019     MAD      0.0       10      40      143       5.719

020     MAD      0.0       10      160     115       5.030

021     MAD      0.0       50      10      112       3.813

022     MAD      0.0       50      40      -101      4.301

023     MAD      0.0       50      160     -12       4.909

024     MAD      0.0       90      10      76        4.059

025     MAD      0.0       90      40      11        4.029

026     MAD      0.0       90      160     97        6.323

027     MAD      0.5       10      10      40        3.974

028     MAD      0.5       10      40      -4        4.658

029     MAD      0.5       10      160     126       6.272

030     MAD      0.5       50      10      -42       4.472

031     MAD      0.5       50      40      25        4.089

032     MAD      0.5       50      160     -135      5.598

033     MAD      0.5       90      10      108       5.251

034     MAD      0.5       90      40      171       5.432

035     MAD      0.5       90      160     140       5.015

036     MAD      0.94                  10      10      -54       5.679

037     MAD      0.94                  10      40      1         4.200

038     MAD      0.94                  10      160     -92       5.915

039     MAD      0.94                  50      10      -2        6.282

040     MAD      0.94                  50      40      161       5.131

041     MAD      0.94                  50      160     -89       5.699

042     MAD      0.94                  90      10      -48       5.025

043     MAD      0.94                  90      40      -99       4.587

044     MAD      0.94                  90      160     54        4.467

045     SANE     -0.94                             10      10      97        5.880

046     SANE     -0.94                             10      40      127       6.388

047     SANE     -0.94                             10      160     -142      6.267

048     SANE     -0.94                             50      10      -59       6.016

049     SANE     -0.94                             50      40      144       3.652

050     SANE     -0.94                             50      160     -110      4.411

051     SANE     -0.94                             90      10      46        3.783

052     SANE     -0.94                             90      40      -129      5.075

053     SANE     -0.94                             90      160     68        6.106

054     SANE     -0.5                  10      10      53        5.276

055     SANE     -0.5                  10      40      76        4.281

056     SANE     -0.5                  10      160     144       4.854

057     SANE     -0.5                  50      10      -173      6.418

058     SANE     -0.5                  50      40      55        5.070

059     SANE     -0.5                  50      160     66        5.000

060     SANE     -0.5                  90      10      154       5.080

061     SANE     -0.5                  90      40      -154      5.236

062     SANE     -0.5                  90      160     42        6.363

063     SANE     0.0       10      10      -31       4.321

064     SANE     0.0       10      40      -145      4.990

065     SANE     0.0       10      160     -13       3.612

066     SANE     0.0       50      10      70        5.684

067     SANE     0.0       50      40      83        4.773

068     SANE     0.0       50      160     -175      6.262

069     SANE     0.0       90      10      104       4.019

070     SANE     0.0       90      40      -68       5.367

071     SANE     0.0       90      160     -49       3.994

072     SANE     0.5       10      10      -116      6.142

073     SANE     0.5       10      40      -76       5.287

074     SANE     0.5       10      160     -17       4.763

075     SANE     0.5       50      10      138       5.613

076     SANE     0.5       50      40      -80       3.717

077     SANE     0.5       50      160     -109      6.333

078     SANE     0.5       90      10      128       6.403

079     SANE     0.5       90      40      -162      5.271

080     SANE     0.5       90      160     -7        3.586

081     SANE     0.94                  10      10      -46       4.768

082     SANE     0.94                  10      40      38        5.548

083     SANE     0.94                  10      160     122       5.875

084     SANE     0.94                  50      10      -36       5.920

085     SANE     0.94                  50      40      142       4.311

086     SANE     0.94                  50      160     177       6.096

087     SANE     0.94                  90      10      84        4.592

088     SANE     0.94                  90      40      138       6.142

089     SANE     0.94                  90      160     165       3.521

090     MAD      0.0       150     160     19        5.694

091     MAD      0.5       70      160     -32       5.271

092     MAD      -0.5                  30      160     46        4.497

093     MAD      0.94                  30      10      86        4.069

094     MAD      0.5       150     160     11        5.437

095     SANE     -0.94                             70      10      -51       4.778

096     SANE     0.5       110     40      176       5.835

097     SANE     -0.94                             130     160     -118      3.984

098     SANE     0.0       150     40      70        5.065

099     SANE     0.5       110     10      -71       5.020

Simulation parameters for each of the GRMHD-based synthetic data sets used for $\theta_g$ calibration (top; indices 000–089) and validation (bottom; indices 090-099).  The sign of the spin follows the convention of , where negative values indicate that the angular momentum of the accretion flow is antialigned with that of the black hole.  The inclination angle is given in degrees, with $90^{\circ}$ indicating an edge-on system and $0^{\circ}$ indicating a system whose spin vector is pointed toward us. The position angle is given in degrees east of north, and refers to the orientation of the black hole spin vector.  For each simulation, the input value of $\theta_g$ is given in .

Example average images from some of the GRMHD movies selected for synthetic data generation; each movie has been light-curve-normalized prior to averaging.
The data set indices are labeled in the upper left corner of each panel, and the corresponding GRMHD parameters are listed in <ref>.  All of the images share a common brightness color scale; the absolute brightness scale is arbitrary because each image has been normalized to have unit total flux density, and a modest amount of saturation has been permitted in the brightest regions to enhance the visibility of low-brightness features.  We note that these average images tend to be have much smoother structure than the individual frames of the movies that were averaged to produce them (see, e.g., <ref> for several example frames from one movie); the synthetic data sets themselves are produced from the movies and not from the average images.

After selecting the GRMHD models and prior to generating synthetic data, we first modify their orientations and angular sizes from their default simulation values.  We rotate each simulated movie by a position angle that is a uniformly chosen integer in the range $[-180^{\circ},180^{\circ}]$.  Each of the simulations from  was produced assuming a mass of $M = 4.14 \times 10^6$ M$_{\odot}$ and a distance of 8.127 kpc <cit.>, corresponding to an angular gravitational radius of $\theta_g = 5.03$ .  To avoid biasing our calibration in favor of any one value of $\theta_g$, we modify the overall spatial scale of each simulated movie by a random factor that is uniformly drawn from the range $[0.7,1.3]$.  The input position angles and gravitational radii for each movie are listed in <ref>.

Once the GRMHD models are selected and their movies rotated and scaled, we generate synthetic data sets in the same manner as for the synthetic data sets used in .  We use the software to first apply artificial scattering to the source structure per the scattering model from <cit.> and then sample the Fourier transform of each movie at a cadence and at $(u,v)$ locations identical to those of the EHT observations of .  The resulting visibilities are then corrupted with thermal noise and station-based gain and leakage effects at a level that is consistent with the data .  Eight synthetic data sets are generated for each GRMHD model, corresponding to the $(u,v)$-coverage on four observing nights – 2017 April 5, 6, 7, and 10 – and two frequency bands (see <ref>).

§ REPRESENTATIVE  VALUES FOR EACH ANALYSIS METHOD

In this section we provide some example representative $\chi^2$ values and associated quantities for each of the analysis methods used in this paper.  Specifically, we report values corresponding to the example fits shown in <ref>, <ref>, and <ref>.

For any fitted data quantity $q$ with modeled counterpart $\hat{q}$ and associated measurement uncertainty $\sigma$, we determine the $\chi^2$ as

\begin{equation}
\chi^2 = \sum_i^{N_{\text{data}}} \left( \frac{q_i - \hat{q}_i}{\sigma_i} \right)^2 ,
\end{equation}

where the sum is understood to be taken over all $N_{\text{data}}$ fitted data points.  We also define a reduced-$\chi^2$ value,

\begin{equation}
\chi_{\text{red}}^2 = \frac{\chi^2}{N_{\text{dof}}} ,
\end{equation}

where $N_{\text{dof}}$ is the number of degrees of freedom remaining in the data after accounting for the free parameters in the model.  However, we note that despite their familiarity, the interpretation of either of these $\chi^2$ statistics is complicated by several aspects of the analyses presented in this paper.

The first complication is that the number of degrees of freedom is generically unknown, rendering $\chi_{\text{red}}^2$ difficult to define in practice.  The use of informative priors and the presence of correlations among model parameters mean that $N_{\text{dof}}$ cannot be determined as simply the difference between $N_{\text{data}}$ and the number of free parameters in the model.  For example, in the RML imaging methods (and ), the number of effective image parameters is implicitly reduced – relative to the number that would be assumed by simply counting the total number of image pixels – by a large factor by the imposition of regularization terms (e.g., smoothness, sparseness) in the objective function .  Additionally, for all methods that fit to complex visibilities or visibility amplitudes, station gains are reconstructed as part of the fitting process; for those methods that simultaneously reconstruct high- and low-band data, the gains for the two bands are necessarily strongly correlated.  Furthermore, strong priors are imposed by network calibration , further reducing the number of effective model parameters and growing the effective number of degrees of freedom.

The second complication is that, with the exception of the snapshot modeling presented in <ref> (and ), all models make use of an added uncertainty budget to account for source variability.  For the imaging and full-track modeling analyses, the parameters describing this excess variability noise are simultaneously fit alongside those describing the image structure.  In both cases, the impact is to drive the $\chi_{\text{red}}^2$ value toward unity, rendering the resulting $\chi_{\text{red}}^2$ value not particularly meaningful as a metric of fit quality.

Nevertheless, in <ref> we present the $\chi^2$ values and relevant properties of the data sets and models used for representative examples from each analysis pathway.  In lieu of a well-defined number of degrees of freedom, we consider two limits.  An optimistic estimate is given by the procedure adopted in , in which $N_{\rm dof}\approx N_{\rm data}$, appropriate when the total number of effective model parameters is much less than the number of data points.  In this limit, $\chi_{\text{red}}^2$ ranges from 0.35 to 0.79 for the fits listed in <ref>.  A more pessimistic accounting is given by $N_{\rm dof}\approx N_{\rm data}-N_{\rm gains}-N_{\rm params}$.  This quantity is negative for some analyses, a consequence of the strong correlations that limit the effective number of model parameters in practice (e.g., for RML imaging and for certain individual snapshot models).  Among those analysis methods that do not exhibit this pathology in the $N_{\rm dof}$, the pessimistic $\chi_{\text{red}}^2$ estimates range from 0.9 to 1.4.

Representative $\chi^2$ values and ancillary model information
imaging                   Apr 7       HI+LO     4174.5     12082       6358      6400

imaging                   Apr 7       HI+LO     2192.9      6806       3562      5625

imaging                  Apr 7       HI+LO     4194.5     12082       6358      2262

imaging                  Apr 6+7     HI+LO     1444.5      2990       1323        58

full-track modeling      Apr 6+7     LO        1237.0      1562        666        19

snapshot modeling       Apr 6+7     LO        1987.2      4977       1345      3328

snapshot modeling           Apr 6+7     LO        1752.7      3176     …     1664

The $\chi^2$ values, number of data points fit ($N_{\text{data}}$), number of station gain parameters ($N_{\text{gains}}$), and number of parameters describing the image structure ($N_{\text{params}}$) for representative example fits from each of the analysis methods used in this paper.  For the and analyses, $N_{\text{params}}$ is equal to the number of pixels in the image ; for the analyses, $N_{\text{params}}$ is determined by the number of CLEAN components, with each component contributing three parameters.  For the imaging and full-track modeling analyses, $N_{\text{params}}$ is equal to the number of model parameters used to specify the image structure (see  and <ref>, respectively).  For the and snapshot modeling analyses, each of the $N_s$ snapshots contributes 13 and 8 parameters, respectively, to the total $N_{\text{params}}$ count (see <ref>).  Note that because the analysis fits only to closure quantities, it does not contain any station gain parameters.
†Note that the listed numbers of gains and model parameters are solely those necessary to forward-compute the model data values, i.e., these values characterize a property of the model specification.  Importantly,
due to model nonlinearity and strong correlations between parameters,
these numbers are not generally suitable for determining $N_{\text{dof}}$ (see <ref>).

§ SNAPSHOT MODELING LIKELIHOOD FUNCTIONS

In this appendix we provide specific expressions for the likelihood functions used during the snapshot geometric modeling analyses described in <ref>.  We assume the high signal-to-noise ratio limit for all data products, which is not strictly satisfied for the relatively short integration times (120 s) employed in the snapshot modeling analyses, but which has the benefit of reducing all likelihood functions to Gaussians.

§.§ Visibility amplitude likelihood

For each snapshot and baseline, the visibility amplitudes are distributed according to a Rice distribution, which in the high signal-to-noise ratio limit becomes Gaussian <cit.>.  We can thus write the visibility amplitude likelihood function as

\begin{equation}\label{eq:ssmod:va}
\mathcal{L}_{|V|, sb} = \frac{1}{\sqrt{2\pi\sigma_{sb}^2}}
\exp\left[-\frac{(|\hat{V}_{sb}| - |g_i||g_j||\hat{\mathcal{V}}_{sb}|)^2}{2\sigma_{sb}^2}\right] ,
\end{equation}

where $b$ is a baseline index that runs over all station pairs $\{i,j\}$ in snapshot $s$.  Here $|\hat{\mathcal{V}}_{sb}|$ is the model visibility amplitude and $|g_i|$ and $|g_j|$ are the individual station gain amplitudes (see <ref>). We use this Gaussian approximation to the Rice distribution for all of the snapshot geometric modeling analyses. For a snapshot $s$, the joint visibility amplitude likelihood across all baselines is then given by

\begin{equation}
\mathcal{L}_{|V|,s} = \prod_{b}\mathcal{L}_{|V|, sb} ,
\end{equation}

where the product is taken over all baselines $b$.

§.§ Closure phase likelihood

The visibility phases in EHT data sets are heavily corrupted by atmospheric fluctuations , so all of our snapshot geometric modeling analyses work instead with closure phases $\psi$ (see <ref>).  In the high signal-to-noise ratio limit, the variance in the closure phase on the triangle containing stations $i$, $j$, and $k$ is given by

\begin{equation}\label{eq:cphaseerr}
\sigma^2_{\psi, ijk} =
\sigma^2_{\ln |V|, ij} +
\sigma^2_{\ln |V|, jk} +
\sigma^2_{\ln |V|, ki} ,
\end{equation}

\begin{equation}
\sigma_{\ln |V|,ij} = \frac{\sigma_{ij}}{|V_{ij}|},
\end{equation}

is the uncertainty in the log visibility amplitude in the same limit and $\sigma_{ij}$ is the uncertainty in $V_{ij}$. Hereafter we replace the triangle indices $ijk$ with a single multi-index $t$ for clarity. We approximate the closure phase likelihood for a single triangle $t$ and snapshot $s$ by a von Mises distribution,

\begin{equation}
\mathcal{L}_{\psi, st} = \frac{1}{2\pi I_0(\sigma_{\psi, t}^{-2})}\exp\left[\frac{\cos(\psi_{t} - \hat{\psi}_{t})}{\sigma^2_{\psi,t}}\right] ,
\end{equation}

where $\hat{\psi}$ denotes a measured closure phase and $I_0(x)$ is a modified Bessel function of the first kind of order 0.  In the high signal-to-noise ratio limit, the von Mises becomes a Gaussian distribution with mean $\hat{\psi}_{t}$ and standard deviation $\sigma_{\psi, t}$. Note that up to a normalization $\ln \mathcal{L}_{\psi}$ is equal to the closure phase $\chi^2$ defined in .

Because a full set of closure phases is highly redundant <cit.>, we use <cit.> to construct a minimal nonredundant set within every snapshot.  The minimal set of closure phases is constructed by selecting all triangles that contain the most sensitive station in the array (typically ALMA), which ensures that the resulting closure phases are minimally covariant.  We ignore the remaining covariances and approximate the joint closure phase likelihood for a snapshot $s$ as

\begin{equation}
\mathcal{L}_{\psi, s} = \prod_{t} \mathcal{L}_{\psi, st},
\end{equation}

where the product is taken over all triangles $t$ in the minimal set.

§.§ Log closure amplitude likelihood

Some of our snapshot modeling fits use log closure amplitudes in place of visibility amplitudes.  In the high signal-to-noise ratio limit, the variance in the log closure amplitude on the quadrangle containing stations $i$, $j$, $k$, and $\ell$ is given by

\begin{equation}\label{eq:logclosureerr}
\sigma^2_{\ln A, ijk\ell} =
\sigma^2_{\ln |V|, ij} +
\sigma^2_{\ln |V|, k\ell} +
\sigma^2_{\ln |V|, ik} +
\sigma^2_{\ln |V|, j\ell} .
\end{equation}

Hereafter we replace the quadrangle indices $ijk\ell$ with a single multi-index $q$ for clarity.
For a single log closure amplitude on quadrangle $q$, we thus have the Gaussian likelihood

\begin{equation}
\mathcal{L}_{\ln A, sq} =\frac{1}{\sqrt{2\pi\sigma_{\ln A, sq}^2}}
\exp\left[-\frac{(\ln A_{sq} - \ln\hat{A}_{sq})^2}{2\sigma_{\ln A,sq}^2}\right],
\end{equation}

in the same limit.  We use to construct a minimal nonredundant set of log closure amplitudes within every snapshot, and for the joint likelihood over all of the linearly independent quadrangles within that snapshot we again ignore covariance and treat the data products as being statistically independent.  The joint likelihood is then given by

\begin{equation}
\mathcal{L}_{A, s} = \prod_{q} \mathcal{L}_{A, sq},
\end{equation}

where the product is taken over all quadrangles $q$ in the minimal set.

§ SNAPSHOT MODELING PRIOR DISTRIBUTION CONSISTENCY

In specifying the hypermodel for combining geometric modeling results from individual snapshots (see <ref>), we have selected a set of snapshot priors $\pi(\bm{\theta}_s)$, an average model prior $\pi(\bar{\bm{\theta}})$, and a hypermodel $\pi(\bm{\theta}_s | \bar{\bm{\theta}})$.  A self-consistent set of hypermodel and snapshot priors should satisfy the relation

\begin{equation}
\pi(\bm{\theta}_s) = \int \pi(\bm{\theta}_s | \bm{\mu}, \bm{\sigma}) \pi(\bm{\mu})\pi(\bm{\sigma}) {\rm d}\bm{\mu}{\rm d}\bm{\sigma}. \label{eqn:PriorConsistency}
\end{equation}

However, because we have carried out the hierarchical modeling from <ref> in two stages – i.e., first fitting the model to individual snapshots, then combining the results from all snapshots – the priors on the model parameters during each stage are selected for local convenience, and the relation in <ref> does not strictly hold.

As a result, the final posteriors for each parameter behave as though the individual snapshot fits were carried out using effective priors that differ from the priors we actually imposed (see <ref>).  <ref> illustrates how the original priors differ from the effective priors induced by the snapshot-combining procedure.  The net result is a modest ($\lesssim$20% maximum probability density deviation) bias for parameter values to fall toward the center of the original (flat) prior range and away from the edges of that range.  Because the widths of the individual snapshot posteriors are typically much narrower than the full prior range (see, e.g., <ref>), we expect the impact of this centralizing bias on our reported values to be negligible.