Detection of Intrinsic Source Structure at ∼3 Schwarzschild Radii with Millimeter-VLBI Observations of SAGITTARIUS A*

The a priori calibration of the visibility amplitudes was done using the system equivalent flux density (SEFD) measurements, which were determined by the antenna gain (K Jy⁻¹) and the opacity-corrected system temperature. Following previous EHT work (e.g., Fish et al. 2011; Doeleman et al. 2012; Lu et al. 2012, 2013; Akiyama et al. 2015), we performed a time-dependent station gain correction on the a priori calibrated amplitudes.²⁷

Our approach has three steps. First, we identified and adopted a total flux density of 3.1 Jy for Sgr A*. This is the average of the reported CARMA local interferometer measurements during the time of VLBI observations, which is consistent with measurements done in parallel at SMA.²⁸ The daily average total flux measurements suggest a change at the ∼10% level, which is the same order of the a priori calibration accuracy.

Second, with an intra-site baseline (e.g., DF) and a third station (e.g., S), one can accurately calculate the gain correction factors for the two co-site stations by assuming that the parallel hand visibilities on the intra-site baselines measure the total flux density of Sgr A* and the gain calibrated amplitudes on the other two baselines (e.g., SD/SF) are identical. The lengths of the intra-site baselines (92 m between the CARMA reference antenna (D/E) and the phased array (F/G) and 156 m between JCMT and phased SMA) are long enough to resolve out the diffuse thermal emission around Sgr A* on arcsecond scales (with projected baseline length >26 m at 1.3 mm, Fish et al. 2016), but short enough to not resolve the compact source itself at all. We refer to a triangle with such a short baseline as a pseudo closure amplitude triangle (Kawaguchi et al. 2015). For the data presented here, we have effectively two pseudo closure amplitude triangles and this calibration can be applied to stations including F/G and D/E (only at low band), and Q and J for the 30 s scans.

Third, we transfer the CARMA station gain correction factors obtained at low band to the high band and each 30 s Q/J correction factor to their following long scan for stations P and J. When transferring the gain correction factors, we have assumed a constant scaling factor during the time of our observations and we set the scaling factors such that the data at both bands from both parallel hands after calibration are self-consistent with each other. Given the known amplitude biases on baselines to station P between the two bands seen in earlier data (Lu et al. 2013), we used a different P/Q scaling factor for the low and high bands.

We note that the station gain correction factors for SMT and APEX cannot be calibrated in an "absolute" manner with the pseudo closure amplitude method. We have assumed that the gain correction factors for SMT and APEX are within 10%–20% from unity. Details of the a priori calibration for APEX can be found in Wagner et al. (2015). We added 10% systematic errors in quadrature to all the calibrated amplitude data to reflect the uncertainties of the a priori gain calibration, although there might still be remaining unaccounted systematic gain offsets for SMT and APEX, which we estimate to be <20%. Table 2 shows the amplitudes of Sgr A* after this calibration.

Our approach differs slightly from the network calibration procedure described in Johnson et al. (2015) due to, e.g., our improved understanding of the SMT gain that allowed us to directly use the SMT measurements. Nevertheless, our calibrated amplitudes are statistically consistent with the amplitudes reported by Johnson et al. (2015).

In the ensemble-averaged limit, the angular broadening of an image due to interstellar scattering is described by the convolution of the unscattered image with a scattering kernel, or equivalently by a multiplication in the Fourier domain. Following Fish et al. (2014), we corrected the visibility amplitudes before fitting models of the source structure by employing the scattering model determined by Bower et al. (2006). In this model, the major axis of the scattering kernel is oriented at 78° (east of north), with the associated full width at half maximum (FWHM) for the major and minor axes of ${\theta }_{\mathrm{FWHM}}^{\mathrm{maj}}=1.309$ (λ/1 cm)² mas = 22 μas and ${\theta }_{\mathrm{FWHM}}^{\min }=0.64$ (λ/1 cm)² mas = 11 μas, respectively. The correcting factors, which follow an elliptical Gaussian distribution in the Fourier domain, decrease monotonically in all directions from unity for the intra-site baselines to ∼0.37 for the longest baseline between APEX and phased SMA (see Pearson 1991, for the formula for calculating the correcting factors). We show the amplitudes of Sgr A* after this correction in Figures 2 and 3.

Zoom In Zoom Out Reset image size

