The Event Horizon Telescope: exploring strong gravity and accretion physics
Abstract
The Event Horizon Telescope (EHT), a global submillimeter wavelength very long baseline interferometry array, is now resolving the innermost regions around the supermassive black holes Sgr A* and M87. Using black hole images from both simple geometric models and relativistic magnetohydrodynamical accretion flow simulations, we perform a variety of experiments to assess the promise of the EHT for studying strong gravity and accretion physics during the stages of its development. We find that (1) the addition of the LMT and ALMA along with upgraded instrumentation in the “Complete” stage of the EHT allow detection of the photon ring, a signature of Kerr strong gravity, for predicted values of its total flux; (2) the inclusion of coherently averaged closure phases in our analysis dramatically improves the precision of even the current array, allowing (3) significantly tighter constraints on plausible accretion models and (4) detections of structural variability at the levels predicted by the models. While observations at 345 GHz circumvent problems due to interstellar electron scattering in lineofsight to the galactic center, short baselines provided by CARMA and/or the LMT could be required in order to constrain the overall shape of the accretion flow. Given the systematic uncertainties in the underlying models, using the full complement of two observing frequencies (230 and 345 GHz) and sources (Sgr A* and M87) may be critical for achieving transformative science with the EHT experiment.
keywords:
accretion — accretion discs — relativistic processes — black hole physics — galaxy: centre — submillimetre — techniques: interferometric1 Introduction
The Event Horizon Telescope (EHT) is a very long baseline interferometry (VLBI) array consisting of radio telescopes around the world.^{1}^{1}1http://www.eventhorizontelescope.org Once complete, it will be capable of achieving a resolution of up to 23 microarcseconds at 230 GHz and 15 microarcseconds at 345 GHz. This resolution is required, as its name suggests, to resolve black hole structure on event horizon scales. Both the supermassive black hole at the center of the Milky Way, Sagittarius A*, and that of the elliptical galaxy M87 are massive enough and close enough (subtending tens of microarcseconds) for the EHT to resolve. Sgr A* is predicted to subtend a larger angle on the sky, and with a declination of , it remains permanently visible in the Southern Hemisphere. However, radio images of Sgr A* are blurred due to interstellar scattering by turbulent electrons (Backer, 1978; Bower, 2006), which is consistent with a thin scattering screen located near the Scutum spiral arm (Bower et al., 2014). Supergiant elliptical galaxy M87, known for its striking jet in the optical, provides an alternate source–it has a black hole nearly as large in angular size, and without any indication of strong interstellar scattering at 7mm (e.g., Hada et al., 2011). Initial EHT observations have made the first detections of event horizon scale structure using both sources (Doeleman et al., 2008; Fish et al., 2011; Doeleman et al., 2012). The array capabilities are rapidly improving both in sensitivity and coverage due to increased bandwidth and the inclusion of additional radio telescopes (e.g., the Atacama Large Millimeter/submillimeter Array (ALMA), Fish et al., 2013).
Here we assess the potential for future EHT observations at various milestones in its development, taking into account the locations and sensitivities of the EHT’s constituent telescopes. We attempt to answer a variety of observational questions. Can the EHT detect signatures of strong field general relativity? How well can we constrain models for black hole images? Should we preferentially observe Sgr A* or M87? What are the tradeoffs between observing at 345 GHz versus 230 GHz? We perform a set of experiments using geometric and theoretical models for black hole images (§2) and our simulated arrays (§3). We statistically quantify the performance of each simulated observing run (§4), and discuss the implications and limitations of this study (§5) before summarizing our results (§6).
2 Modeling
The sparse uvcoverage of the current EHT array prevents imaging (but see Lu et al., 2014, for an example of image formation with the future EHT array). Maximizing the science return on the data or simulating future observations requires fitting models in the Fourier domain. Both geometric shapes and physical emission models have been used to interpret the observations. The mm emission from Sgr A* and M87 is thought to be synchrotron radiation (e.g., Reynolds & McKee, 1980; Melia, 1994; Falcke et al., 1998), either from a radiatively inefficient accretion flow (Narayan, Yi & Mahadevan, 1995; Yuan, Quataert & Narayan, 2003; Reynolds et al., 1996) or a jet (Falcke & Markoff, 2000; Broderick & Loeb, 2009a). Recently, there has been considerable interest in making predictions for the observed event horizon scale structure, both using semianalytic (Broderick & Loeb, 2009b; Yuan et al., 2009, e.g.,) and MHD (e.g., Mościbrodzka et al., 2009; Dexter, Agol & Fragile, 2009; Shcherbakov, Penna & McKinney, 2012) accretion flow and jet models. Thus far both geometric and physical models can produce satisfactory fits (e.g., Broderick & Loeb, 2009a; Dexter et al., 2010; Fish et al., 2011), but physicallymotivated ones are statistically favored (Broderick et al., 2011b; Kamruddin & Dexter, 2013). In this work, we use both physical models from numerical simulations of the innermost regions of accretion flows and jets and a parameterized geometric description of their resulting “crescent” images.
Name  Source  a  i ()  Type  Field  References (simulation, image) 

