SUSTAINED MAGNETOROTATIONAL TURBULENCE IN LOCAL SIMULATIONS OF STRATIFIED DISKS WITH ZERO NET MAGNETIC FLUX

We utilize a number of diagnostic tools based on Fourier analysis to study the simulation output. This is straightforward in a periodic domain, but in a shearing periodic system the basis functions are shearing waves with a time-dependent wavevector. This complication can be handled with simple remappings before and after Fourier transforming, as outlined in Hawley et al. (1995).

The quantities of principal interest will be the power spectral density (PSD) of the magnetic and kinetic energies. The PSD is highly anisotropic in Fourier space, and contours of constant Fourier amplitude are (approximately) ellipsoidal and triaxial in shape. In order to plot one-dimensional quantities, we adopt two averaging procedures. First, we define the shell-integrated power spectrum of the magnetic field as $B^2_k \equiv 4\pi k^2 |\tilde{B}^2(k)|^2$, where $|\tilde{B}(k)|^2$ denotes the average of $|\tilde{B}({\mbox{\boldmath {$k$}}})|^2$ over shells of constant $|{\mbox{\boldmath {$k$}}}|$, and $\tilde{B}({\mbox{\boldmath {$k$}}}) = \int B({\mbox{\boldmath {$x$}}}) \exp {(-i {\mbox{\boldmath {$k$}}} \cdot {\mbox{\boldmath {$x$}}})} d^3 {\mbox{\boldmath {$x$}}}$ is the Fourier transform of B.⁴ To facilitate comparison with the previous work (e.g., Fromang & Papaloizou 2007), we also consider the PSD integrated over planes of constant k_z defined as $B^2_k(k_z) \equiv \int \int \tilde{B}^2({\mbox{\boldmath {$k$}}}) d k_x d k_y$. Note that this is equivalent (S. Fromang 2009, private communication) to the quantity plotted in Figure 6 of Fromang & Papaloizou (2007), although the description in their text implies that $\tilde{B}^2(k_x=k_y=0,k_z)$ (i.e., the power along the z-axis in k-space) is plotted.

It is also instructive to calculate the Fourier transform of the induction equation. Taking Fourier transforms of the x and y components of Equation (3) we find

$$\begin{eqnarray} &&\frac{1}{2}\frac{\partial |\tilde{B}_x({\mbox{\boldmath {$k$}}})|^2}{\partial t} = A_x + E_{z,y} + E_{y,z} \end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray} &&\frac{1}{2}\frac{\partial |\tilde{B}_y({\mbox{\boldmath {$k$}}})|^2}{\partial t} = S + A_y + E_{z,x} + E_{x,z}. \end{eqnarray} \tag{ 10 }$$

We will focus on these two components as they appear to be the most important for understanding the disk dynamo.

The definitions of the terms on the right-hand sides of Equations (9) and (10) are

$$\begin{eqnarray} &&S({\mbox{\boldmath {$k$}}}) = Re \left[\tilde{B}^{\ast }_{y}({\mbox{\boldmath {$k$}}}) \cdot \int B_x \frac{\partial v_{\rm sh}}{\partial x} \exp {(-i {\mbox{\boldmath {$k$}}} \cdot {\mbox{\boldmath {$x$}}})} d^3 {\mbox{\boldmath {$x$}}} \right],\qquad \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &&A_i({\mbox{\boldmath {$k$}}}) = -Re \left[\tilde{B}^{\ast }_{i}({\mbox{\boldmath {$k$}}}) \cdot \int v_{\rm sh} \frac{\partial B_i}{\partial y} \exp {(-i {\mbox{\boldmath {$k$}}} \cdot {\mbox{\boldmath {$x$}}})} d^3 {\mbox{\boldmath {$x$}}} \right],\qquad \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&E_{z,y}({\mbox{\boldmath {$k$}}}) = Re \left[\tilde{B}_x^{\ast }({\mbox{\boldmath {$k$}}})\cdot \int \frac{\partial {\mathcal {E}_z}}{\partial y}\exp {(-i {\mbox{\boldmath {$k$}}} \cdot {\mbox{\boldmath {$x$}}})} d^3 {\mbox{\boldmath {$x$}}} \right],\qquad \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&E_{y,z}({\mbox{\boldmath {$k$}}}) = - Re \left[\tilde{B}_x^{\ast }({\mbox{\boldmath {$k$}}})\cdot \int \frac{\partial {\mathcal {E}_y}}{\partial z}\exp {(-i {\mbox{\boldmath {$k$}}} \cdot {\mbox{\boldmath {$x$}}})} d^3 {\mbox{\boldmath {$x$}}} \right],\qquad \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&E_{z,x}({\mbox{\boldmath {$k$}}}) = - Re \left[\tilde{B}_y^{\ast }({\mbox{\boldmath {$k$}}})\cdot \int \frac{\partial {\mathcal {E}_z}}{\partial x}\exp {(-i {\mbox{\boldmath {$k$}}} \cdot {\mbox{\boldmath {$x$}}})} d^3 {\mbox{\boldmath {$x$}}} \right],\qquad \end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray} &&E_{x,z}({\mbox{\boldmath {$k$}}}) = Re \left[\tilde{B}_y^{\ast }({\mbox{\boldmath {$k$}}})\cdot \int \frac{\partial {\mathcal {E}_x}}{\partial z} \exp {(-i {\mbox{\boldmath {$k$}}} \cdot {\mbox{\boldmath {$x$}}})} d^3 {\mbox{\boldmath {$x$}}} \right],\qquad \end{eqnarray} \tag{ 16 }$$

where the subscript i in A_i refers to either the x or y coordinate. The EMF is defined as ${\mbox{\boldmath {$\mathcal {E}$}}} = {\mbox{\boldmath {$v$}}_t} \times$ B, with ${\mbox{\boldmath {$v$}}_t} = {\mbox{\boldmath {$v$}}} - {\mbox{\boldmath {$v$}}_{\rm sh}}$ and

