MODELING SEVEN YEARS OF EVENT HORIZON TELESCOPE OBSERVATIONS WITH RADIATIVELY INEFFICIENT ACCRETION FLOW MODELS

The amplitudes of the complex visibilities are weakly dependent on the potentially large phase delays that occur during propagation through the atmosphere. While they are impacted by atmospheric absorption and gain uncertainties within the telescopes, this can be addressed via careful calibration. Encoded in the amplitudes is primarily information on characteristic scales within the emission region. More importantly, these can be constructed from observations by pairs of telescopes, permitting probes of horizon-scale structures without multiple long baselines. Hence, the initial mm-VLBI observations exclusively reported visibility magnitudes. Here, we employ the observations reported in Doeleman et al. (2008) and Fish et al. (2011), in which can be found the full details of the observations, calibration, and data processing.

Doeleman et al. (2008) describes mm-VLBI observations conducted over two nights in 2007 April, using the JCMT, SMT, and a single CARMA antenna. Nineteen visibility amplitudes were measured on the CARMA-SMT and JCMT-SMT baselines, with only an upper limit obtained on the JCMT-CARMA baseline, which we ignore in favor of more recent detections due to the weak constraint it applies. The signal-to-noise ratios of the incoherently averaged visibility magnitudes of 4 and 8 were typical on the short and long baselines, respectively. The full CARMA array was operated concurrently as a stand-alone array, and measured an effective zero-baseline flux of $2.4\pm 0.25\;\mathrm{Jy}$. This is similar to the visibility magnitudes measured on the shorter CARMA-SMT baseline. It is, however, anomalously low in comparison to the more typical 3 Jy flux at 1.3 mm, implying that Sgr A* was in a quiescent state.

Fish et al. (2011) present subsequent observations over three nights in 2009 April that employed the JCMT, SMT, and multiple CARMA antennas operated independently. A total of 54 visibility magnitudes were obtained on JCMT-SMT, CARMA-SMT, and both JCMT-CARMA baselines on two of the three days. The signal-to-noise ratios of the incoherently averaged visibility amplitudes were considerably higher, reaching 17 and 5 on the short and long baselines, respectively. In addition, the two independent CARMA antennae formed a very short baseline, accessing angular scales on the order of 10'' and finding substantially more correlated flux than implied by detections on the CARMA-SMT baselines. As in Broderick et al. (2011a), we will assume that this arises from a separate large-scale component not present in the 2007 observations. This is supported by the similarities in the structures inferred using only the 2007 and 2009 visibility magnitudes despite significant variations in their overall normalizations (Fish et al. 2011). Hence, we do not consider the inter-CARMA visibilities further here.

Due to the challenges in calibrating visibility magnitudes and the novel nature of the closure phases (see the following section), we leave the consideration of the visibility magnitudes obtained in subsequent observing runs for future work.

The measured phase of a complex visibility, ${{\rm{\Phi }}}_{E}$, contains information on the structure of the observed source but is corrupted by variations introduced by the propagation of radiation through the atmosphere. At longer wavelengths, these delays can be removed by rapid nodding between a target source and a nearby calibrator, but at 1.3 mm the timescale over which the atmospheric phase contribution changes by a radian can be too rapid—as short as a few seconds depending on weather conditions—to permit phase transfer from a nearby calibrator.

Fortunately, it is possible to construct VLBI observables that depend on the visibility phases of a source while being robust against atmospheric corruption. The simplest such quantity, closure phase, is constructed by taking the directed sum of the visibility phases along a closed triangle of baselines. The phase introduced by the unknown atmospheric delay at a station on one baseline is exactly canceled by the phase on the other baseline including that station. In fact, closure phases are insensitive to almost all station-based phase effects, whether atmospheric or instrumental in origin.

The closure-phase data set used in this work consists of 181 measurements on the California–Hawaii–Arizona triangle taken from 2009 March through 2013 March. The data are described in greater detail in Fish et al. (2016). These are shown explicitly in Figure 1 as a function of index, which indicates (roughly) the time ordering. As with the amplitude data, errors on the closure-phase data are smaller in later epochs due to increases in sensitivity of the array, especially the use of phased arrays at SMA and CARMA. An analysis of the closure-phase data determined that the error distributions are very well approximated by Gaussians with a standard deviation in radians equal to the inverse of the signal-to-noise ratio.

Zoom In Zoom Out Reset image size

We begin by assessing if our RIAF model class provides a statistically adequate description of the EHT observations. We do this via the χ² statistic, both for the entire data set and for contributions from sub-regions within it, listed in Tables 2 and 3.

Table 2.  Fit Results by Epoch