90h  Sgr A*  0.9  60  Disk  Poloidal  Fragile et al. (2007); Dexter, Agol & Fragile (2009); Dexter et al. (2010) 
MBD  Sgr A*  0.92  60  Disk  Dipolar  McKinney & Blandford (2009); Dexter et al. (2010) 
MBQ  Sgr A*  0.94  50  Disk  Quadrupolar  McKinney & Blandford (2009); Dexter et al. (2010) 
515h  Sgr A*  0.5  40  Tilted Disk  Poloidal  Fragile & Meier (2009); Dexter & Fragile (2013) 
DJ1  M87  0.92  25  Disk/jet  Dipolar  McKinney & Blandford (2009); Dexter, McKinney & Agol (2012) 
J2  M87  0.92  25  Jet  Dipolar  McKinney & Blandford (2009); Dexter, McKinney & Agol (2012) 
2.1 Simulated Black Hole Images
The theoretical black hole images used here are the result of relativistic radiative transfer calculations from MHD simulations of accretion flows (e.g., Dexter et al., 2012). We use images corresponding to the global best fit to submm Sgr A* data (spectral and EHT) both from aligned simulations (Fragile et al., 2007; McKinney & Blandford, 2009) and from those where the angular momentum of the accreting material was misaligned from the black hole spin axis by (Fragile et al., 2007; Fragile & Meier, 2009). More details about the calculations and fitting procedure can be found in Dexter et al. (2010, 2012); Dexter & Fragile (2013). Some basic properties of the models we use are given in Table 1, including the black hole spin , inclination , the initial magnetic field configuration, and the portion of the gas distribution which dominates the observed submm emission. The images are shown in Figure 1. These simulated images are necessarily a small subset of those in the literature. However, as discussed below the predicted images are strongly shaped by relativistic effects, and so we believe the results of simulated EHT observations using these models will be similar to those using images calculated by other groups (e.g., Bromley, Melia & Liu, 2001; Noble et al., 2007; Mościbrodzka et al., 2009; Broderick et al., 2011b; Shcherbakov, Penna & McKinney, 2012). These are taken as a range of plausible models for Sgr A* and M87, but are not intended to span the true parameter space of possibilities for the outcome of EHT observations.
2.2 Geometric Parameterization
The theoretical models described thus far to calculate event horizon scale black hole images either from semianalytic models or MHD simulations suffer from systematic uncertainties due to incomplete physics and uncertain initial conditions. For this reason, we also employ a geometric crescent (+ ring) model for the image of an accretion flow. The crescent is the difference between two constant intensity disks (Kamruddin & Dexter, 2013). This model aims to capture the dominant relativistic effects near black holes: Doppler beaming and gravitational light bending. The accretion flow moving toward our line of sight is “beamed” to a higher flux, leading to asymmetry in the image, while light from the back side of the accretion flow is deflected above and below the black hole to the observer. These combined effects lead to a crescent image morphology (e.g., Bromley, Melia & Liu, 2001, and see Figure 1).
A geometric crescent model is defined by five parameters (Kamruddin & Dexter, 2013): , the radius of the outer edge of the crescent; , the fractional radius of the hole in the crescent; , the position angle, defined as the angle between the uaxis and the line between the center of the crescent’s hole and the center of the crescent; , the total flux of the crescent; and , the displacement of the crescent’s hole from the center, parameterized as a fraction of . The asymmetry parameter can be thought of as a proxy for the strength of Doppler beaming, and in the context of standard accretion flow models as a proxy for the inclination of the accretion disk with respect to our line of sight. A faceon disk should have symmetric flux (, at least when averaged over long timescales), while an edgeon disk can be highly asymmetric from Doppler beaming ().
In addition to the accretion flow, captured by the crescent, general relativity predicts that black holes should cast a “shadow” corresponding to the region where photon orbits reaching the observer are bound to the black hole (Bardeen, Carter & Hawking, 1973; Falcke, Melia & Agol, 2000). Outside of this location, optically thin images show a bright ring corresponding to the unstable circular photon orbit, located at a physical radius of . Detection of this photon ring is tantamount to the discovery of the shadow. In addition, the size of the ring allows a measurement of the mass and distance of the black hole (Johannsen & Psaltis, 2011), while its shape can in principle be used to test the “nohair theorem” of general relativity (Johannsen & Psaltis, 2010). Simulated models (Dexter, Agol & Fragile, 2009; Dexter et al., 2010; Dexter & Fragile, 2013) predict that this ring could contribute of the total flux, estimated by extracting the fraction of the flux concentrated just outside the location of the predicted location of the lensed photon orbit.
The ring component, when added, has only one parameter: , its total flux. It is parameterized as a delta function fixed at a radius of , the lensed orbit for the “photon ring” of a nonspinning black hole. Although the photon ring becomes highly asymmetric at high values of black hole spin (Bardeen, Carter & Hawking, 1973), its size only varies by . This component is predicted by general relativity, but its existence has not yet been confirmed via observation. Given its physical extent, we determine the angular size of this ring from the mass and distance of each of our target objects, whose basic parameters are listed in table 2. Note that while recent gasdynamical studies of M87 suggest a black hole mass of (Walsh et al., 2013), we use the larger value determined by stellardynamical studies.
Parameter  Sgr A*  M87 

()  
(kpc)  
Diameter of Shadow (as)  53.0  37.8 
Right Ascension  17h45m40.03599s  12h30m49.42338s 
Declination  2900’28.1699”  +1223’28.0439” 
The complete Fourier Transform of the crescent model, , where and are in units of , can be expressed analytically as:
(1) 
where and are the Fourier transforms for a crescent and ring,
(2)  
(3) 
is the Fourier transform of a constant intensity disk of size R at radial spatial frequency ,
(4) 
and are Bessel functions of the first kind. Finally, is the Fourier transform of the Gaussian, which accounts for the blurring effect of strong interstellar scattering towards Sgr A*, using empirical values from Bower (2006).
2.3 Making Measurements
The coordinates sampled by a given EHT array depend on the location of the telescopes, the time of year, and the position of the source. The calculations to determine appropriate coordinates are carried out using the MIT Array Performance Simulator^{2}^{2}2http://www.haystack.mit.edu/ast/arrays/maps/ (MAPS), a program designed to simulate radio interferometry observations. Each observing run is a 24 hour time period set on the April 5th of 2009, with 10 minute scans onsource beginning every 20 minutes. (The actual date is unimportant, since we simulate observations for the full 24 hour period and do not consider the effects of the Sun or other celestial objects in our analysis.) We assume 10 second integrations, which are limited by the atmosphere coherence time. Each telescope is assumed to have elevation limits of between 15 and 85 degrees, with the exception of PdBI, which can observe down to 10 degrees. Gaussian noise is applied to the real and imaginary components of these measurements, characterized by a standard deviation of , where is the System Equivalent Flux Density of each telescope in the baseline, is the bandwidth, and is the integration time (Thompson, Moran & Swenson, 2001). The observables we consider, visibility amplitudes and closure phases, are insensitive to atmospheric phase delays which we ignore.
2.3.1 Real, Imaginary, and Amplitudes
EHT data thus far have almost exclusively consisted of visibility amplitudes. However, most of the information in an image is encoded in its phase. One peculiarity of the absence of phase information is a degeneracy due to the fact that . Unfortunately, atmospheric phase delays make it challenging to record Fourier phases using VLBI. Hence, we often simulate measurements comprised only of visibility amplitudes by taking the magnitude of simulated real and imaginary visibilities. Such a measurement is assigned a Gaussian error bar equal to the same error quoted above in the context of visibilities, which follows from linear error propagation theory. Amplitude measurements have approximately Gaussian errors, so normal statistics can be applied to their measurements. The appropriate errors in closure phase measurements are described below.
2.3.2 Closure Phases
While atmospheric phase delays make measurements of Fourier phases nearly impossible, closure phases of Sgr A*, the sum of phases between three different telescopes, have been successfully measured (Fish et al., 2011) and will constitute a crucial EHT observable. Although closure phases do not include all of the full Fourier phase information, they are insensitive to all sitespecific errors such as atmospheric and hardware delays (Rogers, 1976). Because of this, closure phases can be coherently averaged for much longer than the atmospheric coherence time.
A straightforward linear Gaussian error propagation exercise would imply that the error on each phase should be , where is the signaltonoise ratio, and that these errors could be added in quadrature when measuring a closure phase. Yet in the limit where is small, phase noise is nonGaussian, and we notice discrepancies from Gaussian statistics. For simulated observations of amplitudes and closure phases, we calculate amplitudes and their errors as above, then assign an error on the phase measurement based on the probability distribution for phase noise (Rogers et al., 1984).
(5)  
In the limit of large , this converges to
(7) 
a normal distribution with a standard deviation from linear Gaussian error propagation formulas. Note that we have corrected errors found in both of these equations in our source, Rogers et al. (1984).
To minimize computational time, we calculate only a minimal set of the independent closure phases, selecting those which can measured with the greatest precision. Hence, for each array, we designate the telescope with the minimum SEFD as the fixed telescope with which all closure triangles are made. Then, for each simultaneous measurement, phases are summed for each triangle containing the fixed telescope. The effective signaltonoise ratios of the resultant closure phases are calculated using equation 8. This is equal to the inverse of the fractional closure phase uncertainty.
(8) 
where is given by
(9) 
and and are hyperbolic Bessel functions.
In the high signaltonoise limit, equation 8 asymptotes to
(10) 
which would be the familiar result from linear Gaussian error propagation theory.
Because closure phases can be coherently averaged for much longer than the atmospheric coherence time ( s), we allow them to be averaged for the the full integration time, minutes, as in Broderick et al. (2011a). This has the effect of shrinking the error bars on closure phases by a factor of under our assumptions. In general, is limited by the inherent variability of the source (unknown a priori).
2.4 Fitting
Using this modeling procedure, we can generate best fit static ringed and ringless crescents to 230 GHz observations of Sgr A* (Doeleman et al., 2008; Fish et al., 2011) and M87 (Doeleman et al., 2012). These fits provide us a template for the truth images used in several of our experiments. Best fit parameters and uncertainties are calculated using the MCMC fitting algorithm described in Goodman & Weare (2010). Our MCMC ensemble is comprised of 100 walkers and is run through 1000 trials ( total samples). The initial configuration of the walkers is a Gaussian distribution in parameter space, instantiated with realistic values and widths corresponding roughly to the uncertainties in each parameter (as determined by previous runs). We verify that the our results are insensitive to these initial conditions. The tunable acceptance parameter, which Goodman & Weare (2010) denote as , is set to 3 in order to force the acceptance ratio below 50%.
Convergence is determined by watching the means of the parameter distributions, as well as the average asymptote to a stable value. While convergence is typically seen within the first % of the MCMC, we conservatively discard the first 50% to ensure that our samples are distributed according to the true posterior PDF. We organize the remainder of the models into bins in each parameter, which span between standard deviations from the parameter’s mean. The number of bins is given by , where is the number of models included in the histogram. This is the optimal number of bins spanning this width estimated from biasvariance theory, assuming an inherently normal distribution.^{3}^{3}3The BiasVariance Tradeoff applied to histogram bins is given by , where is the bin width, is the number of data points, and is the first derivative of the PDF. The formula we apply is obtained by assuming , where is the peak of the PDF, and is its standard deviation, estimated in the usual manner. The histogram created this way defines each parameter’s marginalized posterior probability distribution. We define the “best fit” parameters to be the modes of these distributions (i.e., the values with the highest probability) when binned as above. In order to obtain uncertainties in these values, we determine the width in parameter space which encloses 68.3% (1 ) of the models on either side of the mode. For parameters without strict boundary conditions, these distributions tend to be roughly Gaussian.
Object  Model  (as)  ()  (Jy)  (Jy)  