$$\begin{eqnarray} &&{\mbox{\boldmath {$v$}}}_{\rm sh} = \frac{\hat{\mbox{\boldmath {$j$}}}}{L_y L_z} \int \int v_y dy dz, \end{eqnarray} \tag{ 17 }$$

where L_y and L_z are the sizes of the computational domain in the y- and z-directions. The terms A_x and A_y are included for completeness, but they are generally much smaller than the other terms so we will not discuss them further.

These relations are similar to the transfer functions used in Fromang & Papaloizou (2007) and Simon et al. (2009). In fact, our definition of S is identical and if we sum A_i over all three spatial dimensions, it would be equivalent to their quantity A. These authors expand

$$\begin{equation*} {\mbox{\boldmath {$\nabla $}}} \times ({\mbox{\boldmath {$v$}}_t}\times {\mbox{\boldmath {$B$}}})= ({\mbox{\boldmath {$B$}} \cdot \nabla }){\mbox{\boldmath {$v$}}_t} -({\mbox{\boldmath {$v$}}_t \cdot \nabla }){\mbox{\boldmath {$B$}}}-({\nabla \cdot \mbox{\boldmath {$v$}}_t}){\mbox{\boldmath {$B$}}}, \end{equation*}$$

and perform similar Fourier analysis on the three right-hand side terms individually, labeling them T_bv, T_bb, and T_div, respectively. One drawback of this expansion is that terms such as B_x∂v_x/∂x are present, even though they do not contribute to the evolution of ${\mbox{\boldmath {$B$}}}$, because they appear with opposite signs in both T_div and T_bv. Such terms can be quite large, complicating the interpretation of T_bv, T_bb, and T_div. We therefore prefer to write the right-hand sides in terms of the EMFs.

For plotting purposes, we find it useful to normalize the quantities on the right-hand sides of Equations (9) and (10) with the power spectrum. To differentiate them from the unnormalized quantities, we will use lowercase letters. For example, $e_{y,z}(k) \equiv 2 E_{y,z}(k)/(|\tilde{B}_x(k)|^2 \Omega)$. This then constitutes the shell-averaged Fourier amplitude of the normalized rate of field production of B_x due to the vertical variation of $\mathcal {E}_y$. The factor Ω has been introduced to make the quantities dimensionless rates.

Our primary goal is to test the robustness of sustained turbulence and angular momentum transport in stratified shearing boxes with zero net flux, and in addition, to characterize the properties of the turbulence in this case. Figure 1 is an image showing the three-dimensional structure of the density at late time (250 orbits) in a typical simulation, computed with a resolution of 64 grid zones per scale height in an 4H × 4H × 4H domain. Spiral density waves characteristic of all (compressible) simulations of the MRI in shearing boxes (Heinemann & Papaloizou 2009) are evident in the density isosurfaces.

Zoom In Zoom Out Reset image size

To investigate the convergence of the stress with numerical resolution, we have performed a resolution study at 32, 64, and 128 grid zones per scale height in an H × 4H × 4H domain (hereafter simulations S32R1Z4, S64R1Z4, and S128R1Z4, respectively). Since interest in MRI turbulence is driven primarily by its role in angular momentum transport, we first focus on stress as a diagnostic. Table 1 and Figure 2 summarize the behavior of the stress as we vary the resolution.

Zoom In Zoom Out Reset image size

Table 1. Simulation Summary

Simulation	Domaina	Resolution	Orbits	−〈〈BxBy〉〉/P0b	〈〈ρvxδvy〉〉/P0b	$\langle \langle B_y \rangle \rangle /\sqrt{P_0}$c
S32R1Z4	H × 4H × 4H	32/H	300	0.012	0.0029	0.066
S64R1Z4	H × 4H × 4H	64/H	300	0.0075	0.0018	0.029
S128R1Z4	H × 4H × 4H	128/H	300	0.0076	0.0016	0.034
S32R1Z6	H × 4H × 6H	32/H	165	0.010	0.0024	0.074
S32R4Z4	4H × 4H × 4H	32/H	250	0.0082	0.0022	0.035
S64R4Z4	4H × 4H × 4H	64/H	160	0.0076	0.0019	0.040

Notes. ^aL_x × L_y × L_z. ^bBrackets denote temporal averages taken from 50 orbits onward and volume averages over whole domain. ^cBrackets denote temporal averages taken from 50 orbits onward and volume averages over innermost two scale heights.

Download table as:  ASCIITypeset image

The MRI turbulence contributes to the stress through the Maxwell stress −B_xB_y and the Reynolds stress ρv_xδv_y, where δv_y is the y component of the velocity with the background shear removed. The domain and time average of these quantities are listed in Table 1. We have normalized them by the initial midplane pressure P₀. Since the initial condition is in hydrostatic equilibrium and magnetic pressure is always significantly lower than gas pressure at the midplane (see, e.g., Figure 6), this value of the midplane pressure does not evolve significantly. These normalized stresses are roughly equivalent to the α parameter of Shakura & Sunyaev (1973). The time average is carried out from 50 to 300 orbits to exclude the transient period of enhanced stress during and immediately after the linear growth phase of the MRI. Alternatively, we could have normalized the stress with the domain-averaged pressure 〈P〉, which is a factor of a few smaller, and therefore yields a slightly larger effective α. For the initial pressure profile, 〈P〉 = 0.441P₀ and 0.295P₀ for the L_z = 4H and L_z = 6H domains, respectively.

Both contributions to the stress decrease as we increase resolution, but the changes are much smaller when going from 64/H to 128/H than from 32/H to 64/H, indicating convergence. This should be compared with the behavior observed in unstratified boxes (e.g., Sano et al. 2004; Pessah et al. 2007; Fromang & Papaloizou 2007) in which the total stress decreases by factors of ∼2 as resolution is increased by a factor of 2. The normalized total stress in the S128R1Z4 run is 0.0095, comparable to previous results for stratified domains with zero net flux (Brandenburg et al. 1995; Stone et al. 1996; Johansen et al. 2009; Suzuki & Inutsuka 2009). The Maxwell stress is 4–5 times greater than Reynolds stress, which is slightly higher than, but roughly consistent with previous results for stratified domains (e.g., Stone et al. 1996). Similar values are also observed in unstratified runs, where the result appears to be independent of field geometry or flux, and depend mainly on the rate of shear (Pessah et al. 2006).

