THE POWER OF IMAGING: CONSTRAINING THE PLASMA PROPERTIES OF GRMHD SIMULATIONS USING EHT OBSERVATIONS OF Sgr A*

General relativistic magnetohydrodynamic (GRMHD) simulations of accretion flows onto black holes have significantly expanded and improved our understanding of accretion physics (see Abramowicz & Fragile 2013; Yuan & Narayan 2014 for recent reviews). Multiple numerical algorithms have been developed, which allow for accurate simulations of the turbulent magnetohydrodynamic (MHD) flow in the curved spacetime near a black hole, the plunging of matter at the location of the innermost stable circular orbit, the accretion through the event horizon, and the generation of outflows and jets.

An important limitation of all current GRMHD simulations is their inability to follow the thermodynamic properties of the electrons in the flow. This is a crucial ingredient in predicting the observational characteristics of accretion flows onto black-holes, because the emission from such systems is dominated by the radiative properties of the electrons. There are several reasons that hamper rapid progress in this aspect of global simulations. First, the heating and, in general, the acceleration of electrons by the dissipation of the MHD turbulence occurs at sub-grid scales (see, e.g., Riquelme et al. 2012). Second, the transport of heat by conduction is highly anisotropic and takes place in a regime that is poorly understood (see, e.g., Sharma et al. 2008). Finally, radiative energy losses are characterized by cross sections that have very strong dependence on photon energy and render approximate radiative transfer schemes inadequate (see, e.g., Davis et al. 2012).

Improving our understanding of the relevant physical process will be stimulated in the near future by technological advances in spatially resolved imaging observations of black holes. The Event Horizon Telescope (EHT) will perform millimeter very long baseline interferometry (VLBI) observations of a number of supermassive black holes and will achieve horizon-scale resolution for at least two well studied systems, Sgr A* and M87. Early EHT observations have confirmed that the size of the emitting region at 1.3 mm in both systems is at most equal to a few Schwarzschild radii (Doeleman et al. 2008, 2012). This simple piece of evidence, in combination with the broadband spectra of the sources, can already place stringent constraints on the physical conditions and geometries of the accretion flows (see, e.g., Broderick et al. 2009, 2011a; Dexter et al. 2009, 2010, 2012; Dexter & Fragile 2011; Shcherbakov et al. 2012; Psaltis et al. 2014). In the near future, incorporating ALMA and the South Pole Telescope to the array of telescopes that comprise the EHT will generate polarization dependent images of Sgr A* and M87 at two wavelengths (1.3 mm and 0.8 mm) and at different epochs (see, e.g., Doeleman et al. 2009; Broderick et al. 2011b). This wealth of data will provide the observational foundation against which the results of GRMHD simulations will be calibrated.

The power of comparing theoretical models to spatially resolved observations of accreting black holes becomes apparent when exploring the relative importance of the shearing flow (or, for simplicity, the disk) and of the bulk outflow or jet in a system. Albeit typically absent from simple analytic models of accretion disks, outflows, such as highly relativistic jets and winds are common features of GRMHD simulations that often form spontaneously (see, e.g., McKinney & Gammie 2004). The plasma density in an outflow is typically much lower than in the disk but the magnetic field and relativistic Lorentz factor are much larger, potentially dominating the emerging radiation from the system. Indeed, semi-analytic models dominated by a radiatively inefficient accretion flow (see, e.g., Narayan et al. 1998) or by a relativistic jet (see, e.g., Falcke & Markoff 2000) have both been used to explain the spectra and inferred image sizes of Sgr A*. Such conceptually different geometries will be readily distinguishable using the combination of horizon-scale resolution, sensitivity to different polarizations, and ability to follow the variability of emission that will become available with the complete EHT array.

Our aim is to investigate the ability of future EHT observations to distinguish between different emission geometries, black-hole properties, and thermodynamic conditions in the accreting plasma around Sgr A*. In this first paper, we explore a large suite of GRMHD simulations onto black holes (Narayan et al. 2012; Sa¸dowski et al. 2013a) with different black-hole spins, different prescriptions regarding the thermodynamic properties of the electrons, and different magnetic field topologies of the saturated state of the turbulent flow. Among the very large range of possible configurations, we select those that agree with the broadband spectrum of Sgr A* as well as with the initial measurement of the size of its image at 1.3 mm. In follow-up articles, we will use this suite of models to make detailed predictions for EHT observations and to develop observing strategies that will maximize the scientific return of the EHT.