Sgr A*  Gaussian  
Sgr A*  Crescent  
Sgr A*  Crescent + Ring  
M87  Gaussian  
M87  Crescent 
Best fit parameters for the crescent and Gaussian models are listed in Table 3, and figure 3 displays the images corresponding to these parameters after accounting for electron scattering. Our ringless crescent and Gaussian model parameters are consistent with previous work (Broderick et al., 2011b; Kamruddin & Dexter, 2013). For Sgr A*, the inclusion of the photon ring causes only a slight shift in the best fit crescent parameters. The inferred flux in the ring component is poorly constrained, but is marginally nonzero and consistent with the fraction found in theoretical models of Sgr A* (see below). The inclusion of the photon ring is not statistically supported by M87 data (see 4.1).
Current EHT data from both sources are insufficient to distinguish between these different models, and this is further illustrated in figure 2, a visualization of our fits. Here, EHT visibility amplitudes are overplotted on the visibility amplitudes of our best fit models given in 3. The xaxis represents radius in the Fourier plane, while each color from purple to red represents a slice at a different angle, sampled evenly in the interval , where is the position angle of the model. (Apart from the subtle asymmetry of scattering at these values, these models are symmetric about the line , in addition to the usual rotational symmetry in the (u,v) plane.) Data points are color coded according to the section of the surface they are meant to fit.
These fits reveal the dangers of overfitting data poorlysampled in Fourier space, as qualitatively different models can produce similar visibility amplitudes at the uv locations measured by current stations. The total flux of M87 can vary noticeably between models, since our M87 data lack an independent flux constraint. Note further that the M87 bestfit crescent has a much larger size than its bestfit Gaussian, made possible because the long baselines in the crescent fit measure flux from a higherorder peak in the power spectrum. In general, such subtleties make it challenging to find global bestfit parameters, as many completely different models can produce similar fits. These degeneracies will be broken with the inclusion of phase information, which discriminate between adjacent peaks of the power spectrum, and additional telescopes, which improve (u,v) coverage and more precisely measure its shape.
3 Arrays
In the context of VLBI, an “array” corresponds to an instantiation of the the entire EHT, while a “telescope” corresponds to a single element of the EHT, which is sometimes an interferometer that has been phased together to function as a single dish. To represent the EHT at different stages of its development, we devise two sets of increasingly complete arrays. One of these sets operates at 230 GHz, for which there is existing data and a greater number of telescopes. The other set operates at 345 GHz, a frequency which is less affected by the astrophysical noise noted in section 2. The telescopes and their System Equivalent Flux Densities (SEFD) for each of these arrays (Doeleman et al., 2009) are displayed in Table 4. The current 230 GHz array uses a bandwidth of 1024 MHz, while all other arrays use a bandwidth of 4096 MHz, a bandwidth which is soon to be achieved.
230 GHz  345 GHz  

