GLOBAL GENERAL RELATIVISTIC MAGNETOHYDRODYNAMIC SIMULATIONS OF BLACK HOLE ACCRETION FLOWS: A CONVERGENCE STUDY

The numerical study of black hole accretion flows has advanced significantly in the last decade. The advent of techniques for numerically solving the equations of general relativistic magnetohydrodynamics (GRMHD) has enabled self-consistent global modeling of accretion driven by the magnetorotational instability (MRI; Balbus & Hawley 1991; Gammie 2004) onto rotating black holes. Qualitative aspects of these simulations are code independent (e.g., De Villiers & Hawley 2003; Gammie et al. 2003; Anninos et al. 2005), but quantitative variations raise the question of numerical convergence. Recent work has shifted focus from dynamical properties of the accretion flow to simulated observations that can potentially constrain parameters for particular sources such as Sgr A* (Dolence et al. 2009; Mościbrodzka et al. 2009; Dexter et al. 2009, 2010), including polarized radiative transfer (Shcherbakov et al. 2010). To assess the credibility of these radiative models, it is necessary to assess quantitative convergence of the underlying GRMHD simulations.

Convergence studies of global accretion models are computationally expensive. An alternative is to use a local model with shearing box boundary conditions to study the dynamics of MRI-driven turbulence. These are simpler in the sense that there are fewer free parameters, and cheaper in that numerical resolution can be focused on a few correlation volumes ∼H³, where H is the disk scale height. The local model has for decades been a key theoretical tool for probing astrophysical disks (e.g., Goldreich & Lynden-Bell 1965; Goldreich & Tremaine 1978; Narayan et al. 1987) coupled to the shearing box boundary conditions and has been widely used for the study of magnetized disks (e.g., Hawley & Balbus 1991, 1992; Hawley et al. 1995, 1996; Stone et al. 1996; Sano et al. 2004; Hirose et al. 2006; Fromang & Papaloizou 2007; Fromang et al. 2007; Guan et al. 2009; Davis et al. 2010; Fromang 2010; Guan & Gammie 2011; Simon et al. 2011).

Shearing box models have been integrated (1) with or without a mean magnetic field; (2) with or without stratification; (3) with or without explicit dissipation; and (4) with or without explicit treatment of energy transport. There are now dozens of shearing box studies that treat aspects of this problem. The only models that clearly do not converge are unstratified, zero-net field models without explicit dissipation (Fromang & Papaloizou 2007). These models have a magnetic field correlation length that decreases proportional to the grid scale (Guan et al. 2009). But with explicit dissipation (Lesur & Longaretti 2007; Fromang 2010), a mean field (Hawley et al. 1995; Guan et al. 2009), or stratification (Davis et al. 2010; Simon et al. 2011), the models do converge. One of the best resolved studies is Davis et al. (2010), who convincingly demonstrate convergence of a stratified, isothermal, zero explicit dissipation model with the athena code at a physical resolution of up to 128 zones per scale height H. These stratified local models are physically closest to global simulations (e.g., Hirose et al. 2004), which are dominated by toroidal magnetic field. Local studies have shown, therefore, that with sufficient resolution numerical studies of MRI-driven turbulence can converge.

Local models can focus on a few H³, while global simulations must contain many H³. Are any of the dozen or so global disk models (e.g., Brandenburg 1996; Matsumoto et al. 1996; De Villiers & Hawley 2003; Gammie et al. 2003, 2004; De Villiers et al. 2003; McKinney & Gammie 2004; McKinney 2006; Fragile et al. 2007, 2009; Beckwith et al. 2008, 2009, 2011; Shafee et al. 2008; Fragile & Meier 2009; Noble et al. 2009, 2010; Penna et al. 2010; Flock et al. 2011; Hawley et al. 2011, and many others) converged? And are synthetic observations based on global models (e.g., Dexter & Fragile 2011; Hilburn et al. 2010; Mościbrodzka et al. 2009; Dexter et al. 2009; Noble et al. 2007; Schnittman et al. 2006; J. C. Dolence et al. 2011, in preparation) sensitive to resolution? While some authors have included limited resolution studies (e.g., Shafee et al. 2008; Noble et al. 2010; Penna et al. 2010), the answer is not yet clear.