In earlier work, GRMHD simulations have been used to explore the effect of changing the black-hole spin magnitude and orientation with respect to the observer (Dexter et al. 2010; Mościbrodzka et al. 2009), the tilt of the black-hole spin with respect to the angular momentum of the accreting flow (Dexter & Fragile 2013), and the thermodynamic properties of the electrons in the outflow/jet (Mościbrodzka & Falcke 2013; Mościbrodzka et al. 2014). Our work improves on three aspects of these earlier studies. First, we consider the influence of the large-scale magnetic field in the saturated state of the flow, by contrasting the simulations with Magnetically Arrested Disks (MAD) to those with Standard And Normal Evolution (SANE) of Narayan et al. (2012) and Sa¸dowski et al. (2013a). These two types of simulations have very different dynamical behaviors, as well as very different magnetic field topologies near the black-hole horizon. Second, the GRMHD simulations that we are using have been evolved long enough for the flow to reach a dynamical steady state out to about a hundred gravitational radii. This is important because a very large volume of the accretion flow contributes to the low-frequency synchrotron emission as well as to the X-ray bremsstrahlung emission. Finally, owing to our use of a very fast radiative transport algorithm, instead of calculating the time-averaged properties of the simulated flows before constructing images and spectra like in many previous studies, we are able to calculate images and spectra for each snapshot of the flow, before averaging them together. This approach was also taken by Dexter et al. (2009, 2010) and Dexter & Fragile (2013). Albeit computationally very expensive, it mimics more closely the averaging that will inevitably occur during the EHT observations and produces results that differ considerably from the alternate procedure.

In order to achieve the high efficiency of ray tracing required to calculate images and broadband spectra of every snapshot of a GRMHD simulation, we employ our algorithm GRay (Chan et al. 2013). Unlike standard central processing unit (CPU)-based ray tracing algorithms, GRay uses graphics processing units (GPUs) to accelerate the computationally intensive geodesic integration. With careful handling of the access to the memory, GRay achieves an order of magnitude speed up and allows us to compute images during 12 hr of wall time using 32 nVidia Tesla K20X GPUs on the El Gato cluster at the University of Arizona.

In the next section, we describe the GRMHD simulations, the plasma methods, our implementation of the radiative processes, and the ray tracing algorithm that we are using in this paper. In Section 3, we study the dependence of the calculated spectra and images on the various model parameters. In Section 4, we compile the current spectral and imaging observations of Sgr A*, which we then use in Section 5 in order to constrain the range of model parameters. Finally, we discuss the implications of our results in Section 6.

We carry out ray tracing calculations on six sets of three-dimensional GRMHD simulation reported in Narayan et al. (2012) and Sa¸dowski et al. (2013a), which consist of two different classes of accretion flows. The first class, called Standard And Normal Evolution (SANE), uses multi-loop initial magnetic fields, while the second class, called Magnetically Arrested Disk (MAD), uses single-loop initial fields. The accretion flows were evolved for exceptionally long times, e.g., ∼200, 000 GM  c⁻³, until the flows reached steady state up to ∼100 GM  c⁻². Three values of spin parameter, a = 0, 0.7, and 0.9 were used in each class. We summarize the setup of these simulations in Table 1.

2.1. Characteristics of the MHD Simulations

The differences in the initial magnetic fields and black-hole spins affect the dynamic and thermodynamic properties of the accretion flows. Because the plasma-β and the gas temperature are important quantities in the emission models and affect the observables, we first calculate, for each simulation, their temporal and azimuthal averages, which we denote by

and

and plot them in Figure 1. In the above definitions, P_gas and P_mag denote gas and magnetic pressure, while u and ρ denote the internal energy and density.

Zoom In Zoom Out Reset image size

All simulations show a funnel region with low plasma-β and high temperature, with only exception the zero spin SANE simulation a0SANE. This makes a0SANE a special configuration in our subsequent studies. The transition from the funnel region to the disk occurs over a rather narrow sheath, across which the plasma-β changes rapidly (note the logarithmic scale for the colors). This is consistent with the findings of Mościbrodzka & Falcke (2013) and allow us to separate the funnel from the disk by setting a threshold on the plasma-β, which here we fix to β_threshold = 0.2 (see Section 5). The saturated red region around the equators of all the SANE simulations indicates that the accretion disks are strongly gas dominated. The MAD simulations, on the other hand, tend to have stronger magnetic fields even in the disk regions. For more details of the simulations, we refer to the original papers by Narayan et al. (2012) and Sa¸dowski et al. (2013a).

We note that the coordinate singularities along the poles can cause numerical difficulties and can lead to a few low-density hot cells. Although these artificial hot cells do not affect the flow dynamics, they may change the predicted spectra by over producing X-rays and γ-rays. In addition, the simulations use a special coordinate transformation developed by Tchekhovskoy et al. (2011) to "cylindrificate" the grid near the poles. While this technique significantly speeds up the GRMHD simulations by allowing a larger time step, it expands the problematic hot cells and enlarges the error. To overcome these numerical artifacts, we excise the cells around the poles by setting their emissivity to zero when they are unphysically hot compared to their neighbors.

The six numerical models described above are scale free in length and density (their time, velocity, and energy scales can be obtained by scaling length and density with different powers of the speed of light c). In order to compute images for Sgr A*, we fix the length scale by setting the mass of the central black hole to that of Sgr A*, M = 4.3 × 10⁶ M_☉ (Gillessen et al. 2009). We introduce the density scale n_e as a free parameter to convert the gas density ρ into the electron number density, which at the same time determines the physical accretion rate. We also define the observer's inclination angle i with respect to the spin axis of the black hole and treat it as a free parameter.

2.2. Thermodynamics of the Accretion Flow