In Figure 3 (right), we show the amplitudes on the APEX baselines. Due to scheduling and technical difficulties, the number of observed VLBI scans to APEX was lower than for other stations. On APEX-SMT and APEX-CARMA baselines, Sgr A* is detected with S/N in the range of 5–12. However, on the longer APEX-Hawaii baseline, Sgr A* is detected only in one scan. The S/N of the Sgr A* detection on this baseline is not very high (5.7 and 7.9 for the low and high band), but the low and high bands have very similar delay and delay rates and both are close to the values for detections of nearby AGN (e.g., OT 081 and PKS B1921−293) scheduled in adjacent scans to Sgr A*, which make a false detection highly unlikely. A fringe fitting over the combined low and high band data gives $\sqrt{2}$ sensitivity improvement and results in a firm detection of this scan. In Figure 3, the upper limits at other times for this baseline are also shown. In addition to Sgr A*, a few other sources (M 87, 3C 273, 3C 279, Centaurus A, 4C +38.41, OT 081, PKS B1921−293, and BL Lac) have been robustly detected on the APEX to Hawaii baseline during the same observing session. The data for these sources will be presented in forthcoming papers.

Zoom In Zoom Out Reset image size

We show the measured closure phases for Sgr A* in Table 3 and Figure 4. The "trivial" closure phases, which are formed on a triangle including an intra-site baseline, are consistent with zero (median: 02 ± 07 and mean: 01 ± 06), as expected, indicating a point-like structure at the arcsecond scale resolution provided by the intra-site baselines. On the CARMA–Hawaii–Arizona triangle, we found a median closure phase of 67 ± 15 and a mean of 78 ± 12 (consistent with an earlier analysis by Fish et al. 2016), with a trend of increase toward later times during an observing night (Figure 4). Following Fish et al. (2016), we can also rule out closure phase errors larger than ∼02 due to bandpass effects.

Zoom In Zoom Out Reset image size

We detected Sgr A* on the SMT–CARMA–APEX triangle, as well, with an S/N in the range of 4–8. The quality and uncertainties in the measured closure phases for this triangle are not sufficient to quantify their dependence, e.g., on sidereal time. However, we can infer their statistical properties by describing the measurements in terms of a Gaussian mixture model. We find that all 11 closure phase measurements in the SMT–CARMA–APEX triangle are consistent with having the same underlying value of ${5.0}_{-4.6}^{+12.9}$ degrees, where the uncertainties correspond to a 99.7% credibility interval (3σ).²⁹ Compared to the CARMA–Hawaii–Arizona triangle, which is open, small, and oriented mostly along the E–W orientation, the SMT–CARMA–APEX triangle is skinny, larger, and oriented mostly along the N–S orientation. The fact that the closure phases in the SMT–CARMA–APEX triangle are positive and comparable to those measured in the CARMA–Hawaii–Arizona triangle provides additional constraints on the properties of the structure probed by the various baselines of the array we use here, as we will show in the following section.

In principle, the VLBI structure of Sgr A* may be variable on short timescales (Fish et al. 2011; Lu et al. 2016; Medeiros et al. 2016; Rauch et al. 2016). However, we note that the uv coverage on the separate observing days is insufficient for a detailed imaging/modeling on a per day basis. We therefore, combined all data into a single data set based on the assumption that the VLBI structure does not change significantly during our observing campaign. This assumption is supported by an unpaired t-test, which shows an insignificant difference in closure phases between the different observing nights.³⁰ Additionally, we also calculated the compactness ratio $R={S}_{\mathrm{long}}/{S}_{\mathrm{short}}$ of the average correlated flux density on long (2.9–3.1 Gλ) and short baselines (400–700 Mλ) on a per day basis. We then performed a χ²-test to check for possible variability of that ratio. The resulting p = 0.81 excludes that R varies significantly. This is consistent with earlier observations of Fish et al. (2011), who showed stationarity of the source size over timescales of a few days, despite total flux density variations at a level of ≤20%, the latter being of the same order as the amplitude calibration accuracy and the level of total flux density variability observed in this experiment. Owing to the limited number of detections on each observing day, we therefore cannot make a strong statement about a possible underlying structural variability of Sgr A*. New data with better and more uniform uv coverage will be needed for such a statement.