Telescopes  Current (0)  Precision (1)  Coverage (2)  Complete (3)  Basic (0)  Coverage (1)  Future (2) 
JCMT  9400  
SMTO  11900  11900  11900  11900  23100  23100  23100 
CARMA  17700  3500  3500  3500  4900  
SMA  4900  4900  4900  8100  8100  8100  
LMT  10000  560  
APEX  6500  6500  12200  12200  12200  
ASTE  14000  14000  14000  
ALMA  110  1000  1000  
PV  2900  
PdBI  1600  3400  3400  
SPT  7300  
Bandwidth (MHz)  1024  4096  4096  4096  4096  4096  4096 
The Current 230 GHz array consists of only the James Clark Maxwell Telescope (JCMT) in Hawaii, the Submillimeter Telescope Observatory (SMTO) in Arizona, and 2 dishes of the Combined Array for Research in Millimeterwave Astronomy (CARMA) in California. This array, which has been used in EHT campaigns for the last three years, has relatively poor sensitivity and little coverage in the NorthSouth direction. The next array, titled Precision, preserves the same baselines while increasing sensitivity. The Submillimeter Array (SMA) replaces the JCMT, both of which are located in Hawaii, and the number of CARMA dishes is increased to 8. Additionally, the bandwidth on all telescopes is increased from 1024 MHz to 4096 MHz, directly improving sensitivity by a factor of two. In the Coverage array, the Large Millimeter Telescope (LMT) located in Mexico and the Atacama Pathfinder EXperiment (APEX) are added to provide better (u,v) coverage, especially in the NorthSouth direction. APEX is currently being used in EHT campaigns, while efforts are underway to use the LMT at 230 GHz. Finally, the Complete array includes 50 dishes of the Atacama Large Millimeter Array (ALMA) in Chile, the Plateau de Bure Interferometer (PdBI) in France, the 30meter telescope located atop the Pico Veleta (PV) in Spain, and the South Pole Telescope (SPT). Instrumentation is also upgraded at the LMT, greatly increasing its sensitivity. This Phase should be reached in the near future, once the ALMA phasing project (Fish et al., 2013) is complete.
So far, no data have been taken at 345 GHz. Our Basic 345 GHz array consists of SMTO, SMA, APEX and the Atacama Submillimeter Telescope Experiment (ASTE). The Coverage array then includes both 10 ALMA dishes, which provides much better sensitivity but no additional baselines, and PdBI, which provides an additional baseline at the edge of the resolutionlimit due to turbulence. Finally, the Future array includes 8 CARMA dishes, providing previouslylacking coverage on small baselines. Such a short baseline as provided by CARMASMTO is essential to fitting models to 345 GHz data: before its introduction, there are no short baselines to determine the largescale structure of the target. The LMT can be added alternatively/in addition for a similar effect. This is discussed further in section 4.3.
Figure 3 illustrates the (u,v) coverage of the EHT at both targets and frequencies, as well as the image of the underlying bestfit crescent model of each target at 230 GHz. Observing Sgr A* at 230 GHz, very long baselines drawn from the south pole to the Northern Hemisphere are readily apparent, but the scales they probe are largely blurred away by electron scattering. Moving to 345 GHz, one can observe the effects of diminished scattering, as well as the widening of all of the baselines at higher frequency. Until CARMA is added to the 345 GHz array, there are no short baselines to gain information about the largescale structure of either target. The EHT has balanced (u,v) coverage in lineofsight to M87, and there is no evidence for strong interstellar scattering along this line of sight.
Figure 4 provides a visualization of the increased coverage and precision of the completed EHT compared to our current dataset. Here, real or simulated data points are overplotted upon our blurred bestfit crescent + ring model to Sgr A*. The colors from purple to red represent rather than as in figure 2, since the asymmetry from the scattering Gaussian becomes noticeable at longer baselines. For clarity, only errors bars are shown. The dashed and dotteddashed lines represent the scattering Gaussian along its minor and major axes respectively. From this figure, it is clear that the Complete array will be capable of making precise detections at even its longest baselines, in spite of the presence of interstellar scattering. These long baselines will therefore be useful in probing smallscale structure in our targets, which may arise e.g. from turbulence driven by accretion disk instabilities.
4 Simulated Experiments
We perform a series of simulated EHT observations using models from §2 and mmVLBI arrays from §3, in order to i) assess its potential for detecting signatures of strong gravity and test black hole accretion theory; and to ii) determine which technical advances are most important for achieving these goals.
4.1 Detecting the Photon Ring
We do this by simulating observations of a model including a photon ring component, and asking when fitting models without this component is no longer statistically acceptable, e.g., when the fit would improve significantly from the inclusion of a photon ring. We first generate 7 crescent models with photon ring flux contributions that range from 1% to 40% of the total flux, sampled evenly in logarithmic space. At 230 GHz, the total flux and geometric parameters are fixed at the best fit ringless crescent parameters to current data of the target source, either Sgr A* or M87. At 345 GHz, where there have been no observations, these parameters are fixed instead to the best fit ringless crescent parameters from simulation images, listed in table 5. For Sgr A*, we use simulation MBQ, while for M87 we use model DJ1, both of which resemble crescents by eye and were simulated for each of the respective targets (Table 1 and Figure 1).
Object  Model  (as)  ()  (Jy)  (Jy)  

