GLOBAL SIMULATIONS OF ACCRETION DISKS. I. CONVERGENCE AND COMPARISONS WITH LOCAL MODELS

As turbulence involves fluctuation quantities in both space and time, diagnostic quantities will invariably involve some combination of spatial and temporal averaging. For simplicity, we define here the quantities we will use below. The most common quantity is simple volume integration for a variable X(R, ϕ, z, t) and is defined as 〈X〉(t) = ∫_VXdV, where V represents the full simulation domain. Related to this is volume averaging, defined as $\bar{X}(t) = \langle X \rangle /|V|$.

The fiducial timescale of these simulations is the orbital period at the inner edge of the disk, P_o = 2π, and for brevity is simply designated as an orbit without further qualification. Often we would like to compute scalar values of quantities that are representative of a simulation as a whole. In contrast to local simulations in which conserved quantities cannot change due to the closed nature of the simulation, global simulations involve significant evolution. The initial growth of certain physical quantities, particularly the stress and magnetic energy, is exponential and driven by the linear phase of the MRI. These quantities reach a peak value and quickly decrease as the linear phase of the MRI transitions into a fully turbulent state in which there is significantly reduced secular variation. We will refer to this stage as the quasi-steady state (QSS), which we define for simulations initialized with a vertical field as between 50 orbits and the end of the simulation. Due to the slower growth rate of the toroidal MRI, we quantify the QSS as between 100 orbits and the end of the simulation. Scalar representative values, like those in Table 1, are denoted as 〈X〉_QSS (or $\bar{X}_{{\rm QSS}}$) and are defined as the temporal average over the time frame defined as the QSS of the volume integral (or volume average) of the quantity. An alternative diagnostic to quantify disk evolution is the use of the measurement T_1/2, which refers to the number of orbits it takes for half of the total mass in the disk to leave the simulation domain. The values of this quantity for all simulations considered, save for simulation BzZ32W, which did not accrete half its mass in the time frame considered, are given in Table 1.

The dynamics of accretion disk turbulence are of particular astrophysical interest due to the anisotropic structure of the resultant turbulent magnetic field and its ability to drive angular momentum radially outward. To study the efficacy of angular momentum transport, many of our diagnostics will focus on the stress that allows this transport to happen. This stress is of two types: the Reynolds stress, R_Rϕ = ρv_Rδv_ϕ, and the Maxwell stress, M_Rϕ = −B_RB_ϕ, with the Maxwell stress generally dominating the Reynolds component. While these two quantities are of the most direct physical relevance, it is common to scale the stress by the gas pressure, resulting in the following diagnostics:

$$\begin{eqnarray} &&\alpha = \frac{\displaystyle {R_{R\phi } + M_{R\phi }}}{\displaystyle {P}},\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray} &&\alpha _{M} = \frac{\displaystyle {M_{R\phi }}}{\displaystyle {P}}. \end{eqnarray} \tag{ 5b }$$

Also of interest is the strength of the magnetic field in relation to the thermal energy of the gas, given by β = P/P_b, the ratio of gas and magnetic pressure. Although an abuse of notation, for brevity we will use the notation 〈α〉 and 〈β〉 to refer to the volume integral of the numerator divided by the volume integral of the denominator.

In an effort to study how well resolved the MRI is, we consider the following diagnostics:

$$\begin{eqnarray} &&F_{z} = |V|^{-1} \int _{V} (\lambda _{{\rm MRI}} \ge 8\Delta z)dV, \end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray} &&F_{\phi } = |V|^{-1} \int _{V} (\lambda _{C} \ge 8 R \Delta \phi) dV, \end{eqnarray} \tag{ 6b }$$

where in the above equations logical statements refer to indicator functions that assume the value of unity or zero based on the truth or falsity of the statement. The characteristic wavelength of the toroidal field, λ_C, is defined analogously to λ_MRI save for the use of the toroidal Alfvén speed in place of the vertical. These scaled integrals represent the fraction of the disk where the fastest-growing modes of the vertical and critical toroidal modes of the MRI are resolvable, using the often employed eight-zone criterion. These are related to the quality factors, first used by Noble et al. (2010), given as

$$\begin{eqnarray} &&Q_{z} = \lambda _{{\rm MRI}} / \Delta z, \end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray} &&Q_{\phi } = \lambda _{c} / (R \Delta \phi), \end{eqnarray} \tag{ 7b }$$

where these values are in general functions of space and time. Use of the quality factor as a diagnostic represents an important step in appreciating the numerical resolvability of the MRI, but we feel that the resolvability fraction diagnostic we employ is an even more stringent resolvability requirement.

It is well known that the anisotropy of the magnetic field is key to angular momentum transport, and indeed the correlation between the radial and azimuthal component of the magnetic field is an important diagnostic of angular momentum transport in disks. An alternative measure of the anisotropy, called the magnetic tilt angle, was first discussed by Guan et al. (2009). This measure is defined as

$$\begin{equation} \theta _{B} = \arcsin (\alpha _{M}\beta)/2. \end{equation} \tag{ 8 }$$

Physically, this can be thought of as an approximation to the angle between the planar magnetic field and the azimuthal axis. It should be pointed out that the equivalence between the quantity defined in Equation (8) and the physical interpretation relies on the assumption that the contribution of the vertical field to the total magnetic pressure is negligible. The growth of the MRI is characterized by strong toroidal field amplification, and thus this assumption will generally be realized; however, for the simulations initialized with vertical magnetic field this assumption will not be satisfied at very early times before the initial growth phase of the MRI. An estimate of the tilt angle is derived as θ_B ≈ 15°, in Guan et al. (2009), by noting that all the local models considered satisfy the relationship αβ ≈ 1/2. Table 1 shows for the various global simulations presented here the value of 〈α〉_QSS〈β〉_QSS, which is found to be in good agreement with the expected scaling αβ ≈ 1/2. Strictly speaking, generally α_M < α, and thus we expect this estimate to be an upper bound.

Studying the spectral structure of turbulent flows is often more natural than alternative diagnostics in physical space. Our primary interest here will be studying the azimuthal structure of power spectra of physical quantities, in particular the density and magnetic pressure. To this end, we begin by defining the subdomain