To calculate spectra, light curves, and images of the accretion flow around Sgr A*, we need to specify the electron temperature in the disk and in the funnel regions, which depend on the details of the heating due to the dissipation of turbulence, the radiative cooling, as well as the energy exchange between the protons and the electrons in the flow. In low density regions, where the collisional timescale is long compared to the dynamical timescale, the electron temperature T_e is expected to be lower than the ion temperature T_i, leading to a two-temperature plasma (Shapiro et al. 1976; Narayan et al. 1995, 1998; see also the recent review by Yuan & Narayan 2014 and references therein). In addition, in the funnel region, where the heating rate due to the dissipation of the MHD turbulence is expected to be low, thermal conduction can be efficient in bringing the electrons to a constant temperature (Mościbrodzka & Falcke 2013; Mościbrodzka et al. 2014). We treat both of these possibilities in the calculation of the radiation from the accretion flow and use the following parameterizations to capture the resulting electron temperatures.

Constant electron–ion temperature ratio model.In the first model for the plasma, we assume that the electron-to-ion temperature ratio is a fixed function of the plasma-β. This is motivated by the fact that the temperature of the electrons can differ significantly from that of the protons as the gas density gets lower and the magnetic field strength, which influences the cooling time for the electrons via synchrotron emission, gets higher. We consider a simple step-function model where the electron temperature takes the following form:

In this equation, the symbols θ_disk and θ_funnel denote the electron temperatures in the disk and in the funnel, in units of the ion temperature. Both quantities are expected to be less than unity under the assumptions that (1) dissipation of turbulence mainly heats up the ions, (2) the electrons cool faster than the ions, and (3) the electron cooling rate is faster than the energy exchange rate between electrons and ions (see discussion in Narayan & Yi 1994, 1995; Quataert & Gruzinov 1999). We specify the electron-to-ion temperature ratio in the disk, θ_disk, and the electron-to-ion temperature ratio in the funnel, θ_funnel, and constrain these parameters in the next section using observations of Sgr A*. Note that this model, which we will refer to as the "constant ratio model" in the remainder of the paper, reduces to the standard two-temperature plasma model used in earlier studies (e.g., Dexter et al. 2009; Mościbrodzka et al. 2009), if we fix θ_disk = θ_funnel.

Constant funnel electron temperature model.The second model accounts for the possible effects of electron conduction in the funnel region. We follow the parameterization in Mościbrodzka & Falcke (2013) and assume that the electron temperature in the funnel is constant. We retain the constant temperature ratio parameterization for the disk. Because the plasma-β is a good indicator of the separation between the highly magnetized, low density funnel region and the higher density disk, we use a threshold value of this parameter to specify the change in the electron temperature, such that

Note that T_{e, funnel} may acquire values larger than unity for ultrarelativistic electrons according to Equation (2). We vary the values of the threshold plasma-β, β_threshold, the electron-to-proton temperature ratio in the disk, θ_disk, and the electron temperature in the funnel, T_{e, funnel}, to investigate the effects of these three parameters on the observables and when fitting the observed spectrum and images of Sgr A*. We refer to this as the "constant temperature model" in the rest of the paper.

2.3. Radiative Processes and Transfer

We calculate the radiation emitted from the accretion flow by solving the radiative transfer equation along null-geodesics through the domain of the GRMHD simulations, using the GRay code. GRay integrates the radiative equation backward from the image plane to the source. This approach has at least two advantages. First, it allows us to solve only for light rays that are normal to the image plane, drastically reducing the number of photon trajectories that need to be integrated. Second, integrating the radiative transfer equation backward allows us to stop the integration when the optical depth is sufficiently large (we use a cutoff at ln 1000 ≈ 6.9) or when the outgoing ray is sufficiently far away (∼1000 GM  c⁻²) from the black hole. ⁴ This significantly speeds up the image calculation for optically thick media. Both of these advantages make ray tracing algorithms much faster than Monte Carlo techniques for our application.

To ensure numerical stability, we follow Younsi et al. (2012) to express the radiative transfer equation in two coupled differential equations:

where λ is the affine parameter, γ⁻¹ ≡ ν₀/ν is the relative energy shift, τ and are the optical depth and Lorentz invariant intensity at frequency ν, and α_{0, ν} and j_{0, ν} are comoving absorption and emission coefficients.

In hot, magnetized accretion flows, synchrotron and bremsstrahlung radiation are the two major radiative processes that contribute to the emission and absorption (see Narayan et al. 1998). In this paper, we assume thermal synchrotron and bremsstrahlung emission and do not treat the possible contribution from non-thermal electrons, which can contribute to both radio and X-ray fluxes (Mahadevan 1998; Özel et al. 2000) and can help explain the variability observed in the X-rays (Chan et al. 2009). For synchrotron emissivities, we use the approximate expression derived by Leung et al. (2011). For thermal bremsstrahlung emission, we use the expression derived in Rybicki & Lightman (1979), with a Gaunt factor taken from Novikov & Thorne (1973). ⁵ Although Compton scattering is generally important in stellar-mass black holes, its contribution to the Sgr A* spectrum is mainly in the optical (see, e.g., Narayan et al. 1998), for which we do not have data to impose any constraints on (see Section 4). Moreover, Mościbrodzka et al. (2009, see their Figure 4) showed that the contribution of Compton scattering to the X-ray flux predicted in GRMHD simulations is typically much smaller than the observed flux, unless the black hole is rapidly spinning and observed from a nearly edge-on orientation. For these two reasons, we neglect here the effects of Compton scattering on the spectrum.