The first systematic convergence test of a global black hole accretion simulation was done by Hawley et al. (2011, hereafter HGK), using a zeus type code to simulate an H/R ≈ 0.1 disk in a pseudo-Newtonian potential. HGK simulate a π/2 wedge in azimuth, varying resolution around a fiducial 256 × 288 × 64 (r, z, ϕ in cylindrical coordinates). After reviewing local model simulations and global non-relativistic models, HGK find that a minimum of 10 cells per vertical characteristic MRI wavelength is required for convergence (HGK's Q_z; e.g., Sano et al. 2004), and 20 cells per azimuthal MRI wavelength (HGK's Q_ϕ). They conclude that most global simulations to date are far from resolved, except Noble et al. (2010) which used barely adequate poloidal resolution.

In this paper we study the same convergence problem considered by HGK, but (1) in relativistic MHD and (2) using slightly different diagnostics. We ask what resolution is required for convergence (if convergence can be achieved), and how the global resolution requirements are related to local models. We are also particularly interested in whether resolution influences the spectra calculated from the models in the weakly radiative limit. This requires a fully relativistic simulation since in weakly radiative accretion flows much of the emission arises from plasma near or even inside the innermost stable circular orbit (ISCO) of a spinning black hole. At these radii, the relativistic models incorporate the dynamics of the plunging region and strong lensing effects on the radiation field.

A third contrast with HGK is that we simulate a full 2π in azimuth rather than π/2. Our experience suggests that there is structure in the disk in all azimuthal Fourier components, with the most power in the m = 1 component. Models with small azimuthal extent have reduced field strength and therefore require higher physical resolution by the HGK Q criteria.

We proceed as follows. Section 2 describes the code and initial and boundary conditions. Section 3 describes convergence of radial profiles of non-dimensional variables. Section 4 describes convergence of azimuthal correlation lengths. Section 5 describes convergence of simulated spectra calculated with a Monte Carlo code. Section 6 gives a brief summary.

Throughout the paper, we follow the standard notation of Misner et al. (1973) and set GM = c = 1. We consider a test fluid (no self-gravity) in the Kerr metric with dimensionless spin a* = 1 − 2⁻⁴ ≈ 0.94. The governing GRMHD equations express conservation of particle number

$$\begin{equation} (\rho u^{\mu })_{;\mu }=0 , \end{equation} \tag{ 1 }$$

and conservation of energy-momentum

$$\begin{equation} {T^{\mu }}_{\nu ;\mu }=0 , \end{equation} \tag{ 2 }$$

together with the source-free Maxwell equations

$$\begin{equation} ^{*}\!{F^{\mu \nu }}_{;\nu }=0 , \end{equation} \tag{ 3 }$$

where u^μ, ρ, T^μν, and *F^μν are the fluid's four-velocity, rest-mass density, GRMHD stress-energy tensor, and dual of the electromagnetic field (EMF) tensor, respectively. The equation of state is

$$\begin{equation} p = (\gamma - 1) u, \end{equation} \tag{ 4 }$$

where γ = 13/9, appropriate for a collisionless plasma with relativistic electrons and non-relativistic protons.

We evolve the GRMHD equations using the HARM3D code (Noble et al. 2006, 2009; Gammie et al. 2003). HARM3D is a conservative high-resolution shock-capturing scheme demonstrated to have second-order convergence in space and time for smooth flows. The zone-centered magnetic field is updated with flux-interpolated constrained transport (flux-ct: Gammie et al. 2003; Tóth 2000) which preserves a particular numerical representation of $\nabla \cdot \mathbf {B}=0$. For this study, we use piecewise parabolic interpolation for both fluxes and EMFs.

The numerical grid is uniform in modified Kerr–Schild coordinates x₁, x₂, and x₃ (Gammie et al. 2003):

$$\begin{eqnarray} x_1 &= \ln r \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \theta &= \pi x_2 + h \sin (2\pi x_2) \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} x_3 &= \phi , \end{eqnarray} \tag{ 7 }$$

where r, θ, and ϕ are the Kerr–Schild radius, colatitude, and azimuth, respectively. We set h = 0.35 to concentrate the grid near the equatorial plane. The grid extends from below the horizon to r = 40, [0.017π, 0.983π] in colatitude, and [0, 2π) in azimuth. HARM3D sets a "floor" for density and internal energy to avoid numerical problems that arise when those values are low: ρ_min = 10⁻⁴r^−3/2 and u_min = 10⁻⁶r^−5/2.

The initial condition is an equilibrium, prograde torus (Fishbone & Moncrief 1976) with inner edge at r = 6, pressure maximum at 12, and outer edge at 40. To make the torus unstable to MRI, it is seeded with weak poloidal magnetic field whose vector potential is

$$\begin{equation} A_{\phi } = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}C (\rho /\rho _{{\rm max}} - 0.2)& {\rm if}\; A_{\phi } > 0 \\\\ 0& {\rm if}\; A_{\phi } \le 0, \end{array}\right. \end{equation} \tag{ 8 }$$

