LOCALITY OF MHD TURBULENCE IN ISOTHERMAL DISKS

Most of our models are evolved using ZEUS (Stone & Norman 1992). ZEUS is an operator-split, finite difference scheme on a staggered mesh. It uses artificial viscosity (not an anomalous viscosity!) to capture shocks. For the magnetic field evolution ZEUS uses the Method of Characteristics-Constrained Transport (MOC-CT) scheme, which is designed to accurately evolve Alfvén waves (MOC) and also to preserve the $\mbox{\boldmath $\nabla $}\cdot \mbox{\boldmath $B$}= 0$ constraint to machine precision (CT).

We have modified ZEUS to include "orbital advection" (Masset 2000; Gammie 2001; Johnson & Gammie 2005) with a magnetic field (Johnson et al. 2008). Advection by the orbital component of the velocity $\mbox{\boldmath $v$}_{{\rm orb}}$ (which may be supersonic with respect to the grid) is done using interpolation. With this modification the time step condition $\Delta t < {\it C} \Delta x/(|\delta \mbox{\boldmath $v$}| + c_{\mbox{max}})$ (c_max≡ maximum wave speed and C≡ Courant number) depends only on the perturbed velocity $\delta \mbox{\boldmath $v$}= \mbox{\boldmath $v$}- \mbox{\boldmath $v$}_{{\rm orb}}$ rather than $\mbox{\boldmath $v$}$. So when $|\mbox{\boldmath $v$}_{\mbox{orb}}| \gtrsim c_{\mbox{max}}$ (for shearing box models with v²_A/c²_s ≲ 1, when L ≳ H) the time step can be larger with orbital advection, and computational efficiency is improved.

Orbital advection also improves accuracy. ZEUS, like most Eulerian schemes, has a truncation error that increases as the speed of the fluid increases in the grid frame. In the shearing box without orbital advection the truncation error would then increase monotonically with |x|. Orbital advection reduces the amplitude of the truncation error and also makes it more nearly uniform in |x| (Johnson et al. 2008).

Do our results depend on the algorithm used to integrate the MHD equations? To find out, we have also evolved a subset of models using ATHENA, a second-order accurate Godunov scheme that solves the equations of ideal MHD in conservative form. The algorithm couples the dimensionally unsplit corner transport upwind (CTU) method of Colella (1990) with the third-order in space piecewise parabolic method (PPM) of Colella & Woodward (1984), and a constrained transport (CT) algorithm for preserving the $\mbox{\boldmath $\nabla $}\cdot \mbox{\boldmath $B$}= 0$ constraint. Details of the algorithm and test problems are described in Stone et al. (2008). The specific application of ATHENA to the shearing box (as used in this work) is described in Simon et al. (2009).⁵

We now consider a set of zero net field shearing box models to introduce our correlation function analysis. We use the same model parameters as FP07 so that these models also serve, by comparison with FP07, as a nonlinear code test.

The models in this section have size (L_x, L_y, L_z) = (1, π, 1)H. The initial magnetic field is B_z = B_z0 × sin(2πx/H), where B_z0 satisfies β ≡ 8πP₀/B²_z0 = 400. Noise is introduced in the initial velocity field to stimulate the growth of the unstable modes. The models are evolved to t_f = 600 Ω⁻¹.

We consider four resolutions: (N_x, N_y, N_z) = N(32, 50, 32), where N = 1, 2, 4, 8. The last three models correspond to runs std32, std64, and std128 in FP07, respectively. The evolution of 〈E_B〉 ≡ 〈B²/(8πρ_oc²_s)〉 is shown in Figure 1 (〈〉≡ volume average). The saturation 〈E_B〉 decreases as the resolution increases.

Zoom In Zoom Out Reset image size

The dimensionless shear stress

$$\begin{equation} \alpha \equiv {\left\langle \rho v_x \delta v_y - {B_x B_y\over {4\pi }} \right\rangle \over {\langle \rho \rangle c_s^2}}. \end{equation} \tag{ 3 }$$

We measured the time average of α, denoted as $\overline{\alpha }$, (from tΩ = 250 to tΩ = 600) for each run, and recorded the results in Table 1.⁶ These averages are nearly identical to those obtained by FP07. This consistency enhances confidence in both sets of results.

Table 1. Shearing Box Runs with a Zero Net Vertical Field