Epoch	Na	kb	min χ2c	fit χ2d	V00e	${\phi }_{E}^{M}$ f	${\bar{\phi }}_{E}$ g
1	19	4	5.77	7.50	2.44	⋯	⋯
2	12	4	10.8	11.0	2.16	⋯	⋯
3	19	4	18.8	28.3	2.12	⋯	⋯
4	20	4	12.4	13.4	2.95	⋯	⋯
5	11	4	9.68	22.3	⋯	1.90	1.21
6	3	4	1.63	4.19	⋯	9.10	2.51
7	10	4	2.90	3.96	⋯	4.02	2.42
8	7	4	0.883	8.41	⋯	−5.92	−2.37
9	2	4	1.74 × 10–12	0.185	⋯	5.05	0.831
10	5	4	1.45	3.73	⋯	8.34	4.53
11	17	4	10.0	11.4	⋯	13.6	10.9
12	25	4	14.5	23.1	⋯	8.30	7.24
13	28	4	27.5	44.7	⋯	0.73	0.65
14	10	4	3.77	9.05	⋯	2.91	1.76
15	15	4	8.62	14.2	⋯	3.54	2.55
16	32	4	37.8	46.9	⋯	5.61	5.23
17	16	4	4.83	6.00	⋯	7.91	5.99

Notes.

^aNumber of data points, including detections only. ^bNumber of fit parameters. ^cEpoch specific minimum $\chi 2$. ^dEpoch specific $\chi 2$ of the most probable model. ^eVisibility magnitude normalization. ^fMost likely closure-phase offset. ^gMarginalized closure-phase offset.

Download table as:  ASCIITypeset image

Table 3.  Model Comparison

Model Class	Na	kb	min χ2	${\rm{\Delta }}\mathrm{AIC}$ c	Odds Ratioc
$\theta \geqslant 90^\circ$ with ${\phi }_{E}=0$	251	7	336.4	57.0	9.7×10−9
$\theta \leqslant 90^\circ$ with ${\phi }_{E}=0$	251	7	338.8	59.5	6.2×10−10
$\theta \geqslant 90^\circ$ with ${\phi }_{E}\ne 0$	251	20	251.1	0.94	0.93
$\theta \leqslant 90^\circ$ with ${\phi }_{E}\ne 0$	251	20	250.2	0	1

Notes.

^aNumber of data points, including detections only. ^bNumber of fit parameters. ^cThe $\theta \leqslant 90^\circ$, ${\phi }_{E}\ne 0$ is used as the reference.

Download table as:  ASCIITypeset image

For the entire data set, the minimum χ² is related to the maximum likelihood in the standard way and has the advantage of being easily interpreted. We find a value of 250.2, consistent with the expected value of 231 at the 2σ level, i.e., at a p-value of $18.4\%$.¹² Thus, as found in earlier analyses, RIAF models continue to provide an excellent fit to the horizon-scale structure of Sgr A*. Critically, we note that the phenomenological models considered in Broderick et al. (2011a), as well as annuli and symmetric double point sources, are unable to produce the small but non-zero closure-phase evolutions observed, excluding these with overwhelming confidence.¹³

A similar analysis can be performed for each epoch individually, assessing whether the RIAF model provides a satisfactory fit for each. This is complicated by the large variation in the number of data points within each epoch; some epochs have too few measured values to permit an independent fit, i.e., ${N}_{e}\leqslant 4$. Excluding these from consideration, the remaining epochs all admit good fits, with the smallest p-value being 10%, for epoch 16. We also show the epoch-specific χ² values after including an isotropic prior on the orientation and priors on ${\phi }_{E}$. While marginally worse than the best-fit models at each epoch, it is also a good fit for all epochs individually (minimum p-value of $2.2\%$, assuming the number of degrees of freedom is simply N_e). Therefore, we conclude that there is no evidence that the RIAF model is insufficient to describe the quiescent emission of Sgr A* for any epoch individually or for all of the epochs combined.

Included in Table 3 are fit results for when the inter-epoch closure-phase variations are not modeled. As anticipated, the quality of the resulting fits is much worse—the minimum χ² is 336.4, with an associated p-value of $4.8\times {10}^{-6}$ or excluded at a level of more than $4.4\sigma$. However, the smaller set of fit parameters, 7 instead of 20, motivates a careful comparison of the fit quality. We do this in two ways, both of which are shown in Table 3. In all cases, we take the $\theta \lt 90^\circ$, with closure-phase fluctuations modeled as our default case, to which the others are compared.

The first employs the Akaike information criterion ($\mathrm{AIC}$), defined by

$$\begin{eqnarray}&&\mathrm{AIC}={\chi }_{\mathrm{min}}^{2}+2k+\displaystyle \frac{2k(k+1)}{N-k-1}.\end{eqnarray} \tag{ 1 }$$

where k is the number of fit parameters, N is the total number of data points, and ${\chi }_{\mathrm{min}}^{2}$ is the minimum χ² achieved (Akaike 1974; Takeuchi 2000; Burnham & Anderson 2002; Liddle 2007). The $\mathrm{AIC}$ is a penalized χ² statistic and presents an approximate measure of the difference between the true and modeled data distributions (Takeuchi 2000). Only differences in the $\mathrm{AIC}$ are statistically relevant, measured in the Jeffrey's scale, for which ${\rm{\Delta }}\mathrm{AIC}\gt 10$ are "decisive" evidence against the model with the higher $\mathrm{AIC}$. The model which ignores the closure-phase fluctuations is therefore decisively excluded by this criterion, exhibiting ${\rm{\Delta }}\mathrm{AIC}\gt 57$ for both potential inclinations relative to the models in which the inter-epoch closure-phase variations are modeled.