where C is a constant and ρ_max is the maximum initial density. This gives dipole field line loops that run parallel to density contours. The field strength is normalized so that the ratio of the maximum gas pressure to the maximum magnetic pressure β is 100. Small perturbations are introduced into the initial conditions to seed the MRI. The density and magnetic field lines are shown in Figure 1 for the initial conditions and for a later snapshot of the turbulent accretion flow.

Zoom In Zoom Out Reset image size

The models have outflow boundary conditions at the inner and outer radial (x₁) boundaries and periodic boundary conditions in the azimuthal (x₃) direction. The remaining (x₂) boundaries are offset slightly from the pole, so the grid excludes a narrow cone around each pole. This avoids having the last polar zone control the timestep via the Courant condition because the polar zones become narrow in x₃ (the computational expense is proportional to N⁵_x if poles are included). While this treatment is essential for a convergence study, it is difficult to implement an appropriate boundary condition on the cone. We consider two different polar-boundary conditions.

The first, "hard" boundary is a solid reflective wall. We manually set the flux through the boundary to zero, and adjust the EMF in the flux-ct routine to make the cutout completely opaque to the magnetic field, since the field vectors are modified in the routine after setting the boundary condition. This boundary condition produces an unphysical relativistic flow in the grids along the polar cone, so in addition we force the poloidal velocity in the zones along the boundary to be zero.

The second, "soft" boundary also models a reflective wall. The variables in the ghost zones are all copied from the first physical zone. The x₂ components of the velocity and magnetic field are inverted across the boundary (as usual for reflecting boundaries), but this only zeros fluxes on the boundary to within truncation error. This version of the polar-boundary condition permits some leakage of magnetic flux through the polar boundaries, but does not produce unphysical flows along the boundary.

We ran a low-resolution simulation with no polar cutout to evaluate both boundary conditions. The results suggest that the difference between the boundary conditions does affect the evolution of the high latitude "funnel" region. The soft boundary condition, in particular, causes a steady drop in the funnel region magnetic flux. On the other hand, all three cases (hard, soft, and no cutout) exhibit remarkably similar disk evolution.

Our runs have numerical resolutions $(N_{x_1}, N_{x_2}, N_{x_3})$ = (96, 96, 64), (144, 144, 96), (192, 192, 128), and (384, 384, 256). The runs last until t_f = 16, 000 for 96 × 96 × 64, 12, 000 for 144 × 144 × 96, 10, 000 for 192 × 192 × 128, and 6000 for 384 × 384 × 256. Each resolution is run for both the soft and hard polar-boundary conditions except the highest resolution case which is run only for the soft polar boundary due to numerical expense. A list of runs is shown in Table 1. The runs required ${\approx} 10^6 (N_{x_1}/384)^4 (t_f/{\rm 6000})$ cpu hours on TACC ranger.

Table 1. List of Runs

Resolution	Duration ($\frac{GM_{{\rm BH}}}{c^3}$)	Polar Boundary Type
96 × 96 × 64	16,000	Soft
144 × 144 × 96	12,000	Soft
192 × 192 × 128	10,000	Soft
384 × 384 × 256	6,000	Soft
96 × 96 × 64	16,000	Hard
144 × 144 × 96	12,000	Hard
192 × 192 × 128	10,000	Hard

Download table as:  ASCIITypeset image

Each simulation's initial data contain noise inserted in each zone with a random number generator. This noise seeds the growth of instabilities in the torus. Each run will therefore differ in the details of the evolution, but over long enough periods one expects the differences to average away. Nevertheless, because our runs have finite duration, we expect some "cosmic variance," and this noise from run-to-run variations is present in every measurement we use to evaluate convergence.

To evaluate run-to-run variation, we have repeated each of the $N_{x_1} = 96$ and $N_{x_1} = 144$ runs three times, and have used the variance of these runs to attach error bars to our measurements. We find that large run-to-run variations are caused by "events" that last a non-negligible fraction of the simulation time. For example, the lowest resolution runs sometimes gather a large mass of plasma near the ISCO, then dump it suddenly into the black hole. We have also observed a bundle of magnetic field directed opposite to the field in the funnel merged into the funnel, leading to a large fluctuation in the run with resolution 144 × 144 × 96 and hard polar boundary. While the nature, frequency, and origin of these events are still unclear (we have only a handful of runs), it appears that run-to-run variation decreases at higher resolution.