Table 1 also includes the rms toroidal field, volume averaged over the central two scale heights and time averaged from 50 orbits onward. The rms field strengths are relatively weak, with 〈B_y〉² ∼ 0.01〈B²_y〉, but may still be dynamically important since the presence of a mean toroidal field in unstratified simulations has been shown to enhance turbulent stresses and energy densities (Hawley et al. 1995). In fact, the rms toroidal field strength correlates well with the stress, although this may be the by-product of stronger turbulence rather than a cause.

Figure 2 shows the temporal variation of the total stress. There is considerable variability, with relatively long-lived (≳50 orbits) departures from the mean. In the S32R1Z4 run, the ∼50 orbit periods of enhanced stress contribute almost as much to the average as the longer periods of "quiescent" stress. There are similar long-term fluctuations in the higher resolution runs, but these are generally smaller in amplitude and less important for determining the overall mean. Nevertheless, it is clear that one must average over relatively long baselines to obtain a representative value, making higher resolutions prohibitively expensive at the present time.

A power spectral analysis of the magnetic field also indicates convergence, at least for the large scales where most of the power resides. Figure 3 shows the PSD for the S32R1Z4, S64R1Z4, and S128R1Z4 simulations. In the top and bottom panels we integrate over spherical shells and planes of constant k_z (respectively, see Section 2.1) in k-space. We also average in time, excluding the first 50 orbits to avoid the initial transients. As in the spacetime plots, there is substantial variability in the power spectra on 50 orbit timescales, and the variance is significant. However, the uncertainties in the time-averaged mean PSD are small due to the long integration times. Therefore, we do not plot error bars, which are only slightly larger than the thickness of the curves used for plotting. Since L_x = H, the largest k_x = 2π/H, and so we plot shells with k ⩾ 2π/H. The plane-integrated PSD extends to the smallest k_z = π/(2H). There is a significant drop in power in going from the 32/H to 64/H resolution runs, combined with a shift in the peak of PSD to larger k. However, when going from 64/H to 128/H, there is a significantly smaller drop in amplitude at most scales and a much smaller shift in the peak wavenumber, also indicating convergence.

Zoom In Zoom Out Reset image size

We reach the same conclusions as Shi et al. (2009), who also find convergence in simulations with stratified domains. Their results are complimentary to ours in that they considered larger domains (in terms of physical dimensions, particularly height), used outflow boundary conditions, and considered more realistic thermodynamics, including radiation. However, this comes at the cost of significantly limiting the computational size (or resolution) of their simulation domains. Their standard run is computed on a domain with 32 × 64 × 256 zones, which has the same total number of zones as our S32R1Z4 simulation. Convergence is then evaluated by doubling the resolution once in each spatial dimension while keeping the other two fixed. Since we simultaneously quadruple the resolution in all three spatial dimensions, our largest domain (S128R1Z4) is a factor of 32 larger than their largest domain.

Although we focus on the magnetic energy density, a similar convergence is observed in the PSD of the kinetic energy density. We show a comparison of the perturbed kinetic and magnetic energy density PSDs for run S128R1Z4 in Figure 4. Specifically, we are plotting the PSD of the perturbed kinetic energy, meaning that the background shear has been subtraced from the azimuthal velocity before taking the Fourier transform. It is notable that power in the magnetic field fluctuations exceeds that of the kinetic energy at all scales. This result differs from simulations of helically driven turbulence where the kinetic and magnetic energy have comparable amplitude at all but the lowest k where magnetic energy dominates (see, e.g., Brandenburg 2001).

Zoom In Zoom Out Reset image size

The convergent behavior of the stratified runs should be contrasted with that of the unstratified simulations shown in Figure 5. This plot shows the PSD for four unstratified runs with resolutions of 32/H, 64/H, 128/H, and 256/H in H × 4H × H shearing boxes. Each factor of 2 increase in resolution results in a decrease by nearly a factor of 2 in the integrated power. There is also a shift in the peak wavenumber to larger k as resolution increases. This lack of convergence is almost identical to that found by previous authors (Fromang & Papaloizou 2007; Guan et al. 2009; Simon et al. 2009). It seems that both the power and characteristic scale of the turbulence are set by the domain resolution. Adding stratification appears to fundamentally change the dynamics and provides a characteristic scale and amplitude of the turbulence which is independent of the resolution. This could be related to the different mechanisms that lead to the saturation of the MRI in the presence of stratification, perhaps associated with the development, and subsequent buoyant rise, of large-scale magnetic fields (Pessah et al. 2007; Vishniac 2009, see also Section 4).

Zoom In Zoom Out Reset image size

In Figure 6, we show spacetime diagrams for several variables associated with the magnetic field. For the sake of brevity, we focus on magnetic quantities, which appear to play the dominant role, as suggested by the PSD analysis above. The spacetime plots are generated by averaging over the x and y coordinates at each height in the domain every quarter of an orbit. The top panel of Figure 6 shows β = 2〈P〉_xy/〈B²〉_xy, where the angle brackets denote horizontal averages. As noted previously, magnetic pressure remains relatively weak near the midplane, but dominates in the surface regions (|z| ≳ 1.5H).

Zoom In Zoom Out Reset image size

The middle panels show the normalized x and y components of the magnetic field. The B_y component is considerably larger than B_x and they are negatively correlated. The symmetry of the boundary conditions and the initial conditions require the vertical average of these quantities to be zero. Nevertheless, there are localized regions with net field that trace out curved trajectories in the spacetime plot, indicative of buoyantly rising regions of the strong magnetic field. This is similar to other "butterfly diagrams" seen in previous shearing box calculations of stratified domains (Brandenburg et al. 1995; Stone et al. 1996; Turner 2004; Johansen et al. 2009; Suzuki & Inutsuka 2009). Consistent with previous simulations, the polarity is usually even about the midplane and, at fixed height, oscillates quasi-periodically on timescales of ≲10 orbits.