The second measure we employ is the Bayesian Odds Ratio, defined by the ratio of the model likelihoods, marginalized over all of the fit parameters (see, e.g., Gregory 2005). That is, given two models, M_I and M_II, with parameters ${{\boldsymbol{p}}}_{I,{II}}$, parameter priors ${{\rm{\Pi }}}_{I,{II}}({{\boldsymbol{p}}}_{I,{II}})$, and associated likelihoods ${L}_{I,{II}}({{\boldsymbol{p}}}_{I,{II}})$, the Odds Ratio is given by

$$\begin{eqnarray}&&\mathrm{Odds}\;\mathrm{Ratio}\equiv \displaystyle \frac{\int {d}^{m}{p}_{I}\;{L}_{I}({{\boldsymbol{p}}}_{I})\;{{\rm{\Pi }}}_{I}({{\boldsymbol{p}}}_{I})}{\int {d}^{n}{p}_{{II}}\;{L}_{{II}}({{\boldsymbol{p}}}_{{II}})\;{{\rm{\Pi }}}_{{II}}({{\boldsymbol{p}}}_{{II}})}.\end{eqnarray} \tag{ 2 }$$

Note that up to an overall prior on the model as a whole, this is simply the ratio of the posterior probabilities of the two models given the data. As such, in the absence of a clear prior preference, the Odds Ratio provides a straightforward assessment of their relative success. Since this necessarily incorporates the assumed priors, and therefore penalizes the effort to model the inter-epoch closure-phase variability in two ways: first, by the inclusion of a larger parameter space volume and, second, by a prior penalty on large ${\phi }_{E}$. Nevertheless, the Odds Ratio conclusively disfavors the analysis in which the closure-phase variability is ignored.

The magnitudes of ${\phi }_{E}$ required are listed in Table 2 both for the most likely model and after imposing the prior on ${\phi }_{E}$ described in Section 4. The marginalized shifts, ${\bar{\phi }}_{E}$, are roughly normally distributed with a mean of 334 and standard deviation of 341, implying a small but significant net shift in the closure-phase measurements, possibly implying systematic deviations from the underlying model. Such a shift could in principle be caused by turbulent structures within the accretion flow that are not modeled here; typically, orbiting compact emission features produce a net bias in closure phases on the California–Hawaii–Arizona triangle that depends on the orientation of the disk (see, e.g., left column of Figure 4 in Doeleman et al. 2009b), although future work is required to assess the magnitude and direction of these biases. The largest closure-phase shifts occur for epochs 11, 12, and 17, all of which are clearly discrepant with the others in Figure 2. In the latter two cases, however, the shift is less than $2{\sigma }_{\phi }=7\_\_AMP\_\_fdg;7$. The inferred mean closure phase during epoch 11 is sensitive to the fringe search method employed, varying between −03 and 35. As a result, it is consistent with being less discrepent than epochs 12 and 17. For this reason, we do not consider the apparently large shifts required for epoch 11 to be significant at this time, and thus all shifts are within $2{\sigma }_{\phi }$ of zero.

Having established that our RIAF model class provides an adequate description of the EHT data, we now turn to assessing the consistency among epochs. With the exception of flare events, there is no reason for the underlying image morphology of an aligned RIAF to evolve across epochs. Misaligned accretion flows can precess on a variety of timescales, leading to secular changes in the reconstructed source parameters. However, within the context of the RIAF model class considered here, we expect that the black hole parameters that characterize the images will be fixed. Thus, while necessary, the existence of good fits for each epoch is insufficient—we also require that the reconstructed parameters be consistent among all epochs.

Posterior probability distributions are obtained by first marginalizing over the epoch-specific flux renormalizations (assuming a flat prior, which is well justified within the small range of fluxes permitted) and closure-phase shifts (assuming a Gaussian prior). This may be done analytically, and results in the same flux normalization estimate (Broderick et al. 2014, see also Appendix B). We further assume a flat prior on the spin magnitude and an isotropic prior on the spin orientation. The resulting three-dimensional posterior probability distribution is a function only of a, θ, and ξ.

Example posterior probability distributions are shown in Figure 3 for the combination of the visibility magnitude data (reproducing Broderick et al. 2011a) and three epochs of indicative closure-phase data with samples from each year: epochs 5, 8, and 16. All of the epochs have solutions consistent with the $\xi =156^\circ$ solution to the visibility magnitude fits. Thus, as anticipated by the reasonable χ² found for the combined fit in the previous section, a consistent set of spin parameters exists among all of the epochs. This is shown for each fit parameter as a function of epoch in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This correspondence within epochs 1–4 has already been discussed in Broderick et al. (2011a), and to a significant degree is a natural consequence of the consistency of the visibility amplitudes with each other. A similar conclusion largely holds for epochs 5–17, where the consistency between the evolution of the closure-phase evolutions implies that any satisfactory model will be similar on all epochs. This is explicitly exhibited by the comparison of the best-fit model with the closure-phase data in Figure 1, which shows the same characteristically rising evolution. Nevertheless, the data consistency is not absolute, as seen in Figure 4, where small variations do appear in the epoch-specific parameter reconstructions.