$$\begin{equation} \mathcal {S} = [R_{0}-2H_{0},R_{0}+2H_{0}] \times [0,2\pi ] \times [-2H_{0},2H_{0}]. \end{equation} \tag{ 9 }$$

For simulations that do not model the full 2π radian azimuthal domain the subdomain definition is modified to include the full azimuthal range modeled. We then consider, for a physical quantity X(R, ϕ, z), the azimuthal Fourier decomposition X(R, m, z). Further, we define X(m) as the volume-weighted radial and vertical mean over the subdomain $\mathcal {S}$. This removes spatial transients, but to remove temporal transients we also average between orbits 50 and 100, for simulations seeded with a vertical field, and between orbits 100 and 150, for simulations seeded with a toroidal field. We refer to this reduced temporal domain as RQSS. To ensure that these results are not adversely affected by secular trends within the subdomain related to the evolution of the disk, we scale, at each time step, by $\langle X \rangle ^{\mathcal {S}}$, the volume integral of the quantity over the subdomain. Formally,

$$\begin{equation} \hat{X}(m) = \left< \frac{\displaystyle {X(m)}}{\displaystyle {\langle X\rangle ^{\mathcal {S}}}} \right>_{{\rm RQSS}}. \end{equation} \tag{ 10 }$$

To more directly study important azimuthal scales, we will employ the transformation of the wavenumber, m, to an effective azimuthal wavevector given by

$$\begin{equation} k_{\phi } = \frac{\displaystyle {m}}{\displaystyle {2\pi R_{0}}}. \end{equation} \tag{ 11 }$$

One final tool we employ when studying the spectral structure of disk turbulence is the use of a fiducial power spectrum. To study the relative importance of large-scale versus small-scale features, we consider a fiducial power spectrum given by P_m = m⁻¹. This fiducial spectrum has the property that the inner and outer scales are equally important to the total integrated power over the accessible wavenumbers. We diagnose the dominant azimuthal scale by considering the location of the peak value of the quantity defined in Equation (10) scaled by this fiducial power spectrum, namely, $m|\hat{X}(m)|^{2}$. This approach has the added benefit of allowing a simpler visual interpretation of the power spectrum when the horizontal axis is logarithmic.

The above formalism will be employed to study the structure of density and magnetic energy, and while it is common to study the structure of stress in a similar manner, we choose a slight variant. Astrophysical interest in stress as a diagnostic is based on its ability to drive radial angular momentum transport. We note, however, that when projecting this quantity onto the Fourier basis all but the m = 0 mode will correspond to an azimuthal average of zero. Instead, we choose to analyze the structure of contributions to stress that result in net radial transport. Formally, we define

$$\begin{equation} \tilde{\alpha _{M}}(m) = \left< \frac{\displaystyle {-B_{R}(m)B_{\phi }(m)}}{\displaystyle {\langle P\rangle ^{\mathcal {S}}}} \right>_{{\rm RQSS}}. \end{equation} \tag{ 12 }$$

One of the goals of this paper will be to explicitly compare local and global models; toward this end we proceed in a similar fashion to Sorathia et al. (2010) and decompose our global simulation into a set of subvolumes through which we can calculate "local" statistics. This is accomplished by considering the subvolume of the physical domain $(R,\phi ,z) \in \mathcal {S}$ and decomposing it into wedges of size [H₀, 2πH₀, H₀]. This yields, at each time step, a set of 160 subvolumes in which relevant physical quantities can be volume-averaged. To calculate statistics, we simply take the mean, denoted as [X](t), and when relevant the standard deviation of these quantities at each time step.

We also consider the instantaneous correlation of flux and stress both as a convergence criterion and as a means of comparing local and global models. The flux–stress connection was first explored by Hawley et al. (1995) in the context of saturation predictors and further refined to take into account numerical effects by Pessah et al. (2007). Saturation predictors map the initial flux threading a local model to the predicted stress after the initialization and saturation of MRI-driven turbulence. It is noted in Hawley et al. (1995) that saturated stress increases with the strength of the vertical field threading a simulation domain, specifically α∝λ_MRI. This observation is augmented in Pessah et al. (2007) to include the numerical constraint that a sufficiently weak net field will be unresolvable and thus behave in a manner identical to a zero-net field. Formally, the saturation predictor presented in Pessah et al. (2007) is given by

$$\begin{equation} \alpha _{M}\Big (\frac{\displaystyle {H}}{\displaystyle {L}}\Big)^{(5/3)} = 0.61 \times \left\lbrace \begin{array}{@{}ll}\Delta / L & : \lambda _{{\rm MRI}} \le \Delta ,\\ \lambda _{{\rm MRI}}/L & : \Delta < \lambda _{{\rm MRI}} \le L, \\ 0 & : L < \lambda _{{\rm MRI}}, \end{array} \right. \end{equation} \tag{ 13 }$$

where L and Δ are the box size and grid scale, respectively. Their saturation predictor includes three regions: an unresolved region in which λ_MRI is unresolvable and consequently α∝Δ, the resolved region studied in Hawley et al. (1995), and a stable region in which λ_MRI exceeds the vertical domain of the simulation and turbulence is absent.

To study the global analog to these local saturation predictors, we proceed in a manner similar to Sorathia et al. (2010). We calculate the instantaneous flux and stress in each subvolume at each time step after orbit 50 to remove possible noise in the data from the linear growth rate of the MRI. The resulting flux–stress pairs are logarithmically binned according to flux in order to diagnose trends. In contrast to the local model and Sorathia et al. (2010), where flux–stress pairs are scaled by L and (H/L)^5/3, respectively, we proceed in a manner more appropriate to global simulations in which the size of the subdomain has no physical meaning. The flux–stress pairs calculated here are of the form (λ_MRI/Δz,α_M), as in Beckwith et al. (2011). Scaling the local flux to the grid scale is more directly meaningful and facilitates a closer inspection of the transition point discussed in Sorathia et al. (2010).

