Imaging an Event Horizon: Mitigation of Source Variability of Sagittarius A*
The black hole in the center of the Galaxy, associated with the compact source Sagittarius A* (Sgr A*), is predicted to cast a shadow upon the emission of the surrounding plasma flow, which encodes the influence of general relativity in the strong-field regime. The Event Horizon Telescope (EHT) is a Very Long Baseline Interferometry (VLBI) network with a goal of imaging nearby supermassive black holes (in particular Sgr A* and M87) with angular resolution sufficient to observe strong gravity effects near the event horizon. General relativistic magnetohydrodynamic (GRMHD) simulations show that radio emission from Sgr A* exhibits variability on timescales of minutes, much shorter than the duration of a typical VLBI imaging experiment, which usually takes several hours. A changing source structure during the observations, however, violates one of the basic assumptions needed for aperture synthesis in radio interferometry imaging to work. By simulating realistic EHT observations of a model movie of Sgr A*, we demonstrate that an image of the average quiescent emission, featuring the characteristic black hole shadow and photon ring predicted by general relativity, can nonetheless be obtained by observing over multiple days and subsequent processing of the visibilities (scaling, averaging, and smoothing) before imaging. Moreover, it is shown that this procedure can be combined with an existing method to mitigate the effects of interstellar scattering. Taken together, these techniques allow the black hole shadow in the Galactic center to be recovered on the reconstructed image.
Subject headings:black hole physics – galaxies: individual (Sgr A*) – Galaxy: center – techniques: image processing – techniques: interferometric
The compact source at the Galactic center (Sgr A*) makes a very strong case that it is linked with a supermassive black hole, which due to its proximity (8 kpc) spans the largest angle on the sky among all known black holes (Melia & Falcke, 2001; Genzel et al., 2010; Falcke & Markoff, 2013). For Sgr A*, one Schwarzschild radius, , is 0.1 AU that subtends an angle of 10 as to us. According to general relativity (GR), a lensed image of the event horizon of Sgr A* (known as the “black hole shadow”) will appear (Bardeen, 1973; Luminet, 1979; Falcke et al., 2000; Takahashi, 2004) and can now be resolved by the Event Horizon Telescope (EHT), a project to assemble a VLBI network of millimeter wavelength dishes that aims to resolve general relativistic signatures in the vicinity of nearby supermassive black holes (Doeleman et al., 2008, 2009a; Broderick et al., 2009; Johannsen & Psaltis, 2010; Broderick et al., 2011a, b; Fish et al., 2011; Doeleman et al., 2012).
Horizon scale imaging promises to test basic predictions of GR and improves our understanding of the physics responsible for accretion and emission in a strong gravitational field. In particular, imaging a black hole shadow has been a long-standing goal of black hole astronomy. However, imaging the black hole shadow feature in Sgr A* has been inherently challenged by two known effects. First, the scattering by interstellar medium blurs the strong GR features near the black hole. In a recent work, it has been shown that this effect can be mitigated based on the fact that the scattering is well understood over the relative range of baseline lengths provided by the EHT (Fish et al., 2014). Second, while the predicted shadow feature is nearly independent of the spin or orientation of the black hole to within 10 % (Bardeen, 1973; Takahashi, 2004), the emission region surrounding the black hole depends on the details of the underlying accretion process and is intrinsically time variable primarily due to the stochastic nature of magnetorotational-instability-driven turbulence and magnetic reconnection in the accretion flow. Magnetorotational instability (MRI, Balbus & Hawley, 1991, 1998) is believed to be the leading mechanism driving turbulence in accretion disks and develops on orbital timescales. The timescale for the Keplerian motion at the innermost stable circular orbit around the black hole in Sgr A* ranges from 30 minutes for a non-rotating black hole to 4 minutes for prograde orbits around a maximally rotating black hole (Doeleman et al., 2009b). These time scales are much less than the typical duration of a Very Long Baseline Interferometry (VLBI) experiment, which violates one of the basic requirements for VLBI Earth-rotation aperture synthesis imaging. In contrast, the corresponding timescales in the nearby giant elliptical galaxy M87, which has the second largest apparent event horizon, are much larger (a minimal time scale of a few days).
In this paper, we show that the short-time scale structural variability of Sgr A* does not prevent construction of time-averaged images that contain distinguishable features like the black hole shadow. Section 2 describes the models we employed in this analysis and Section 3 details the data simulation, imaging analysis and quality assessment metrics. Our results on imaging strategy are presented in Section 4 and are discussed in Section 5. We summarize our conclusions in Section 6.
2. GRMHD Simulations of Sgr A*
We performed time-dependent simulations of black hole accretion using fully conservative 3D GRMHD code HARM (Gammie et al., 2003; Noble et al., 2006). The simulation assumes the flow is radiatively inefficient, as in time-independent, phenomenological models (e.g. Rees et al., 1982; Narayan & Yi, 1995; Yuan & Narayan, 2014), which is appropriate for low-luminosity galactic nuclei such as Sgr A* (see also Drappeau et al., 2013). The simulation starts from a geometrically thick hot disk with its pressure maximum at 24 surrounding a black hole with dimensionless spin . The simulation uses modified spherical-polar coordinates with logarithmically spaced radial grids spanning the range 1.2 to 240 and an azimuthal range of rad. The spatial resolution of the simulation is 260192128 cells in radial, poloidal, and azimuthal directions, respectively. The disk is initially seeded with a weak poloidal magnetic field that makes the disk unstable to the MRI. The simulation is run for 14,545 , which is long enough for saturation of the turbulence to be attained at the pressure maximum (this occurs at about 6,000 ).
To generate images, we perform general relativistic radiative transfer on the result of the GRMHD simulation using the bothros ray-tracing code (Noble et al., 2007). The main emission source at radio wavelengths is synchrotron radiation from the tenuous magnetized gas. The source function is integrated along geodesics that leads to each pixel of a “camera” which is placed 8 kpc from the model. A thermal distribution function is assumed for the electrons. There are several parameters that control the radiative properties of the model: black hole spin (), proton-to-electron temperature ratio (), and viewing angle (). Here we adopted the model with , , and for simulating observations of Sgr A*, which are consistent with existing mm VLBI and spectral measurements, but we also consider simulations with different parameters in the following sections. The length and timescales in the GRMHD model are set by the mass of the black hole, but the density of the accretion flow (equivalently: the accretion rate) is a free parameter. We adjust this free parameter so that the time averaged flux at 230 GHz after the MRI saturation (t = 6,000–14,545 ) gives 3.4 Jy, as observed by Marrone (2006). The dimensions of the camera frame of the movie are 210 210 as with a resolution of 256256 pixels.The interval between frames is 221.3 seconds, adding up to a total movie length of 53 hours. Figure 1 shows sample images from this simulation.
3.1. Data Simulation
VLBI observations were simulated using the MIT Array Performance Simulator (MAPS) software, following Lu et al. (2014, hereafter L14). Data were simulated at 230 GHz with a total bandwidth of 16 GHz assuming that the model images represent the mean flux over the entire bandwidth. The assumed 16 GHz bandwidth is consistent with the targeted recording bandwidth of near future EHT observations. As in L14, the full EHT array was used for simulation, which included the following sites: Submillimeter Array and James Clerk Maxwell Telescope on Mauna Kea, the Arizona Radio Observatory Submillimeter Telescope, the Combined Array for Research in Millimeter-wave Astronomy, the Large Millimeter Telescope (LMT), the Atacama Large Millimeter/submillimeter Array (ALMA), the Institut de Radioastronomie Millimétrique (IRAM) 30-m telescope on Pico Veleta, the IRAM Plateau de Bure Interferometer (PdBI), and the South Pole Telescope (SPT). For the simulations reported here, telescope elevations were restricted to be above 15 degrees (except 10 degrees for the PdBI), where calibration issues are expected to be reduced.
As a reference to the underlying quiescent structure, a static image was generated by averaging on a pixel basis all the frames in the GRMHD simulation. Then we generated synthetic VLBI data sets by calculating the complex visibilities and errors on each EHT baseline during a typical night of observing (with 12-hr long total time coverage). These data sets were generated using the static image, and also using the time evolving GRMHD movie with time resolution that matched the movie frame cadence. To simulate multi-epoch data sets, each consecutive block of 12 hrs of the movie (corresponding to 192 frames) was sampled by the array as one epoch with identical uv-coverage. The duration of the simulation allows 4 epochs of 12-hr long observations without overlap in frames. In order to further increase the number of observing epochs, we also considered a case where the input block of frames for each epoch were sampled with a halfway overlap, leading to a total of 8 epochs. Figure 2 illustrates how the movie was sampled over time in this eight-epoch case.
3.2. Imaging and quality assessment
Images were reconstructed with the BiSpectrum Maximum Entropy Method (BSMEM, Buscher, 1994) software, due primarily to its user-friendliness and speediness. We refer the reader to [and references therein]2014ApJ…788..120L for details concerning imaging reconstruction algorithms. Compared to the widely used deconvolution-based imaging techniques (e.g., CLEAN), forward imaging techniques, like the BSMEM, are well suited for mm-VLBI (L14;Fish et al. 2014).
As in L14, two image comparison metrics, i.e., mean square error (MSE) and structural dissimilarity (DSSIM) index, were applied to quantify the quality of the reconstructions. MSE compares the two images on a pixel-by-pixel basis and it is good for comparing all pixel intensities (Equation 1 in L14). Unlike MSE, DSSIM is derived from the human visual perception metric, structural similarity (SSIM, Wang et al., 2004) by DSSIM=(1/—SSIM —)-1. SSIM attempts to measure the change in the structural information between the two images by taking into account the change in luminance, contrast, and structure (see Equations 2–6 in L14 for details). In spite of this, these widely used metrics may not be perfectly suited to assessing the reconstruction quality when preserving the visual quality of salient features is crucial. Instead, visual inspection can often do a better job. In any case, future dedicated algorithms being able to detect and characterize specific image features and visual perception experiments appears very important. Here we use the average of all frames of the employed movie as a reference image. For both metrics, lower values indicate better reconstruction quality.
Figure 3 shows the average image (a) and the reconstruction of this static structure from a 12-hr observation with the assumed array (b). The reconstruction produces a high-fidelity image with critical features such as the photon ring and shadow preserved. However, a straightforward reconstruction of simulated data from a 12-hr movie is very unsuccessful (Figure 4, a), indicating that structural variability is a major hurdle to successful horizon-scale imaging of Sgr A*. Averaging of simulated data up to 8-epoch observations improves the reconstruction quality, but noisy features in the reconstruction can make recognition of critical features (e.g., the black hole shadow) difficult in practice (Figure 4, b).
Figure 5 (left, upper panel) compares the visibility amplitudes of the average image and the averaged amplitudes of the 8 epochs. The average was done in the complex plane for the visibilities with identical u-v coordinates. The data for the static image look smooth, while for the movie reconstruction there are “wiggles” in the visibilities. Closure phases, which are phases of triple products of the complex visibilities around closed triangles (Jennison, 1958; Rogers et al., 1974) and are calculated from the averaged visibilities of the 8 epochs, also show similar effects (Figure 5, right, upper panel). This is not surprising because the averaged visibilities are the average of Fourier components of different images. In addition, different movie frames were sampled at different u-v points. The visibilities at the various u-v points are thus inconsistent with both each other and themselves. Furthermore, all visibilities of the static reconstruction contain information about all movie frames, as the image under observation is the average of all frames. On the other hand, each visibility of the movie reconstruction only contains information of the seven or eight frames being averaged.
In order to reduce this effect, two subsequent data processing steps were applied: scaling and smoothing. The scaling is motivated by the observation that the brightness of pixels in the model images fluctuates in a highly correlated way. That is, there is a component of the variability that can be thought of as scaling the entire image up and down. The total flux information can thus be used to partially remove the variability. The remaining variable component of the structure can be treated as a high-frequency noise on top of the underlying quiescent image and therefore a Fourier smoothing algorithm can be used as a denoising technique.
For the first step, all visibility amplitudes were normalized by dividing each visibility by the total (zero spacing) flux density of the then observed frame (Figure 4, c). In practice, the total flux density can be obtained by continuously observing the source with a connected interferometer. In cases where a mismatch between the measured total flux and the zero-spacing flux density of the horizon-scale structure exists, the scaling factors could possibly be determined by how well the scaling works on short baselines. Scaling significantly reduces the irregularities in the amplitudes on short baselines ( 2 G), but does not change the closure phases (Figure 5, middle panels). After the normalization, significant deviations from the static reconstruction in the amplitudes still exist on baselines longer than G ( G for closure phases in Figure 5). A smoothing algorithm was then applied to make the visibilities on baselines longer than G resemble the reconstruction of the averaged image. This smoothing algorithm is a moving average: each new data point is the average of all old data points within a certain time window, centered on the timestamp of the current data point. This was done in the complex plane for each baseline separately. In order to be able to set a large enough window without losing too much information on large timescales, the data was convolved with a Gaussian weighting function. The smoothing is thus in fact a low pass filter. High frequency structures are averaged out, while longer existing structures are preserved.
The “cutoff frequency” of the filter is determined by the standard deviation of the Gaussian. If the cutoff frequency is too high, the wiggles in the data will still be followed closely by the smoothed data. On the other hand, if the cutoff frequency is too low the smoothing outcome will be a flat line corresponding to the time-averaged visibility on a particular baseline. A standard deviation of 100 data points and a window size twice as large gave the best correspondence to the static reconstruction data. With an integration time of 20s ( ), this standard deviation corresponds to 2000s or 100 , which corresponds to a few ISCO (Innermost Stable Circular Orbit) rotation periods for a maximally rotating black hole. The smoothing algorithm thus filters out variability on shorter timescales. It is worth pointing out that time-averaging of the visibility data smears out the response of a point source located away from the center of the field of view (time smearing effect). Employing an averaging time of 2000 seconds would lead to a fall off of 10 % in the response to the flux 0.12 mas from the phase center for a 10000 km baseline. Since the emission of Sgr A* is known to be very compact (within 0.04 mas, Doeleman et al., 2008; Fish et al., 2011), this effect is very small.
After smoothing, the visibility amplitudes and closure phases of the movie reconstruction show much more overlap with those of the static reconstruction than before (Figure 5, lower panels). The resulting reconstruction is shown in Figure 4 (d), which shows the same black hole features as the static reconstruction (Figure 3), but slightly less prominent, as indicated by the increase in MSE and DSSIM.
As shown by Fish et al. (2014), the effects of interstellar scattering can be mitigated by correcting the visibilities for the scattering kernel before imaging. As expected, this technique works well for a static structure (Figure 6, a–c). We applied this technique to the movie reconstruction by convolving all movie frames with the scattering kernel from Bower et al. (2006) before the MAPS simulation. After the observation, all visibilities were divided by the Fourier transform of the scattering kernel at their particular u-v coordinates. Because de-blurring amplifies the thermal noise on long baselines, the degree to which de-blurring works depends on the noise at each telescope, and also on the baseline length. The order in which averaging, scaling and deblurring are performed is irrelevant as they are linear operations, but smoothing was always performed as the last step before imaging. In Figure 6 (d), the reconstruction from the movie is able to clearly recover the characteristic signature of the black hole: the photon ring. From the MSE and DSSIM values, it is clear that scattering and deblurring only marginally decrease the image quality (Figure 4 (d) and 6).
The above strategies are indispensable to each other for imaging the horizon-scale signatures of Sgr A*. Scaling mainly affects the large-scale structure (short baselines), whereas smoothing and deblurring mainly affects the small-scale structure (long baselines). When both techniques are applied to the data, the more visibilities are averaged, the better the quiescent structure is approached (Figure 7 and Table 1). In addition, averaging visibilities will also mitigate refractive noise (i.e., deviations from ensemble-average scattering, Johnson & Gwinn, 2015). In practice, the source can simply be observed for multiple days to fulfill this requirement. Since observing the scattered movie and deblurring the visibilities makes little difference in the final image, the characteristic shadow and photon ring of the black hole can indeed be recovered in Sgr A* with the short timescale source variability and interstellar scattering present.
In the current work we have restricted ourselves to the reconstruction of a static structure out of a GRMHD simulation and have thus used the averaged image as a reference. By averaging in time, however, some of the strong GR effects presented in the GRMHD simulation are smeared out and therefore future horizon-scale imaging should extend the present work to the reconstruction of a variability image, featuring not only the black hole shadow, but also details of the turbulent accretion flow.
The appearance of the horizon-scale image depends on the black hole properties (mass and spin) and details of the accretion structure and process. Most of these are currently still very uncertain. In this work, we have only considered one time-dependent model to explore the variability mitigation strategy for horizon-scale imaging. Models with different black hole spin, proton-to-electron temperature ratio, and inclination angle may be used to explore the effect of these parameters on the observations, with the knowledge of black hole mass and mass accretion rate. A parameter survey by Mościbrodzka et al. (2009) favors 0.94, , and , respectively, though the parameters will be model-dependent (e.g., Dexter et al., 2010; Broderick et al., 2011a).
Our procedure, however, is not limited by the uncertainties in the parameter space. As an example, Figure 8 shows a reconstruction of the model movie with 0.94, , and assuming 8-day observations. With a higher inclination angle, the approaching side of the accretion flow becomes brighter and the receding side becomes darker (almost invisible) due to Doppler effects. In this case, the model movie can still be fairly well reconstructed. This suggests that comparison of future EHT observations with simulated observations may tightly constrain model parameters, but to obtain a quantitative comparison new algorithms would have to be developed to detect and characterize features in the image such as the shadow size and position.
Our simulations have implicitly assumed that the visibility amplitude and phase can be measured and calibrated. Accurate amplitude calibration has traditionally been challenging at (sub)millimeter wavelengths, especially when the array is small and limited in sensitivity. However, with more stations being added to the array and planned increases in the data recording rate, the effects of calibration errors should diminish in the near future (see Fish et al., 2014, for more discussions on the potential improvement on amplitude calibration).
On the other hand, fringe phase could also possibly be corrupted by the fluctuations in the atmospheric path lengths. With an array of telescopes of three or more, however, the closure phase, which is inherently robust against station-based phase errors, can be used to retrieve phase information. The fraction of phase information retained by the closure phase monotonically increases with the number of telescopes in the array as , where is the number of telescopes. Thus, with the anticipated EHT array of eight stations, 75 % phase information will be recovered and can be used with measured amplitudes to generate visibilities. In practice, visibility amplitude and closure phase information can be measured in terms of incoherently averaged quantities with well established algorithms to overcome coherence losses (Rogers et al., 1995). Recent EHT observations have shown that both amplitude and closure phase can indeed be successfully measured on Sgr A* (Doeleman et al., 2008; Fish et al., 2011, 2015).
It is worth pointing out, however, that the triple products from the visibilities averaged over multiple days are not the same as averaging the triple products in the complex plane. Figure 9 compares the triple product amplitudes and closure phases for these two cases. The average triple amplitudes and triple amplitudes from the averaged visibilities are similar everywhere. The closure phases, however, are consistent only on small triangles and begin to differ on large triangles (triangle longest leg 4G), where the triple products start rotating rapidly in the complex plane. As a result, using the averaged triple products, the black hole features are less well-preserved (Figure 10), indicating sufficient visibility phase information is critical for proper imaging of black hole features.
The source of Sgr A*’s submillimeter variability is not well understood. Although there will almost certainly be turbulent variability, as in our GRMHD model, other mechanisms can cause variability on length and timescales detectable by EHT. Variability may be caused by orbiting hot spots (Doeleman et al., 2009b, and references therein), jets (Mościbrodzka et al., 2014), tilted disks (Dexter & Fragile, 2013) or episodic particle acceleration. Our proposed imaging technique may not function equally well in all these cases. If the variability is very rapid in the uv domain, for example, the width of the Gaussian smoothing kernel will need to be adapted accordingly. It will therefore be important to test our technique on as large a universe of theoretical models as possible in future studies.
Our ability to image the horizon-scale signatures of the black hole depends on the properties of the observing array. In Figure 11 we show the reconstruction degradation compared to what is obtainable with the full array when a given site is unavailable, e.g., the phased CARMA, Pico Veleta and PdBI, or phased ALMA. A visual inspection indicates that the most severe degradation happens when the phased ALMA is missing (panel d). This is because all the longest and most sensitive baselines are provided by the phased ALMA. However, both the MSE and DSSIM statistics (Table 2) do not confirm the same assessment as perceived by human observers, indicating that these pixel-based metrics provide little understanding on how the morphology of black hole features differs from image to image. Future development of feature-based metrics (i.e., metrics that characterize the morphological properties of black hole features) can potentially provide a more unbiased way for black hole image comparison. Figure 12 shows the recording bandwidth impact on the reconstruction fidelity from 4 to 16 GHz by powers of two. As suggested by the MSE, the image quality gets better with wider recording bandwidths (i.e., higher sensitivity), although the DSSIM does not follow the same trend. Since the bandwidth () and integration time () equivalently improve the signal-to-noise ratio of a coherently integrated signal as , it is preferable to record data with maximum bandwidth, as the atmospheric coherence time is usually short (10s).
We have shown that the variability of Sgr A* at 1.3 mm can be significantly mitigated and the general relativistic black hole features, such as the shadow and photon ring, predicted by a GRMHD model movie of Sgr A*, can in principle be imaged by the EHT. To get a high-quality image, it is essential to observe Sgr A* for multiple days and to average visibilities in the complex plane before imaging. Normalizing the visibilities with respect to the zero-spacing flux density, which can be measured in practice, is an important tool to obtain high image quality, especially for large-scale structures. Applying a smoothing algorithm and increasing the observation time will further increase the image quality. If the properties of the scattering kernel are well known, the reconstructions can be corrected for interstellar scattering.
The inclusion of phased ALMA in future observations and recording with wide bandwidth will be critical for imaging the black hole shadow. In order to detect and characterize the black hole shadow features, development of dedicated algorithms is needed in near future studies. Given the currently limited understanding of the origin of flaring structures in Sgr A*, it is also important to explore a wider range of time-dependent source models to improve the capabilities of the proposed imaging techniques.
- Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
- Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1
- Bardeen (1973) Bardeen, J. M. 1973, Black Holes, ed. C. DeWitt & B. S. DeWitt (New York: Gordon and Breach), 215
- Bower et al. (2006) Bower, G. C., Goss, W. M., Falcke, H., Backer, D. C., & Lithwick, Y. 2006, ApJL, 648, L127
- Broderick et al. (2009) Broderick, A. E., Fish, V. L., Doeleman, S. S., & Loeb, A. 2009, ApJ, 697, 45
- Broderick et al. (2009) Broderick, A. E., Loeb, A., & Narayan, R. 2009, ApJ, 701, 1357
- Broderick et al. (2011a) Broderick, A. E., Fish, V. L., Doeleman, S. S., & Loeb, A. 2011a, ApJ, 735, 110
- Broderick et al. (2011b) Broderick, A. E., Fish, V. L., Doeleman, S. S., & Loeb, A. 2011b, ApJ, 738, 38
- Buscher (1994) Buscher, D. F. 1994, in IAU Symp. 158, Very High Angular Resolution Imaging, ed. J. G. Robertson & W. J. Tango (Dordrecht: Kluwer), 91
- Dexter et al. (2010) Dexter, J., Agol, E., Fragile, P. C., & McKinney, J. C. 2010, ApJ, 717, 1092
- Dexter & Fragile (2013) Dexter, J., & Fragile, P. C. 2013, MNRAS, 432, 2252
- Doeleman et al. (2008) Doeleman, S. S., Weintroub, J., Rogers, A. E. E., et al. 2008, Nature, 455, 78
- Doeleman et al. (2009a) Doeleman, S., Agol, E., Backer, D., et al. 2009, astro2010: The Astronomy and Astrophysics Decadal Survey, 2010, 68
- Doeleman et al. (2009b) Doeleman, S. S., Fish, V. L., Broderick, A. E., Loeb, A., & Rogers, A. E. E. 2009, ApJ, 695, 59
- Doeleman et al. (2012) Doeleman, S. S., Fish, V. L., Schenck, D. E., et al. 2012, Sci, 338, 355
- Drappeau et al. (2013) Drappeau, S., Dibi, S., Dexter, J., Markoff, S., & Fragile, P. C. 2013, MNRAS, 431, 2872
- Falcke & Markoff (2013) Falcke, H., & Markoff, S. B. 2013, Classical and Quantum Gravity, 30, 244003
- Falcke et al. (2000) Falcke, H., Melia, F., & Agol, E. 2000, ApJ, 528, L13
- Fish et al. (2011) Fish, V. L., Doeleman, S. S., Beaudoin, C., et al. 2011, ApJ, 727, L36
- Fish et al. (2014) Fish, V. L., Johnson, M. D., Lu, R.-S., et al. 2014, ApJ, 795, 134
- Fish et al. (2015) Fish, V. L., Johnson, M. D., Doeleman, S. S., et al. 2015, ApJ, submitted.
- Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444
- Genzel et al. (2010) Genzel, R., Eisenhauer, F., & Gillessen, S. 2010, Reviews of Modern Physics, 82, 3121
- Jennison (1958) Jennison, R. C. 1958, MNRAS, 118, 276
- Johannsen & Psaltis (2010) Johannsen, T., & Psaltis, D. 2010, ApJ, 718, 446
- Johnson & Gwinn (2015) Johnson, M. D., & Gwinn, C. R. 2015, arXiv:1502.05722
- Lu et al. (2014) Lu, R.-S., Broderick, A. E., Baron, F., et al. 2014, ApJ, 788, 120 (L14)
- Luminet (1979) Luminet, J.-P. 1979, A&A, 75, 228
- Marrone (2006) Marrone, D. P. 2006, Ph.D. Thesis,
- Melia & Falcke (2001) Melia, F., & Falcke, H. 2001, ARA&A, 39, 309
- Mościbrodzka et al. (2009) Mościbrodzka, M., Gammie, C. F., Dolence, J. C., Shiokawa, H., & Leung, P. K. 2009, ApJ, 706, 497
- Mościbrodzka et al. (2012) Mościbrodzka, M., Shiokawa, H., Gammie, C. F., & Dolence, J. C. 2012, ApJ, 752, L1
- Mościbrodzka et al. (2014) Mościbrodzka, M., Falcke, H., Shiokawa, H., & Gammie, C. F. 2014, A&A, 570, A7
- Narayan & Yi (1995) Narayan, R., & Yi, I. 1995, ApJ, 452, 710
- Noble et al. (2006) Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626
- Noble et al. (2007) Noble, S. C., Leung, P. K., Gammie, C. F., & Book, L. G. 2007, Classical and Quantum Gravity, 24, 259
- Rees et al. (1982) Rees, M. J., Begelman, M. C., Blandford, R. D., & Phinney, E. S. 1982, Nature, 295, 17
- Rogers et al. (1974) Rogers, A. E. E., Hinteregger, H. F., Whitney, A. R., et al. 1974, ApJ, 193, 293
- Takahashi (2004) Takahashi, R. 2004, ApJ, 611, 996
- Wang et al. (2004) Wang, Z., Bovik, A. C., Sheikh, H. R., et al. 2004, IEEE Trans. Image Process., 13, PP. 600, 612
- Rogers et al. (1995) Rogers, A. E. E., Doeleman, S. S., & Moran, J. M. 1995, AJ, 109, 1391
- Yuan & Narayan (2014) Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529