More remarkable is the consistency among qualitatively distinct classes of intereferometric data. Earlier analyses relied on the visibility magnitudes alone (resulting in the 180° position angle degeneracy evident in Figures 3 and 4). While it may be reasonably expected that high-quality fits to subsequent, consistent visibility magnitude measurements will produce similar parameter estimates, this is certainly not true for visibility phases, and therefore closure phases. A clear counterexample is the statistically successful (though comparatively disfavored) Gaussian spot model for the emission region which predicts identically vanishing closure phases (e.g., Broderick et al. 2011a). Thus, the recent set of closure-phase measurements presents an a priori test of the original visibility-magnitude selected RIAF models. As a consequence, the consistency among epochs, and therefore the consistency among classes of interferometric observables, is an impressive success of the RIAF picture.

The impact of black hole spin on the millimeter-wavelength images of RIAF is nontrivial and arises from a number of sources, including modifications to the null geodesics traversed by the millimeter photons and the underlying dynamics of the emitting material. Importantly, for RIAFs, it is not possible to assign the impact of spin to the modification of the ISCO, as is commonly the case for thin disks (this is evident in efforts to constrain perturbations to Kerr with similar analyses, e.g., Broderick et al. 2014). For the set of spectrally fit, aligned RIAF models considered here, the magnitude of the spin is determined primarily by the size of the emission region: higher spins result in both more rapidly orbiting material in the inner regions of the accretion flow and a greater degree of alignment between the emitting gas and the millimeter-photon directions, leading to more severe Doppler boosting and beaming, all of which at fixed total flux produces a more compact image. Because we consider an aligned class of RIAF models, we necessarily mix the impact of the black hole spin and the accretion flow inclination. The degree to which one dominates the constraint over the other depends on the orientation and magnitude of the black hole spin; at vanishing spins the orientation reconstruction is essentially dominated by that of the accretion flow.

Having verified that high-quality fits exist and are consistent among epochs, we now turn to the estimation of the black hole spin magnitude and orientation. The combined epoch fits are shown in Figure 5 with the two-dimensional probability distribution marginalized over position angle shown in Figure 6 and the one-dimensional marginalized probability distributions for each parameter shown in Figure 7. As found in Broderick et al. (2011a), the allowed parameter space is restricted to a narrow region, being primarily degenerate in a and θ. Unlike Broderick et al. (2011a), the allowed region is now becoming constrained in all directions, with the result that the three-dimensional and one-dimensional parameter estimates are becoming similar. That is, the large-scale correlations that have complicated the characterization of the allowed spin parameters are disappearing.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.4.1. Spin Orientation

For the quiescent image models considered here, an approximate reflection symmetry across the equatorial plane of the black hole, broken only by absorption by the black hole itself and the nontrivial optical depth of the accretion flow, prevents distinguishing between inclinations with the same $| \mathrm{cos}\theta |$, i.e., between θ and $180^\circ -\theta$. The images produced by the two cases are related by an image reflection, and thus we independently fit model libraries for both $\theta \leqslant 90^\circ$ and $\theta \geqslant 90^\circ$. At present, the two model classes are indistinguishable, both having nearly equally good fits and having an Odds Ratio of unity (see Table 3). Figures 3–6 show the case when $\theta \leqslant 90^\circ$, illustrative of both. Figure 7 shows results marginalized over both potential inclinations. Note that dynamical observations within the context of shearing turbulent flows can explicitly break this degeneracy (Johnson et al. 2015).

Nevertheless, despite the degeneracy, the quantitative estimates for the inclination are significantly improved, with constraints that are nearly a factor of two better than those in Broderick et al. (2011a): $\theta ={60^\circ }_{-{8}^{^\circ }-{13}^{^\circ }}^{+{5}^{^\circ }+{10}^{^\circ }}$.¹⁴ When marginalized over all of the other parameters, this results in an inclination estimate of $\theta ={60^\circ }_{-{4}^{^\circ }-{8}^{^\circ }}^{+{2}^{^\circ }+{5}^{^\circ }}$, which is again a more than a factor of two improvement over the marginalized estimates of Broderick et al. (2011a). In both cases, these are consistent at the 1σ level with the previous constraint, as shown explicitly in the center panel of Figure 7.

Unlike the visibility magnitudes, closure phases are able to conclusively select a single position angle, marking both a qualitative and quantitative improvement in its determination. The value selected is $\xi ={156^\circ }_{-{17}^{^\circ }-{27}^{^\circ }}^{+{10}^{^\circ }+{14}^{^\circ }}$ east of north. After marginalizing over the spin magnitude and inclination, the position angle estimate is $\xi ={156^\circ }_{-{3}^{^\circ }-{18}^{^\circ }}^{+{10}^{^\circ }+{16}^{^\circ }}$, which is consistent with the values reported in Broderick et al. (2011a) between the 1σ and 2σ level.