Earlier 1.3 mm VLBI observations of Sgr A* with baselines primarily oriented in the E–W orientation justified no more than a circular Gaussian model fit, which results in an approximate source size of ${37}_{-10}^{+16}$ μas (3σ errors) for the brightness distribution of Sgr A* (Doeleman et al. 2008; Fish et al. 2011). The addition of the APEX telescope to the VLBI array allows us to constrain the source size also in the N–S direction and almost doubles the angular resolution. In fact, the measurement of a nonzero closure phase for the triangles SMT–CARMA–Hawaii and CARMA–SMT–APEX (both Figure 4 and Fish et al. 2016) and the observed visibility amplitudes (Figure 2) are both inconsistent with a simple elliptical Gaussian model, which is characterized by the different size in the two orthogonal directions (see also Johnson et al. 2015).

In order to demonstrate this, we show in Table 4 the results for a model fit with one elliptical Gaussian. The inferred size of the major axis is 52 ± 1 μas, the axial ratio is 0.4 ± 0.1, and the position angle of the major axis is 81 ± 3° east of north. However, this model provides a very poor fit to the data, as measured by the difference in the Bayesian information criterion (BIC, Liddle 2007) for this model, in comparison with the more asymmetric models, which are described in Section 4.3 (see Table 4). For the purposes of our work, we compute the BIC as

$$\begin{eqnarray}&&{\rm{BIC}}=k\mathrm{ln}(N)-2\mathrm{ln}P(\mathrm{data}| {\boldsymbol{w}}),\end{eqnarray} \tag{ 4 }$$

where k and N are the number of model parameters and data points, respectively. Of course, the measurement of nonzero closure phases is also inconsistent with other point-symmetric models, like a ring or a symmetric annulus. It requires more complex, asymmetric structures.

Table 4.  Model-fitting Results

	Flux Densitya	xb	yb	Sizec	Ratio	P.A.d	ke	ΔBICf
	(Jy)	(μas)	(μas)	(μas)		(degree)
Single elliptical Gaussian	3.0 ± 0.2	0	0	52 ± 1	0.4 ± 0.1	81 ± 3	4	980
Model A (S1)	0.9 ± 0.1	0	0	20 ± 1	...	...	6	14
(S2)	2.1 ± 0.2	31 ± 1	31 ± 3	51 ± 2
Model B (S1)	4.9 ± 0.6	0	0	52 ± 2	...	...	6	0
(S2)	−1.7 ± 0.4	1 ± 1	2 ± 1	25 ± 2

Notes.

^aFlux density for each component. ^bRelative position of each component in R.A. and decl. The total displacement of the two components is $d={({x}^{2}+{y}^{2})}^{1/2}$. ^cThe full width at half maximum (FWHM) of each Gaussian component. ^dThe position angle of the major axis for the elliptical Gaussian model in degrees east of north. ^eNumber of model parameters. ^fΔBIC = BIC-682.

Download table as:  ASCIITypeset image

It is important to emphasize here that the above result does not depend on the precise model of scattering induced blurring that we use. In deblurring the visibilities, we have assumed a Gaussian kernel at 1.3 mm with major and minor axis sizes extrapolated from the longer-wavelength measurements following a ${\lambda }^{2}$ law. The extrapolated sizes using the slightly different inferences by Bower et al. (2006) and Shen et al. (2005) differ by only 1 μas at 1.3 mm, which has little impact on our results. The scattering correction does depend weakly on the assumption of a Gaussian kernel, of isotropy, and of the quadratic λ² slope of the scattering law. There is evidence that, at millimeter-wavelengths, the diffractive scale becomes comparable to the inner scale of turbulence and, hence, the angular broadening scaling becomes steeper than λ² and that the shape of the kernel becomes non-Gaussian on long baselines (Gwinn et al. 2014; Johnson & Gwinn 2015). To assess the impact of these effects, we deblurred the calibrated data by considering extreme cases of reasonable inner scales of turbulence (Johnson 2016) and found that the model parameters do not change significantly.

In principle, the visibility amplitudes that we have measured on the various APEX baselines may not be caused by intrinsic source structure alone but rather affected by refractive scattering, which is caused by the interstellar medium, if the APEX baselines already resolve the underlying source structure. Refractive scattering effects on long VLBI baselines have been detected in Sgr A* at 1.3 cm by Gwinn et al. (2014). The properties of refractive substructure have been calculated by, e.g., Johnson & Gwinn (2015) and Johnson & Narayan (2016). In this subsection, we follow the work of Johnson & Narayan (2016) to demonstrate that the visibility amplitudes observed at the APEX baselines are unlikely to have been caused by refractive scattering effects, but they must represent (at least partially) evidence for intrinsic source substructure.