Within the inner three scale heights, the horizontally averaged net field is a sum over a fluctuating $\mbox{\boldmath {$B$}}$ field so that 〈B_y〉²_xy is much less than 〈B²_y〉_xy and similarly for B_x. As these regions of net field rise, the ratio of 〈B_y〉²_xy/〈B²_y〉_xy increases. Near the vertical boundaries magnetic dominated regions of rather uniform ${\mbox{\boldmath {$B$}}} \sim B_y {\hat{\mbox{\boldmath {$j$}}}}$ develop. Their presence is very likely related to our assumption of periodicity in the vertical boundary condition, so they are likely not physically relevant to real accretion flows. The degree to which these regions affect the dynamics is discussed further in Section 3.3.

The bottom panel of Figure 6 shows the Maxwell stress. In addition to the temporal fluctuations, there is also considerable variation with height. The bulk of the stress is contributed by the middle three scale heights and the regions of greatest stress tend to be found off the midplane. The stress is weakest in the magnetically dominated regions very near the boundary and is even slightly negative in some places (although the colorbar only goes to zero). A striking result is the correlation of regions of stronger than average 〈B_y〉_xy with regions of larger than average stress. This lends support to the idea that the presence of a mean toroidal field leads to enhanced turbulent stresses and energy densities. Although not shown, we note that the spacetime plot of the Reynolds stress is very similar to the Maxwell stress and the two are well correlated in both z and time.

Figure 6 clearly displays a significant amount of large-scale and long timescale structure in space and time coordinates, respectively. To better understand and quantify this structure, we perform a complimentary Fourier analysis in both space and time. Since we are interested in horizontally averaged quantities, we focus on vertical k-vectors with ${\mbox{\boldmath {$k$}}}=k_z \hat{\mbox{\boldmath {$z$}}}$. Every one-quarter of an orbit, we compute the discrete Fourier transform $B^2(k_{z}) \equiv 4\pi k_{z}^2 |\tilde{B} (k_{z})|^2$. Note that this quantity corresponds to the k_z-axis (k_x = k_y = 0) in k-space, and therefore differs from the quantity plotted in the bottom panel of Figures 3 and 5, which is the integral over planes of constant k_z in k-space. It is equivalent to taking the vertical Fourier transform of the horizontally averaged domain. For each k_z, we Fourier transform in time to obtain B²(k_z, f) where f is the time domain frequency. We divide the data into five 50 orbit time series (between 50 and 300 orbits), Fourier transform each separately, and then average. Although we lose access to the longest timescales, the averaging reduces the variance in the resulting power spectra.

In Figure 7, we plot both log [fk_zB²_x(k_z, f)/(2P₀)] (top panel) and log [fk_zB²_y(k_z, f)/(2P₀)] (bottom panel). Note that the horizontal average 〈B_z〉_xy is conserved (at zero) to round-off error, so there is no physical information in the z component for vertical wavevectors. We have multiplied by both k_z and f before taking the logarithm. Since we use a logarithmic scale for both the k_z and f axes, this weights each pixel so that its color scales linearly with its contribution to the overall power (i.e., in the same sense that νF_ν is commonly used in astrophysics).

Zoom In Zoom Out Reset image size

There are significant differences in the morphologies of the B_x and B_y power spectra. For the largest spatial scales (smallest k_z), both B_x and B_y show a peak in power at ∼5–10 orbit timescales. As we go to larger k_z > 4π/H the B_x power is widely distributed across the whole frequency range. In the case of B_y, power is also broadly distributed in frequency, but there is a drop in power at high frequencies for small k_z. The locus of maximum power shifts to higher k_z as f increases.

Due to the low densities and high magnetic field strengths near the vertical boundary that arise from stratification, increasing the vertical extent of the domain is particularly computationally expensive. Therefore, we performed our resolutions study in boxes with four vertical scale heights, which seemed suitable to get a reasonable separation between the midplane dynamics and the vertical boundary. In order to confirm that our results are not strongly dependent on this choice, we have repeated our 32/H resolution run with six vertical scale heights (hereafter S32R1Z6). As one can see in Table 1, the resulting volume average of the stress in the two 32/H simulations is in reasonable agreement, although slightly smaller in the S32R1Z6 run.

One downside to choosing vertically periodic boundary conditions is the buildup of strongly magnetized regions with uniform B_y near the boundaries. Although these regions do not contribute significantly to the angular momentum transport, they are likely unphysical so one might worry that they feedback on the dynamics closer to the midplane. In order to get a better sense of their effect on the flow, we have plotted spacetime diagrams comparing the S32R1Z4 and S32R1Z6 simulations in Figures 8 and 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows the y component of the magnetic field. In a larger domain, there are still regions of rather uniform B_y near the vertical boundaries. They are somewhat more extended in height but with a slightly weaker net field. The regions of net B_y generated near the midplane rise to larger heights before interacting with the boundary region. Since the horizontally averaged field tends to increase as the fluid rises, this leads to further enhancement of the field strength over those found in the smaller domain.

Figure 9 shows the Maxwell stress in the two runs. In the central four scale heights the two plots look very similar, suggesting that the four scale height runs are yielding a fairly robust estimate for the angular momentum transport. Regions of enhanced Maxwell stress are again correlated with regions of strong net B_y. In the S32R1Z6 run, the Maxwell stress tends to peak somewhat off the midplane at about |z| ≃ 1 − 1.5H, and then drops nearly to zero for |z| ≳ 2H. Therefore, the stress does not seem to be as impacted by the vertical boundary as in the L_z = 4H runs which (on average) only drop to 30%–50% of their peak values at the boundary.

The volume-weighted average stresses in Table 1 are lower for S32R1Z6 than for S32R1Z4 in part because we are averaging the stress over the whole box, including the regions of weak stress near the boundaries, but dividing by the constant midplane pressure rather than vertically averaged pressure. Clearly the vertically averaged pressure is lower in the box with larger vertical extent, simply because the average includes the lower pressure regions at |z|>2H. If we instead normalize the Reynolds and Maxwell stresses with the vertically averaged pressure (to be consistent with our vertical average of the stress), the S32R1Z6 run yields ∼30% larger normalized stresses than the S32R1Z4 run.