In addition to studying the relationship between vertical flux and stress, we also consider the analogous connection in the presence of toroidal flux. Arguably, the presence of strong net toroidal field is more astrophysically relevant (van Ballegooijen 1989). The evolution of the toroidal MRI is more complex than its vertical counterpart, and there is no simple form for its most unstable mode. We proceed, as in Hawley et al. (1995), by considering the toroidal flux defined by $\lambda _{c} = 2\pi \sqrt{16/15} \:v_{A,\phi }/\Omega$. The saturation predictor there is given as

$$\begin{equation} \alpha _{M} \propto \frac{\displaystyle {L_{y}}}{\displaystyle {H^{2}}} \lambda _{c}, \end{equation} \tag{ 14 }$$

where L_y is the effective azimuthal length of the shearing box. To test this in the global context, we consider flux–stress pairs of the form (λ_c/RΔϕ, α_M) and perform the same logarithmic binning as described for the vertical flux–stress pairs.

The simplest and most astrophysically relevant convergence criterion is accretion efficacy as measured by the dimensionless stress α_M. Figure 5(a) shows the evolution of the Maxwell stress over the course of the simulation for all three standard ZF models. The initial peaks in the stress associated with the linear growth of the vertical MRI are resolution dependent. While the initial magnetic field is constructed so that the most unstable mode, λ_MRI, is resolvable in each simulation, the largest growing modes of the non-axisymmetric MRI are associated with large k_z (Balbus & Hawley 1992) and thus depend on resolution. Thus, initial peaks in the stress associated with the early growth of the MRI are resolution dependent. The true test of convergence is the behavior of the stress in the saturated QSS. The lowest resolution simulation, BzZ8, exhibits 〈α_M〉_QSS ≈ 0.01, while both higher resolution simulations exhibit a comparable stress to each other that is roughly 50% greater than BzZ8.

Zoom In Zoom Out Reset image size

While this fundamental criterion of convergence in stress is satisfied, other diagnostics paint a more subtle picture. Figure 5(b) illustrates the evolution of the scaled magnetic energy, β⁻¹. Again we see that the initial field amplification is monotonic with resolution and dominated by the growth of the toroidal magnetic field; however, the saturation and transition into the fully nonlinear state are more complex than expected. During the linear growth phase of the MRI simulation, BzZ32 exhibits the largest scaled magnetic energy; however, during the QSS it is the lower resolution simulations, BzZ8 and BzZ16, that maintain a stronger magnetic field. These lower resolution simulations also lack the steep drop-off in magnetic field energy often associated with the saturation of the MRI. The exact nature of this resolution dependence is unclear and likely depends on the details of the saturation of the MRI and transition to nonlinearity. Also of note is the late-time field growth associated with BzZ8. The nature of this growth is unclear but may be suggestive of a very low frequency temporal behavior.

Our analysis of the physical metrics of models initialized with a seed field possessing net flux proceeds in much the same way as our analysis of the ZF models. As above, we begin by considering the time evolution of the global quantities, α_M and β⁻¹. The results are shown in Figures 5(c) and (d), respectively. The evolution of the linear MRI displays a significant resolution dependence, as would be expected, but in all cases the values of α_M in the saturated state are comparable for the same initial field topologies. Due to the significantly higher values of α_M in the net field simulations, we see a corresponding significant increase in mass loss of the disk (Table 1) compared to the ZF runs. Field amplification for the runs initialized with a vertical field is monotonic with resolution. However, regarding field amplification the net toroidal runs behave more analogously to the net-zero field runs discussed above in which lower-resolution simulations seem to exhibit greater field amplification in the QSS. The presence of a net field, regardless of topology, results in order-of-magnitude increases in the peak values of α_M and β⁻¹ compared to the ZF simulations considered above. Insofar as angular momentum transport is concerned, there is only a weak dependence on resolution in the saturated state even for the most poorly resolved simulations. The evidence from these global metrics suggests that net-flux simulations converge more quickly; however, we will see through the consideration of the other metrics that the picture is more subtle.

Overall, we see that while the behavior of α_M paints a simple picture and is suggestive of convergence, the magnetic energy is far more volatile and does not follow a clear pattern in its dependence on resolution. Further, as a convergence criterion α_M suffers from its dependence on initial field topology.

The use of orbital advection allows us a unique opportunity of exploring isotropic resolutions in a cost-effective manner. Without orbital advection, the time step constraint is set by v_K/Δ_ϕ, and thus doubling azimuthal resolution results in a steep price, specifically a quadrupling of the necessary operations from the doubled number of cells to time advance and the half-time-step being used. The result is that in the majority of the literature, azimuthal resolution is compromised for computational expediency. The goal of runs BzZ32R and BzZ32RR is to study the importance of azimuthal resolution in a controlled way in comparison to run BzZ32 and deduce a maximum aspect ratio from which converged turbulence is guaranteed. Figures 6(a) and (b) show the evolution of α_M and the scaled magnetic energy for these runs. We note that in both cases we see non-monotonic behavior; specifically reducing the resolution by a factor of two results in increased accretion efficacy and toroidal field amplification. Further reduction results in a significant drop-off in both of these quantities. While increased field amplification with reduced resolution is expected from Figure 5(b), that a reduction in resolution could increase angular momentum transport is unexpected. Assessing convergence from these global quantities is murky, at best. For now, we merely note the importance of azimuthal resolution by pointing out that reducing azimuthal resolution can severely impact accretion efficiency. The convergence of these reduced azimuthal runs is returned to in Section 4.4, where it is demonstrated that run BzZ32R is converged, whereas BzZ32RR is not.

Zoom In Zoom Out Reset image size

Next, we consider a comparison of the resolvability fractions defined in Equation (6) and shown in Figures 7(a) and (b). The initial radial profile of the vertical seed field is constructed so that λ_MRI/Δz ⩾ 8 for simulation BzZ8 in a small neighborhood around each maximal value of the sinusoid. The variations in the initial values of F_z are due to the size of the neighborhood about each maximal value of the initial sinusoid in which the resolvability criterion is satisfied. Of all the diagnostics considered, the most startling behavior is seen when considering F_z. The two lower resolution simulations, BzZ8 and BzZ16, both show a significant initial drop in the resolvability of the vertical MRI correlated with the linear growth and saturation phase of the evolution. Following the transition, both of these simulations show a slight increase in the resolvability fraction, but clearly in most of the disk the vertical MRI is not adequately resolved. In contrast to this, BzZ32 seems to exhibit substantially differing behavior and shows a value of F_z that is roughly constant during the full evolution of the simulation. The evolution of F_ϕ, shown in Figure 7(b), is by comparison much simpler and monotonic in nature. As the toroidal field grows due to shear amplification driven by the vertical MRI, the toroidal MRI becomes resolvable throughout a significant portion of the disk.