MBQ  Gaussian  
MBQ  Crescent  
DJ1  Gaussian  
DJ1  Crescent 
For each fraction of ring flux, each array, each frequency, and each target we find the best fit ringless crescent using the MCMC methods described in §2 and record the minimum value. To minimize computational time, we use only 200 trials, which we found to reliably converge to the minimum value. To represent the ability to run multiple observing campaigns, we simulate observations and model fitting three times and find the maximum value of from each of the three runs. The maximum represents the case in which a ringless crescent fits the visibilities of a ringed crescent most poorly. These values are used in order to determine the probability of a false negative: the probability of concluding that a ringless crescent is consistent with the ringed crescent data. We do not account for possible false positives, which could be caused by degeneracies between different models.
In the case of the Gaussian distributed amplitude measurements, we can calculate the false negative probability using the CDF of the distribution:
(11) 
where and is the number of observed data points.
Closure phases are not Gaussian distributed, and so the values derived from the usual formula do not follow the chisquared distribution. In fact, the reduced value for the best fit model when measuring amplitudes and closure phases is not always 1, sometimes varying from . This is most evident in the case of few closure triangles and low signaltonoise. Since the true PDF of values can be characteristically wider or thinner than its standard deviation would imply, this is not surprising. Since our analysis hinges on correctly representing the tail of the distribution of , we make a separate estimate of the CDF of values. We simply “bruteforce” sample the distribution by repeatedly calculating values of models fit to themselves. We generate a large sample by simulating observations of the “true” image and calculating when it is fit with its own model. This is computationally intensive, since the distribution is a function of both the instantiation of the EHT array and the model image. The array determines which closure phases are measured, while the model image affects the signaltonoise ratio, an input to equation 8. We calculate 20000 values of for each image and array for which a closure phase analysis is done, giving us the true distribution of . This number of samples allows us to determine up to with some accuracy. Rejection probabilities for test models are then calculated by determining where its value ranks in this probability distribution of ; we compute how likely it would be for the testmodel’s value of to be drawn from the distribution of truth models.
Figure 5 displays the false negative probabilities as a function of array and ring contribution for a variety of experiments. The most pessimistic results originate from Sgr A* at 230 GHz, where in all cases, the photon ring is undetectable even with the most advanced arrays and the greatest ring contributions. However, this particular case suffers three problems: (1) a crescent model can be degenerate with a ring model, since both have Bessel functions as their Fourier transforms; (2) the best fit ringless crescent to Sgr A* happens to subtend the same angular size as the photon ring ( as), placing them on top of one another; and (3) interstellar scattering suppresses flux on long baselines where the Bessel functions differ significantly. By repeating this experiment for Sgr A* without accounting for scattering, we determine that it is indeed the intrinsic geometry of the observed image that creates these pessimistic results (problems 1 and 2), rather than the presence of scattering (problem 3). This is consistent with our result that moving to 345 GHz does little to rectify this problem, although the Future array manages to detect the photon ring when its flux is the brightest. This geometry therefore represents the worstcase scenario.
The most significant difference between a ring and a crescent in Fourier space is a phase asymmetry that corresponds to the offset of the crescent’s hole from the image’s center. In addition, closure phases can be measured with much greater precision than amplitudes. Hence, we see a marked improvement in our results by including closure phase information. The most optimistic scenario for the Sgr A* crescents is one in which amplitudes and closure phases are measured at 345 GHz.
These results are model dependent: The similarity of a crescent to a ring makes ring detection particularly difficult. We include for comparison results where the crescent in the observed image is substituted with the best fit Gaussians to Sgr A* and M87 (or their respective simulations at 345 GHz). Here, even without phase information, all arrays aside from the current 230 GHz array are capable of detecting a ring with some amount of flux. In fact many arrays fall off the plot, immediately distinguishing even a 1% ring flux beyond 4 significance. This is because Gaussian models provide little power on long baselines. In this setup, all of the flux on those baselines must originate from the photon ring.
Results are also more optimistic if the source is changed to M87 (Figure 5). With a significantly larger best fit crescent size than photon ring size, model degeneracies are less severe for the M87 crescent. This, combined with the absence of astrophysical noise, allows for ring detection at high significance even with amplitudes alone. Advancing to closure phases, the arrays reach our limit quite easily in almost all cases. However, M87 Gaussian results happen to be more pessimistic than those for Sgr A*, due to the fact that the Gaussian’s smaller intrinsic size does provide power at some of the longer baselines.
If the models used represent the true images of M87 and Sgr A*, a photon ring is more likely to be detected either by observing at 345 GHz or by observing M87. Yet in practice, both the model parameters and general shapes are not certain enough to warrant this conclusion. Using both sources and both observing frequencies provides the best opportunity for the EHT to detect a photon ring given possible degeneracies. In addition, a timevariability test, such as that developed in section 4.4, could be more effective in this regard, since the photon ring stays fixed in radius while the accretion flow can undergo significant structural variations.
4.2 Differentiating Single Accretion Flow Images
There are a variety of dynamical models that all provide good fits to existing mmVLBI Sgr A* data at 230 GHz. As the EHT develops, we will be able to distinguish these models. We illustrate this using some of the best fit simulation images from §2. We conduct an experiment to determine whether the EHT at various stages can distinguish between the best fit frames of timedependent images calculated from different magnetohydrodynamical models of Sgr A* (Table 1).
For each of these images, we calculate its Fourier Transform and interpolate onto the coordinates observed by a given array to obtain observed visibilities, with noise and errors applied as in §2. The fluxes of the M87 images are normalized to the values of our bestfit crescent model listed in table 3, but we do not perform this step for Sgr A* bestfit images. To obtain , we then interpolate each model’s Fourier Transform onto the model visibilities of the truth image. Then, we calculate the number of with which each observation can rule out the fitted model. For amplitude measurements, we compare to the distribution and the builtin IDL function gauss_cvf
to obtain the number of with which we can determine that these models are dissimilar. For amplitude + closure phase measurements, we compute a sample of 20000 measurements and compare to this distribution rather than the canonical distribution as in §2.
Table 6 displays the arrays which are capable of distinguishing two models by or greater given different modes of observation. We find that all models can be distinguished from one another with 4 significance prior to the inclusion of ALMA into the EHT. For Sgr A*, the extra sensitivity provided by the Precision array is sufficient to make two of the distinctions, while the most difficult distinction is accomplished with the additional telescopes of the Coverage array. For M87, the disk and jet models (DJ1 and J2 respectively) are distinguishable even with the Current array, implying that the current data may already be capable of constraining properties of its accretion flow. Finally, the inclusion of closure phases allows even the Current array to make any distinctions.
Sgr A*  M87  

Measurements  90h vs. MBD  MBD vs. 915h  MBD vs. MBQ  MBD vs. 915h (345 GHz)  MBD vs. MBQ (345 GHz)  DJ1 vs. J2 
Amplitudes  Coverage  Precision  Precision  Basic  Basic  Current 
Closure Phases  Current  Current  Current  Basic  Basic  Current 
These results imply that additional data with the current array may already be able to provide new constraints on the accretion flows of Sgr A* and M87. In practice, however, the inherent time variability of these models combined with parameter degeneracies make these comparisons less straightforward.
4.3 Constraints on Crescent Parameters
The previous subsection showed that the EHT can already begin to distinguish between viable theoretical models. In order to quantitatively measure the ability of the different instantiations of the EHT to constrain accretion models, we perform an experiment to determine how well each of the arrays could constrain crescent model parameters. We select the crescent model with the best fit parameters to Sgr A* as the observed image. To isolate array performance, we use these same parameters at 345 GHz as well as 230 GHz. As in §2.4, we fit crescent models using MCMC, where the number of walkers increased to 300 in order to ensure that the posterior distribution was wellsampled. We then calculate the uncertainties on these parameters by enclosing the 99.7% of the models () on either side of the mode of the posterior distribution, and averaging the uncertainty of either side.
Figure 6 displays the results of this test. For the parameter , uncertainties were normalized by dividing by the mean value of each of those parameters. Uncertainties in are in units of degrees. Solid black lines represent runs where measurements are made using amplitudes only, while dashed blacklines represent runs involving both amplitudes and closure phases. For the 230 GHz plots, the red solid and dashed lines correspond to the RMS variations of the best fit parameters to simulations 515h and 90h, one of the most and least varying simulations respectively. Fluctuations between different frames of 515h are just at the threshold necessary to detect timevariability with the Current array, and this is later developed in section 4.4. The blue solid and dashed lines represent the differences in bestfit crescent model parameters for MBD vs. MBQ and MBD vs. 915h, two of the most similar and least similar models respectively. Interestingly, while the differences between MBD and MBQ are difficult to determine at 230 GHz, they become less similar at 345 GHz.
The progression from amplitude to closure phase measurements is more pronounced at 345 GHz than at 230 GHz, due to increased signaltonoise on long baselines in the absence of electron scattering. The plateau in closure phase constraints in the last two 230 GHz arrays, most evident in and , suggest that a limit is reached due to electron scattering.
4.3.1 Possible Large Scale Problems at 345 GHz
A subtle issue that we encounter is that the 345 GHz arrays poorly constrain largescale (as) structure without the inclusion of short baselines. This is not evident in the previous test (and those to come), since crescent and ring models both have power at high spatial frequencies, and can therefore be constrained using only large baselines. Figure 7 demonstrates this shortcoming. Here, we simulate observations of a symmetric Gaussian 50 as in radius, with a flux of 1.7 Jy, a reasonable flux taken from simulation MBQ. We compare data points obtained using the Coverage array with data points obtained using the Future array, with the addition of the LMT, another telescope capable of providing shorter baselines. Data points from baselines which include CARMA or the LMT are plotted in red, while those from baselines without either are plotted in black. To calculate error bars, the SEFD of the LMT is taken to be 13700, as in Lu et al. (2014). Model visibilities of Gaussians with the true radius, half the radius, and twice the radius are overplotted. It is only with the inclusion of CARMA that the true size can be constrained–all other baselines either measure the total flux, or resolve out the source.
This problem should be exacerbated on even larger scales, and it is not uncommon for simulated accretion flows to subtend areas 23 times larger. Hence, caution should be used interpreting future data obtained with the EHT’s 345 GHz arrays if shorter baselines are not acquired. It should be noted that even our Future array does not constrain largescale structure in the NorthSouth direction. Obtaining 345 GHz capabilities at additional telescopes such as the LMT help could rectify this problem. However, even without short baselines, 345 GHz data may still be quite valuable for testing models in conjunction with 230 GHz observations. For example, current viable models make significantly different predictions for even the total flux at 345 GHz (e.g., Dexter et al., 2010).
4.4 Structural Variability On Event Horizon Scales
Both Sgr A* and M87 exhibit significant () variability on relevant timescales for EHT observations: hr for Sgr A*, yr for M87, and the simulated images we use are also timevariable at consistent amplitudes. Fish et al. (2011) found no significant evidence of structural variability in Sgr A* data even as the total flux changed by between observations, but it should be readily detected with future data. Here we assess the sensitivity of the EHT to structural variability as predicted by the models. We answer this question using both timedependent accretion flow and jet images calculated from MHD simulations, and geometric crescents with rings, which make fewer assumptions about the physical environment around the black hole. For each test, we assume three different days of observation. The model crescent or simulation image is static during each individual day, but no two days’ images are alike. This assumption is justified for M87, where the lightcrossing time of the event horizon is day. Sgr A* varies on hour timescales (e.g., Dexter et al., 2014). Our simulated days of observation could also represent separate segments within a single day for the case of Sgr A*.
4.4.1 Crescent Models
Geometric crescent models are computationally cheap to generate, and their variability is easier to quantify. We consider variability between days in each model parameter except total flux of 1%, 2%, and 5% through 25% with 5% steps, levels consistent with predictions of MHD simulations (see Table 8). We begin with the best fit crescent parameters to Sgr A* or M87 given by table 3 at 230 GHz, or the best fit crescents to simulations MBQ and DJ1 given by table 5. Then, we multiply each geometric parameter by for each day, where is a random variable drawn from a Gaussian distribution with a mean of 0 and a standard deviation equal to the value of the variability in question. Since variations in the total flux of the image are trivial to deduce from data, we keep the fluxes of the crescent fixed, so that we only explore variations in the spatial parameters of the crescent model.
We then use MCMC to fit to a combined dataset of three days of observation in two different ways: First, a static fit in which on model must fit all three days of observation, and secondly a timevarying fit in which parameters are allowed to vary between different days. In order to determine whether a static or timevarying fit is statistically favored, we record the best value of in each case and determine whether it is consistent with the expected distribution of values, given the number of data points and free parameters. When we include closure phase measurements, nonGaussian statistics impel us to “bruteforce” sample the distribution directly. As in §4.1, we calculate 20000 samples of the distribution of values with which to compare our data.
Source  Frequency  Measurements  1%  2%  5%  10%  15%  20%  25% 