The small shift in the reconstructed position angle may indicate a mild tension between the orientations inferred from the visibility magnitudes alone and those from the closure phases alone. The origin of this tension remains unclear, though it may derive from a variety of sources. Much of the angular information in the visibility magnitude analyses is found in the early and late visibility measurements along the long Hawaii-CARMA and Hawaii-SMT baselines, where Sgr A*'s elevation at all stations is at its lowest, and therefore the flux calibration most challenging. Errors in the reconstructed station fluxes would mimic structure in the north–south direction, modifying the reconstructed image orientation. These calibration difficulties do not impact the inferred structure from closure-phase measurements, which instead derive from their non-zero value at all times. Alternatively, the short time variability associated with compact features in the accretion flow can induce significant dynamical deviations in the closure phases that are not accounted for in our attempts to address the inter-epoch closure-phase variability (Doeleman et al. 2009b). As with calibration uncertainties, this will modify the position angle most strongly as a result of the small north–south baselines at present. Therefore, unsurprisingly, the imminent improvements in the north–south coverage of the EHT will be critical in addressing the position angle of Sgr A* and the attendant implications for accretion modeling.

The reconstructed orientation of Sgr A* is in remarkable agreement with a variety of features within the Galactic center (Psaltis et al. 2015), shown in Figure 8. Most notably, one of the inclination solutions is aligned at the 1σ level with the inferred orbit of the infrared-luminous gas clouds G1 and G2 (Pfuhl et al. 2015). Equally suggestive is the near-alignment with the clockwise disk of young stars (CWD), believed to be responsible for feeding the accretion flow onto Sgr A* (Lu et al. 2009; Genzel et al. 2010). These are natural for two reasons, relating either to the structure of the accretion flow or the spin of the black hole.

Zoom In Zoom Out Reset image size

First, the strong viscous coupling within RIAFs, arising from the large-scale magnetic fields that result from their large scale heights, prevents the Bardeen–Petterson effect from efficiently aligning the disk with the black hole spin (Bardeen & Petterson 1975; Fragile et al. 2007). Thus, for small black hole spins the orientation of an RIAF is determined by the original angular momentum of the accreting gas.¹⁵

Second, if these stars formed in situ, then the accretion of the associated gas disk is sufficient to reorient the black hole spin, naturally accounting for the inferred alignment. To overcome the local tidal forces requires a disk mass of ${M}_{\mathrm{disk}}\approx {10}^{4}$–${10}^{5}\;{M}_{\odot }$ extending to a distance of ${R}_{\mathrm{disk}}\approx 0.4\;\mathrm{pc}$, of which only roughly ${10}^{3}\;{M}_{\odot }$ would have resulted in stars, with the remainder accreting onto Sgr A* over the past ${10}^{6}\;\mathrm{year}$ (Levin 2007; Bartko et al. 2009). The orbital angular momentum in this gas is roughly

$$\begin{eqnarray}&&{J}_{\mathrm{gas}}\approx \sqrt{{{GMR}}_{\mathrm{disk}}}{M}_{\mathrm{disk}}=\displaystyle \frac{{J}_{\bullet }}{a}\sqrt{\displaystyle \frac{{R}_{\mathrm{disk}}{c}^{2}}{{GM}}}\displaystyle \frac{{M}_{\mathrm{disk}}}{M}\gtrsim 3\displaystyle \frac{{J}_{\bullet }}{a}\end{eqnarray} \tag{ 3 }$$

where J_• is the angular momentum of the black hole. Thus, generally, there is enough angular momentum in the associated accretion disk to align the black hole if the two can be efficiently coupled. This depends, in turn, on the details of the disk accretion.

The accretion rate at the black hole is highly uncertain, depending critically on the impact of disk winds. For wind models similar to those in the models we have considered here, $\dot{M}(r)\propto {r}^{0.45}$, and thus the accretion rate at the horizon is only $0.15\%$ of that at the star-forming region. Nevertheless, even in the absence of wind losses, from Levin (2007) the accretion rate near $r={R}_{\mathrm{disk}}$ is ${\dot{M}}_{\mathrm{disk}}\approx 0.02{\dot{M}}_{\mathrm{Edd}}$ assuming typical values. Thus, in the presence of even marginal mass loss the accretion will proceed within the RIAF regime. As a consequence, magnetic torques are expected to strongly couple the accretion flow over the large distances required to align the black hole spin.

It is noteworthy that the implied total energy output, $3\times {10}^{53}\mbox{--}2\times {10}^{56}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$ (assuming radiative efficiency of 1%) is comparable to the energy budget required to inflate the Fermi bubbles, kiloparsec-scale gamma-ray features above and below the Galactic plane believed to have occurred contemporaneously with the formation of the CWD. This is consistent with the emerging picture of Sgr A*'s recent accretion history, characterized by the recent end of a more active phase roughly ${10}^{6}\;\mathrm{year}$ ago (see, e.g., Ponti et al. 2014 and references therein).

The reconstructed orientation of Sgr A* is potentially in alignment with the X-ray jet feature reported by Li et al. (2013). In addition, it is consistent with the orientation of the large-scale accretion flow as inferred from dynamical drag required to place G1 and G2 on similar orbits (Madigan et al. 2016; McCourt & Madigan 2016). However, either would necessarily preclude alignment with the CWD. Notably, no solution for the orientation of Sgr A* aligns with the more distant counter-CWD of stars (CCWD), circumnuclear disk at 3 pc (CND), Galactic rotation axis (and thus kiloparsec-scale Fermi bubbles, Su et al. 2010; Ackermann et al. 2014), or the nearby S-stars (Gillessen et al. 2009b; Genzel et al. 2010).