Model	Resolution	$\overline{\alpha }$	λB,min	λB,maj	λB,z	θB,tilt	λv,min	λv,maj	λv,z	θv,tilt	λρ,min	λρ,maj	λρ,z	θρ,tilt
z32	32 × 50 × 32	3.8 × 10−3	0.090	0.62	0.080	11	0.078	0.57	0.13	5.5	0.076	0.82	0.39	8.6
z64	64 × 100 × 64	4.2 × 10−3	0.059	0.38	0.050	13	0.074	0.45	0.18	7.1	0.056	0.60	0.33	6.6
z128	128 × 200 × 128	2.1 × 10−3	0.037	0.22	0.032	14	0.053	0.32	0.11	6.7	0.043	0.62	0.33	6.4
z256	256 × 400 × 256	1.1 × 10−3	0.024	0.17	0.024	14	0.035	0.24	0.10	5.9	0.019	0.34	0.18	5.8

Download table as:  ASCIITypeset image

For N_x ⩾ 64

$$\begin{equation} \overline{\alpha } \simeq 0.0021 \left({{N_x}\over {128}}\right)^{-1} \end{equation} \tag{ 4 }$$

is a good fit to the numerical results. The magnetic field energy density, α, and the kinetic energy density are almost inversely proportional to N_x, and so, like FP07, we conclude that the zero net field models do not converge.

All shearing box models considered in this paper have $\overline{\alpha } \simeq \overline{\langle E_B\rangle }/2 = 1/(2\overline{\beta })$. This implies a characteristic orientation of the field, since (neglecting the Reynolds stress ρv_xv_y ∼ 0.25 × [ − B_xB_y]/[4π]) $\overline{\alpha } \approx - \overline{\langle B_x B_y\rangle }/(4\pi \langle \rho \rangle c_s^2) \equiv \overline{\langle B^2 \cos \theta _B \sin \theta _B \rangle }/(4\pi \langle \rho \rangle c_s^2)$, where θ_B is the angle between the field and the y-axis. Then $\overline{\alpha } = \overline{\langle E_B\rangle } 2 \cos \overline{\theta }_B \sin \overline{\theta }_B = \overline{\langle E_B\rangle }/2$. So $\overline{\theta }_B = \pi /12$ (15°) is the characteristic angle between the magnetic field and the y-axis.

Next, we turn to the structure of the zero net field turbulence. Consider the two-point correlation function for the density fluctuations

$$\begin{equation} \xi _{\rho } \equiv \langle \delta \rho (\mbox{\boldmath $x$}) \delta \rho (\mbox{\boldmath $x$}+ \Delta \mbox{\boldmath $x$}) \rangle \end{equation} \tag{ 5 }$$

(δρ ≡ ρ − 〈ρ〉), for the trace of the velocity fluctuation correlation tensor

$$\begin{equation} \xi _{v} = \langle \delta v_i(\mbox{\boldmath $x$}) \delta v_i(\mbox{\boldmath $x$}+ \Delta \mbox{\boldmath $x$})\rangle \end{equation} \tag{ 6 }$$

(δv_i ≡ v_i − v_i,orb − 〈v_i〉), and there is an implied summation over i and for the trace of the magnetic field correlation tensor

$$\begin{equation} \xi _B = \langle \delta B_i(\mbox{\boldmath $x$}) \delta B_i(\mbox{\boldmath $x$}+ \Delta \mbox{\boldmath $x$}) \rangle \end{equation} \tag{ 7 }$$

(δB_i ≡ B_i − 〈B_i〉). All correlation functions are calculated for fluctuating dynamical variables with zero mean. Figure 2 shows slices of correlation functions through ξ at z = 0 for run z128. The cores of the correlation functions are ellipsoidal, with three principal axes, and concentrated at $|\Delta \mbox{\boldmath $x$}| < H$. The correlations are localized.

Zoom In Zoom Out Reset image size

We measure four features of the correlation functions: the angle θ_tilt between the correlation function major axis and the y-axis, and the correlation lengths along the major, minor, and z-axes (λ_maj, λ_min, λ_z), where the correlation length λ_i is defined ⁷ by

$$\begin{equation} \lambda _i \equiv {1\over {\xi (0)}} \int _{0}^{\infty } \xi (\hat{\mbox{\boldmath $x$}}_i l) dl. \end{equation} \tag{ 8 }$$