We will compare poloidally, azimuthally, and time-averaged radial profiles of the flow variables for the different resolution runs. We take a density-weighted average to focus on the accretion flow within ∼ H of the equatorial plane. The explicit expression for the averaged radial profile F(x₁) for a variable f is

$$\begin{equation} F(x_1) = \frac{\int _{t_1}^{t_2} \bar{f}(t,x_1) dt}{t_2-t_1}, \end{equation} \tag{ 9 }$$

where

$$\begin{equation} \bar{f}(t,x_1) = \frac{ \int _{(x_2)_1}^{(x_2)_2} \int _{(x_3)_1}^{(x_3)_2} \sqrt{-g} \rho (t,\vec{x}) f(t,\vec{x}) dx_2 dx_3 }{ \int _{(x_2)_1}^{(x_2)_2} \int _{(x_3)_1}^{(x_3)_2} \sqrt{-g} \rho (t,\vec{x}) dx_2 dx_3 } \end{equation} \tag{ 10 }$$

is the density-weighted poloidally and azimuthally averaged radial profile of the variable f and $g=g(\vec{x})$ is the determinant of the metric. For our case, ((x₂)₁, (x₂)₂) = (0.01, 0.99) and ((x₃)₁, (x₃)₂) = (0, 2π).

We compare only non-dimensional variables since dimensional variables depend on the accretion rate, which decreases in time as the initial torus is accreted by the black hole. Our choice of the non-dimensional variables is scaled electron temperature θ_e = kT_e/m_e (=m_pp_g/(2m_eρ) if T_e = T_p) and β ≡ p_g/p_B = (Γ − 1)u/(b²/2), where b² ≡ b^μb_μ,

$$\begin{eqnarray} b^{\mu } &\equiv \frac{1}{\gamma }({g^{\mu }}_{\nu }+u^{\mu }u_{\nu })\mathcal {B}^{\nu } \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \mathcal {B}^{\mu } &\equiv -n_{\nu }{^* \! F}^{\mu \nu } \,\, {\rm where} \,\, n_{\mu }=(-\sqrt{-1/g^{tt}},0,0,0) \, , \end{eqnarray} \tag{ 12 }$$

γ is the Lorentz factor of the flow measured in the normal observer's frame, and Γ, k, m_p and m_e, and T_p and T_e are the adiabatic index, Boltzmann constant, proton and electron mass, and proton and electron temperature, respectively. When calculating β we average p_g and p_b separately using Equation (10) and take the ratio of the averages. This prevents zones with near-zero magnetic energy from dominating the average.

Figure 2 shows the radial profile of β and θ_e calculated using Equation (9) for all the runs. All time averages run from t = 4000 to the end of the run; at t = 4000 the disk at r ≲ 10 is in a steady state for all runs except for the lowest resolution model, which shows a clear upward trend in β over the entire run. The lowest (96 × 96 × 64) and medium (144 × 144 × 96) resolution runs are averaged over three runs with different initial seeds to reduce run-to-run variation. The figure shows profiles for both the hard and the soft polar-boundary conditions described in Section 2.

Zoom In Zoom Out Reset image size

Figure 3 shows β and θ_e plotted against radial resolution $N_{x_1}$ for r = 2.04 (ISCO) and 8. The soft and hard polar-boundary results are shown as solid black and red lines, respectively. Most quantities vary sharply from $N_{x_1}=96$ to 144 and then far less at higher resolution. For example, the soft polar-boundary models have β(ISCO) = (11.6, 7.3, 7.8, 6.6) and θ_e(ISCO) = (31, 47, 48, 57) at the four resolutions.

Zoom In Zoom Out Reset image size

Note that at resolutions greater than 144 × 144 × 96 there are only small quantitative differences between the hard and soft polar-boundary conditions, as seen in Figures 2 and 3. We conclude that the effect of the polar-boundary conditions on the main, equatorial flow is small for these dimensionless variables.

What part of the variations at $N_{x_1} \ge 144$ is real variation with resolution, and what part is run-to-run noise? The error bars in Figure 3 show standard deviation of the three runs performed for the lowest (96  ×  96  ×  64) and medium (144  ×  144  ×  96) data points with different initial seeds. Error bars are not available for the higher resolution data points due to computational expense. The size of the error bars is comparable to the differences between models run with different resolution. One might hope to gain additional information by measuring, e.g., β at several radii and averaging the trend with resolution, but, interestingly, the entire radial profile varies in a correlated way. Nevertheless, Figures 2 and 3 show a clear trend of decreasing β and θ_e with increasing resolution. It seems likely, therefore, that there is a genuine but weak trend with resolution.