Zoom In Zoom Out Reset image size

The behavior of the resolvability fractions paints an interesting picture. Runs BzZ16 and BzZ32 exhibit similar saturated values of stress, but their ability to resolve the vertical MRI seems to be significantly different. This suggests the intriguing possibility that run BzZ16, while initially seeded with a vertical field, is actually reliant on the resolvability of the toroidal field to reach a comparable stress to run BzZ32. However, if these two simulations are actually taking differing routes to turbulence, there then must be some additional mechanism that accounts for the similar values of stress achieved in the QSS. The resolvability fractions suggest that run BzZ32 has reached a categorically different state in its resolvability of the vertical MRI, and the toroidal MRI is resolvable throughout almost all of the disk.

Next we consider the resolvability fractions of the NF runs, given in Figures 7(c) and (d). The general resolvability of the vertical MRI is, at first glance, generally somewhat poor, as demonstrated in Figure 7(c). Of all the simulations considered, only BpN32 resolves the vertical MRI in a majority of the disk over the full evolution. As was the case with the ZF runs considered, all of the simulations resolve the toroidal MRI in the majority of the disk (Figure 7(d)). From these resolvability fractions, only BpN32 is clearly well resolved. Of interest is that the resolvability fraction in the saturated state appears to have stronger dependence on resolution than initial field topology.

The resolvability fractions employ the often-used eight-zone resolvability criterion; however, this is somewhat arbitrary. They allow us to get a sense of how much of the disk is "well" resolved. Now we employ the quality factors to understand how "well" resolved the disk is. Figures 8(a) and (b) illustrate the evolution of the vertical and toroidal quality factors, respectively, for the ZF runs. The vertical quality factors roughly double with corresponding resolution doublings, starting at Q_z ≈ 2 for the lowest resolution simulation. This results in a situation in which run BzZ16 is crudely resolved with an average of four zones per λ_MRI and BzZ32 is demonstrably resolved. Consideration of the toroidal quality factor results in a reiteration of the point made regarding the toroidal resolvability fraction, F_ϕ, specifically that the two highest resolution simulations clearly are resolving the toroidal MRI. In particular, run BzZ32 has an average of 40 zones per critical wavelength.

Zoom In Zoom Out Reset image size

Consideration of the quality factors for the NF runs (Figures 8(c) and (d)) reveals a similar resolution dependence as that seen in the resolvability fractions. With the exception of BpN32, all of the net-flux runs evolve to a point with an average of a vertical quality factor less than eight, the fiducial criterion. The toroidal quality factors, in contrast, are quite large, and the simulations with all but the lowest resolution simulations exhibit toroidal quality factors above 20.

Overall, consideration of the resolvability fractions and quality factors is taken as evidence that runs BzZ32 and BpN32 are resolved whereas BzZ16 is barely resolved.

While integrated global quantities are undoubtedly important when diagnosing turbulence, it is often the case that the spectral structure of turbulence is more amenable to study than is the physical structure. Convergence in this domain is again inherently ill-defined due to the non-dissipative nature of these simulations. Increasing the resolution of a simulation increases the number of modes in which power can reside, and it would be unphysical to expect that these newly opened modes would remain free of power. We adopt as a definition of convergence that the mode at which power peaks remains constant with increasing resolution. To probe the spectral structure at the smallest scales, we include in consideration the wedge runs, BzZ32W and BzZ64W, where the function of the former is primarily as a control to ensure that the small-scale behavior is not altered by the reduction of azimuthal domain.

Figure 9 shows the time-averaged azimuthal power spectra of the density, magnetic pressure, and stress scaled to remove the secular evolution and temporally averaged for 50 orbits beginning at orbit 50. It is clear from these figures that reducing the azimuthal domain of the simulation, as done in BzZ32W, does not significantly change the small-scale distribution of power. The importance of large-scale structure to the overall density profile of the disk is clear from Figure 9(a). Comparing the structure of the magnetic pressure and stress, Figures 9(b) and (c), respectively, suggest that the magnetic pressure is clearly dominated by intermediate scales while the stress is more equally distributed at large and intermediate scales. From a visual inspection, it is clear that all of the power spectra are self-similar and, for resolutions above 32 zones per scale height, all peak at approximately the same scale. At small scales the power spectra look quite like those of the non-converged, zero-net-flux models presented by Fromang & Papaloizou (2007). The large-scale behavior, however, is quite different as the power at large scales is at most weakly dependent on resolution, in stark contrast to the non-converged models of Fromang & Papaloizou (2007).

Zoom In Zoom Out Reset image size

Next, we consider the azimuthal spectral structure of the turbulence for the NF runs. As for the ZF runs, we consider the density (Figure 10(a)), the magnetic pressure (Figure 10(b)), and the stress (Figure 10(c)). Broadly, we note again that the structure is dependent largely on resolution and only minimally on field topology. Simulations of lower resolution are more dependent on large-scale features, and as resolution increases we see a more even distribution of power between the intermediate and small scales. While the density and pressure have a clear peak at an intermediate scale, the scale distribution of stress is more evenly distributed at large to intermediate scales.

Zoom In Zoom Out Reset image size

To quantify the dominant scale, we consider in Table 2 the wavenumber at which the wavenumber-scaled power peaks for each simulation and for each of the variables considered. For the eventual comparison of these results with local models, it is useful to consider this measurement both as wavenumber, appropriate for global simulations, and as an effective azimuthal wavevector, k_ϕH₀, which is more directly analogous to local models. The dominant density scale is much larger than the corresponding dominant magnetic scales. All of the simulations considered exhibit a dominant density scale of approximately 4H₀, and in this regard resolution does not appear to play a significant role save for the absolute lowest resolution simulation. Reducing the azimuthal domain results in an increase of the effective scale associated with the magnetic quantities. While the magnetic quantities do not clearly converge with increasing resolution, the peak of the highest resolution simulation, BzZ64W, is almost precisely displaced in wavenumber by eight. The reduced azimuthal domain simulations use only an eighth of the full 2π domain, which means that the accessible wavenumbers are limited to modes that are integer multiples of eight. In this regard, the peak of runs BzZ32W and BzZ64W is within one accessible wavenumber of the peak of BzZ32. With this in mind, we again take this as evidence of the convergence of BzZ32.