Using GRay for radiation transport and the plasma models described in the previous section, we calculate the broadband spectra emerging from the accretion flow in each GRMHD simulation. We compute a radiation spectrum from a series of snapshots taken from the GRMHD simulations. In order to take into account the emission from near the event horizon as well as at the outer parts of the computational domain, we compute the images using three different fields of view: the innermost images cover a 32r_S × 32r_S plane, the middle images cover a 128r_S × 128r_S plane, and the outermost images cover a 512r_S × 512r_S plane, where r_S ≡ 2 GM  c⁻² is the Schwarzschild radius. The outermost image size is chosen so that the X-ray fluxes computed within the simulation volume converge. We resolve each domain using 512 × 512 pixels and obtain the total flux by appropriately combining and summing the contributions from the three images.

In order to explore the effect of the parameters describing the thermodynamics of the plasma, as well as of the inclination of the observer with respect to the spin axis of the accretion flow, on the spectra and images of the flows, we first calculated a large suite of simulations for a range of values of the parameters. Each of the two plasma models described in the previous section requires five parameters to be fully specified. For the constant ratio model, these are: the normalization of the electron density n_e, the inclination of the observer i, the threshold plasma-β, β_threshold, that separates the funnel from the disk, the electron-to-ion temperature ratio in the funnel, θ_funnel, and the electron-to-in temperature ratio in the disk, θ_disk. For the constant temperature model, these are: the normalization of the electron density n_e, the inclination of the observer i, the threshold plasma-β, β_threshold, the funnel electron temperature T_{e, funnel}, and the electron-to-ion temperature ratio, θ_disk in the disk.

In Figure 2, we show the effect of independently varying each of these parameters on the simulated radio-to-X-ray spectrum for Sgr A*. Each column corresponds to one of the GRMHD simulations listed in Table 1, while each row shows how the spectra change in response to varying one of the five model parameters, while keeping the others fixed. The observational data points shown in the leftmost panels were collected by Broderick et al. (2011a).

Zoom In Zoom Out Reset image size

As mentioned at the beginning of this section, the X-ray fluxes computed within the simulation volume converge at a radius of ∼512r_S. However, this radius is much smaller than the Bondi radius of Sgr A*, ∼10⁵ r_S, at which the accretion flows is expected to still have significant contribution to the X-ray emission (see, e.g., Quataert 2004). Our models do not include this large volume of X-ray emission because the central black-hole is fed by a torus that lies at a few hundred Schwarzschild radii. To consider the contribution of the large scale flow, we computed a set of spectra up to the Bondi radius by using the extrapolation scheme developed in Sa¸dowski et al. (2013b). The results are shown as green dotted curves in the first row in Figure 2. The X-ray fluxes of the extrapolated flows are one-to-two orders of magnitude larger, as expected. Note that even though we show the results from this extrapolation here, we will opt to limit our domain to the inner accretion flow and consider the fraction of X-rays that originate from the inner flow when actually fitting the X-ray data to the models (see Section 4).

Fiducial model. The resulting broadband spectrum for each GRMHD simulation and a "fiducial" plasma model with a representative set of parameters (n_e = 3.2 × 10⁸, i = 60°, β_threshold = 0.2, θ_disk = 0.01, and T_{e, funnel} = 10) is depicted in the first row of Figure 2. Note that this model is not a specific fit to Sgr A* observations. We use the fiducial setup, before embarking on fitting the observations, in order to study several theoretical aspects of the spectra, such as the predicted variability across the spectrum and the different emission regions that give rise to the structures of the images at different wavelengths.

The strong turbulence in the accretion flow naturally causes the resulting spectra and images to be time dependent. This complicates the comparison of the simulations to the observations, as both need to be averaged properly in order to avoid comparing a particular realization of the turbulent flow in the simulations to a different realization in the observed flow. In order to assess the variability of the simulated spectra and images we perform our ray tracing calculation for 100 different snapshots in each of the simulations, as listed in Table 1. Each snapshot of the simulation was obtained at regular time intervals of 10 GMc⁻³. For the mass of Sgr A*, 100 snapshots correspond to ≈5.9 hr, which is similar to the time interval over which EHT observations will take place.

The black solid curves on the panels of the first row in Figure 2 show the mean spectra obtained by averaging the simulated spectra in the 100 snapshots of each simulation. The shaded gray area around them marks the maxima and minima of the spectra emerging from these 100 snapshots and is representative of the expected spectral variability. The red dotted curves on the same panels show the average spectra computed using only the image of size 32r_S × 32r_S, which we refer to as the emission from the inner accretion flow. As expected, when the emission originates from the inner accretion flow, the spectra acquire their largest degree of variability while the opposite is true when the emission originates in a much larger volume. Indeed, in the optically thick, low-frequency region of the spectrum, the emission is weakly variable because the radius of the photosphere at these frequencies is equal to tens to hundreds of Schwarzschild radii and the dynamical timescales there are very long. The variability is also weak at the optically thin, high-frequency region of the spectrum, where the emission is generated by thermal bremsstrahlung over a very large volume around the black hole. In the millimeter to IR range (i.e., ν  ∼  10¹¹–10¹⁴ Hz for the a0SANE simulation), however, the emission is optically thin and originates very close to the event horizon. The characteristic timescales there are very short and both the spectra and images show significant variability.