For the sake of computational expediency, we have performed our resolution study on domains with only one scale height in the radial direction. This has traditionally been the box size employed in most shearing box computations, mostly due to CFL constraints on the timestep which are imposed by the background shear. Using Athena's orbital advection scheme (discussed in Section 2), we can consider larger domains to examine the effect of this choice on our results. We have computed shearing boxes with 4H × 4H × 4H domains at 32/H and 64/H resolutions (hereafter S32R4Z4 and S64R4Z4, respectively).

We plot the evolution of the total stress in these two simulations in Figure 10 and the normalized, time- and volume-averaged stresses are listed in Table 1. The mean values of the stress are very similar to each other and also to those found at higher resolutions runs in the smaller boxes (S64R1Z4 and S128R1Z4). This suggests that convergence is occurring at even lower resolution than in the smaller domain computations. The time evolution differs from that seen in Figure 2 in that the amplitude of fluctuations is much lower.

Zoom In Zoom Out Reset image size

In Figure 11, we compare the PSDs from simulations with different radial extent but the same resolution (S64R1Z4 and S64R4Z4). The PSDs are very similar at all but the lowest k. At low k the differences are accounted for in part by our shell averaging scheme. Since the larger box is a cube with L_x = L_y = L_z = 4H, we can isotropically average shells all the way to k = π/(2H). Since we cannot do this average isotropically with smaller box, some of this low k < 2π/H power is included in the k = 2π/H bin. Overall, the PSDs seem to be rather independent of the aspect ratio, consistent with nearly identical values for the box-integrated stresses.

Zoom In Zoom Out Reset image size

Figure 12 shows the spacetime diagram for the S64R4Z4 run. Overall, it is rather similar to S128R1Z4 spacetime diagram in Figure 6. There are long timescale (≳50 orbit) fluctuations in the Maxwell stress, as well as ∼10 orbit quasi-periodic variations in B_y and the Maxwell stress which are qualitatively similar to those in Figure 6. However, the amplitude of fluctuations is generally smaller in the larger box, consistent with the lower amplitude fluctuations in the total stress found in Figure 10.

Zoom In Zoom Out Reset image size

In order to better understand the lower fluctuation amplitudes, we have split the larger domain into four sub-domains of one scale height each in the radial direction. We find that the volume-averaged stresses in each subdomain are highly correlated with each other. Therefore, the lower amplitude is not simply the result of "averaging" over several independent boxes, but appears to be a global property of a correlated domain. This behavior is somewhat surprising in light of the results of Guan et al. (2009), which show that the turbulence decorrelates on scales ≳H in unstratified boxes.

To investigate the correlation in stratified boxes, we have followed Guan et al. (2009) and calculated the trace of the two-point correlation of the magnetic field

$$\begin{eqnarray} &&\xi _B=\langle \delta B_i({\mbox{\boldmath {$x$}}}) \delta B_i({\mbox{\boldmath {$x$}}}+\Delta {\mbox{\boldmath {$x$}}})\rangle, \end{eqnarray} \tag{ 18 }$$

where there is an implied summation over i and δB_i = B_i − 〈B_i〉. We computed the average in two ways: using the full domain and only the innermost two scale heights. The two different procedures give different results for the correlations in the x–y plane at large separations, since full domain average is dominated by the magnetically dominated regions near the vertical boundary. Since these are likely unphysical, we report only the two scale height average.

We compute the correlation for 21 and 18 evenly spaced snapshots for the S64R1Z4 (50–300 orbits) and S64R4Z4 (50–150 orbits) simulations, respectively. Although correlation at small separations is relatively consistent from snapshot to snapshot, there is significant variation at larger scales. Therefore, we average over all snapshots to get a mean correlation for each run. We plot the results of the two scale height average for the x–y plane in Figure 13. For both simulations we normalize ξ_B by their maximum values, ξ₀, which agree to within 10%. Our results are qualitatively consistent with those of Guan et al. (2009) in that we find comparable tilt angles θ_t, which is the angle between the major axis of the correlation function and the azimuthal axis. The correlation is very similar in both simulations, although the tilt angle from the S64R4Z4 calculation is slightly smaller (θ_t ∼ 15° rather than 18°).

Zoom In Zoom Out Reset image size

We also plot the correlation along the minor axis (defined by $\hat{\mbox{\boldmath {$x$}}} \cos \theta _t + \hat{\mbox{\boldmath {$y$}}} \sin \theta _t$) in Figure 14. Again, the core of the correlation at small separations is nearly identical in the two simulations and consistent with the results of Guan et al. (2009). Near Δl ∼ 0.1H the slope of the curve from the S64R1Z4 run flattens and the large-scale correlation plateaus at a value of ξ_B ≃ 0.04ξ₀. The S64R4Z4 curve declines further, also flattens with ξ_B ≃ 0.01ξ₀ until Δl ∼ 1.3H where it drops nearly to zero.

Zoom In Zoom Out Reset image size

As Equation (18) requires, we have subtracted the mean field which is generally non-zero when only the inner two scale heights are considered (see Table 1), but is identically zero when using the whole domain. These mean fields can be substantial and would dominate the large-scale correlation if included. It seems that these fields are sufficient to enforce relative uniformity in the magnetic energy and Maxwell stress throughout the four scale height domains. This uniformity may disappear for sufficiently large domains, and we see some suggestion of this in calculations at 32/H resolution where we have compared 8H × 8H × 4H and 4H × 4H × 4H domains. Although we find that in both cases the Maxwell stress and magnetic energy density in one scale height wide subdomains remain correlated, they show greater variance in the eight subdomains of the large box than in the four subdomains of the smaller box.