Sgr A*  230 GHz  Amplitudes  Complete  Coverage  Precision  Current  Precision  Precision  Current 
Sgr A*  230 GHz  Closure Phases  Precision  Current  Current  Current  Current  Current  Current 
Sgr A*  345 GHz  Amplitudes  Never  Coverage  Coverage  Basic  Basic  Basic  Basic 
Sgr A*  345 GHz  Closure Phases  Basic  Basic  Basic  Basic  Basic  Basic  Basic 
M87  230 GHz  Amplitudes  Complete  Precision  Coverage  Precision  Current  Current  Current 
M87  230 GHz  Closure Phases  Current  Current  Current  Current  Current  Current  Current 
M87  345 GHz  Amplitudes  Future  Basic  Coverage  Basic  Basic  Basic  Basic 
M87  345 GHz  Closure Phases  Basic  Basic  Basic  Basic  Basic  Basic  Basic 
Table 7 shows the results of this experiment, displaying the first array which detects time variability at the level for different levels of fluctuation. We notice no preferred source for this case. By performing a test in which we repeat this experiment for M87 with its predicted flux doubled (to a 230 GHz flux more similar to Sgr A*), we determine that these results are not limited by the signaltonoise ratio. After all, we expect these largescale fluctuations to be detectable by the short baselines with the highest signaltonoise. Neither do we see greater success rates at 345 GHz versus 230 GHz. One would hope to perform this test at both frequencies, however. First, the sources may have different fluxes at each frequency; neither source has been resolved at 345 GHz, so this remains an open question. Second, observations at 345 GHz plumb a deeper optical, and therefore physical, depth and may therefore undergo larger amplitude structural variations. The factor which has the greatest effect on detecting timevariability is the inclusion of closure phase information. As discussed in §2.3.2, longer integration times will allow more precise measurements of closure phases than visibility amplitudes.
For comparison with physical models, we fit crescent models to each frame of various simulation images, and record the best fit parameters. A mesh in (u,v) space is generated, linearly sampling with 30 steps in each parameter from 0 to 5 G, and from 0 to . Physical models and crescent models are then fit using MCMC with visibility amplitudes obtained at these data points. Data points are obtained by this method rather than by simulating observations with a real array because (i) this ensures an even sampling of (u,v) space, with equallyweighted points and (ii) crescent models are not expected to fit data points where , where power is dominated by small scale structure. Given these fits, we compute the standard deviation in each parameter and divide by the average to obtain fractional time variability for each model. These are listed in table 8. Since most of the model’s crescents vary much more than the threshold needed to detect time variability for the Current array, we predict that geometric time variability should be detectable in EHT data once coherently averaged closure phases become available. In the following section, we perform a more rigorous experiment to lead to this conclusion.
Model  (as)  

515h  16%  27%  58%  14% 
90h  4%  6%  10%  6% 
MBQ  8%  15%  19%  10% 
DJ1  16%  19%  22%  27% 
J2  17%  15%  19%  10% 
4.4.2 Simulated Models
We perform a similar test to compare our results with timevarying crescents to physical simulated images. For a given simulation, we first calculate the obtained by fitting each of its images to each of its other images when observed with each of the different arrays. Variations in flux alone from frame to frame are enough to make this experiment trivial for even the most basic arrays. Hence, to isolate geometric fluctuations, we normalize the flux in each frame to our current bestfit fluxes of each target.
We determine the fraction of the time that each of the arrays can detect timevariability by simulating 20000 observing campaigns, each of which contains observations for three different epochs. For each of these epochs, we assume the source remains fixed at a random frame from its simulation. The bestfit timevarying model is simply given the three different frames themselves, while the bestfit static model is the single frame which minimizes the total when fit to all three epochs simultaneously. Saving each of these values for each campaign, we obtain the distribution of that we expect for both timevarying and static fits.
Given these distributions, we then calculate the percentage of the time that the static fits are more than inconsistent with the timevarying fits. In table 9, we display this fraction for each image, each array, and each mode of observation. We find that it is only the current 230 GHz array that has an appreciable difficulty in detecting timevariability. The Precision array, which has the same baselines as the Current array but a more sensitive Hawaii telescope, additional CARMA dishes, and higher bandwidth, performs dramatically better than the Current array. This suggests that precision, rather than (u,v) coverage, is the limiting factor in determining timevariability, as it occurs at all spatial scales plumbed by the EHT. Once again, the inclusion of closure phases also boosts the chances of distinguishing models significantly.
230 GHz  345 GHz  
Measurements  Simulation  Current  Precision  Coverage  Complete  Basic  Coverage  Future 
Amplitudes 
515h  72%  100%  100%  100%  
90h  6%  94%  100%  100%  
MBQ  49%  97%  100%  100%  
DJ1  3%  83%  93%  100%  84%  100%  100%  
J2  44%  99%  100%  100%  42%  100%  100%  
Closure Phases 
515h  100%  100%  100%  100%  
90h  98%  100%  100%  100%  
MBQ  99%  100%  100%  100%  
DJ1  90%  100%  100%  100%  100%  100%  100%  
J2  100%  100%  100%  100%  100%  100%  100%  