The dominant density scale is roughly halved when going from run BpN8 to BpN16; however, the transition from BpN16 to BpN32 results in only a minor (∼ 10%) decrease and suggests that the latter run is near convergence. This is contrasted with the dominant scales of the magnetic energy and stress, which decrease much more significantly in the transition from BpN16 to BpN32. While the dominant density scale of run BzZ32 is twice that of the dominant stress scale, the toroidal field runs display a much more comparable value of these two scales. Overall, we take this evidence as suggesting that run BpN32 is approaching convergence but not fully so.

The most striking feature of the net vertical field simulations is the significantly larger dominant stress scale (∼ 10H₀). Also, larger than expected is the dominant magnetic energy scale, which is again in excess of the other topologies considered. These large scales may suggest the existence of a memory of the initial field configuration. In particular, this large structure may be a remnant of the linear shearing phase in the initial growth of the MRI. This linear shearing phase is particularly strong in the simulations initialized with a net vertical field and exhibits magnetic pressures in excess of the gas pressure.

The metrics discussed thus far are useful tools to measure convergence, but they suffer from the limitation that their behavior must be compared against other simulations using identical field topologies. Meaningful convergence studies can be computationally quite expensive, and therefore a more robust indicator of convergence that is independent of field topology would be useful. The evidence we will present here suggests that the magnetic tilt angle, defined in Equation (8), may indeed be this diagnostic. As a precursor to this, we consider the evolution of the tilt angle in all of the ZF simulations (Figure 11(a)). The magnetic tilt angle reaches an approximately steady-state value that is almost identical for all of the higher resolution simulations (above 32 zones/H₀). This value, θ_B ≈ 13°, is remarkably close to the estimated value of 15° (Guan et al. 2009). Also of note is that while the initial peaks in both stress and magnetic energy are strongly resolution dependent, the tilt angle, a function of the product of the two, has an initial peak that is seemingly independent of resolution. This may suggest that the transition from the linear growth of the MRI into saturated turbulence depends on a precise relationship between the quantities α_M and β.

Zoom In Zoom Out Reset image size

Next, we consider the behavior of the magnetic tilt angle in the NF runs, given in Figure 11(b). As resolution increases, we see a corresponding increase in the tilt angle. Regarding the toroidal field simulations, we note the strong similarity between runs BpN16 and BpN32. This suggests that a further doubling in resolution may indeed prove convergence. The vertical field simulations display a similar behavior and ambiguity. Unfortunately, with the simulations available we are not able to conclusively demonstrate convergence for these net field runs; however, we do argue that BpN32 is likely converged.

At this point, we recall the difficulty of defining convergence using the physical metrics in the context of the reduced azimuthal resolution runs, BzZ32R and BzZ32RR. The physical metrics are ambiguous, but when we consider the tilt angle (Figure 11(c)), we see a clear monotonicity with resolution. A reduction in the azimuthal resolution by a factor of two results in a very minor change in the tilt angle; however, further reduction of the azimuthal resolution results in a much more significant alteration. We take this as evidence that run BzZ32R is converged, whereas BzZ32RR is not. This demonstrates the importance of azimuthal resolution and indeed suggests that treating azimuthal resolution on nearly equal footing with vertical resolution is of vital importance toward ensuring simulations that are numerically converged.

The results we have presented indicate that the magnetic tilt angle is a powerful diagnostic tool toward demonstrating convergence of MRI-driven accretion disk simulations. The fact that it is monotonic with resolution, exhibits minimal variation in the saturated state, and appears to converge suggests that fiducial values of converged tilt angle can be computed and compared against simulations. Additionally, as demonstrated in Figure 11(d), the value in the saturated state appears to be almost independent of initial field topology, and there may exist a single fiducial value of tilt angle. We conjecture that higher resolution simulations initialized with a net field will indeed definitively demonstrate the existence of a single scalar value of the magnetic tilt angle that characterizes accretion disk turbulence. Assuming that this can be done, this would prove a powerful tool toward verifying convergence without resorting to the significant expense of running higher resolution control simulations. Further, the existence of a single fiducial tilt angle would imply the ability to unify the often-disparate phenomenology associated with varying initial field topologies.

The natural question raised by the existence of a single scalar that appears to characterize the saturated, fully nonlinear state of MRI-driven MHD turbulence is whether this can be derived from theory. The relationship, αβ ≈ 1/2, was first observed by Hawley et al. (1995) and discussed at length by Blackman et al. (2008). The latter work includes a heuristic argument that α and β should be related by a constant; however, predicting that constant is beyond the means of the dimensional analysis used. Recent work by Pessah (2010) provides the potential for a more quantitative understanding of this relationship. Through an analysis of the growth and subsequent saturation of a single vertical MRI mode by parasitic instabilities, as initially studied by Goodman & Xu (1994), they find that when Kelvin–Helmholtz parasitic modes dominate, the saturation of the primary MRI mode occurs when αβ ≈ 0.4 (θ_B ≈ 12°). The regime in which Kelvin–Helmholtz parasites dominate is when the Elsasser number, Λ_η = v²_A/ηΩ_K, is larger than unity, and it is this regime that we expect our ideal simulations to correspond to. This work is encouraging as it represents an important point of connection between our numerical results and theory. The analysis of Pessah (2010) finds that the magnetic tilt angle will depend on dissipative coefficients, and finding agreement between this dependence and future numerical simulations would bolster the importance of this quantity. Additionally, the fact that simulations initialized with a toroidal field are characterized by the same tilt angle as simulations initialized with a vertical field may suggest that the formation and subsequent parasitic destabilization of channel solutions may play an important role in the saturation of the non-axisymmetric MRI.