It has been shown that the behavior of MRI in the linear regime (e.g., Sano & Miyama 1999; Pessah & Chan 2008), and the saturation of turbulence in unstratified shearing boxes depends on the nature of the dissipation (Sano et al. 1998; Fleming et al. 2000; Fromang et al. 2007; Lesur & Longaretti 2007; Simon & Hawley 2009). Simulations with explicit dissipation yield different results than those with only numerical dissipation, and the strength and evolution of the turbulence depends on both the viscosity ν and resistivity η, or equivalently the Reynolds number Re ≡ c_sH/ν and the magnetic Reynolds number Rm ≡ c_sH/η. Specifically, it has been shown that Re and Rm, or alternatively the magnetic Prandtl number Pm = Rm/Re, determine whether turbulence is sustained in zero net flux, unstratified simulations (Fromang et al. 2007; Simon & Hawley 2009).

To see if these results hold in stratified domains, we perform three simulations with differing values of viscosity and resistivity: Re = 800 with Pm = 4 (hereafter Re800Pm4), Re = 800 with Pm = 2 (Re800Pm2), and Re = 1600 with Pm = 2 (Re1600Pm2). Examination of Figure 11 of Fromang et al. (2007) or Table 1 of Simon & Hawley (2009) show that none of these would sustain turbulence in an unstratified domains with no net field, regardless of whether Zeus or Athena is used for the simulations. We have confirmed these results for unstratified domains with our own Athena calculations.

Figure 15 shows the stress as a function of time for the three simulations with a logarithmic vertical scale. The behavior of Re800Pm4 and Re1600Pm2 is rather different from the unstratified domains where the simulations decay rapidly to zero on timescales of 10–20 orbits after the initial linear growth of the MRI. Even Re800Pm2, which drops rapidly to a dimensionless stress of ∼10⁻³ and then ∼10⁻⁴ in the stratified domain, decays much more rapidly and continues to even lower values in the unstratified domain. The behavior of the stratified domains is also considerably more complicated. In two of the three simulations, turbulence never decays away completely, but vigorous turbulence is not sustained in any of the calculations for longer than 100 orbits. The amplitude of variability is large, and turbulence decays slowly on timescales of hundreds of orbits. The Re1600Pm2 and Re800Pm2 runs both show recoveries to nearly peak values after spending over 100 orbital periods in stagnation or slow decay. Only the Re800Pm4 simulation completely decays, but it takes over 200 orbits to do so.

Zoom In Zoom Out Reset image size

Using the criteria of Fromang et al. (2007), we would probably have labeled the Re800Pm4 run as having sustained turbulence (over the first 100 orbits, which is the baseline used there), the Re1600Pm2 run as marginal, and the Re800Pm4 as either marginal or not having sustained turbulence, although the complex variability of the stratified runs makes this somewhat subjective. Figure 11 of Fromang & Papaloizou (2007) maps out a locus of sustained turbulence in Re–Pm space, and it is notable that Re800Pm4 and Re1600Pm2 simulations are on the cusp of showing sustained turbulence while Re800Pm2 more firmly in the non-turbulent regime. Therefore, it would seem that stratification slightly increases the parameter space for which sustained turbulence is possible, but does not qualitatively alter the conclusion that turbulence dies out for sufficiently low Pm.

It is suggestive that all three sets of dissipation terms show sustained turbulence in unstratified boxes once a net toroidal field is imposed (Simon & Hawley 2009). As noted previously, it is conceivable that the main impact of stratification is the production of the (local) mean toroidal magnetic field, which then leads to enhanced turbulence. The turbulence in the stratified runs is significantly less vigorous, but this may be consistent with the rms field strengths in the stratified simulations being much weaker than the toroidal fields considered by Simon & Hawley (2009).

In Figure 16, we plot the time-averaged power spectra of the magnetic energy density from Re800Pm4 and Re1600Pm2, including S64R1Z4 for comparison. Since the magnetic energy in the Re1600Pm2 drops rapidly to a low amplitude, we have elected to exclude it. Of course, the amplitude of the power spectrum depends on the interval used in the time average, which is 25–100 orbits. The Re800Pm4 and Re1600Pm2 power spectra are similar in shape, falling off somewhat more rapidly than S64R1Z4 as k increases.

Zoom In Zoom Out Reset image size

The power in Re800Pm4 exceeds that in Re1600Pm2 at all k. Note that these two calculations have the same resistivity but Re800Pm4 has a higher viscosity, so it has higher amplitude despite having larger overall dissipation (although this conclusion depends on the interval used for the comparison). Similar behavior is also observed in two-dimensional simulations of MRI driven turbulence with a vertical magnetic field and non-zero viscosity described in Masada & Sano (2008). In some of their runs saturation is not achieved since the turbulent stresses increase with time until the end of the runs. In this case, the fact that the magnetic field generated by the MRI can reach high amplitudes could be due to the viscous quenching of Kelvin–Helmholtz parasitic modes (Goodman & Xu 1994; Pessah & Goodman 2009). Although plausible, it is less obvious that this process is responsible for the similar behavior observed in the simulations with stratification and non-mean magnetic flux presented here.

First, we compare our results with the discussion of Brandenburg et al. (1995), who noted the similarities of their stratified results with an α–Ω dynamo model. They showed that the azimuthal EMF $\mathcal {E}_y =(\delta {\mbox{\boldmath {$v$}}} \times {\mbox{\boldmath {$B$}}})_y$ was correlated with the azimuthal field B_y in a manner that leads to growth in the radial B_x. Coupled to the shear, which regenerates B_y from B_x, this describes a simple dynamo. Our results are in good agreement with their Equation (21). For the S128R1Z4 simulation, we find

$$\begin{eqnarray} &&\langle \mathcal {E}_y \rangle \simeq (-3, 5) \times 10^{-3} \langle B_y\rangle \langle {\delta v}^2 \rangle ^{1/2}, \end{eqnarray} \tag{ 19 }$$

where the values in parentheses correspond to volume averages one scale height above and below the midplane, respectively. They also report a correlation of $\langle \mathcal {E}_x \rangle$ with 〈B_y〉 in their Equation (22). Again, our results are qualitatively similar to theirs with