l is the distance from $\Delta \mbox{\boldmath $x$}= 0$ along the principal axis defined by the unit vector $\hat{\mbox{\boldmath $x$}}_i$, and $\hat{\mbox{\boldmath $x$}}_{\rm maj} = \hat{\mbox{\boldmath $x$}} \sin \theta _{\rm tilt} - \hat{\mbox{\boldmath $y$}} \cos \theta _{\rm tilt}$, $\hat{\mbox{\boldmath $x$}}_{\rm min} = \hat{\mbox{\boldmath $x$}} \cos \theta _{\rm tilt} + \hat{\mbox{\boldmath $y$}} \sin \theta _{\rm tilt}$, and $\hat{\mbox{\boldmath $x$}}_z = \hat{\mbox{\boldmath $z$}}$.^8,9 Figure 3 shows ξ_B along each of the principal axes in run z128. The dotted lines show ξ(0)exp(−l/λ_i). The correlated regions in the magnetic field are narrow (λ_min ≃ λ_z ≃ λ_maj/6) filaments with a trailing spiral orientation.

Zoom In Zoom Out Reset image size

What do the correlation functions mean? For the magnetic field, there is a characteristic orientation of the field $\overline{\theta }_B$ obtained through our measurement of the shear stress. The major axis of the correlation function is very nearly parallel to this, $\theta _{\rm tilt} \simeq \overline{\theta }_B$. It is reasonable to view the magnetic field correlation function, then, as tracing out a characteristic, filamentary structure in the magnetic field.

It is worth recalling briefly what we might expect for ξ_v in isotropic, homogeneous turbulence with an outer scale L₀ = 2π/k₀ and velocity dispersion σ_v. In the inertial range v²_k ∝ k^−11/3, so a reasonable functional form for the power spectrum is

$$\begin{equation} v_k^2 = N \left(\sigma _v^2 \over {k_0^3}\right) {1\over {(1 + (k/k_0)^2)^{11/6}}}, \end{equation} \tag{ 9 }$$

where N is a nondimensional normalization constant. The corresponding autocorrelation function is

$$\begin{equation} \xi _v = \sigma _v^2 {\sqrt{3} \Gamma \big(\frac{2}{3}\big)\over {2^{1/3}\pi }} (k_0 r)^{1/3} {K}_{1/3}(k_0 r), \end{equation} \tag{ 10 }$$

where K is the modified Bessel function and r ≡ |Δx|. For k₀r ≪ 1,

$$\begin{equation} \xi _v \approx \sigma _v^2 (1 - 0.955 \, (k_0 r)^{2/3} + O(r^2)). \end{equation} \tag{ 11 }$$

This is the usual r^2/3 Kolmogorov dependence at small separation. For k₀r ≫ 1,

$$\begin{equation} \xi _v \approx \sigma _v^2 \,\, 0.743 \,\, (k_0 r)^{-1/6} \, e^{-k_0 r}, \end{equation} \tag{ 12 }$$

which yields the expected decorrelation for k₀r ≫ 1. The correlation length defined in Equation (8) is 0.838/k₀. Note that the power spectrum does not have zero power as k → 0; rather it asymptotes to Nσ²_vk⁻³₀.

For MHD turbulence in disks, however, the correlation function is not isotropic, and its structure is not anticipated by any predictive theory. In the absence of such a theory it may be useful to have a convenient analytical representation of the numerical results. This can be obtained by stretching Equation (10) along each of the principal axes, that is, by replacing k₀r by $u \equiv ((\Delta \mbox{\boldmath $x$}\cdot \hat{\mbox{\boldmath $x$}}_{\rm min}/\lambda _{\rm min})^2 + (\Delta \mbox{\boldmath $x$}\cdot \hat{\mbox{\boldmath $x$}}_{\rm maj}/\lambda _{\rm maj})^2 + (\Delta \mbox{\boldmath $x$}\cdot \hat{\mbox{\boldmath $x$}}_{\rm z}/\lambda _{\rm z})^2)^{1/2}$.

Do the correlation lengths converge? We find that the correlation lengths are resolution dependent, with

$$\begin{equation} (\lambda _{\rm min}, \lambda _{\rm maj}, \lambda _{z}) \simeq (0.04, 0.24, 0.03) \, H \, \left({{N_x}\over {128}}\right)^{-2/3}; \end{equation} \tag{ 13 }$$

λ_z and λ_min are at most six zones. The scaling of λ_min with zone size is clearly not linear, but the 2/3 power-law scaling is just a fit to the data and should not be taken too seriously. The nonlinear scaling does hint at the possibility that, as resolution is increased and λ_min and λ_z are better resolved, there could be a transition in the outcome.