To visualize this point in a different manner, we show in Figure 3 the averaged snapshot images of the a9MAD simulation, for the fiducial parameters of the plasma model. The different rows, from top to bottom, depict images at frequencies ν = 10¹¹ Hz, 10¹³ Hz, and 10¹⁸ Hz. The different columns, from left to right, are for image sizes equal to 512r_S × 512r_S, 128r_S × 128r_S, and 32r_S × 32r_S. As discussed above, the emission at the lowest and highest frequencies shown originates primarily from large distances away from the black hole. On the other hand, the emission at infrared frequencies originates in a very small region, close to the horizon. The black-hole shadow is obscured at the optically thick radio frequencies but is visible in both infrared and X-ray, which are optically thin. For the X-ray images, although the region around the black-hole shadow has a larger surface brightness than all other regions, most of the (integrated) flux actually originates from a few hundred Schwarzschild radii in our models.

Zoom In Zoom Out Reset image size

In calculating the average spectra and images discussed above, we computed individual spectra and images for each snapshot of the simulations and then averaged together the resulting surface brightness and fluxes on the image plane of the observer. This procedure generates results that can be very different compared to calculating the average hydrodynamic and thermodynamic properties of each simulation and then computing a single image and spectrum for this mean flow. This is because of the fact that the plasma properties are substantially variable and the radiative transfer equation, which we solve along geodesics to calculate the image brightness, is a highly non-linear function of the plasma properties. In the top row of Figure 2, we show as blue dashed curves the spectra computed using the mean properties of the flows. The relative difference between the two averaging procedures is largest in the case of the non-spinning black holes.

In Figure 4, we compare the images computed by the two averaging approaches for a frequency of 10¹¹ Hz, using the a0SANE simulation. The left panel shows the image calculated using the mean properties of the flow and the right panel shows the average of the images in each snapshot. There is a striking difference between the two images. The image of the mean flow (left panel) is very dim and almost invisible in the plot. On the other hand, the average image of the snapshots is dominated by magnetic filaments, which are short lived (and hence do not contribute significantly to the average properties of the flow) but are very bright (and hence dominate the average emission).

Zoom In Zoom Out Reset image size

Finally, comparing the spectra calculated for the different GRMHD simulations but for the same, fiducial plasma model, we find that the flux at the thermal peak increases monotonically as we move from the leftmost to the rightmost columns. In other words, there is more low-frequency emission in the MAD and in the high-spin simulations. This is a direct consequence of the fact that MAD and high-spin simulations are characterized by significant relativistic jets. In contrast, the MAD simulations generate less X-ray radiation than the SANE simulations because there is less flux from the accretion flows at large radii (i.e., r ≳ 128r_S). This is also a direct consequence of the fact that the density in the MAD simulations is more centrally concentrated than in the SANE simulations (see Narayan et al. 2012).

Parameter study. The remaining rows in Figure 2 show how the spectra are affected when we vary one of the parameters of the plasma model, i.e., n_e, i, β_threshold, T_{e, funnel}, θ_disk, and θ_funnel, while holding all other parameters fixed to the fiducial model.

We first vary the density scale in the range n_e = 10⁸ to n_e = 5.6 × 10⁸ and show the result for each GRMHD simulation in the second row of Figure 2. At frequencies ν ≳ 10¹² Hz, increasing the density scales causes an increase in the overall flux that is proportional to . This dependence is trivial to understand for the X-rays, which are generated by optically thin bremsstrahlung emission, since the emissivity of bremsstrahlung scales as the square of the electron density. In the infrared, near the peak of the thermal bump, the emission is primarily due to optically thin synchrotron processes, with an emissivity that scales as j_synch∝n_e B². However, the natural scaling of the GRMHD equations causes the dimensional magnetic field strength in the flow to be always proportional to the square root of the electron density. In other words, , and the optically thin synchrotron emissivity also scales as .

At lower frequencies, i.e., when ν ≲ 10¹² Hz, the emerging flux has a weaker dependence on the density scale. In this part of the spectrum, the inner accretion flow becomes optically thick to synchrotron self absorption. Increasing the density scale causes the cooler outer disk to partially block the hotter inner disk and, therefore, to flatten the dependence of the flux on the electron density. To demonstrate this point, we show in Figure 5 the images at a frequency of ν = 10¹¹ Hz for the a7SANE simulation and for the four different values of the electron density scale. As we increase the density scale, the inner accretion disk and the funnel become brighter as expected. An increasing fraction of that bright emission, however, is obscured by the colder outer disk, reducing the strong dependence of the flux on the scale of the electron density.

Zoom In Zoom Out Reset image size