For an isotropic scattering medium, the amplitude of refractive noise (in units of the zero-baseline flux density) at a baseline that resolves the image is given by Equation (18) in Johnson & Narayan (2016)—see also Figure 7 and Equation (19) of Johnson & Gwinn (2015)—i.e.,

$$\begin{eqnarray}\sigma (b) & \simeq & \sqrt{\displaystyle \frac{{\rm{\Gamma }}(4/\alpha ){\rm{\Gamma }}(1+\alpha /2)}{{2}^{2-\alpha }{\rm{\Gamma }}(1-\alpha /2)}}{\left(\displaystyle \frac{{r}_{0}}{{r}_{F}}\right)}^{2-\alpha }{\left[\displaystyle \frac{b}{(1+M){r}_{0}}\right]}^{-\alpha /2}\\ & & \times {\left(\displaystyle \frac{{\theta }_{\mathrm{scat}}}{{\theta }_{\mathrm{img}}}\right)}^{2},\end{eqnarray} \tag{ 5 }$$

where Γ(x) is the gamma function, α is the power-law index of the turbulent power spectrum, r₀ is the phase decoherence length, r_F is the Fresnel scale, $M\equiv D/R$ is the magnification factor, D is the observer-screen distance, R is the screen-source distance,

$$\begin{eqnarray}&&{\theta }_{\mathrm{scat}}\equiv \displaystyle \frac{\sqrt{2\mathrm{ln}2}}{\pi }\displaystyle \frac{\lambda }{(1+M){r}_{0}}\end{eqnarray} \tag{ 6 }$$

is the scattered angular size of a point source, and ${\theta }_{\mathrm{img}}=\sqrt{{\theta }_{\mathrm{src}}^{2}+{\theta }_{\mathrm{scat}}^{2}}$ is the ensemble-average angular size of the source. At 1.3 mm, ${\theta }_{\mathrm{scat}}\simeq 22$ μas (Bower et al. 2006), whereas the present observations constrain the size of a potentially resolved source to ${\theta }_{\mathrm{src}}\simeq 52$ μas (see Table 4). For a screen at 2.7 kpc (as inferred by observations of the galactic center magnetar; Bower et al. 2015), the Fresnel scale is equal to ∼10¹⁰ cm. The APEX baselines have a length of ≃6 Gλ ≃ 8 × 10⁸ cm. Finally, under the hypothesis (which we try to reject) that the correlated flux in the APEX baselines is primarily due to refractive noise, these baselines have resolved the underlying image, i.e., $b\gt (1+M){r}_{0}$. Using this inequality and substituting the above values into Equation (5), we get

$$\begin{eqnarray}\sigma (b) & \leqslant & 0.065\sqrt{\displaystyle \frac{{\rm{\Gamma }}(4/\alpha ){\rm{\Gamma }}(1+\alpha /2)}{{2}^{2-\alpha }{\rm{\Gamma }}(1-\alpha /2)}}{\left(\displaystyle \frac{b}{8\times {10}^{8}\mathrm{cm}}\right)}^{2-\alpha }\\ & & \times {\left(\displaystyle \frac{{r}_{F}}{{10}^{10}\mathrm{cm}}\right)}^{\alpha -2}{\left(\displaystyle \frac{{\theta }_{\mathrm{scat}}}{22\mu \mathrm{as}}\right)}^{2}{\left(\displaystyle \frac{{\theta }_{\mathrm{img}}}{56\mu \mathrm{as}}\right)}^{-2}.\end{eqnarray} \tag{ 7 }$$

For a 3D Kolmogorov spectrum of turbulence, α = 5/3 and σ < 2%. This value depends weakly on α, since for α = 3/2, σ < 3% and for α = 1.9, σ < 1.5%. Note that these upper limits are very conservative since we have assumed that the APEX baselines have barely resolved the ensemble-average image of the source.