The major axis for ξ_B lies ∼15° from the y-axis. For the density and peculiar velocity field the tilt angle is ∼7°. The latter tilt is consistent with the measured Reynolds stress: 〈ρv_xv_y〉/〈ρv²〉 ∼ 1/8, so the average perturbed velocity is tilted at ∼7° to the y-axis.

Finally, note that there are low-amplitude features in ξ_v and ξ_ρ at scales of a few correlation lengths (see particularly in Figure 4(a)). These features may be due to the excitation of rotationally modified sound waves by MHD turbulence. To test this hypothesis note that, for tightly wrapped (k_y ≪ k_x), linear sound waves δρ²_k/ρ²₀ ≃ δv²_k/c²_s. A field composed of these waves would then have ξ_ρ/ρ²₀ = ξ_v/c²_s. Taking the differential correlation function ξ_v/c²_s − ξ_ρ/ρ²₀ should therefore remove those pieces of ξ_v that are due to sound waves. Figure 4(b) shows ξ_v/c²_s − ξ_ρ/ρ²₀ at Δy = 0. Evidently much of the large-scale power is removed. Figure 5 shows another slice through ξ_v/c²_s − ξ_ρ/ρ²₀ at Δz = 0. Again the large-scale power is removed. This is consistent with the hypothesis that the largest scale features in the correlation functions are acoustic waves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To test convergence we use the same size models as the zero net field runs, (L_x, L_y, L_z) = (1, π, 1)H. 〈B_y〉 is set so that β = 400; in the initial conditions all other magnetic field components vanish. The models are evolved to t = 250 Ω⁻¹. We use five different resolutions: (N_x, N_y, N_z) = N(32, 50, 32), where N = 1, 2, 4, 6, 8. In each run we average over the second half of the evolution to measure $\overline{\langle E_B\rangle }$ and $\overline{\alpha }$. We also measure the correlation lengths from ξ_ρ, ξ_v, and ξ_B using an average of the correlation function calculated from eight data dumps in the second half of the run. Parameters and results for these runs are listed in Table 2.

Table 2. Shearing Box Runs with a Net Azimuthal Field

Model	Algorithm	Resolution	$\overline{\alpha }$	$\overline{\langle E_B \rangle }/\rho _{0}c_s^2$	λB,min	λB,maj	λB,z	θB,tilt	λB,min/Δx
y32a	ZEUS	32 × 50 × 32	0.0094	0.019	0.10	0.40	0.084	12	3.2
y64a	ZEUS	64 × 100 × 64	0.014	0.024	0.066	0.36	0.064	15	4.2
y128a	ZEUS	128 × 200 × 128	0.015	0.028	0.053	0.28	0.049	15	6.8
y192a	ZEUS	192 × 300 × 192	0.020	0.037	0.055	0.30	0.049	15	11
y256a	ZEUS	256 × 400 × 256	0.021	0.041	0.049	0.27	0.045	15	13
y32b	ATHENA	32 × 50 × 32	0.018	0.032	0.11	0.63	0.09	15	3.3
y64b	ATHENA	64 × 100 × 64	0.015	0.027	0.070	0.38	0.060	16	4.5
y128b	ATHENA	128 × 200 × 128	0.018	0.035	0.058	0.33	0.051	16	7.4
y192b	ATHENA	192 × 300 × 192	0.025	0.050	0.052	0.30	0.047	17	10
y256b	ATHENA	256 × 400 × 256	0.027	0.055	0.053	0.32	0.049	16	14

Download table as:  ASCIITypeset image

Figure 6 shows 〈E_B〉(t) for various resolutions. In our two highest resolution runs, $\overline{\langle E_B\rangle } = 3.7\times 10^{-2} \rho _{0} c_s^2$ for run y192a and 4.1 × 10⁻²ρ₀c²_s for run y256a, respectively. A consistent fit is $\overline{\langle E_B\rangle } \simeq 0.03 (N_x/128)^{1/3}$ for 32 < N_x < 256. The saturation energy increases with resolution. It is unlikely that this trend continues indefinitely. As we will see below, the magnetic field correlation lengths are unresolved at N_x = 32 but resolved at N_x = 256. This suggests that the increase in $\overline{\langle E_B\rangle }$ is caused by resolution of magnetic structures near the correlation length. If there is little energy in structures much smaller than the correlation length, then $\overline{\langle E_B\rangle }$ should saturate at somewhat higher resolution.