5.4.2. Spin Magnitude

The magnitude of Sgr A*'s spin remains tightly correlated with the inclination, as clearly seen in Figures 4 and 5. This is marginally weaker than was found in Broderick et al. (2011a), with a and θ related by

$$\begin{eqnarray}&&\theta \approx 63^\circ -31^\circ a\pm 1\_\_AMP\_\_fdg;5\pm 3\_\_AMP\_\_fdg;5.\end{eqnarray} \tag{ 4 }$$

The spin magnitude is not significantly correlated with the position angle. Consequently, the improvements of the constraints on the spin magnitude mirror those for the spin inclination. Unlike in Broderick et al. (2011a), the most probable spin magnitude, $a={0.10}_{-0.10-0.10}^{+0.30+0.56}$, is formally non-zero, although it remains consistent with a = 0, in which case the orientation corresponds to that of the accretion flow angular momentum. The most probable model is shown in Figure 9.

Zoom In Zoom Out Reset image size

Marginalizing over orientation produces strict upper limits of $a={0}^{+0.18+0.45}$, ignoring the possibility of anti-alignment (corresponding here to $a\lt 0$), and is shown in Figure 7. As before, the resulting spin magnitude estimates are consistent at the 1σ level with those in Broderick et al. (2011a). This is complicated, however, by the choice of spin magnitude prior. Regardless, the conclusion in Broderick et al. (2011a) that the spin magnitude must be small continues to be supported by the subsequent EHT observations.

For clarity, we have adopted a flat prior. While ideally a prior may be obtained from simulations of supermassive black hole growth (e.g., Volonteri et al. 2005, 2013; Barausse 2012), considerable astrophysical uncertainties persist preventing this in practice. A completely agnostic spin magnitude prior would favor high spins and is ∝a², corresponding to the volume of the angular momentum phase space of a spinning top. However, the value of such an arbitrary prior is questionable, even in light of its emphasis on high spins $a\lt 0.51$ at 2σ and $a\lt 0.73$ at 3σ. Thus, regardless of the prior adopted, within the context of the RIAF models considered here, very high spins are excluded at high significance.

Typical refractive angles are expected to be on the order of $20\;\mu \mathrm{as}$, which is comparable to the scattering kernel width at $1.3\;\mathrm{mm}$. This is similar to the scale of small structures in the image. Nevertheless, it is instructive to consider the weak refraction limit, i.e., where the image distortions due to scattering occur over scales that are small compared to the typical structures in the image. We treat the strong-scattering limit numerically in Section A.2.

In the small-angle scattering limit, a refractive screen modifies the image via a distortion field ${\boldsymbol{\xi }}({\boldsymbol{x}},t)$ that varies as the screen passes across the source. That is, the observed intensity map is given in terms of the fixed intrinsic image $\bar{I}({\boldsymbol{x}})$ by

$$\begin{eqnarray}&&I({\boldsymbol{x}},t)=\bar{I}[{\boldsymbol{x}}+{\boldsymbol{\xi }}({\boldsymbol{x}},t)]\approx \bar{I}({\boldsymbol{x}})\left[1+{\boldsymbol{\xi }}({\boldsymbol{x}},t)\cdot \displaystyle \frac{{\boldsymbol{\nabla }}\bar{I}({\boldsymbol{x}})}{\bar{I}({\boldsymbol{x}})}\right],\end{eqnarray} \tag{ 13 }$$

i.e., ${\rm{\Delta }}=(\xi \cdot {\boldsymbol{\nabla }}\bar{I})/\bar{I}$. The resulting expression for $T({\boldsymbol{u}},t)$ is then given by

$$\begin{eqnarray}&&T({\boldsymbol{u}},t)=1+\displaystyle \frac{1}{\bar{V}({\boldsymbol{u}})}\;\left[\tilde{{\boldsymbol{\xi }}}*\;\left(\displaystyle \frac{2\pi i{\boldsymbol{u}}}{\lambda }\bar{V}\right)\right]({\boldsymbol{u}},t).\end{eqnarray} \tag{ 14 }$$

While this may be immediately inserted into Equation (12) to obtain the variations in the closure phase, it is instructive to consider the limit in which there is a clear separation between the scales of $\bar{I}({\boldsymbol{x}})$ and the fluctuations.

For our purposes here, because the perturbations are expected to be small, the second term in Equation (14) may still be sub-dominant when this is realized in practice. That is, we assume that $u\gg {u}_{0}$, where u₀ is the baseline length above which $\bar{V}({\boldsymbol{u}})$ begins to decrease substantially. In the interest of concreteness, we will assume that