5 Discussion
This work presents a concrete demonstration of the potential of the EHT to do transformative science within the framework of plausible observing parameters and models for the source structure in Sgr A* and M87. However, there are a number of limitations to the methods used here. The arrays we consider (Table 4) correspond roughly to the technical roadmap of the EHT experiment, but are necessarily approximations of the true evolution of its technological capabilities in time. Future array developments which we have not considered in our analysis include the addition of the 12m Greenland Telescope (GLT, Inoue et al., 2014), which will provide long baselines for M87 observations, and the addition of LMT functionality at 345 GHz, which will help constrain largescale structure at this frequency.
Simulating future observations requires assuming models for each source. We have made plausible choices based on best fit images from state of the art simulations, and a modelindependent geometric crescent approximation aiming to capture the important relativistic effects. Both choices are first steps. The former models are subject to systematic uncertainties in the microphysics and initial conditions of the simulations, which should be explored in the future. The geometric crescent model is physically motivated and has more parameters (5) than asymmetric Gaussian (4) or ring (3) shapes, but is still too simplistic to reproduce the theoretical images in any detail. In addition, while we consider an azimuthally symmetric photon ring contribution for simplicity this feature is also asymmetric in theoretical models due to Doppler beaming. Analytic extensions including additional parameters (Benkevitch et al., 2013) are a promising means to extend these models to quantitatively fit future EHT data sets. We have also separately carried out static and timevariable tests. The most effective method for detecting a photon ring may depend on its stationarity relative to the variable emitting region.
While beyond the scope of this paper, the EHT array will also be able to measure spatially resolved dual polarizations, which will provide independent constraints on the structure of magnetic fields and geometry of the emitting region. Advantages of polarimetric interferometry include the immunity of fractional polarization to scatterbroadening, thereby circumventing electron turbulence in line of sight to Sgr A*, and the ability to reach greater astrometric precision than the size scales probed by an array’s baselines. Recent work by Johnson et al. (2014) demonstrates the promise of polarimetric interferometry for achieving precise relative astrometry of a variable polarized image compact.
Similar to previous work (Broderick et al., 2011a), we find that coherently averaged closure phases are a powerful observable for EHT science. Closure phase can now be measured, and increasing their sensitivity is an important future goal for the observations. Even with current capabilities, we predict that structural variability should be detectable with additional data. The amplitude and type of variability observed will provide a powerful constraint on accretion models, possibly including the direct detection of magnetic turbulence in a black hole accretion flow.
Here we only consider the detectability of a Schwarzschild photon ring. In principle, the characterization of the photon ring, particularly its shape and phase offset from the rest of the image, could allow a measurement of black hole spin (Takahashi, 2004) and a test of the nohair theorem (Johannsen & Psaltis, 2010). This can be explored with simulated EHT data in future work, either using static models as we have done here or by match filtering the expected feature in timedependent simulated data.
The models we consider assume general relativity. If the Kerr geometry is modified, it becomes difficult to disentangle signatures of modified gravity from different accretion flow geometries. A computationallyexpensive, full parameter space search is required, as recently explored by Broderick et al. (2014).
Despite these caveats, the techniques that we develop in this paper do not depend on the underlying models. As new simulation models or analytic parameterizations are introduced, they can be used while applying the same methodology.
6 Conclusion
With both geometric crescents and hydrodynamical accretion models, we have devised and simulated a variety of tests that can be used on current and future EHT data and assessed the array’s capabilities during its stages of development. The main results of this paper are summarized as follows.

The use of closure phase information greatly increases the constraining power of the EHT. Due to the ability to coherently average closure phase information over time periods much longer than atmospheric coherence time, the increase in precision gained through measuring closure phases more than makes up for the lack of full phase information.

The EHT should be able to discriminate between many bestfit theoretical models from numerical simulations even with additional data from the current array, provided that closure phases can be robustly measured with high precision (Section 4.2). Model constraints will become far stricter still in the next few years as bandwidth is increased and stations are added to the array (e.g., LMT and ALMA). EHT observations can therefore provide a stringent test of contemporary accretion theory and a precise measurement of black hole and accretion flow parameters.

EHT closure phases are also currently sensitive to structural variability at the level present in the models (Section 4.4). Current data may therefore allow the first detection of structural variability on event horizon scales around a black hole. Timedomain EHT measurements provide could allow the direct study of instabilities in accretion flows both on the scales of the entire image and its fine structure, e.g. as may result from MHD turbulence driven by the MRI (right panel of Figure 4).