$$\begin{eqnarray} &&\langle \mathcal {E}_x\rangle \simeq (-1, 1) \times 10^{-2} \langle B_y\rangle \langle {\delta v}^2\rangle ^{1/2}, \end{eqnarray} \tag{ 20 }$$

above and below the midplane. The first correlation confirms that there is a mechanism for regenerating the poloidal field from a toroidal field. The second relation suggests that $\langle \mathcal {E}_x\rangle$ generally acts to reduce the magnitude of 〈B_y〉 at the midplane, and is at least partially the result of vertical motion of buoyant magnetic structures, as Brandenburg et al. (1995) speculate.

It is instructive to examine this further with the Fourier analysis methods described in Section 2.1. We focus on the large-scale field and consider the smallest vertically oriented vector ${\mbox{\boldmath {$k$}}} = k_1 \hat{\mbox{\boldmath {$z$}}} =2 \pi /L_z \hat{\mbox{\boldmath {$z$}}}$. Since, k_x = k_y = 0, this term represents the Fourier amplitude of "horizontally averaged" quantities on the largest vertical scale. In Figure 17, we plot time variation of magnetic fields and EMFs over a single dynamo cycle for this choice of $\mbox{\boldmath {$k$}}$. The top panel shows the Fourier amplitudes of magnetic energy densities $|\tilde{B}_x(k_1)|^2$ (solid) and $|\tilde{B}_y(k_1)|^2$ (dotted) for this wavevector. We have multiplied $|\tilde{B}_x(k_1)|^2$ by a factor of 400 to plot both on the same scale. There is considerable variation from cycle to cycle. As Figure 7 shows, these oscillations are only quasi-periodic with a broad range of power for periods near ten orbits, depending somewhat on the wavenumber used for the analysis. However, the phase difference between $|\tilde{B}_y(k_1)|^2$ and $|\tilde{B}_x(k_1)|^2$, and the more uniform variation of $|\tilde{B}_y(k_1)|^2$ relative to $|\tilde{B}_x(k_1)|^2$ are generic to all cycles.

Zoom In Zoom Out Reset image size

In Figure 17, the bottom panel shows the corresponding right-hand side quantities in Equations (9) and (10), which, along with numerical dissipation, drive the evolution of the Fourier amplitudes. We normalize these quantities with the power spectra as described in Section 2.1, e.g., $e_{y,z}(k_1) = 2 E_{y,z}(k_1) /(|\tilde{B}_x(k_1)|^2 \Omega)$. We plot e_y,z (solid), e_y,z (dotted), and s (dashed). Note that e_z,y(k₁) and e_z,x(k₁) are zero for k₁ because of the periodic boundaries. The shear term s primarily drives the variation of $|\tilde{B}_y(k_1)|^2$, flipping sign as $|\tilde{B}_x(k_1)|^2$ goes to zero. In contrast, e_x,z is generally negative, acting as turbulent resistivity. The e_y,z(k₁) term is more erratic, frequently flipping sign over a single cycle, but the net effect is an overall oscillation of $|\tilde{B}_x(k_1)|^2$ over ∼7 orbital periods.

In many respects, the behavior we see in the stratified simulations is similar to that observed in the unstratified, zero-net flux calculations of Lesur & Ogilvie (2008). Using an incompressible spectral code, they find dynamo cycles with a ∼5 orbit periodicity. This is similar to the oscillations in our stratified runs where rms power on large scales is broadly distributed on timescales ∼5–10 orbits (see Figure 7). The normalized quantities plotted in Figure 17 are equivalent⁵ to their Equation (16). Comparison of Figure 17 with Figures 4 and 5 in their paper, show that the behavior of the EMFs during oscillations are also quite similar, suggesting that a common (or, at least, related) mechanism may be responsible for these oscillations. This motivates a more detailed comparison between unstratified and stratified runs in future work.

Although it is useful to focus on the vertical wavevectors when trying to understand properties of large-scale fields, an understanding of the overall power spectrum benefits from an analysis of the shell-integrated quantities. We plot the time and shell average EMFs terms in Figure 18 for S32R1Z4 (black), S64R1Z4 (blue), and S128R1Z4 (red). As in Figure 17, these quantities are normalized by the shell-integrated power spectra. Each normalized term is then time averaged from 50 to 300 orbits. Since the amplitudes of the magnetic energy densities seem to be in statistical steady states over this period, we presume the left-hand sides of Equations (9) and (10) are nearly zero. Therefore, the sum of the terms in each panel must be balanced by numerical dissipation terms, as discussed in the previous work (Fromang & Papaloizou 2007; Simon et al. 2009).

Zoom In Zoom Out Reset image size

In the top panel, we plot the terms e_z,y(k) (solid) and e_y,z(k) (dotted) which contribute to the evolution of B_x. At the large scales, we find that the e_z,y term is more important for field generation and its normalized amplitude is nearly independent of resolution. The e_y,z term is smaller in amplitude and slightly negative as large scales. Even though e_y,z tends to oscillate about zero over an individual dynamo cycle while e_z,y usually remains positive, the amplitude of e_z,y is generally larger, so the dominance of e_z,y at large scales is not simply the result of time averaging. As one moves to smaller scales, e_y,z rises and eventually dominates the generation of B_x. The characteristic wavenumber at which the crossing occurs shifts to higher values as the resolution increases.

The bottom panel shows e_z,x(k) (solid), e_x,z(k) (dotted), and s (dashed), the terms which contribute to the growth in B_y. At large scales growth of B_y is dominated by the shear term, while both e_z,x and e_x,z are of comparable magnitude and negative. At small scales, e_z,x and e_x,z both increase, becoming positive and dominating over the shear term. Again, the wavenumber of the crossover increases with resolution.

Previous work has often focused primarily on horizontally average properties of the flow or (equivalently) the power spectral variation only along vertical wavevectors (e.g., Brandenburg et al. 1995; Lesur & Ogilvie 2008, which were discussed above). We note that the behavior of e_y,z and e_z,y we have described differs significantly from what one would infer if only vertical wavevectors were considered. As previously mentioned, the symmetries of the periodic box force e_z,y(k_z) to be zero, and only e_y,z(k_z) contributes. However, it is clear from Figure 18 that the vertical EMF and its toroidal variation is also essential for understanding the mechanism which sustains turbulence in these simulations.