Hawley et al. (2011) survey a series of local simulations to test various proposed convergence metrics and then compare these metrics with a series of pseudo-Newtonian simulations of thick accretion disks. Our comparison to their work begins with their choice of metrics: α_Mag = M_Rϕ/P_B; the quality factors, Q_z and Q_ϕ (though the definitions given there differ from ours by a numerically negligible factor of $\sqrt{16/15}$); and the correlations, 〈B²_R/B_ϕ²〉 and 〈B²_z/B_R²〉. We note that the measure α_Mag used in Hawley et al. (2011) is not the same as our metric α_M; while both involve a pressure scaling of the magnetic stress, the former uses the magnetic pressure whereas the latter uses the gas pressure. We also note that two of the metrics proposed are in an information content sense equivalent to the magnetic tilt angle, θ_B. Due to the nature of the magnetic tilt angle being a measure of anisotropy in the planar field, we can rewrite α_Mag = α_Mβ, which implies sin (2θ_B) = α_Mag. Similarly, the correlation B²_R/B_ϕ² = tan ²(θ_B). Regarding the latter correlation, B²_z/B_R², the authors note that this term does not appear to exhibit a strong trend with increasing resolution, and as such we will not consider it further. To aid in our comparison, we include a summary of convergence metrics in Table 3.

Table 3. Convergence Metrics

Simulation	〈θB〉QSS	$\bar{Q_{z}}_{,{\rm QSS}}$	$\bar{Q_{\phi }}_{,{\rm QSS}}$
BzZ8	861	1.68	8.38
BzZ16	1144	4.57	16.82
BzZ32	1286	9.91	41.22
BzZ32W	1275	9.98	33.95
BzZ64W	1276	24.51	73.77
BzZ32RR	1112	6.64	7.38
BzZ32R	1259	9.08	16.13
BpN8	808	4.46	18.86
BpN16	1098	9.23	31.92
BpN32	1209	22.2	64.22
BzN8	937	3.54	16.24
BzN16	1183	8.97	32.41

Download table as:  ASCIITypeset image

For the local models reviewed in Hawley et al. (2011), simulations of 64 zones per scale height result in quality factors of Q_z ≈ 10 and Q_ϕ ≈ 40. These are fairly comparable to each other and run BzZ32 (Table 3), despite the disparate initial field topologies. A zero-net vertical field is used in Davis et al. (2010), whereas a net toroidal field is used in Simon et al. (2011). Of note is that the net toroidal field model considered here, BpN32, shows significantly higher quality factors at lower resolutions. While Davis et al. (2010) also consider a model utilizing 128 zones per scale height, which results in an increase in the quality factors by approximately 150%, this does not result in a change in α_Mag = 0.36 from the 64-zone-per-scale-height run. This value of α_Mag corresponds to θ_B ≈ 105. The simulations discussed in Simon et al. (2011) result in a largest value of tilt angle of θ_B ≈ 118. These are both in contrast to the values seen in our global runs with a slightly larger value of θ_B ≈ 12°–13°.

In their discussion of stratified local models, Hawley et al. (2011) also note the importance of azimuthal resolution, which is bolstered by the work presented here. In Section 4.4, we present evidence that when the azimuthal resolution is reduced by a factor of four from an aspect ratio of unity, the resultant simulations exhibit significantly lower accretion and magnetic tilt angles. Also discussed is the ability for large toroidal quality factors to compensate for a poorly resolved vertical MRI. Indeed, this likely explains the results for simulation BzZ16, in which a comparable value of α in the steady state is achieved despite the significant discrepancy between the resolvability fractions of simulations BzZ16 and BzZ32.

The global simulations presented by Hawley et al. (2011) are diagnosed to have lower quality factors and smaller tilt angles than the ones presented here and the local models they consider. However, we observe that the stratified local models also seem to involve a smaller tilt angle than what may be expected from our unstratified global models. Beckwith et al. (2011) measure a tilt angle of approximately 9° in a stratified global simulation. This suggests that stratification itself may suppress tilt angle somewhat, and a future resolution study of stratified global disks using orbital advection would be useful. Overall, the work suggests that Q_z ≳ 10–15 for poorly resolved toroidal quality factors (Q_ϕ ≈ 10) and that larger values of toroidal quality factor (Q_ϕ ≳ 25) can alleviate the constraint on the vertical quality factor. Based on these constraints, we find that all of the simulations utilizing a resolution exceeding or equal to 32 zones per scale height are converged as well as the simulations seeded with a net field and above 16 zones per scale height. The criterion presented here of a converged tilt angle would not consider the runs BzN16 and BpN16 converged, but would consider run BzZ32R converged in contrast to the quality factor criterion.

The benefit of treating a global simulation as a local ensemble is that it allows us to study the local dynamics of accretion disk turbulence over a large range of parameter space simultaneously. Additionally, it offers the opportunity to directly test the importance of curvature terms and the degree to which local flux is truly conserved. Local statistics can be used to study not only the average value of a quantity but also the variation in the quantity and how it evolves. Figure 13 illustrates the time evolution of [α_M] and [β⁻¹] for a selection of ZF and NF simulations. Comparing [α_M] and [β⁻¹] in the ZF simulations, Figures 13(a) and (b) illustrate the resolution dependence analogous to the physical metrics considered above; however, the variability of these quantities is inversely related to the resolution. The behavior of run BzZ8 is statistically separated from the behavior of the higher-resolution simulations, with runs BzZ16 and BzZ32 demonstrating more comparable values of these quantities.

Zoom In Zoom Out Reset image size

The behavior and evolution of the accretion efficiency and scaled magnetic energy for a selection of NF runs are given in Figures 13(c) and (d). Most striking is the extreme initial transient associated with run BzN16, characterized by values of α_M well in excess of unity and strongly magnetized regions with Alfvén speeds several times the sound speed. This transient, however, is short-lived, with the long-term evolution comparable to that of runs BpN16 and BpN32. The toroidal runs, like the ZF simulations, display an inverse relationship between variability and resolution. The resolution dependence of the variability is itself interesting. While the global treatments of these quantities discussed above can demonstrate their convergence, that the variability weakens with resolution may suggest a separate resolution requirement.