$$\begin{eqnarray}\bar{V}({\boldsymbol{u}})\approx {V}_{0}\left\{\begin{array}{ll}1 & u\lt {u}_{0}\\ {(u/{u}_{0})}^{-1} & u\gt {u}_{0},\end{array}\right.\end{eqnarray} \tag{ 15 }$$

where the asymptotic power is typical of that arising from the power-law brightness decline in images for RIAFs. As a result,

$$\begin{eqnarray}&&\left[\tilde{{\boldsymbol{\xi }}}*\;\left(\displaystyle \frac{2\pi i{\boldsymbol{u}}}{\lambda }\bar{V}\right)\right]({\boldsymbol{u}},t)\approx \displaystyle \frac{2\pi {{iu}}_{0}}{\lambda u}{V}_{0}\;{\boldsymbol{u}}\cdot \tilde{{\boldsymbol{\xi }}}({\boldsymbol{u}},t),\end{eqnarray} \tag{ 16 }$$

and $T({\boldsymbol{u}},t)\approx 1+[2\pi i{\boldsymbol{u}}\cdot {\boldsymbol{\xi }}({\boldsymbol{u}},t)/\lambda ]$. Inserting this into Equation (12) yields

$$\begin{eqnarray}&&\delta {{\rm{\Phi }}}_{123}\approx {\boldsymbol{\delta }}{\boldsymbol{u}}\cdot {{\boldsymbol{\nabla }}}_{{\boldsymbol{u}}}\displaystyle \frac{2\pi {\boldsymbol{u}}\cdot {\boldsymbol{\xi }}({\boldsymbol{u}},t)}{\lambda }.\end{eqnarray} \tag{ 17 }$$

Between the characteristic scales of the quiescent image and the scattering screen, ${\boldsymbol{\xi }}$ is roughly constant, giving

$$\begin{eqnarray}&&\delta {{\rm{\Phi }}}_{123}\approx \left(\displaystyle \frac{2\pi u}{\lambda }\xi \right)\left(\displaystyle \frac{\delta u}{u}\right).\end{eqnarray} \tag{ 18 }$$

For $\delta u/u\approx 0.2$, which is appropriate for the CARMA-SMT versus Hawaii baselines, reproducing the $\delta {{\rm{\Phi }}}_{123}\approx 0.07$ requires typical refractive distortions of $\xi \approx 0.05\lambda /u\approx 3\;\mu \mathrm{as}$ on scales of $60\;\mu \mathrm{as}$.

The intra-day timescale is consistent with recent models of a "nearby" scattering screen (Bower et al. 2014) motivated by observations of the recently discovered magnetar (Kennea et al. 2013) with a velocity of $\approx 30\;\mathrm{km}\;{{\rm{s}}}^{-1}$, similar to the velocity dispersion of stars in the disk. It is also consistent, however, with a "distant" scattering screen (e.g., Lazio & Cordes 1998) assuming velocities of $\approx 100\;\mathrm{km}\;{{\rm{s}}}^{-1}$, comparable to those expected in the bulge. Thus, both the magnitude and timescale of the closure-phase variations are consistent with an origin in the scattering screen.

In the strong-scattering regime, i.e., where the scattering induces angular rearrangements on scales comparable to the structures in the image, the preceding perturbative analysis is insufficient. Here, we briefly describe an attempt to simulate this scattering and numerically infer the typical variations in the observed closure phases.

As in the weak case, strong refractive scattering in the interstellar medium produces stochastic fluctuations in images, and hence in the visibility magnitudes and closure phases (Goodman & Narayan 1989; Narayan & Goodman 1989; Johnson & Gwinn 2015). Although the strength of these fluctuations depends on the properties of the scattering, the most significant uncertainties arise from the unknown source image: an extended source quenches the fluctuations in a manner that depends on its size and structure. As in the perturbative case, both nearby and distant scattering screen models are consistent with the observed intra-day variability.

Figure 10 shows an example of the effects of refractive scattering on our best-fit model. The scattering kernel was taken from Bower et al. (2006), and we assumed an inner scale of $1.5\times {10}^{4}\;\mathrm{km}$ for the turbulence. Although this inner scale is somewhat larger than expected for the interstellar medium, it simplifies the scattering simulation and has little effect ($\sim 10\%$) on the resulting refractive fluctuations. Following Johnson & Gwinn (2015), we scattered the model image by first generating a random phase screen with ${2}^{13}\times {2}^{13}$ random phases and then shifting the unscattered image by the scaled gradient of the phase screen.

Zoom In Zoom Out Reset image size

Figure 11 shows the estimated root mean square (rms) fluctuations of the closure phase on the SMT-CARMA-SMA triangle, estimated by sampling the visibilities on an ensemble of scattered images. These vary with time due to the time-variable orientation of the participating telescopes. The times at which observations were made extend from 0.5 GST to 3.8 GST, with a median near 1.8 GST, suggesting a typical rms near 35. This is very similar to the 386 standard deviation observed, implying that the bulk of the closure-phase variation may be due to interstellar scattering.

Zoom In Zoom Out Reset image size

The impact of turbulence on the image is complicated by anisotropy and inhomogeneity as well as the opacity within the emission region. Here, we will ignore these complications in the interest of obtaining a qualitative result, assuming an optically thick, homogeneous emission region, appropriate for the brightest component of the image of Sgr A*. In this case, the intensity is proportional to ${n}_{e}{B}^{-1/2}$, and hence fluctuations in the local electron density and magnetic field strength produce corresponding fluctuations in the intensity map. For magnetohydrodynamic turbulence, these are related via $\delta B/B=\alpha \delta {n}_{e}/{n}_{e}$, where α depends on the particular turbulence model under consideration (e.g., for fixed total pressure $\alpha =-1/2$, which we adopt), and thus $\delta I/I=(5/4)\delta n/n$. Particle conservation relates the density variation to the fluid displacement field, i.e., assuming uniform density, $\delta {n}_{e}/{n}_{e}=-{\boldsymbol{\nabla }}\cdot {\boldsymbol{\xi }}$. While the turbulence is manifestly three-dimensional, the observational consequences are dominated by the transverse mass redistribution, and hence we will concern ourselves only with the two-dimensional dimensional structure of the density (and magnetic field) on the sky. Therefore, the resulting intensity map is given by

$$\begin{eqnarray}&&I({\boldsymbol{x}},t)=\bar{I}({\boldsymbol{x}})\left[1+\displaystyle \frac{5}{4}{\boldsymbol{\nabla }}\cdot {\boldsymbol{\xi }}({\boldsymbol{x}},t)\right].\end{eqnarray} \tag{ 19 }$$

This is similar to the weak refraction case, with the exception that now it is the divergence of ${\boldsymbol{\xi }}$ that enters.

There is no compelling reason for a separation of scales between features in the quiescent image and the impact of turbulence. Nevertheless, we will continue to make this assumption in the interest of computational expedience, providing only a qualitative assessment of the impact of turbulence on the closure-phase variability. As a result,

$$\begin{eqnarray}&&T({\boldsymbol{u}},t)=1+\displaystyle \frac{5}{4\bar{V}({\boldsymbol{u}})}\left[\left(\displaystyle \frac{2\pi i{\boldsymbol{u}}}{\lambda }\cdot \tilde{{\boldsymbol{\xi }}}\right)\ast \bar{V}\right]({\boldsymbol{u}},t).\end{eqnarray} \tag{ 20 }$$

As before, we assume that $\bar{V}({\boldsymbol{u}})$ takes the approximate form in Equation (15), in which case

$$\begin{eqnarray}&&\left[\left(\displaystyle \frac{2\pi i{\boldsymbol{u}}}{\lambda }\cdot \tilde{{\boldsymbol{\xi }}}\right)\ast \bar{V}\right]({\boldsymbol{u}},t)\approx \displaystyle \frac{2\pi i}{\lambda }{V}_{0}\;{\boldsymbol{u}}\cdot {\boldsymbol{\xi }}({\boldsymbol{u}},t).\end{eqnarray} \tag{ 21 }$$

Inserting this into Equation (12) then gives

$$\begin{eqnarray}\delta {{\rm{\Phi }}}_{123} & \approx & \displaystyle \frac{5}{4}{\boldsymbol{\delta }}{\boldsymbol{u}}\cdot {{\boldsymbol{\nabla }}}_{{\boldsymbol{u}}}\displaystyle \frac{2\pi u{\boldsymbol{u}}\cdot {\boldsymbol{\xi }}({\boldsymbol{u}},t)}{\lambda {u}_{0}}\\ & \approx & \displaystyle \frac{5}{2}\left(\displaystyle \frac{2\pi u}{\lambda }\xi \right)\left(\displaystyle \frac{u}{{u}_{0}}\right)\left(\displaystyle \frac{\delta u}{u}\right).\end{eqnarray} \tag{ 22 }$$

Again, we insert typical scales, setting $\delta u/u\approx 0.2$, ${u}_{0}/u\approx 0.2$, which require $\xi \approx 0.003\lambda /u$, i.e., turbulence displacements on scales on the order of $0.3\%$ of those probed by the Hawaii-SMT/CARMA baselines or $0.2\;\mu \mathrm{as}$. With strong gravitational lensing, and therefore compression of the features within the disk around the photon ring, this implies turbulent displacements on the order of 5% of the disk scale height, i.e., 5% density variations.

The difference in physical scale from the weak refraction model is due entirely to the different coupling to the background intensity field combined with the assumed scale hierarchy. For the weak refraction case, it is gradients of the quiescent intensity map that enter, which are limited to angular scales of $\lambda /{u}_{0};$ for the accretion turbulence case it is gradients of the displacement field, which exist on scales of $\lambda /u$ and are thus better matched to the interferometric measurements.

The small turbulence amplitude required appears problematic initially. However, this is the projected turbulence amplitude, obtained after integrating over the column, which will generically average down the impact of the turbulent variations. Strong turbulence (i.e., 100% local density variations) requires roughly 400 contributing turbulent eddies to produce the required 5% variations. Since the photosphere volume is roughly $\pi {r}^{3}$, where r is the orbital radius, this implies only $\approx 5$ turbulent eddies per scale height (commensurate with r).

The timescale for strong accretion-driven turbulent driven variability is ostensibly on the order of the orbital periods, which is roughly $0.5\;\mathrm{hr}$ at the ISCO of a non-spinning black hole in Sgr A*. However, this is a reasonably strong function of the radius, scaling as ${r}^{3/2}$. This size is well matched to the projected scales being probed by the Hawaii-SMT baseline, which has a nominal resolution of $60\;\mu \mathrm{as}$, corresponding to a linear scale of $10{GM}/{c}^{2}$, at which the orbital timescale is $1\;\mathrm{hr}$. For even moderately sub-Keplerian orbits, this can grow to the scales needed to produce the observed inter-epoch variations.