In the third row of Figure 2, we show the effect of varying the inclination of the observer in the range i = 20° to i = 80°. At low frequencies, i.e., when ν ≲ 10¹² Hz, the flux has a very weak dependence on the inclination. The dependence becomes stronger, however, at somewhat higher frequencies. In order to understand this behavior, we show in Figure 6 the corresponding images of the accretion flow in the a9SANE simulation. At low frequencies (top row), as the observer's inclination increases, a larger fraction of the hot inner flow is obscured by the colder outer disk. This effect cancels out the fact that both the upper and lower funnel regions become visible, and produces an overall weak dependence of the inclination. At intermediate frequencies (bottom row), the optically thin emission originates very close to the black-hole shadow; increasing the inclination causes the Doppler effect to boost more radiation toward the observer and, therefore, increases the resulting flux. In X-rays (not shown in the figure), which originate in a large, quasi-spherical volume, there is very little dependence of the resulting spectrum on the inclination.

Zoom In Zoom Out Reset image size

In the fourth row of Figure 2, we show the effect of varying the threshold plasma-β that distinguishes the disk from the funnel, in the range β_threshold = 0.2 to β_threshold = 1.0. For the parameters of the fiducial model, the funnel is brighter than the disk in the radio to infrared frequencies. Therefore, increasing β_threshold, which creates a larger funnel, increases the radio and infrared flux for all simulations. The situation is reversed, however, at higher frequencies, where the emission is dominated by the disk as shown in the lower row of Figure 6. Note also that, primarily in the MAD simulations, the γ-ray flux is also affected by the threshold value of the plasma-β. This arises from the few hot zones along the pole near the event horizon, as we described in Section 2.2, and we have simply left it here to demonstrate this numerical artifact.

In the fifth row of Figure 2, we show the effect of varying the electron-to-ion temperature ratio in the disk, in the range θ_disk = 0.01 to θ_disk = 0.06. As expected, because the synchrotron and bremsstrahlung emissivities depend on the electron temperature, the flux in the optically thin part of the radio-to-infrared spectrum increases with θ_disk.

Finally, in the sixth row of Figure 2, we show the effect of varying the thermodynamic properties of the electrons in the funnel. In the case of the constant temperature model, we vary the electron temperature in the funnel in the range T_{e, funnel} = 10 to T_{e, funnel} = 56; in the case of the constant ratio model, we vary the electron-to-ion temperature ratio in the funnel in the range θ_funnel = 0.01 to θ_funnel = 0.06. As discussed above, for the fiducial values of the parameters, the emission from the funnel dominates the radio-to-infrared part of the emission. As a result, changing the electron temperature in the funnel has a very large effect on that part of the spectrum. On the other hand, the X-ray emission comes primarily from the outer disk and, therefore, is not affected by the thermodynamic properties of the electrons in the funnel.

In Section 3, we performed an extensive parameter study aiming to understand the dependence of our results on the model assumptions related to the MHD and thermodynamic properties of the accretion flow. Our next goal is to identify the set of models and parameters that are consistent with the current spectral and imaging observations of Sgr A*. The data points in the leftmost panels of Figure 2 correspond to non-simultaneous observations of Sgr A* at a broad range of frequencies and with different instruments. Several different, partially overlapping data sets exist in the literature. The particular data points we are using here were prepared by Broderick et al. (2011a), while a comparable collection of data is available in the more recent review by Falcke & Markoff (2013).

When comparing the results of our simulations to spectral observations, we make the following choices.

• 1.  
  In the radio, we do not consider the data points at frequencies below 10¹¹ Hz, because even a small fraction of non-thermal electrons in the accretion flow (which we do not include here) can affect significantly the flux at these low frequencies (see Mahadevan 1998; Özel et al. 2000; Chan et al. 2009).
• 2.  
  In the infrared, i.e., at frequencies 10¹³ Hz <ν < 10¹⁵ Hz, we do not consider the individual data points but rather require the simulated spectra to fall within the lowest and highest observed fluxes at ν = 1.38 × 10¹⁴ Hz. This is justified by the fact that both the observational data (Genzel et al. 2003; Ghez et al. 2004; Dodds-Eden et al. 2011) and our simulations shown in Figure 2 indicate that the infrared emission is highly variable and specific observed or simulated fluxes will depend entirely on the specific realization of the turbulent flow.
• 3.  
  In the X-rays, we investigated the option of extrapolating the accretion flow out to the Bondi radius following the prescription of Sa¸dowski et al. (2013b, see also discussion of Figure 2) and computing the X-ray emission using this entire volume. We found, however, that our results depended very strongly on the assumed power-law indices of the density and temperature profiles as well as on the choice of the point at which the extrapolation matched the numerical solutions. Instead of following this approach, we opt to use the result of Shcherbakov & Baganoff (2010, see also the discussion in Neilsen et al. 2013) that 10% of the quiescent X-ray flux from Sgr A* originates in a point source and attribute 10% of the observed flux to the emission from our simulated volume.

The size of the image of Sgr A* has been measured over many wavelengths, from the radio to the millimeter (see Falcke & Markoff 2013 for a recent review). At most wavelengths, the size measurement is dominated by the blurring of the image caused by interstellar scattering. There is strong evidence that, at wavelengths below 1 cm, the intrinsic size of Sgr A* can be discerned (see, e.g., Bower et al. 2014). However, because of the λ² dependence of the size of the scattering ellipse, the most accurate measurements occur at the smallest wavelengths. Early EHT observations of Sgr A* at 1.3 mm measured its size at . Correcting for the blurring using the scattering law of Bower et al. (2006) resulted in an inferred intrinsic size of the source equal to  μas (Doeleman et al. 2008).