In Sorathia et al. (2010), we make a preliminary effort toward connecting local and global models of accretion disks using the method described in Section 2.1 to generate "local" statistics. The simulations considered were useful but suffered due to the added difficulty of separating effects of stratification and resolvability from the desired comparison. However, the fact that the resolvability of the linear MRI in saturated turbulence is important or sufficient is not clear from analytic theory. The simulations considered in this work are designed from first principles to be more directly applicable to the comparisons we would like to make. We have previously considered the instantaneous correlation between vertical flux and stress; here we also study the relationship between toroidal flux and stress. That there is a structure to the relationship between toroidal flux and stress is interesting (Figures 14(b) and (d)). At first glance, one might not expect a transition in the correlation, since the comparison is between toroidal field and stress, a term dominated by the toroidal field. The nature of this correlation is unclear, as in contrast to the vertical MRI, which grows quite rapidly and would be expected to correlate to the presence of field, it is not clear that the toroidal MRI would be associated with timescales small enough to correlate to the presence of field. Assuming that the connection between toroidal flux and stress is causal, as is believed to be the case in the vertical flux–stress correlation, it may suggest that toroidal flux generates stress in the nonlinear regime in a manner much faster than it would in the linear. An alternative possibility is that the correlation between toroidal flux and stress is a consequence of the vertical flux–stress relationship. In this case, the presence of vertical flux would drive stress with toroidal flux amplification as an intermediate step. This possibility would suggest that the relationship between the turning points in the two flux–stress relationships is meaningful. We note that while the fluxes are scaled to the grid resolution, these are not the same as the quality factors presented earlier as the flux is averaged over each wedge as opposed to taken on a cell-by-cell basis. This averaging reduces the strength of the net field, which results in lower values than the quality factors.

Zoom In Zoom Out Reset image size

The calculated flux–stress relationship for net-zero fields is given for both vertical flux (Figure 14(a)) and toroidal flux (Figure 14(b)). It was noted in Sorathia et al. (2010) that the transition between the flat unresolved region and the linear resolved region suggested a seeming super-resolvability, specifically that we saw this transition at λ_MRI ≈ Δz/20. We note similar behavior in run BzZ8 in Figure 14(a), in which the transition occurs at λ_MRI ≈ Δz/10. The remaining simulations exhibit a transition at λ_MRI ≈ Δz, in agreement with the scaling law given by Equation (13) (Pessah et al. 2007) and the simulations of Beckwith et al. (2011). Related to this, we see a qualitatively similar resolvability threshold in the toroidal flux–stress relationship (Figure 14(b)). While the two highest resolution simulations exhibit a transition at approximately 10RΔϕ, the transition for BzZ8 occurs at approximately half of that. In Sorathia et al. (2010), this super-resolvability was attributed to the effect of slower-growing unstable MRI modes, and while this still may explain the placement of the transition point in the higher resolution simulations, it is likely not the reason for the behavior of run BzZ8. The cause of this discrepancy between BzZ8 and the other simulations is unclear; however, it is clear that this super-resolvability of the transition seems to be a characteristic of unresolved turbulence.

Next, we consider the flux–stress relationship for the net-field simulations for both vertical flux (Figure 14(c)) and toroidal flux (Figure 14(d)). We note that in general the net-flux runs exhibit a general structure very similar to the net-zero flux runs. Specifically, the transition points between the flat and linear stress-response regimes are independent of initial field topology. None of the net-field simulations considered exhibit the super-resolvability seen in BzZ8, which suggests that even low-resolution net-field runs are more resolved than their zero-flux counterparts. The general trend is that, with increasing resolution, the stress response to unresolved flux is generally increased and the slope of the linear regime, in which flux is resolved, becomes shallower. An interesting feature is seen in the stress response to vertical flux for the net toroidal field runs. For vertical flux λ_MRI/Δz ≈ 6, we see a flat stress response. This region occurs well below what we would expect for the MRI-stable region in which the most unstable vertical mode exceeds the vertical domain of the simulation, and we also note that the net vertical field runs maintain a linear slope in this region.

A further comparison we can make between local and global models is based not just on the structure of the flux–stress relationships but also on the precise manner in which increased flux generates an increased stress response. While the dependence of local saturation predictors on the box size prevents direct comparisons, we can work backward and use the slope of the resolved region to define an effective local box size. The physical meaning of this length scale is unclear and is presented here as merely a convenient manner in which to quantify the slope of the linear-response regime. Formally, we use the saturation predictors given by Equations (13) and (14) to constrain the value of the box size to fit the instantaneous correlations we find. The values of the box size from the vertical flux–stress and toroidal flux–stress are denoted ℓ_z and ℓ_ϕ, respectively.

The values of this effective box size are given in Table 4. Both the vertical and azimuthal length scales are monotonic with resolution, with the net-field simulations corresponding to larger effective box sizes. Of interest is the approximate relationship, ℓ_z ≈ 2πℓ_ϕ, in particular due to its connection with commonly used local domain sizes. Also of note is that by reducing the azimuthal domain we see a significant drop in ℓ_z, while ℓ_ϕ remains roughly constant.

The manner in which a localized region of the disk responds to the instantaneous presence of flux makes up an important part of the dynamics of the disk; the natural analog to this is the distribution of magnetic flux through the localized regions of the disk. Here, we focus on the latter. We compute a distribution function of the vertical and toroidal magnetic flux through each wedge during the QSS. As in the flux–stress diagrams, these fluxes are scaled to the grid resolution. The distribution is calculated using logarithmic binning in the flux (80 bins per decade), and the result is smoothed using a 10-bin moving window. The results for simulations BzZ32, BpN32, and BzN16 are displayed in Figures 15(a), (b), and (c).

Zoom In Zoom Out Reset image size

Of interest is the fact that for simulation BzZ32, the majority of the disk is below the threshold in which there is a linear stress response to flux. In contrast to this, simulations BpN32 and BzN16 have the majority of the disk well above the resolvability threshold of (Q_z, Q_ϕ) ≈ (1, 10). Related to the initial conditions, BpN32 is biased toward toroidal flux and BzN16 is biased toward vertical. Forming an intermediate case, BzZ32 possesses a more even distribution of flux, albeit weaker than the net-field cases.