We have looked at one-point statistics for non-dimensional variables. What about two-point statistics, which measure the spatial structure of the turbulence, and in particular the correlation length? The correlation length is a natural measure of the outer scale of the turbulence, and should be resolved and independent of resolution in a converged simulation.

We consider only the azimuthal correlation length, as this is most straightforward to compute, and is most often underresolved in global simulations (HGK). The correlation function at radius r on the equatorial plane is

$$\begin{equation} R(\phi) = \int _{0}^{2\pi } \delta f(\phi _0) \delta f(\phi _0 + \phi) d\phi _0 \, , \end{equation} \tag{ 13 }$$

where δf is deviation from the average value of variable f at r. In practice, we average R over a small area rΔrΔθ across the equatorial plane, normalize, and average in time:

$$\begin{eqnarray} \bar{R}(r,\phi ,t) &= \int _{r\pm \frac{\Delta r}{2}, \, \pm \Delta \theta } R(r,\theta ,\phi ,t)rdrd\theta / (r\Delta r \Delta \theta) \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} \bar{R}(r,\phi) &= \int _{t_1}^{t_2} R(r,\phi ,t)/R(r,0,t) dt \, . \end{eqnarray} \tag{ 15 }$$

Note that the correlation function for magnetic field is defined as

$$\begin{equation} R(\phi) = \int _{0}^{2\pi } \delta b^{\mu }(\phi _0) \delta b_{\mu }(\phi _0 + \phi) d\phi _0 \, , \end{equation} \tag{ 16 }$$

where b^μ is defined in Section 2. Then

$$\begin{equation} \lambda : \bar{R}(r,\lambda) = \bar{R}(r,0)/e \, \end{equation} \tag{ 17 }$$

is the correlation length at radius r.

Figure 4 shows the azimuthal correlation length for density ρ, internal energy u, magnetic field b, and θ_e for all runs. Evidently the correlation lengths (angles) are nearly independent of r, except close to the outer boundary where the models are not in a steady state. The correlation length varies from about 0.2π at the lowest resolution to 0.1π at the highest resolution for all variables except b. Since H/r ∼ 0.3 for all models⁵ over a wide range in radius (Figure 5), this corresponds (assuming flat space geometry) to 1–2 vertical scale heights.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The non-dimensional resolution λ/Δϕ ≃ 12(λ/(H/r))  × $(N_{x_1}/384)$, where $\Delta \phi = 2\pi /N_{x_3}$, is marginal even for our highest resolution simulation. For b, the correlation length of the highest resolution is smaller than that for any other variable. The magnetic field structure is underresolved.

Figure 6 plots correlation length against resolution at the ISCO for the same variables as in Figure 4; here, red is the hard polar boundary and black is the soft polar boundary. The dotted lines show how the correlation length would vary if it were fixed at 2, 5, and 10 grid zones.

Zoom In Zoom Out Reset image size

For ρ, u, and θ_e (the non-magnetic variables), the correlation length is ∼5 grid zones for the two lowest resolution simulations. At higher resolution—$N_{x_1} = 192$ and 384—the correlation length increases to >10 grid zones, and the slope of the change in correlation length with resolution decreases. This suggests that for the two highest resolution runs some structures in the turbulence are beginning to be resolved.

For b, on the other hand, the correlation length decreases nearly proportional to the grid scale, with the correlation length fixed at around five grid zones per correlation length. There are small signs of an increase at the highest resolution, but in light of run-to-run variations the significance of this increase is marginal at best. The outer scale for the magnetic field is not resolved.

For all variables the correlation lengths for hard and soft boundary polar conditions are consistent. Evidently, the polar boundary does not influence the structure of turbulence in the equatorial disk.

How do these correlation lengths correspond to those found in local model simulations? Guan et al. (2009) found in their unstratified shearing box model that the three-dimensional correlation function was a triaxial ellipsoid elongated in the azimuthal direction and tilted into trailing orientation. The relationship between our azimuthal correlation length λ_b and the Guan et al. (2009) results is

$$\begin{equation} \lambda = \left(\frac{\cos ^2 \theta _{{\rm tilt}}}{\lambda _{{\rm maj}}^2} + \frac{\sin ^2\theta _{{\rm tilt}}}{\lambda _{{\rm min}}^2} \right)^{-1/2}, \end{equation} \tag{ 18 }$$

where θ_tilt ≈ 15° is the tilt angle of the correlation ellipse, and λ_maj and λ_min are the major and minor axes of magnetic correlation lengths, respectively. For the best resolved net azimuthal field model in Guan et al. (2009) (y256b, which like our global models saturates at β ≃ 20), this implies λ ≃ 0.17H ≃ 0.05 rad, or 0.016π rad. Therefore, it is surprising that correlation length as large as ≃ 0.3 rad  ∼ H is measured in our model for the non-magnetic variables.

