ABSTRACT
The recent very long baseline interferometry (VLBI) observation of the Galactic center black hole candidate Sgr A* at 1.3 mm shows source structure on event-horizon scales. This detection enables a direct comparison of the emission region with models of the accretion flow onto the black hole. We present the first results from time-dependent radiative transfer of general relativistic MHD simulation data, and compare simulated synchrotron images at black hole spin a = 0.9 with the VLBI measurements. After tuning the accretion rate to match the millimeter flux, we find excellent agreement between predicted and observed visibilities, even when viewed face-on (i ≲ 30°). VLBI measurements on 2000–3000 km baselines should constrain the inclination. The data constrain the accretion rate to be (1.0–2.3)×10−9 M☉ yr−1 with 99% confidence, consistent with but independent of prior estimates derived from spectroscopic and polarimetric measurements. Finally, we compute light curves, which show that magnetic turbulence can directly produce flaring events with 0.5 hr rise times, 2–3.5 hr durations, and 40%–50% flux modulation, in agreement with observations of Sgr A* at millimeter wavelengths.
1. INTRODUCTION
Due to its large angular size, the Galactic supermassive black hole candidate Sgr A* is a promising laboratory for precision black hole astrophysics using very long baseline interferometry (VLBI) at millimeter wavelengths. Previous measurements (Shen et al. 2005; Doeleman et al. 2001; Krichbaum et al. 1998) at 7 mm and 3.5 mm found a small intrinsic source size for Sgr A* (100–200 μas FWHM), but at those wavelengths interstellar scattering is the dominant contribution to the observed size. Recent measurements by Doeleman et al. (2008) at 1.3 mm are the first at short enough wavelengths to avoid contamination due to interstellar scattering and at long enough baselines to achieve event horizon scale resolution. These measurements have been used by Broderick et al. (2009) to constrain the black hole and accretion flow parameters in terms of a 1.5D, nonrelativistic, time-steady radiatively inefficient accretion flow model (RIAF; Yuan et al. 2003).
Sgr A* is also known to exhibit multiwavelength, polarized flares (e.g., Yusef-Zadeh et al. 2007; Marrone et al. 2008; Eckart et al. 2008a; Dodds-Eden et al. 2009). This activity has been modeled as hotspots orbiting in the inner radii (Broderick & Loeb 2005); adiabatically expanding blobs (Yusef-Zadeh et al. 2006; Eckart et al. 2008b); and emission from a short, mildly relativistic jet (Falcke et al. 1993).
As a complementary approach to the relatively simple models outlined above, we perform relativistic radiative transfer on data from a recent three-dimensional, general relativistic MHD (GRMHD) simulation of a black hole accretion disk (Fragile et al. 2007), and compare simulated synchrotron images to the VLBI observations. In addition, we produce light curves at 1.3 mm as a first step toward explaining Sgr A*'s flaring activity in terms of magnetic turbulence.
In contrast to the results from RIAF modeling, we find that for the single black hole spin considered (a = 0.9), all inclinations provide excellent fits to the VLBI data. Visibility profiles from low (face-on) inclinations (i ≲ 30°) resemble “ring” models (e.g., Doeleman et al. 2008), whereas those at higher (edge-on) inclinations are more Gaussian. Ongoing work with increased array sensitivity should be able to distinguish between these two types of visibility profiles.
Detection of a ring-like profile would provide an important constraint not only on the geometry of Sgr A*, but also on the necessary complexity of accurate accretion flow models. These results also demonstrate the power of millimeter VLBI to distinguish between accretion models as well as the model dependence of attempts to extract black hole parameters.
2. METHODS
2.1. Simulation Data
Fragile et al. (2007) presented global GRMHD simulations of black hole accretion disks. They used a spherical grid with polar Kerr–Schild coordinates and the main untilted simulation had an effective resolution of 1283 except near the poles which are underresolved. The torus was initialized with the analytic, time-steady, axisymmetric solution of De Villiers & Hawley (2003) and threaded with a weak poloidal magnetic field with minimum Pgas/Pmag = 10 initially. The magnetorotational instability (MRI) arose naturally from the initial conditions and the disk quickly became fully turbulent. The simulation was evolved for ∼8000 M, or ∼40 orbits at r = 10 M in units with G = c = 1. Only data from the final 2/3 of the simulation are used in this Letter once the disk is fully turbulent and transient effects from the initial conditions have died down.
This simulation evolved an internal energy equation, and injected entropy at shocks. Such a formulation does not conserve energy, and produces a more slender, cooler torus than conservative simulations which capture the heat from numerical dissipation of magnetic fields (Fragile & Meier 2009).
2.2. Ray Tracing
We performed relativistic radiative transfer on the simulation data via ray tracing. Starting from an observer's camera, rays are traced backward in time assuming they are null geodesics (geometric optics approximation), using the public code described in Dexter & Agol (2009). In the region where rays intersect the accretion flow, the radiative transfer equation is solved along the geodesic (Broderick 2006) in the form given in Fuerst & Wu (2004), which then represents a pixel of the image. This procedure is repeated for many rays to produce an image, and at many time steps of the simulation to produce time-dependent images (movies). Light curves are computed by integrating over the individual images.
To calculate fluid properties at each point on a ray, the spacetime coordinates of the geodesic are transformed from Boyer–Lindquist to the Kerr–Schild coordinates used in the simulation. Since the accretion flow is dynamic, light travel time delays along the geodesic are taken into account. Data from the 16 nearest zone centers (eight on the simulation grid over two time steps) were interpolated to each point on the geodesic. Between levels of resolution near the poles, data from the higher resolution layer were averaged to create synthetic lower resolution points, which were then interpolated. Little emission originates in the underresolved regions of the simulation.
In this work, we consider synchrotron emission from Sgr A* using the formulae from Mahadevan et al. (1996) and neglecting polarization. The simulation provided mass density, pressure, velocity, and magnetic field in code units. These were converted into cgs units following the procedure described in Schnittman et al. (2006), except that we scaled the torus mass to match the observed flux of Sgr A* at 1.3 mm. The length and timescales were set by the black hole mass, taken to be 4 × 106M☉, at a distance of 8 kpc.
The temperature was taken to be thermal, with the ion and electron temperatures equal. Many authors have used two-temperature models of Sgr A*, with weak coupling between the dynamically important ions and the cooler electrons. However, our single-temperature model fits the observed flux from Sgr A* at reasonable accretion rates ∼10−9M☉ yr−1 and produces a spectrum near millimeter wavelengths that is consistent with observations. In this model, the gas is optically thin except at high inclinations during flares. The gas temperature is ∼1011 K in the regions of peak synchrotron emissivity, which is consistent with peak electron temperatures in RIAF models. We also followed Goldston et al. (2005) and produced images with (Te + Ti)/Te = 10, keeping the total gas pressure fixed. These models became optically thick before reaching the observed 1.3 mm flux, and were inconsistent with the location of the peak in the quiescent spectrum of Sgr A*.
3. RESULTS
3.1. Visibility Fitting
Our fitting procedure closely mirrors that described in Broderick et al. (2009). Simulated images were produced for a range of inclination angles, time-averaged accretion rates, and observer times (sets of simulation time steps). The images were then averaged over 10 minute observation periods (eight simulation time steps) to match the integration times used in Doeleman et al. (2008), padded sufficiently to resolve the shortest baselines, and Fourier transformed to two-dimensional visibilities. The visibilities were rotated to a range of different orientations, ξ, the projection of the black hole spin axis on the sky measured E of N. Interstellar scattering was then applied by multiplication by an elliptical Gaussian as described in Fish et al. (2009) using the fits from Bower et al. (2006). Sample images and scatter-broadened visibilities are shown in Figure 1, and animations are available in the online version of the journal. At low inclinations, the emission region is essentially a ring around the shadow, which causes a minimum in the two-dimensional visibilities. For higher inclinations ≳30°, Doppler beaming causes the emission to be concentrated on one side of the black hole. The visibilities are then approximately elliptical, with the minor (major) axis corresponding to the major (minor) axis of the emission region.
Figure 1. Best-fit ray-traced images and corresponding two-dimensional, scatter-broadened, and rotated visibilities for i = 10° (left) and i = 70° (right). The intensities are scaled linearly to the maximum of each panel (from blue to red to yellow), and the UV plane locations of the VLBI observations from Doeleman et al. (2008) are overplotted on the visibilities as triangles. The ring of peak brightness at a radius of ∼25 μas in the i = 10° image is due to the transition between bound and unbound photons. Despite the complex structure of the i = 70° image, the visibility is similar in shape to that in Broderick et al. (2009), Figure 2.(Animations of this figure are available in the online journal.)
Download figure:
Standard image High-resolution imageVisibilities were interpolated to the detections in Doeleman et al. (2008), and an effective χ2 was computed as defined in Broderick et al. (2009). Grids of minimum effective, reduced χ2 are shown in Figure 2 over inclination angle and sky orientation. At inclinations i ⩾ 30°, our results are in at least qualitative agreement with Broderick et al. (2009). The orientation is well constrained at high inclinations, since the general image shape is nearly static and asymmetric, and the long baseline VLBI measurements pick out a specific orientation for the visibility ellipse. The two distinct bands of good fits are due to the (approximate) up–down symmetry of the image. However, RIAF models ruled out inclinations i ≲ 30°, whereas for our GRMHD models, good fits (reduced χ2 ≲ 1.2) are possible at all inclinations. This is especially evident from the curve of p(i) versus i overplotted in Figure 2. Inclinations i ≲ 20° are less probable due to the sin i prior and large variation of probability density with observer time. Although the best-fit χ2 is roughly the same for low inclinations at all ξ, at any given observer time good fits are restricted to less than half the range of sky orientations.
Figure 2. Grid of best-fit reduced effective χ2 values vs. inclination and sky orientation ξ. The scale is from χ2 = 1 (white) to 2.5 and greater (black). We find excellent fits at low inclinations, which are roughly independent of sky orientation. At high inclinations, our results show a similar shape to that found by Broderick et al. (2009). Overplotted is the probability density vs. inclination, marginalized over observer time, accretion rate, and sky orientation.
Download figure:
Standard image High-resolution imageVisibility profiles from the best fits at many inclinations are shown in Figure 3. The profiles represent the two-dimensional visibility plotted along the line representing the average angle in the UV plane of the long baseline measurements from Doeleman et al. (2008). Also plotted are the data from that work. At low inclinations, the profile reaches a minimum as described above, whereas at i ⩾ 30°, the profile decreases monotonically with baseline length. Our best-fit face-on visibility profile is almost identical to that of the ring model from Doeleman et al. (2008).
Figure 3. Best-fit visibility profiles for low and high inclinations, plotted along the line in the UV plane corresponding to the average location of the long baseline measurements from Doeleman et al. (2008). The visibilities naturally divide into two types. At low inclinations, the nearly circular shadow leads to a minimum in the visibility profile, similar to the ring model in Doeleman et al. (2008). At inclinations ≳30°, the profiles monotonically decrease with the baseline length. A detection in place of the current upper limit (open circle with arrow) should favor one set of profiles.
Download figure:
Standard image High-resolution imageOur images are strongly peaked at the circular photon orbit, as can be seen in Figure 1. Low inclination (i ≲ 30°) images produced from a single time step of the energy-conserving simulation by McKinney & Blandford (2009) are similar to that in Figure 1, suggesting that our results are not due to the nonconservative nature of the simulation. It is not surprising that the VLBI observations are well fit by face-on disks. In addition to the ring model, a delta function intensity profile at the circular photon orbit roughly matches the observations, as does the spherical accretion model from Dexter & Agol (2009) scaled to the mass of Sgr A*. In fact, any model sharply peaked at the circular photon orbit should approximately match the visibility data. This is because the ratio of the zeroth to first maximum in the visibility for such a model is roughly the same as that of the zeroth-order Bessel function, and agrees with the ratio of the visibility amplitude at ∼3500 Mλ to the total flux.
If this GRMHD simulation approximates the correct model for Sgr A*, we can also constrain the accretion rate for the single spin value used. Upper and lower limits for 68% and 99% confidence from probability distributions marginalized over observer time and sky orientation are shown in Figure 4. Doppler beaming leads to larger fluxes at high inclinations, and hence a lower accretion rate is necessary to match the observed flux. The limits are also more stringent at high inclinations, where the flux is less variable (see Figure 5). Further marginalizing over inclination, we constrain the accretion rate to be (1.0–2.3)×10−9 M☉ yr−1 at 99% confidence. This constraint on the accretion rate of Sgr A*, while model-dependent, only uses the VLBI visibility measurements.
Figure 4. Upper and lower limits to the average accretion rate as a function of inclination at 68% (triangles) and 99% (plus signs) confidence from probability distributions marginalized over observer time and sky orientation for the a = 0.9 simulation. The narrow allowed regions demonstrate the ability of VLBI observations in conjunction with GRMHD models to constrain the accretion rate of Sgr A*. Other simulations may produce different ranges of allowed accretion rates.
Download figure:
Standard image High-resolution imageFigure 5. Light curves of the last 2/3 of the untilted GRMHD simulation from Fragile et al. (2007) at inclinations of 10° (top) and 70° (bottom). The i = 70° light curve has been shifted downwards by 1.5 Jy for clarity. Both light curves exhibit consecutive flares starting at t ∼ 2 hr, which are consistent with those observed from Sgr A* at millimeter wavelengths. They are more prominent at i = 10° due to higher optical depth between the observer and the flaring region at high inclinations. The arrows denote the 10 minute intervals corresponding to the best-fit images shown in Figure 1.
Download figure:
Standard image High-resolution image3.2. Variability
A major advantage of using GRMHD simulations in place of RIAF models is the ability to probe variability from the same dynamical model used to produce spectra and visibilities. Previous work has used hotspots (Broderick & Loeb 2005; Doeleman et al. 2009) to study variability, but this introduces extra free parameters, complicating attempts to extract black hole parameters. The light curves of our Sgr A* models at 1.3 mm shown in Figure 5 are the first from simulation data to consider synchrotron emission, and the first to include the time delay along geodesics.
The variability on short timescales is slightly more noticeable for the i = 70° case. Interestingly, though, the longer timescale flaring behavior is more pronounced in the i = 10° light curve. The power spectra are well described by red noise spectra with power-law index Γ = 2.4, 1.7 for i = 10°, 70°. These are both steeper than the observed power-law index of Γ = 1 in Mauerhan et al. (2005) at 3 mm wavelength.
The twin flares occur simultaneously in both light curves, rise over half an hour, and last 2–3.5 hr. The flux modulation is 50% (40%) at i = 10° (70°), measured from the peak of the flare to the average of the light curve immediately preceding it. All of these features are consistent with millimeter flare observations of Sgr A* (Eckart et al. 2008b; Marrone et al. 2008; Zhao et al. 2003; Li et al. 2009). Since the flares are seen at low inclinations, they are not caused by Doppler shifts from hotspots. The peak intensity is attenuated at large inclination due to rising optical depth. The flares are caused by a rise in magnetic field strength near the midplane in the inner radii (r ∼ 2–6 M). They are not due to heat from magnetic dissipation since this is not possible in the simulation used here. Thus, heating from magnetic reconnection is not necessary to produce the millimeter flares of Sgr A*.
4. DISCUSSION
In this Letter, simulated images of a single-temperature synchrotron emission model from three-dimensional, GRMHD simulations by Fragile et al. (2007) have been fit to recent 1.3 mm VLBI observations of Sgr A*. For the single black hole spin considered, the high inclination results are in at least qualitative agreement with those presented for 1.5D, stationary, nonrelativistic models used by Broderick et al. (2009). However, we find that, in contrast to previous work, images at low inclination provide excellent fits to the data, which are relatively insensitive to sky orientation and present a qualitatively different visibility profile due to approximate circular symmetry.
Radial intensity profiles from face-on images of recent RIAF models (Yuan et al. 2009) appear more extended than those from our GRMHD simulations. The minimum of the visibility profile for these images occurs at too short a baseline to match the VLBI data. The model in Broderick et al. (2009), on the other hand, has increasing number density up to the event horizon, in stark contrast to the results from the simulations, where the number density peaks at r ∼ 5M, well outside of the event horizon at r = 1.4M and the marginally stable orbit at r = 2.3M. Large number densities in the inner radii of the model from Broderick et al. (2009) cause the “shadow” to be much smaller than expected (Bardeen 1973; Falcke et al. 2000), weakening the power on small scales (large baselines), and preventing a good fit to the VLBI data at low inclinations.
Increased 1.3 mm array sensitivity should favor either the low or high inclination profiles in Figure 3. Our predicted profiles would also be sensitive to measurements at 1.3 mm along baselines of length 2000–3000 km. A future detection favoring a ring-like profile would constitute the first measurement of a black hole shadow.
The visibility analysis constrains the average accretion rate of the a = 0.9 simulation to be (1.0–2.3)×10−9 M☉ yr−1 with 99% confidence. While this constraint is model-dependent, it is independent of and consistent with constraints from spectral fitting (Yuan et al. 2003) and polarization measurements (Agol 2000; Quataert & Gruzinov 2000; Macquart et al. 2006; Marrone et al. 2006; Sharma et al. 2007), and in good agreement with the previous spectral fits of single-temperature GRMHD synchrotron images by Noble et al. (2007).
We have also produced the first light curves of synchrotron emission from GRMHD modeling of Sgr A*. At 1.3 mm, flares arise naturally from magnetic turbulence, and are consistent with observations. Unfortunately, our analysis is poorly suited for studying variability in other wavebands such as the radio, NIR, or X-ray. In the radio, most of the emission originates outside of the initial torus pressure maximum at r = 25M, while the NIR and X-ray luminosity is produced by some combination of synchrotron emission from nonthermal electrons and Compton scattering of synchrotron photons. Modeling of both is left for future work. In addition, polarization measurements from flares are an important observational clue. Adding this capability to our relativistic radiative transfer code is in progress.
There are many uncertainties in the analysis presented in this Letter. While GRMHD allows for a self-consistent solution of the primitive variables, it is unknown whether the MHD approximation is valid for the low densities at the Galactic center. In addition, there is a range of possible solutions, depending on the energy prescription used. Current simulations also neglect the interplay between ions and electrons, which should be handled self-consistently.
If Sgr A* is radiatively inefficient, simulations which conserve total energy are probably more appropriate. We estimate the time-averaged energy dissipation rate from numerical losses in the simulation as ∼1037 erg s−1 for an accretion rate of 1.6 × 10−9 M☉ yr−1. This is larger than the bolometric luminosity of Sgr A*, and gives an effective radiative efficiency of ∼0.1. The accretion flow is advection dominated, in that an order of magnitude more energy is lost to the black hole than to numerical reconnection. Face-on synchrotron images from a single time step of the three-dimensional, conservative code in McKinney & Blandford (2009) were similar to those found here, suggesting that our results may not be specific to nonconservative simulations.
Despite the uncertainties, relativistic radiative transfer is an important tool for connecting state of the art, three-dimensional, global GRMHD simulations with observations. Our analysis shows that ongoing millimeter VLBI observations will not only constrain black hole parameters, but accretion models as well. Expanding this analysis to the full range of spins will help determine the model dependence of such future constraints. If Sgr A* is at a low inclination, millimeter VLBI may soon provide the first direct evidence of a black hole shadow.
In addition to the visibility fits and light curves considered here, the variability of Sgr A* can be probed through measurements of closure phases along triangles of baselines (e.g., Doeleman et al. 2009). Identifying salient features of these closure phases is a goal of future work. If the flares seen here are robust features of millimeter synchrotron light curves of GRMHD simulations, multiwavelength and polarization measurements can be used to constrain emission mechanisms and the nonthermal electron distribution.
J.D. thanks Shep Doeleman for the VLBI data, and Omer Blaes for useful discussions. This work was partially supported by NASA grants 05-ATP05-96 and NNX08AX59H, and a graduate fellowship at the Kavli Institute of Theoretical Physics at the University of California, Santa Barbara under NSF grant PHY05-51164.