An alternative manner by which to consider the distribution of magnetic flux is through the consideration of the local statistics as representative of the trajectory of an ensemble. That is, consider the path ([Q_z](t), [Q_ϕ](t)) and its temporal evolution. Figure 16 illustrates the trajectory for simulations BzZ32, BpN32, and BzN16, where markers are placed on the trajectory approximately every eight orbits. For the net-field simulations the evolution of the flux during the initial transient is quite stark, resulting in migration from the axes to large values of both azimuthal and vertical fluxes, specifically large enough to exceed the transition point to linear stress response. After the initial transient, the net-flux simulations are characterized by a general decrease in the mean flux. In contrast to the net-flux simulations, the trajectory of BzZ32 is constrained to a much smaller region of the flux space. The trajectory only briefly crosses into the region associated with a linear stress response, unlike the net-flux simulations, whose trajectories spend the majority of the simulation in this region. Compared to truly local simulations, in which the boundary conditions prevent migration in the Q_z × Q_ϕ space, global simulations are characterized by significant evolution, particularly those initialized with a net flux. This suggests the possibility that local simulations initialized with a net flux could result in an artificially maintained accretion efficiency due to the boundary conditions preventing the advection of flux off the grid. The trajectory of simulations BzZ32 and BzN16 evolves to possess comparable fluxes despite exhibiting significantly different behavior during much of the simulation. The limited temporal domain for which simulation BpN32 was evolved prevents making a similar comparison; however, the behavior is consistent with the possibility of evolving to a state of magnetic flux comparable to simulations BzZ32 and BzN16. Run BpN16, while run for 50 orbits more than BpN32, displays similar behavior that is consistent with this possibility but does not achieve it during the time domain simulated.

Zoom In Zoom Out Reset image size

Here we consider differentially rotating systems with Ω(R) = Ω₀R^−q. Rayleigh's criterion states that systems of this form will be stable for q < 2 and unstable otherwise. Simulations of differentially rotating systems at the cusp of stability are used to demonstrate that the momentum source terms are being correctly implemented. Additionally, the Keplerian case (q = 1.5) is demonstrated to be stable due to its particular importance. To this end, the values of q considered are q ∈ {1.5, 1.95, 1.99, 2.01, 2.05}, using both orbital advection and the standard cylindrical integrator.

The physical domain simulated is given by (R, ϕ, z) ∈ [3, 7] × [0, π/2] × [ − 0.5, 0.5], using a resolution of (N_R, N_ϕ, N_z) = (128, 200, 32). The simulation is initialized using a constant density ρ₀ = 200, Ω₀ = 2π, and an isothermal sound speed of c_s = 0.1, which results in a flow that is highly rotationally dominated. Note that this is just a three-dimensional extension of the Rayleigh stability test presented in Skinner & Ostriker (2010). As in Skinner & Ostriker (2010), the instability is seeded using perturbations to the azimuthal velocity of the form v_ϕ = v_K(1 + Δ), where Δ ∈ [ − 10⁻⁴, 10⁻⁴] and is uniformly distributed. The simulations are evolved to t = 300 the number of cycles taken and the wall-time to complete each simulation are given in Table 5. Taking as an example the case q = 3/2, corresponding to a Newtonian disk, the simulation without orbital advection takes approximately 24 times as long to run. Calculating the cycles/minute, the full operator-split method takes roughly 1.5 times as long for each cycle. Even in the case of an unstable disk, e.g., q = 2.05, where mass is driven off the grid, there is an almost order-of-magnitude performance increase.

To measure the instability of these systems, we consider a radial scaling of the Reynolds stress normalized to the gas pressure given by

$$\begin{equation} \frac{\displaystyle {\langle R\rho v_{R} v_{\phi }^{\prime }\rangle }}{\displaystyle {\langle RP\rangle }} = \frac{\displaystyle { \int R\rho v_{R}(v_{\phi }-v_{K})dV}}{\displaystyle {c_{s}^{2} \int R \rho dV }}, \end{equation} \tag{ A7 }$$

where the integrals are performed over the entire domain. The time evolution of this quantity for all simulations is shown in Figure 17(a). As expected, there is a clear difference in the evolution of the stress depending on Rayleigh's criterion. For q > 2 (unstable) we see that the scaled stress grows by many orders of magnitude, whereas in the case of q < 2 (stable) we see that the scaled stress remains fairly flat.

Zoom In Zoom Out Reset image size

The final test we consider is perhaps the most relevant. Here we compare the behavior of two identical Newtonian MRI simulations using both orbital advection and the standard cylindrical integrator. We compare the results of run BzZ8 with and without orbital advection. Figure 17(b) illustrates the time evolution of α_M in both simulations. We note that the initial growth and saturation of the MRI are in excellent agreement. Further, the agreement is quite strong well into the highly nonlinear regime. At late times there is a discrepancy between the cylindrical and orbital advection runs, with the cylindrical run exhibiting slightly larger values of stress. This resembles the discrepancy in the behavior of runs BzZ32 and BzZ32R, in which a reduction of azimuthal resolution results in an increased stress. This resemblance is not entirely surprising as, in many ways, the use of orbital advection allows a more precise treatment of azimuthal structure.

Finally, and perhaps most importantly, we discuss the quantitative benefit of orbital advection. To address this point, we compare the time steps used in the evolution of each simulation. The time steps are sampled 10 times per orbit and denoted Δt_C and Δt_F for the cylindrical and orbital advection runs, respectively. We define the speedup as

$$\begin{equation} S(t) = \frac{\Delta t_{F}}{\Delta t_{C}}. \end{equation} \tag{ A8 }$$

The time dependence of the speedup is shown in Figure 17(c). The time step associated with orbital advection is significantly more volatile than without, as the time step is associated with turbulent quantities as opposed to the time-invariant Keplerian rotation. In particular, the time step drops significantly during the initial linear growth of the MRI due to the impulsive accretion associated with exponential growth phase. Upon the saturation and breakup into turbulence, the time step increases again steadily and overall maintains an order-of-magnitude increase in the allowable time step.