Davis et al. (2010) have computed correlation lengths in stratified, isothermal models with zero net flux. In a model run with athena at a resolution of 64 zones per scale height, the implied azimuthal correlation length (averaged over −H < z < H) for the magnetic field is slightly larger than in the unstratified models of Guan et al. (2009), about 0.23H, or 0.02π rad. Guan & Gammie (2011) have also run stratified, isothermal models at lower resolution with a zeus type code. They find an implied azimuthal midplane correlation length (similarly averaged) for the magnetic field that is even larger, about 0.9H, or 0.09π rad. Since correlation length decreases with increasing resolution, it is possible that Guan & Gammie (2011) are not resolving the correlation length, and that at higher resolution the correlation length would be closer to that measured by Davis et al. (2010).

The correlation length of our highest resolution run spans 0.6(H/r) to 0.4(H/r) from ISCO to r ∼ 10 where the corresponding β is 7 and 16, respectively. This is larger than the stratified shearing box results of Davis et al. (2010) but smaller than that of Guan & Gammie (2011). To resolve the correlation length found in Davis et al. (2010), we would need another factor of two in linear resolution. Note that recently Beckwith et al. (2011) found in their global thin disk MHD simulation that azimuthal correlation length to be about 1.3(H/r) by averaging |z| < H and 5 < r < 11. This is larger than our result but also falls between Davis et al. (2010) and Guan & Gammie (2011).

An interesting application of GRMHD models is to simulate observations of sources such as Sgr A* (Dolence et al. 2009; Mościbrodzka et al. 2009; Hilburn et al. 2010; Dexter et al. 2009, 2010; Dexter & Fragile 2011). Are the simulated spectra converged?

The dynamical models underlying the spectral models are run with zero cooling, and the spectra are produced in a post-processing step. This is self-consistent as long as the flows are advection dominated: the accretion timescale is much shorter than the cooling timescale. We calculate the emergent radiation using grmonty, a general relativistic Monte Carlo radiative transfer code (Dolence et al. 2009).

grmonty makes no symmetry assumptions and includes synchrotron emission, absorption, and Compton scattering. Using the rest-frame emissivity for a hot, thermal plasma (Leung et al. 2011), the code produces Monte Carlo samples of the emitted photons—"superphotons" that carry a "weight" representing the number of photons per superphoton. The superphotons follow geodesics, with weight varying continuously due to synchrotron absorption. They also Compton scatter and produce new, scattered superphotons with weight proportional to the scattering probability. We use a "fast light" approximation, where for each snapshot of simulation data a spectrum is created by treating the fluid variables as if they were time-independent. This approximation is excellent for the time-averaged spectra we consider here. Superphotons that reach large radius are collected in poloidally and azimuthally distributed bins, and each bin produces a spectrum. A complete description of the code is given in Dolence et al. (2009).

To compare runs we generate spectra for 200–1200 time slices (depending on the length of the run) and time-average them. The spectrum of each time slice is produced from azimuthally averaged bins that extend from 0.12π < |θ − π/2| < 0.18π rad with respect to the equatorial plane.

We modify the simulation-provided data in one respect before calculating the spectrum. The quality of the non-magnetic fluid variable integration in the funnel region is poor due to truncation error. In particular, the temperature can be high (θ_e > 10⁴) and the particle density is determined entirely by a density floor in HARM3D. We therefore zero the emissivity if b²/ρ > 1 to avoid contaminating the spectrum with possibly unphysical emission.

It is necessary to fix a mass, length, and time unit to generate a radiative model. The combination GM_BH sets a length and time scale but not a mass scale because the mass of the accretion flow is negligible in comparison to the black hole. We set M_BH = 4.5 × 10⁶M_☉, comparable to the mass of Sgr A*. The mass unit for the torus $\mathcal {M}$ is still free; we set it so that the 1.3 mm flux matches the observed flux from Sgr A* of ≃ 3.4 Jy (Marrone et al. 2006).

We want to model emission from a statistically stationary accretion flow. Because we start with a finite-mass torus and it accretes over time, however, there is a steady decrease in density, field strength, accretion rate, etc., as the simulation progresses. We scale away this long-term evolution using a smooth model, as follows. We set the mass unit $\mathcal {M} = M_{0} s(t)$, where M₀ is a constant and s(t) is a two-parameter scaling function. Then