The question remains as to why the addition of stratification leads to convergence in the turbulent stresses and energy densities. One possibility is that development of the local toroidal field is a key to sustaining turbulence in both stratified and unstratified domains. It is possible that the strength of the toroidal field is entirely set by the resolution in unstratified domains, while stratified domains offer a characteristic scale which is independent of resolution, due to the action of the large-scale dynamo. Indeed, it has already been demonstrated (e.g., Hawley et al. 1995; Simon & Hawley 2009) that an imposed net toroidal field leads to enhanced turbulent energy densities and stresses, and leads to convergence in the stress as resolution increases (Guan et al. 2009). Furthermore, our simulations show a correlation between the stress and the strength of the mean toroidal field, both globally in the two scale height averages (Table 1) and locally in the spacetime plots (Figures 6 and 12).

Recently, Vishniac (2009) discussed possible differences in the mechanism for saturation of the MRI between stratified and unstratified simulation domains. In both sets of simulations he postulates that the saturation of the turbulence is set by a competition between an MRI dynamo driven by a vertical helicity flux and vertical mixing due to some secondary instability. Balancing the growth rate of the dynamo with the dissipation rate due to the mixing gives a saturation amplitude which depends linearly on the size of a "typical" eddy. The main difficulty with applying this condition is defining what is meant by a typical eddy, since the MRI appears to excite power on a large range of scales. Motivated by the power spectral analysis of Fromang & Papaloizou (2007), Vishniac assumes a roughly constant Fourier amplitude, but with a high k_z cutoff. Although relatively large scales are expected to dominate the helicity flux, these only account for a fraction of the total energy, giving an effective wavenumber which depends on the cutoff. If this cutoff scales linearly with resolution, the model predicts a saturation amplitude proportional to the resolution to the 2/3 power. In stratified domains, Vishniac argues that buoyancy will set a minimum eddy thickness that corresponds to a high k_z cutoff that is independent of the resolution. The resulting cutoff is strongly dependent on the ratio of the characteristic size of the large-scale field L_B to the pressure scale height L_P.

Our results in unstratified domains find a scaling for the magnetic energy density that varies linearly with resolution, in agreement with Fromang & Papaloizou (2007). This variation is somewhat stronger than predicted by Vishniac (2009), although he offers some brief speculation why this might be the case. For the stratified domains, the Vishniac's model does predict a saturation amplitude that is independent of resolution, consistent with our results for the highest resolution domains. However, uncertainties in defining the characteristic sizes of a typical eddy and specifying the size of large-scale field make it difficult to test the model more quantitatively. Examination of the bottom panel of Figure 3 shows that power is distributed roughly evenly over a range of k_z, which two broad "peaks." Recall that our preferred convention for log–log plots is to multiply the Fourier amplitude by k_z so that the plot conveys how power is distributed per logarithmic interval. Therefore, a constant value in Figure 3 corresponds to the Fourier amplitude declining as k⁻¹_z. Therefore, Vishniac's assumption of a constant Fourier amplitude, with a high k_z cutoff is not a good description of our simulations. Even if we ignore this discrepancy in the assumptions, the prediction for the saturated magnetic energy density depends on (L_B/L_P)⁴. In the stratified domains this ratio will vary with height and, although we expect L_B ∼ L_P ∼ H, L_B is not defined precisely enough to make a quantitatively testable prediction. Ultimately, a better test of the ideas described in Vishniac (2009) would involve examination of the vertical helicity flux and its dependence on k_z, but such calculations are beyond the scope of the current work.

It is clear from the above discussions that buoyancy plays some role in determining the dynamics of the large-scale field. Therefore, it is worth clarifying the role of the Parker instability in our simulations. Consistent with Stone et al. (1996), we find that the temporal and vertical average of our L_z = 4H domains are stable (at least marginally) to both modes (undulatory and interchange) of the Parker instability. Since our equation of state is isothermal, this is equivalent to noting that our time average profiles of magnetic energy density are roughly constant with height. In the outer 1.5H of the L_z = 6 domains, the time and horizontally averaged profiles of magnetic energy decrease outward, and are therefore Parker unstable. These results are largely consistent with those of Shi et al. (2009), who also find that the inner ±2H are stable, but that Parker instability has a significant impact on the dynamics farther from the midplane. Since they consider domains with larger vertical extent and more realistic thermodynamics, they are able to explore this physics in much more detail.

Due to our choice of periodic vertical boundaries, our use of simplified thermodynamics, and the inherent limitations for the shearing box approximation, we have largely avoided detailed discussion of observational implications. Such questions are generally better addressed by studies which include more physically realistic vertical boundary conditions (e.g., Miller & Stone 2000) and more realistic thermodynamics, including the treatment of radiation (e.g., Turner 2004; Hirose et al. 2006; Shi et al. 2009). However, it is worth briefly noting that our work confirms some important results of earlier studies (see, e.g., Brandenburg et al. 1995; Stone et al. 1996; Miller & Stone 2000). Figures 3 and 6 show that a significant fraction of the magnetic energy in these simulations resides in large-scale magnetic fields that rise buoyantly to the low density regions above the disk midplane. Blackman & Pessah (2009) argue that the magnetic field structures that power accretion disk coronae must be associated with characteristic lengths that are large compared to the typical turbulent eddies. If this were not the case, the timescales associated with turbulent diffusion would be smaller than the corresponding buoyant rise time, making it difficult to transport significant magnetic energy to the coronae. In other words, if the corona is a consequence of magnetic field structures that originate within the turbulent disk via the MRI (or other magnetic instabilities), but that dissipate above the disk midplane, then these structures must be of large enough scale to survive the buoyant rise without being shredded by the turbulence within the disk. Thus, the results presented in this paper provide support to the prevailing paradigm for X-ray emission in accreting systems which involves an optically thin, hot corona powered by the dissipation of magnetic fields (e.g., Haardt & Maraschi 1993; Field & Rogers 1993).