Zoom In Zoom Out Reset image size

To check the algorithm dependence of the results we ran an identical set of models using ATHENA. The results are listed in Table 2. Both sets of models show a weak upward trend in $\overline{\langle E_B\rangle }$ with resolution. If one corrects for the approximately 2 times higher effective resolution of ATHENA then the ATHENA and ZEUS results are quantitatively consistent with each other.

The correlation lengths for the zero net field runs do not converge. What about the net azimuthal field models? The magnetic field correlation lengths are listed in Table 2 (other correlation lengths are omitted for brevity, but they behave similarly). For N ⩾ 4 both ZEUS and ATHENA find

$$\begin{equation} (\lambda _{\rm min}, \lambda _{\rm maj}, \lambda _{z}) \simeq (0.05, 0.32, 0.05) \, H. \end{equation} \tag{ 14 }$$

For the N ⩾ 4 models λ_min ∼ λ_z > 8Δx. In the highest resolution (N = 8) ATHENA model, λ_min/Δx ≃ 14. The correlation lengths are both converged and resolved.

If MHD turbulence in disks has a forward energy cascade, as do three-dimensional hydrodynamic turbulence and the Goldreich–Sridhar model for strong MHD turbulence in a homogeneous medium, then it is natural to identify the correlation lengths with the outer, or energy injection, scale. Since λ_min ≃ 15 grid zones at our highest resolution there is no resolved inertial range. We anticipate that future, higher-resolution numerical experiments with a mean azimuthal field will show the development of an inertial range.

At what scale is magnetic energy generated in MRI-driven turbulence? To investigate this, we have studied the correlation function for each term driving the evolution of the volume-averaged magnetic energy:

$$\begin{equation} \dot{E_B} = - \left\langle \mbox{\boldmath $\nabla $}\cdot \left(\frac{1}{2}B^2{\mbox{\boldmath $v$}}\right)\right\rangle - \left\langle \frac{1}{2} B^2 \mbox{\boldmath $\nabla $}\cdot {\mbox{\boldmath $v$}}\right\rangle + \langle {\mbox{\boldmath $B$}} \cdot ({\mbox{\boldmath $B$}} \cdot \mbox{\boldmath $\nabla $}{\mbox{\boldmath $v$}})\rangle - D, \end{equation} \tag{ 15 }$$

where D is the volume-averaged numerical dissipation rate. The terms on the right-hand side can be interpreted as describing the effects of advection, compression and expansion, field-line stretching, and numerical dissipation. On average the first term is small (it should vanish exactly for shearing box boundary conditions, but roundoff and truncation error make it nonzero), and the third term dominates the second by a factor of 20 in run y128b. In a time- and volume-averaged sense the right-hand side must be zero, so numerical dissipation must approximately balance energy injection by field-line stretching.

Previous studies (FP07; Simon et al. 2009) analyzed a version of this equation in the Fourier domain to study turbulent energy flow as a function of length scale in the shearing box simulations. Both studies found that the magnetic energy is generated on all scales by the background shear. Here, we perform a complementary analysis in the spatial domain.¹⁰

We have computed the autocorrelation function for the field-line stretching term from run y128b (without subtracting the mean). Figure 7 shows the autocorrelation function at Δz = 0. From this figure we conclude that: (1) the scale and shape of the correlation function is similar to that of the fundamental variables ($\mbox{\boldmath $B$}$, $\mbox{\boldmath $v$}$; recall that in y128b (λ_B,min, λ_B,maj, λ_B,z) = (0.058, 0.33, 0.051)H); (2) energy is injected at scales comparable to the correlation length of the fundamental variables; (3) a superposition of magnetic structures similar to ξ_B that are distributed with uniform probability in space would have a power spectrum that is flat (white noise) at low k. It is plausible that terms in the Fourier-transformed magnetic energy equation would also be flat at low k. This would not imply that energy is injected by dynamically meaningful structures at large scales; it would simply be the consequence of an uncorrelated superposition of small, localized features in the turbulence.

Zoom In Zoom Out Reset image size

Does the saturation energy or correlation length depend on the size of the computational domain (box size)? To investigate, we fix the physical resolution at 64 zones/H and vary the model size: (L_x, L_y, L_z) = (1, π, 1)H, (2, π, 1)H, (1, 2π, 1)H, (2, 2π, 1)H, and (4, 4π, 1)H. The model parameters and outcomes are listed in Table 3.