$$\begin{eqnarray} \rho _{{\rm unit}} &=& \mathcal {M}\big/\left(\frac{{\it GM}_{{\rm BH}}}{c^2}\right)^3 \quad u_{{\rm unit}} = \rho _{{\rm unit}}c^2 \nonumber\\ B_{{\rm unit}} &=& c\sqrt{4\pi \rho _{{\rm unit}}} \, , \end{eqnarray} \tag{ 19 }$$

or expressing with s(t),

$$\begin{equation} \rho _{{\rm unit}} = \rho _{0}s(t) \quad \quad u_{{\rm unit}} = u_{0}s(t) \quad \quad B_{{\rm unit}} = B_{0}\sqrt{s(t)} \end{equation} \tag{ 20 }$$

where they are the unit mass density, internal energy, and magnetic field strength, respectively, and ρ₀, u₀, and B₀ are constants. Conversion from the simulation unit to the cgs unit is, e.g., ρ_cgs = ρ_simρ_unit.

The scaling function we employ has a form

$$\begin{equation} \frac{1}{s(t)} = At^{-5/3}\exp {\left(-\frac{t_{\nu }}{t}\right)}, \end{equation} \tag{ 21 }$$

where A and t_ν are free parameters determined by a fit to the numerical evolution. The form comes from fitting one-dimensional relativistic viscous disk models (see J. C. Dolence et al. 2011, in preparation for more complete discussion). Note that without this time-dependent scaling procedure, or with a different scaling procedure, the spectra would vary systematically over the course of the simulation. The spectra would also differ systematically with resolution because the plasma β varies with resolution.

We fit for A and the viscous timescale t_ν from simulation data after a quasi-steady state has been reached, typically from t = 2000 onward. A sample fit to $\dot{M}$, for the 192 × 192 × 128 run, is shown in Figure 7. The variance of the normalized accretion rate decreases with resolution, that is, at higher resolution the fluctuations are smaller and Equation (21) gives an increasingly good fit. The maximum of the normalized accretion rate is nearly independent of resolution, when models with different resolutions are compared over the same time interval.

Zoom In Zoom Out Reset image size

Broadband, time-averaged synthetic spectra are shown in Figure 8. The mass unit of the torus is fixed by the condition that f_ν(230 GHz) = 3.4 Jy for an Sgr A* model measured at the solar circle. The shape of the spectrum is broadly similar at all resolutions for both polar-boundary conditions.

Zoom In Zoom Out Reset image size

Figure 9 shows flux density plotted against resolution in the infrared ($3.8\,\rm \mu \rm m$) and X-ray (integrated from 2 keV to 8 keV) where most of the emission is from direct synchrotron and single Compton scatterings, respectively. Some of the variation is likely due to run-to-run variation, as indicated by the error bars on the $N_{x_1} = 96$ and $N_{x_1} = 144$ models. The flux varies with resolution by less than about 50% at infrared and 30% at X-ray for $N_{x_1}>144$. The spectra therefore appear remarkably consistent and independent of resolution, at least for the M and $\dot{M}$ appropriate to Sgr A*.

Zoom In Zoom Out Reset image size

In a sense this is not surprising, because (1) our normalization procedure removes much of the variation that might arise from the decrease of β with resolution, and (2) the temperature is very well converged. The combined effect of the fixed flux normalization and the variation with resolution is to strengthen the magnetic field slightly and move the synchrotron peak slightly further into the infrared. This is echoed in the first Compton bump in the X-ray, which is forced to slightly higher energy by the increase in infrared input photons. While we have demonstrated this for only a single set of the model parameters (M, f_ν(230 GHz)), exploration of slightly different models with similarly consistent results shows that this is not a unique case.

We have investigated convergence of global GRMHD simulations of hot accretion flows onto a black hole and the emergent spectrum. We have run GRMHD simulations for four different resolutions, 96 × 96 × 64, 144 × 144 × 96, 192 × 192 × 128, and 384 × 384 × 256, in spherical polar coordinates. We have probed convergence using three diagnostics: time-averaged radial profiles of non-dimensional quantities (plasma β and electron temperature θ_e), azimuthal correlation lengths for several variables including the magnetic field, and artificial spectra generated with a Monte Carlo code.

For most of our diagnostics there are substantial differences between the lowest (96 × 96 × 64) and next lowest (144 × 144 × 96) resolution, and relatively minor changes at higher resolution. Run-to-run variations in the lower resolution models tend to be larger than the differences between the higher resolution (192 × 192 × 128 and 384 × 384 × 256) models.