The assumed underlying image model can significantly impact the results achievable by the EHT. Using the full complement of two sources (Sgr A* and M87) and frequencies (230 and 345 GHz) can mitigate possible degeneracies between the accretion flow or jet image and the photon ring (e.g., the difference between the results for 230 GHz M87 or 345 GHz Sgr A* and 230 GHz Sgr A* in Figure 5). Using M87 and 345 GHz are particularly important to avoid the strong interstellar scattering towards Sgr A*, which can still limit the science capabilities of observations at 230 GHz.
Frequency  Array  Science Goals within Reach 
230 GHz  Current  With amplitudes only, detect expected structural time variability across observation epochs for some physical 
models and distinguish jet from disk models. With closure phases, detect structural time variability for all physical  
models, and break degeneracies between current bestfit models to Sgr A*.  
Precision  With amplitudes only, detect structural time variability across observation epochs for all physical models,  
break degeneracies between most current bestfit models to Sgr A*, and distinguish M87 photon ring from crescent.  
Coverage  With amplitudes only, break all degeneracies between current bestfit models to Sgr A* and probe structural  
timevariability at 5% levels.  
Complete  Detect photon ring in Sgr A* and M87 for all levels of predicted flux for a wide range of models.  
With amplitudes alone, probe structural timevariability at % levels.  
345 GHz  Basic  With amplitudes alone, place constraints on timevariability of physical models. With closure phases, 
detect this timevariability with certainty.  
Coverage  With closure phases, distinguish Sgr A* photon ring from crescent.  
Future  With amplitudes only, distinguish Sgr A* photon ring from crescent. Constrain large scale accretion flow. 
acknowledgements
We thank S. Doeleman and E. Quataert for useful discussions. We are grateful to L. Benkevitch, S. Doeleman, V. Fish, and R. Lu for sharing the MAPS software and for assistance with its use.
References
 Backer (1978) Backer D. C., 1978, ApJ, 222, L9
 Bardeen, Carter & Hawking (1973) Bardeen J. M., Carter B., Hawking S. W., 1973, Communications in Mathematical Physics, 31, 161
 Benkevitch et al. (2013) Benkevitch L. et al., 2013, in American Astronomical Society Meeting Abstracts, Vol. 221, American Astronomical Society Meeting Abstracts #221, p. #143.10
 Bower (2006) Bower G. C., 2006, Journal of Physics Conference Series, 54, 370
 Bower et al. (2014) Bower G. C. et al., 2014, ApJ, 780, L2
 Broderick et al. (2011a) Broderick A. E., Fish V. L., Doeleman S. S., Loeb A., 2011a, ApJ, 738, 38
 Broderick et al. (2011b) Broderick A. E., Fish V. L., Doeleman S. S., Loeb A., 2011b, ApJ, 735, 110
 Broderick et al. (2014) Broderick A. E., Johannsen T., Loeb A., Psaltis D., 2014, ApJ, 784, 7
 Broderick & Loeb (2009a) Broderick A. E., Loeb A., 2009a, ApJ, 697, 1164
 Broderick & Loeb (2009b) Broderick A. E., Loeb A., 2009b, ApJ, 697, 1164
 Bromley, Melia & Liu (2001) Bromley B. C., Melia F., Liu S., 2001, ApJ, 555, L83
 Dexter, Agol & Fragile (2009) Dexter J., Agol E., Fragile P. C., 2009, ApJ, 703, L142
 Dexter et al. (2010) Dexter J., Agol E., Fragile P. C., McKinney J. C., 2010, ApJ, 717, 1092
 Dexter et al. (2012) Dexter J., Agol E., Fragile P. C., McKinney J. C., 2012, Journal of Physics Conference Series, 372, 012023
 Dexter & Fragile (2013) Dexter J., Fragile P. C., 2013, MNRAS, 432, 2252
 Dexter et al. (2014) Dexter J., Kelly B., Bower G. C., Marrone D. P., Stone J., Plambeck R., 2014, MNRAS, 442, 2797
 Dexter, McKinney & Agol (2012) Dexter J., McKinney J. C., Agol E., 2012, MNRAS, 421, 1517
 Doeleman et al. (2009) 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. et al., 2012, Science
 Doeleman et al. (2008) Doeleman S. S. et al., 2008, Nature, 455, 78
 Falcke et al. (1998) Falcke H., Goss W. M., Matsuo H., Teuben P., Zhao J.H., Zylka R., 1998, ApJ, 499, 731
 Falcke & Markoff (2000) Falcke H., Markoff S., 2000, A&A, 362, 113
 Falcke, Melia & Agol (2000) Falcke H., Melia F., Agol E., 2000, ApJ, 528, L13
 Fish et al. (2013) Fish V. et al., 2013, arXiv:1309.3519
 Fish et al. (2011) Fish V. L. et al., 2011, ApJ, 727, L36
 Fragile et al. (2007) Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007, ApJ, 668, 417
 Fragile & Meier (2009) Fragile P. C., Meier D. L., 2009, ApJ, 693, 771
 Gebhardt et al. (2011) Gebhardt K., Adams J., Richstone D., Lauer T. R., Faber S. M., Gültekin K., Murphy J., Tremaine S., 2011, ApJ, 729, 119
 Gillessen et al. (2009) Gillessen S., Eisenhauer F., Trippe S., Alexander T., Genzel R., Martins F., Ott T., 2009, The Astrophysical Journal, 692, 1075
 Goodman & Weare (2010) Goodman J., Weare J., 2010, Comm. App. Math. Comp. Sci., 5
 Hada et al. (2011) Hada K., Doi A., Kino M., Nagai H., Hagiwara Y., Kawaguchi N., 2011, Nature, 477, 185
 Inoue et al. (2014) Inoue M. et al., 2014, arXiv:1407.2450
 Johannsen & Psaltis (2010) Johannsen T., Psaltis D., 2010, ApJ, 718, 446
 Johannsen & Psaltis (2011) Johannsen T., Psaltis D., 2011, Advances in Space Research, 47, 528
 Johnson et al. (2014) Johnson M. D., Fish V. L., Doeleman S. S., Broderick A. E., Wardle J. F. C., Marrone D. P., 2014, arXiv:1408.6241
 Kamruddin & Dexter (2013) Kamruddin A. B., Dexter J., 2013, MNRAS
 Lambert & Gontier (2009) Lambert S. B., Gontier A.M., 2009, A&A, 493, 317
 Lu et al. (2014) Lu R.S., Broderick A. E., Baron F., Monnier J. D., Fish V. L., Doeleman S. S., Pankratius V., 2014, arXiv:1404.7095
 McKinney & Blandford (2009) McKinney J. C., Blandford R. D., 2009, MNRAS, 394, L126
 Melia (1994) Melia F., 1994, ApJ, 426, 577
 Mościbrodzka et al. (2009) Mościbrodzka M., Gammie C. F., Dolence J. C., Shiokawa H., Leung P. K., 2009, ApJ, 706, 497
 Narayan, Yi & Mahadevan (1995) Narayan R., Yi I., Mahadevan R., 1995, Nature, 374, 623
 Noble et al. (2007) Noble S. C., Leung P. K., Gammie C. F., Book L. G., 2007, Class. and Quant. Gravity, 24, 259
 Petrov et al. (2011) Petrov L., Kovalev Y. Y., Fomalont E. B., Gordon D., 2011, AJ, 142, 35
 Reynolds et al. (1996) Reynolds C. S., Di Matteo T., Fabian A. C., Hwang U., Canizares C. R., 1996, MNRAS, 283, L111
 Reynolds & McKee (1980) Reynolds S. P., McKee C. F., 1980, ApJ, 239, 893
 Rogers (1976) Rogers A. E. E., 1976, Methods of Experimental Physics, 12, 139
 Rogers et al. (1984) Rogers A. E. E., Moffet A. T., Backer D. C., Moran J. M., 1984, Radio Science, 19, 1552
 Shcherbakov, Penna & McKinney (2012) Shcherbakov R. V., Penna R. F., McKinney J. C., 2012, ApJ, 755, 133
 Takahashi (2004) Takahashi R., 2004, ApJ, 611, 996
 Thompson, Moran & Swenson (2001) Thompson A. R., Moran J. M., Swenson, Jr. G. W., 2001, Interferometry and Synthesis in Radio Astronomy, 2nd Edition
 Walsh et al. (2013) Walsh J. L., Barth A. J., Ho L. C., Sarzi M., 2013, ApJ, 770, 86
 Yuan, Quataert & Narayan (2003) Yuan F., Quataert E., Narayan R., 2003, ApJ, 598, 301
 Yuan et al. (2009) Yuan Y.F., Cao X., Huang L., Shen Z.Q., 2009, ApJ, 699, 722