Table 3. Shearing Box Runs with a Net Azimuthal Field: Effect of the Box Size

Model	Size	$\overline{\langle E_B \rangle }/\rho _{0}c_s^2$	${\lambda _{B,{\rm maj}}} \over {H}$
y64	(1, π, 1)H	0.024	0.36
y64.x2	(2, π, 1)H	0.028	0.49
y64.y2	(1, 2π, 1)H	0.035	0.45
y64.x2y2	(2, 2π, 1)H	0.038	0.49
y64.x4y4	(4, 4π, 1)H	0.038	0.57

Download table as:  ASCIITypeset image

Evidently, there is a weak dependence of $\overline{\langle E_B\rangle }$ on box size; it increases from 0.024ρ₀c²_s for the smallest run to 0.038ρ₀c²_s for the largest run. The magnetic field correlation lengths also increase with box size, with λ_maj = 0.36H for the smallest box (so L_y/(2λ_maj) = 4.4) to λ_maj = 0.57H for the largest box (so L_y/(2λ_maj) = 11). This upward trend in correlation length is probably real, but it is sufficiently small that it is difficult to separate from noise in the correlation length measurements.

The correlation functions ξ_ρ and ξ_v, unlike ξ_B, have low amplitude tails extending out to the box size. This can be seen in Figure 8, which shows ξ_ρ and ξ_B for run y64.x4y4. The tails are likely due to sound waves, and their absence in the differential correlation function ξ_v/c²_s − ξ_ρ/ρ²₀, shown in the middle panel of Figure 8, is consistent with this.

Zoom In Zoom Out Reset image size

We now compare two models that differ only in their initial field strength. Both have a size L_x, L_y, L_z = (1, 2π, 1)H with resolution N_x, N_y, N_z = 64, 200, 64. One model has the same initial field strength as our other models, β₀ = 400, while the other one starts with a stronger field, β₀ = 100.

We found that the saturated magnetic energy for the β₀ = 400 run is $\overline{\langle E_{B}\rangle } = 0.035\rho _{0}c_s^{2}$. For the β₀ = 100 run $\overline{\langle E_{B}\rangle } = 0.079\rho _{0}c_s^{2}$, slightly more than twice the saturation level of the higher β₀ run. Resolution may be playing a role here: the β₀ = 100 run has twice the resolution per most unstable MRI wavelength, and we know that the saturation level depends on resolution when the field strength is constant.

Our results are consistent with the linear relation between $\overline{\langle E_B\rangle }$ and initial field strength for 〈B_y〉 ≠ 0 models reported in Hawley et al. (1995) (hereafter HGB95; although it may also be consistent with a wide range of exponents for this relation). Our results are inconsistent with HGB95's claim that $\overline{\langle E_B\rangle } \propto L_y$, at least if the box size is ≫λ_maj (the good agreement found for HGB95's predictor may be a coincidence).

At the level we can determine from two data points, our results are consistent with $\overline{\langle E_B\rangle } \propto \rho _0 c_s V_{A,y0}$, where V_A,y0 is the initial azimuthal Alfvén speed (i.e., here scale height c_s/Ω replaces L_y in HGB95; there are no other length scales in the problem). This is interesting: it implies that α depends on the gas pressure, a result first reported by Sano et al. (2004) and thus compressibility plays a role in the saturation of the MRI!

Our results suggest that $\overline{\langle E_B\rangle }$ should scale differently in compressible and incompressible models. In the incompressible models the only length scale available is the size of the box, so in incompressible models we must have $\overline{\langle E_B\rangle } \sim \rho _0 (L\Omega)^a V_{A,y0}^b$, where L is some combination of L_x, L_y, and L_z, and the exponents a and b are not determined.

The correlation lengths for the magnetic field in the two runs are (λ_min, λ_maj, λ_z) ≃ (0.08, 0.45, 0.08)H for the β₀ = 400 run and (λ_min, λ_maj, λ_z) ≃ (0.11, 0.58, 0.10)H for the β₀ = 100 run. To sum up, the correlation length increases weakly as the initial field strength and the box size increases. It is not consistent with the scaling ∼B_y,0/(ρ₀Ω) one would expect if the correlation length scaled with the most unstable wavelength of the background field, and it is not consistent with the scaling ∼〈B²_y〉^1/2/(ρ₀Ω) ∼ B^1/2_y,0 one would expect if the correlation length is related to a characteristic MRI length scale for the (larger) fluctuating field.