We find that the magnetic correlation length is not converged. It decreases nearly linearly with resolution, with the number of grid cells per magnetic correlation length fixed at ∼5, although we do see a slight increase as resolution increases. Comparison with local model/shearing box simulations suggests that the turbulence does not change qualitatively at higher resolution. Such comparisons also suggest that another factor of ≈2 in linear resolution (costing about 1.6 × 10⁷ cpu hours) would resolve the azimuthal magnetic correlation length. None of the existing simulations (local or global) resolve scales more than a factor of ≈4 smaller than the correlation length (particularly the minor axis correlation length, which is oriented nearly along the radial unit vector and which we have not investigated here). If we identify the correlation length with the outer scale of MRI-driven turbulence, as seems reasonable, then none of these models have a resolved inertial range.

On the other hand, time-averaged synthetic spectra based on the GRMHD models, with parameters fixed to match Sgr A*, are remarkably reproducible from resolution to resolution. This suggests that simulated observations from existing simulations have some predictive power. We think it likely that the leading source of error in the high-resolution radiative models is now related to the underlying physical model (particularly the fluid model treatment of the plasma, and the absence of conduction) rather than the finite resolution of the models.

A similar convergence study has been conducted by HGK for non-relativistic global models. It is worth asking whether our models are converged according to the dimensionless resolution Q, the ratio of the most unstable MRI wavelength⁶ to the grid cell size in the azimuthal and vertical directions. In the azimuthal direction, ignoring relativistic corrections,

$$\begin{eqnarray} Q_3 = & \frac{\lambda _{{\rm MRI}}}{r \Delta \phi } \simeq 2 \pi \left(\frac{H}{r}\right) \frac{|B_3|}{c_s \sqrt{2 \rho }} \frac{1}{\Delta \phi } \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} \simeq & 2 \pi \left(\frac{H}{r}\right) \beta ^{-1/2} \frac{1}{\Delta \phi } \end{eqnarray} \tag{ 23 }$$

(Q_y or Q_ϕ in HGK's notation), where c_s ∼ HΩ is the sound speed. This gives Q₃ ≳ 22 and ≳ 10 for $N_{x_1}=384$ and 192, respectively, for all radii less than 10. In the vertical direction,

$$\begin{equation} Q_2 \simeq 2 \pi \left(\frac{H}{r}\right) \frac{|B_2|}{c_s \sqrt{2 \rho }} \frac{1}{\Delta \theta } = Q_3 \frac{|B_2|}{|B_3|} \frac{\Delta \phi}{\Delta \theta} \end{equation} \tag{ 24 }$$

(Q_z in HGK's notation), where Δθ is the zone size in Kerr–Schild coordinates at the midplane. Since |B₃/B₂| is usually ∼3–10 and Δϕ/Δθ = 4, this gives Q₂ ≳ 9–29 and 4–13 for $N_{x_1}=384$ and 192, respectively, for all r < 10. The required Q values to resolve the characteristic wavelength are Q₃ ≳ 6 (Sano et al. 2004) and Q₂ ≳ 20–60. Hence, MRI in the toroidal direction is resolved but not in the poloidal direction in these runs according to HGK's Q criterion.

We summarize our findings in the form of guidance for future simulators. (1) The resolution 96 × 96 × 64 is too low. The convergence measurements differ by factors of several from the highest resolution runs, and the magnetic field weakens steadily in a relative sense (β increases) over the course of the run. (2) The resolution 144 × 144 × 96 shows early signs of convergence except for the correlation length of the magnetic field; (3) the resolutions 192 × 192 × 128 and 384 × 384 × 256 differ relatively little from each other and show signs of convergence in the azimuthal correlation lengths, the temperature, and spectra, but not in the correlation length of magnetic field; (4) the observed trends with increasing resolution (to the extent that they are significant at the highest resolution) are that β decreases, θ_e increases, correlation lengths decrease, and IR and X-ray fluxes increase relative to millimeter fluxes, which we use to normalize the spectrum.

This work was supported by the National Science Foundation under grants PHY 02-05155 and AST 07-09246, by NASA under grant NNX10AD03G, through TeraGrid resources provided by NCSA and TACC, and by a Richard and Margaret Romano Professorial scholarship, a NESS fellowship, NNX10AL24H, to J.C.D., and a University Scholar appointment to C.F.G. Part of this work was completed while C.F.G. was a visitor at Max-Planck-Institut für Astrophysik, and C.F.G. thanks Henk Spruit and Rashid Sunyaev for their hospitality. The authors are grateful to Shane Davis and Xiaoyue Guan for providing unpublished data from their simulations.