The visibility amplitudes that we detected on the APEX baselines is ∼4%–13% of the zero-baseline amplitude (see Table 2), while the dispersion of refractive noise at the same baselines is at most 2%–3%. Furthermore, for anisotropic scattering, the amplitude of refractive noise along the N–S direction (i.e., along the minor axis of the scattering kernel) is further suppressed by factors of a few (see, e.g., discussion in Johnson & Narayan 2016), making the estimate of the amplitude of the refractive noise for the APEX baselines at most 1%. Even though we have only measured one particular realization of the refractive noise over a narrow range of baselines, it is unlikely that we detected a distribution of amplitudes in these baselines at the 4%–13% level caused by refractive scattering effects that are expected to be at the ∼1% level. In other words, it is very unlikely that the visibility amplitudes that we measured in the APEX baselines have resulted predominantly from refractive noise, with no contribution from intrinsic source substructure.

In the previous subsections, we argued that our measurements of the correlated flux on the APEX baselines are indicative of source intrinsic emission on spatial scales comparable to ∼3 Schwarzschild radii. We now assess whether physically plausible structures are consistent with the data we report here.

The sparse uv coverage of our data does not warrant the reconstruction of a VLBI image of the inner accretion flow or fitting the parameters of detailed GRMHD models to the data. However, as discussed in Section 1, models of radiatively inefficient accretion flows predict images that are often shaped either like crescents, for disk-dominated models, or like disjoined compact emission regions at the jet footprints, for jet-dominated models. For this reason, we employ here two simple geometric models that were constructed in the past to capture the basic structural characteristics of images that are generated by complex GRMHD simulations (see, e.g., Kamruddin & Dexter 2013; Benkevitch et al. 2016; Medeiros et al. 2016 and references therein for a discussion).

Model A includes two displaced, positive, circular Gaussian components. It provides a simplified generic description of the expected image of the footprints of a jet (e.g., Chan et al. 2015), although a jet may appear more complex on event horizon scales than described by our Model A  (e.g., Mościbrodzka et al. 2014). It may also be used to describe the dominant emission from a tilted disk (Dexter & Fragile 2013) or a compact emission region, which is located off-axis to a more extended emission region, e.g., a hot spot in an accretion flow (Broderick & Loeb 2006). Because we cannot measure absolute phases in millimeter-VLBI, our data are only sensitive to relative positions. For this reason, we fix the center of one of the Gaussian components to the origin. Model A then requires six parameters: the flux density of each component, the R.A. and decl. displacement of the second Gaussian component, and the FWHM sizes of each component.

Model B is meant to provide a simplified generic description of the expected crescent-like emission around a BH for disk-dominated emission models. It is similar to the models of Kamruddin & Dexter (2013) and Benkevitch et al. (2016) in the sense that it is constructed as a difference between two displaced Gaussian components. The Gaussian taper of Model B produces smoother variations of the closure phases in comparison to a similar model consisting of uniform disks (Kamruddin & Dexter 2013). The sharp edges of the latter lead to steep gradients in the closure phase on SMT–CARMA–Hawaii, which are not observed. On the other hand, our model B has fewer parameters than the more complex geometric model of Benkevitch et al. (2016), as warranted by the limited uv coverage of our current data. Since Model B also involves two Gaussian components (albeit one with a negative normalization), it has the same number of parameters as Model A.

Table 4 summarizes the results from the model fitting. The single Gaussian model does not fit the data.³¹ For the two-component models A and B, we list the most likely values of their parameters (uncertainties reflecting formal errors) and the relative difference between their BIC. We note that the difference in BIC between models A and B depends largely on the treatment of the measurement uncertainties and the gain calibration. Although model B formally fits the data better, we cannot rule out model A, owing to the residual uncertainties in the error budget.

Figure 5 shows the relative sizes and orientations of the model components, the resulting u − v maps, and the locations on the maps where our measurements reside. It is important to emphasize here that we have assumed flat priors for all model distributions with no constraints on the space of model parameters. It is, therefore, instructive that the most likely characteristic sizes of both models, as measured by the component separation in Model A or the size of the larger component in Model B, are comparable to 50 μas, i.e., the expected size of the black hole shadow. This is consistent with the predictions of complex GRMHD simulations: either two jet footprints surrounding the black hole shadow for jet-dominated images or a crescent with a size comparable to that of the shadow for disk-dominated images (see references above). Moreover, in both models, the visibility amplitudes at the APEX baselines are generated by physical structures (i.e., jet footprint sizes or crescent widths) that are smaller than the size of the shadow.

Zoom In Zoom Out Reset image size