The inferred image size for Sgr A* was based on fitting sparse visibility data in the interferometric u−v plane with a Gaussian model (alternate models have also been considered; see Doeleman et al. 2008). The images from our simulations, however, have significant asymmetry, either because they are dominated by emission in the funnel or because of Doppler effects in the disk. In principle, in order to compare our simulations to the observed image size, we will need to calculate the predicted scattering broadened visibilities in the u−v plane and compare them directly to the data (as is done, e.g., in Broderick et al. 2011a and in Dexter et al. 2009). However, this introduces additional free parameters in the model, such as the orientation of the black-hole spin vector on the plane of the sky, and the current coverage of the u-v plane is too sparse to allow us to constrain the model parameters significantly better than the simple estimate of the size (see, e.g., the large areas within the confidence contours in Broderick et al. 2011a and in Dexter et al. 2009). For this reason, we follow a more approximate procedure in comparing our simulations to the current estimates of the image size at 1.3 mm.

We take into account the effects of interstellar scattering using the elliptical scattering kernel of Bower et al. (2006) that has a major axis

and a minor axis

with a position angle of the major axis at PA = 78° East of North. Because we do not know a priori the position angle of the black-hole spin axis, we consider the entire range of relative position angles between the major axis of the scattering kernel and the spin axis, ΔPA. For each configuration, we convolve the ray traced image from the simulation with the scattering kernel and fit the resulting blurred image with a single Gaussian profile.

In Figure 7, we show how the relative angle ΔPA between the black-hole spin and the major axis of the scattering kernel affects the blurred image sizes. The top row shows 1.3 mm images for the a0SANE simulation, using the constant ratio model; the bottom row shows 1.3 mm images for the a9SANE simulation using the constant temperature model. Note that, in anticipation of the results presented in the following section, the plasma parameters for these two sets of images correspond to the best-fit values for Sgr A*, which we are going to obtain in the following section.

Zoom In Zoom Out Reset image size

In each row, the first column shows the results of the direct ray tracing simulations. The remaining columns show the same image after taking into account the effects of interstellar scattering, for three values of the relative position angles. The dependence of the inferred image size on this parameter is weak. Moreover, in some cases, because of the asymmetries in the images, the blurred image size is not a monotonic function of the relative position angle. For this reason, for each simulation and for each set of model parameters, we scan the entire range of relative orientation angles and consider the range of predicted sizes as an uncertainty in the model. We then compare the predicted values to the observed image size of .

For each of the GRMHD simulations, which corresponds to different black-hole spins and magnetic field configurations, and for either the constant ratio or the constant temperature models, we use observations described in the last section to constrain the five additional model parameters. In particular, we use, (1) the quiescent flux in the X-rays, (2) the observed flux and spectral shape in the frequency range 10¹¹ Hz ⩽ ν ⩽ 10¹² Hz, (3) the range of observed fluxes at ν ≈ 1.3 × 10¹⁴ Hz, and (4) the image size at 1.3 mm. This is a large parameters space, which is difficult to explore computationally and may lead to large degeneracies between model parameters.

Following Mościbrodzka & Falcke (2013), we fix the threshold plasma-β value that separates the disk from the funnel to β_threshold = 0.2. Changing the value of this parameter would simply lead to correlated changes in the electron-to-ion temperature ratio in the disk, without affecting significantly the overall results. We take advantage of the fact that the spectra are insensitive to the inclination angle at low frequencies (see the third row of Figure 2) and set i = 60° (this value was justified by Psaltis et al. 2014, and will be checked later for consistency). Finally, as we discussed in Section 3, several aspects of the simulated spectra and images depend very weakly on some of the model parameters. This allows us to follow the procedure described below, which leads us to use different aspects of the observations to constrain successive subsets of the remaining model parameters.

Fixing the electron density scale using the X-ray flux. In Figure 2, we showed that, among the five model parameters, the predicted X-ray flux at ν ≈ 10¹⁸ Hz for each simulation is only sensitive to the density normalization n_e when we fix the inclination. This property lets us find a correlation between n_e and θ_disk (with θ_funnel = T_{e, funnel} = 0) such that the simulated X-ray flux agrees with the observed quiescent flux. The result is shown as a set of solid curves in the second row of Figure 8. Note that, because the accretion rates, in code units, are comparable for all simulations (see top panel of Figure 4 in Sa¸dowski et al. 2013a), the density normalizations shown here indicate that a0SANE has the highest physical accretion rate, while a9MAD has the lowest.

Zoom In Zoom Out Reset image size

Fixing the electron temperature in the disk using the millimeter-to-centimeter spectrum. For each value of the electron-to-ion temperature ratio in the disk, the observed quiescent flux in the X-rays sets the electron density scale. Given that we have also fixed the values of two parameters, i.e., i = 60° and β_threshold = 0.2, the predicted millimeter-to-centimeter spectrum of Sgr A* is now only a function of the remaining two parameters, which are the electron-to-ion temperature ratio in the disk, θ_disk, and either the constant electron density in the funnel T_{e, funnel} or the electron-to-ion temperature ratio in the funnel θ_funnel, depending on which plasma model we are considering. For each GRMHD simulation and for each plasma model, we show contours of χ² values obtained by comparing the model prediction to the observed millimeter-to-centimeter spectrum of Sgr A*. In each panel, dark blue/purple colors represent small χ² values (see the color bar) and hence are better fits.

Ruling out bright funnels using the infrared flux. In Figure 2, we showed that the predicted flux at infrared wavelengths is sensitive to the thermodynamic properties of the electrons in the funnel, with a weaker dependence on the electron temperature in the disk. In particular, when the electron temperature or the electron-to-ion temperature ratio in the funnel become large (and the funnel becomes very bright), the predicted infrared flux becomes too large to account for the range of observed infrared fluxes from Sgr A*. The same is true for the electron-to-ion temperature in the disk. In Figure 8, we use forward slanted lines "///" to show the regions of the parameter space that generate infrared fluxes outside of the observed range.

Rejecting models using the 1.3 mm image sizes. The final constraint on the model parameters arises from the comparison of the predicted to the observed 1.3 mm image sizes. We use backward slanted lines "\\\" to indicate the parameter space that is excluded by this constraint. With this final constraint folded in, we find the best-fit set of parameters for each simulation within the allowed regions of the parameter space and mark them with white circles or crosses in Figure 8 and list them in Table 2. The white circles in Figure 8 are for fits with χ² < 1.5, which we consider to be good fits; the white crosses are for χ² > 1.5. Note that for the a0SANE models and the constant temperature a0MAD temperature model, no set of parameters satisfies all of the observed constraints.

In Figure 9, we plot the best-fit model spectra together with the broadband spectral data that we use for the fitting. The gray band between ν ≈ 10¹¹ Hz and ≈10¹² Hz marks the frequency range over which we perform the least-squares fits. The gray line at ν ≈ 10¹⁴ Hz marks the infrared frequency at which we use the range of fluxes observed at different times to impose an upper and lower bound on the models. The gray line at ν ≈ 10¹⁸ Hz marks the X-ray frequency where we used 10% of the observed quiescent flux (the open circle below the bow-tie) to fix the density normalization in the flow. The white dotted line inside the gray band marks λ = 1.3 mm, where we use the EHT measurement of the image size. The five curves correspond to the best-fit models with parameters marked as white circles in Figure 8 and listed in Table 2. For these fits, the two SANE simulations have spectral peaks near ν ≈ 10¹³ Hz with , and are indistinguishable from each other. The remaining three MAD simulations have spectral peaks near ν ≈ 3 × 10¹² Hz with half an order of magnitude lower fluxes at these frequencies.

Zoom In Zoom Out Reset image size

As a final consistency check, we plot in Figure 10 the range of image sizes for the best-fit models as functions of the observer's inclination. As discussed in Psaltis et al. (2014), in models that are dominated by disk emission, the inclination affects primarily the size of the 1.3 mm image, because of the Doppler effect. This is clearly seen in the first two panels that correspond to the SANE simulations. However, the dependence of the 1.3 mm image size on inclination for the a0MAD simulation is very weak, with a size that actually increases with inclination. In all cases, the value of i = 60° that we have adopted for the inclination of the observer is consistent with observations.

Zoom In Zoom Out Reset image size

In this paper, we analyzed the predicted spectra and images of various GRMHD simulations and plasma models for the accretion flow around Sgr A*, and narrowed the model parameter range by imposing the requirement that our predictions be consistent with the observed spectra and the 1.3 mm image size of Sgr A*. As seen in Figure 9, the best-fit model spectra agree with observational data very well. The models are most variable at infrared wavelengths, which is also consistent with observations. It is important to emphasize here that, if the spectral properties of Sgr A* were considered without any image size constraints, the range of allowed model parameter would be significantly wider. This demonstrates the power of using images with horizon-scale resolution to distinguish between models that would otherwise make seemingly similar predictions.

Future imaging observations from the completed EHT is expected to more easily distinguish disk-dominated from all funnel-dominated models. In Figure 11, we show the predicted 1.3 mm images for Sgr A* for the GRMHD simulations and plasma model parameters that are consistent with all current spectra and imaging observations. With the completed EHT, it will be straightforward to distinguish between the first two pairs of images from the remaining three pairs of images shown, i.e., the disk-dominated from the funnel models, respectively. However, distinguishing among the various disk-dominated models and measuring the parameters of the black hole and of the plasma, will require additional information from the EHT, including the polarization and scale-dependent variability. We will explore these aspects of our models in future work.

Zoom In Zoom Out Reset image size

This research is conducted at the University of Arizona and is supported by NSF grant AST 1312034. The Arizona team and R. N. acknowledge the support of NASA/NSF TCAN awards NNX14AB48G and NNX14AB47G, respectively. All ray tracing calculations were performed in the El Gato GPU cluster that is funded by NSF award 1228509. We thank Alexander Tchekhovskoy for providing the original derivations of the coordinate transformation used in the HARM simulations. We also thank Eliot Quataert and Jason Dexter for numerous discussions. F. Ö. thanks the Miller Institute for Basic Research in Science at the University of California Berkeley for their support and hospitality.