MAGNETOROTATIONAL TURBULENCE TRANSPORTS ANGULAR MOMENTUM IN STRATIFIED DISKS WITH LOW MAGNETIC PRANDTL NUMBER BUT MAGNETIC REYNOLDS NUMBER ABOVE A CRITICAL VALUE

Figure 2 shows that the dependence of α_SS on Pm is not well described by a single power law at a given Re, as had been previously claimed for unstratified shearing boxes by Lesur & Longaretti (2007) and Fromang et al. (2007). The points are the mean α_SS for all times t > 20 t_orb, and the error bars represent the standard deviation over that range. These runs all use the VF boundary conditions. The figure shows evidence of a cutoff that moves to higher Pm as Re decreases. This is indicative of a critical Rm, most visible for Re = 3200 (center left) and Re = 9600 (lower left). Above Rm_crit, it appears that α_SS is consistent with a constant value. These data show that the behavior of the stratified, zero-net-flux MRI system has two distinct regimes, controlled by Rm_crit. When Rm < Rm_crit, the transport is strongly dependent on Pm; when Rm > Rm_crit, the transport is independent of Pm. Given that astrophysical disks have typical Re ≳ 10¹⁶, this result means that the MRI is capable of robust angular momentum transport even at low Pm, so long as Rm > Rm_crit. We discuss these estimates in more detail in Section 5 and suggest that this criterion is easily met in most disks.

Zoom In Zoom Out Reset image size

We recast these results in the (Rm, Re) plane in Figure 3. This figure makes clear the fact that along with a reduced but still present Pm dependence, there is a fairly clear critical Rm: transport is significantly reduced left of the vertical line near Rm ∼ 3000. As the simulations approach Rm_crit, the α_SS variance increases significantly.

Zoom In Zoom Out Reset image size

Figure 4 shows that the MRI appears to operate in two distinct dynamo states, one in which regular 〈B_y〉 cycles appear with a characteristic period of τ_B ∼ 10 t_orb, and one with irregular variations with timescales on the order of ∼50 t_orb. The key control parameter appears to be Rm, as the second, third, and fourth panels from the top show the regular behavior despite being at different Pm. There is no discernible trend in cycle period with Pm, though the top panel, with Rm = 12, 800, Pm = 2 seems to show an intermediate behavior. However, all runs with Rm ≲ 3200 conclusively show the irregular behavior. There appears to be a secondary Pm effect as well, since the two Rm = 3200 models show different behavior depending on Pm (third and fifth panels from the top).

Zoom In Zoom Out Reset image size

These two drastically different behaviors complicate efforts to ascertain scaling properties of the MRI as a function of dimensionless parameters. The range of available Pm is limited by numerical resolution both from above and below: large values of Pm imply a viscous scale much larger than the resistive one, while small Pm implies the opposite. Since both length scales must fit on (and be resolved by) the grid, and indeed neither can be so large as to stabilize the largest linear MRI modes that fit in the simulation box, our parameter range is necessarily limited. Furthermore, since low Rm runs transition to a very different behavior with much larger cycle period and very different turbulent properties in a discontinuous fashion, our range of available Pm space for the regular dynamo mode is further limited.

Thus, in what follows, we briefly describe the irregular dynamo before focusing on connecting the details of the MRI turbulence in the regular region of parameter space to better-established results on small- and large-scale dynamo action. The MRI is known to produce both small- and large-scale dynamo action, the latter typically requiring stratification or open boundary conditions (though Lesur & Ogilvie 2008 demonstrate a large-scale dynamo with neither). The small-scale dynamo for driven, isotropic, homogeneous, and incompressible turbulence has Rm_crit that increases with decreasing Pm (Schekochihin et al. 2005), and thus it behooves us to understand how MRI turbulence fits into this picture.

When Rm drops below Rm_crit, the mean-field dynamo switches to an irregular regime, showing quasi-periodic magnetic cycles, with longer cycle lengths (one hesitates to call them periods) ∼40 t_orb (see the lower two panels of Figure 4). This is consistent with Figure 15 of Davis et al. (2010), who also showed long cycle lengths in highly resistive simulations. In this regime, the turbulence often appears only in one half-plane, either z > 0 or z < 0, sometimes staying that way for nearly the duration of the simulation. This kind of behavior at first appears unphysical, as the stratified shearing-box system, while having odd parity about the midplane, should not necessarily damp perturbations of one helicity more than the other. Indeed, in all simulations, the kinetic helicity shows erratic fluctuations of sign with respect to the midplane.

However, this is explainable in a simple α − Ω treatment. The only prerequisite for this explanation is that the two half-planes be dynamically decoupled from one another. For the case of VF boundary conditions, B_y(z = ±L_z)↦0. The largest wavenumber modes compatible with this boundary condition are B_y∝sin (2π/L_zz) and B_y∝sin (π/L_zz), corresponding to k = 1 and 1/2, respectively. The α − Ω dynamo operates when the mean induction equation takes the form

$$\begin{equation} \partial _t \langle \mathbf {B} \rangle = \alpha \langle \mathbf {B} \rangle - q \Omega \langle B_x \rangle \mathbf {\hat{y}} + \eta \nabla ^2 \langle \mathbf {B} \rangle. \end{equation} \tag{ 2 }$$

If during a cycle period the MRI turbulence shuts off, the α effect that it produces will cease, leaving a mean induction equation that looks like

$$\begin{equation} \partial _t \langle B_x\rangle = \eta \nabla ^2 \langle B_x\rangle \end{equation} \tag{ 3 }$$

for the x component and

$$\begin{equation} \partial _t \langle B_y \rangle = -q \Omega \langle B_x \rangle +\eta \nabla ^2 \langle B_y \rangle \end{equation} \tag{ 4 }$$

for the y component. Because the decay time for modes is roughly t_decay ≃ 1/k²η, small-scale structure will undergo selective decay, leaving only the largest scale modes, the decay time of which is ∼108 t_orb for Rm = 1600. At this stage, 〈B_y〉 has about this much time to grow linearly via stretching of any residual 〈B_x〉 to amplitudes at which non-axisymmetric MRI can re-establish fluid turbulence and hence an α effect. If the two half-planes are not strongly coupled, it is possible that the turbulence will die out in one half of the box first, leading to an α effect only in the upper or lower midplane. This scenario is essentially the same as the one demonstrated in the midplane by Simon et al. (2011).

Understanding the origin of the scaling of angular momentum transport with magnetic Prandtl number requires understanding the underlying field generation mechanism. Gressel (2010) has demonstrated that the field patterns present in the MRI can be explained in terms of a mean-field dynamo model that includes dynamical α-quenching, a mechanism that modulates α by enforcing magnetic helicity conservation as the mean field grows (see Section 4.2 for more details). Here, we decompose fields from the simulations into their mean and fluctuating components, $\mathbf {B = \bar{B} + b}$, where $\mathbf {\bar{B}}$ is a horizontal (x–y) average and $\mathbf {b}$ is the fluctuating field. We denote full box averages as $\langle\mathbf {B}\rangle$. In the standard α − Ω dynamo mechanism, the α effect of isotropic, helical turbulence generates a poloidal field from toroidal fields, which are in turn sheared out by differential rotation Ω and regenerate a toroidal field, thus leading to exponential amplification. Based on the early work of Brandenburg et al. (1995), the traditional α − Ω scenario can explain the observed periodic dynamo generation and propagation of $\bar{B_y}$ away from the disk midplane if the α term has the opposite sign of that expected from rotating, stratified turbulence (the strong Keplerian shear has no problem stretching poloidal field to generate toroidal field; the issue at hand is generating the former). However, that expected sign was derived from an analysis that only assumed the presence of a kinetic α_K effect from the helicity of the fluid turbulence. Once the magnetic α_M from the Lorentz force is also included (Pouquet et al. 1976; see also Section 4.2), the total α = α_K + α_M is dominated by $\alpha _{M} = 1/3 \tau \mathbf {\langle \nabla \times v_a \cdot v_a \rangle}$, where τ is a typical turbulent correlation time and $\mathbf {v_a = b/\sqrt{4 \pi {\bm \rho} }}$. The total α has the required sign (Gressel 2010). While we similarly find that α_M ≃ 10α_K, our results for the z profile of α_M (and thus the total α) appear to contradict those of Gressel (2010). Figure 5 shows the α_M profile for three simulations with Pm = 1, 4 and Re = 3200, 6400, 12, 800. There is no monotonic trend with Pm, and there are some differences in shape, but the overall profile is negative in the upper plane (z > 0) and positive in the lower, opposite to that found by Gressel (2010). However, it is worth noting that both the Re = 12, 800 and the Re = 3200 runs show a reversal of the α_M(z) profile near the midplane, just as found by Gressel (2010), though with the opposite sign.

Zoom In Zoom Out Reset image size

The only significant differences between his study and ours are the vertical boundary conditions on the fluid, which are outflow in his case, and the fact that his domain covers 6H, while our domains typically cover 4H. We have run one simulation, with Re = 6400 and Pm = 4, with a vertical domain of 6H. Figure 6 shows α_M(z) for both the standard and extended domain for this model. The figure suggests a possible explanation for the discrepancy between our models and Gressel's. When we enlarge our domain, we see regions with |z| ≳ 2 showing a positive α_M effect, while the regions |z| ≲ 2 are roughly the same between the 4H and 6H models, outside of narrow boundary layers in both cases. Thus, our models just appear to show a stronger reversal of the sign of α_M near the midplane than does Gressel (2010). Given that buoyancy and Parker instabilities become more prominent with height, it seems that the most likely explanation for the discrepancy between our results and Gressel's is a combination of a reduced domain for our results and the outflow fluid boundary conditions in Gressel's work. The outflow boundary conditions allow buoyant fluid to escape out of the top of the box, while our closed, stress-free lids force that flow to recirculate. Thus, we do not consider our apparently discrepant results to actually be evidence against the standard α − Ω mechanism acting to drive the observed dynamo phenomenology.

Zoom In Zoom Out Reset image size

In order to place our results in the context of other recent studies, we begin by comparing our results to the unstratified results presented by Käpylä & Korpi (2011). These authors performed two series of unstratified shearing-box simulations of the MRI, both initialized with zero-net magnetic flux in the z-direction. One set of simulations had periodic boundary conditions on the magnetic field on the z boundary, the other had VF boundary conditions. Käpylä & Korpi (2011) show that in the latter case, the angular momentum transport coefficient α_SS is independent of Pm, while with periodic boundary conditions, α_SS∝Pm². We find that even with periodic boundary conditions, a stratified disk with high enough Rm has α_SS basically independent of Pm.

We are not the first to suggest a critical Rm for MRI turbulence. In the unstratified case, Fleming et al. (2000) and Sano & Stone (2002) both noted the existence of a sharp cutoff in α_SS when the resistivity exceeded a certain value. The latter study used the Elsasser number, Λ = v²_A0/ηΩ (though they refer to this as Rm). However, aside from using Λ, their results also measure the α_SS cutoff as a function of the initial magnetic field strength, for both net-flux and zero-net-flux initial fields. In our case, we do not make reference to the initial field strength; instead, by using the instantaneous value of Rm as our parameter, we measure the relative action of shear against dissipation. This is the appropriate measure for a nonlinear, self-sustained system such as the stratified, zero-net-flux MRI that we study here, because any dynamo-generated energy must ultimately come from the free energy in the shear flow.

More recently, a controversy regarding the role of channel modes in saturating the unstratified MRI when a net z field is present has led to similar results regarding the Rm_crit for MRI-driven transport. Pessah (2010) demonstrated the existence of a critical Λ below which the angular momentum transport by the MRI dropped off dramatically. Numerical simulations of the same problem by Longaretti & Lesur (2010), though seemingly disproving the Pessah (2010) saturation mechanism, also suggests a trend similar to the one we report, with two regimes separated by Pm ∼ 1. With Pm < 1, the transport scaling is dominated by Rm, while above it the scaling is mostly with Pm. Their primary aim was to elucidate the role of axisymmetric channel modes in the saturation of MRI turbulence in the presence of a net z field. In our study, there are no channel modes, as there is no net z field, and it is not too surprising that our results are slightly different. Indeed, the physics our control parameter is attempting to represent is not the growth rate and wavenumber of a (quasi-)linear channel mode, but instead a turbulent dynamo that tends to create a net toroidal field (y-field in our simulations), on which any MRI growth is non-axisymmetric and thus transient though with extremely high growth factors (Balbus & Hawley 1992).

We interpret our results in light of the ability of the system to build up a large-scale magnetic flux. Our results suggest that, with stratification, the MRI can act as a mean-field dynamo with any boundary conditions, unlike the unstratified case. In order to clarify terms, we will refer to a mean-field dynamo as any system capable of building magnetic fields at the lowest possible wavenumber in the box, either k = 1 in the periodic case or k = 0 (a true mean field) in the case including a VF boundary condition.

In the modern picture of mean-field dynamo theory, ensuring the conservation of magnetic helicity,

$$\begin{equation} H \equiv \langle \mathrm{A \cdot B} \rangle, \end{equation} \tag{ 5 }$$

provides an important constraint on the evolution of the mean field. Magnetic helicity conservation provides both a compelling physical mechanism for the existence of mean-field dynamos as well as quantitative predictions for saturated field strengths and the timescales necessary to reach those strengths. We briefly sketch out some background related to magnetic helicity conservation and then apply it to our MRI simulations.

The evolution equation for H can be written as

$$\begin{equation} \frac{{\it dH}}{dt} = -2 \eta \langle \mathbf {J \cdot B} \rangle - 2 \oint \left(\mathbf { A \times E} + \phi \mathbf {B} \right) \cdot \mathbf {\hat{n}} dS, \end{equation} \tag{ 6 }$$

where Ohm's law gives $\mathbf {E = -U \times B} + \eta \mathbf {J}$ and ϕ is an arbitrary gauge term in the definition of the vector potential (Brandenburg & Subramanian 2005). In our simulations we use the Kepler gauge, $\phi = \mathrm{\mathbf {U}^{\mathrm{kep}}\cdot A}$; see the Appendix as well as Brandenburg et al. (1995) for more details. The first term is due to resistivity acting on the current helicity C ≡ 〈J · B〉, while the second is a boundary flux term that can entirely determine the behavior of the MRI in unstratified simulations (Käpylä & Korpi 2011). In the stratified case, the vertical boundary conditions are considerably less important. This is evident in the periodic, stratified simulations of Davis et al. (2010), who find sustained turbulence for lower Pm than similar unstratified boxes. We note that our gauge choice eliminates all possible horizontal magnetic helicity fluxes (see Hubbard & Brandenburg 2011 for a detailed explanation), leaving only the possibility of vertical fluxes, which could be driven by a turbulent diffusivity or shear, as in the case of Vishniac & Cho (2001).

By using an eddy-damped, quasi-normal, Markovian closure scheme, Pouquet et al. (1976) demonstrated that current helicity drives a magnetic α_M effect, akin to the fluid α effect of classical mean-field dynamo theory. This α_M in turn leads to an inverse cascade of magnetic energy, and thus the build up of large-scale magnetic field. That inverse cascade can be seen simply as a result of the conservation of magnetic helicity. We outline this process following Brandenburg & Subramanian (2005). Consider a fully helical field, which has its scale-dependent relative helicity,

$$\begin{equation} \mathcal {H}(k) = \frac{k H(k)}{2 M(k)} = 1. \end{equation} \tag{ 7 }$$

H(k) and M(k) are the Fourier transforms of magnetic helicity and magnetic energy, respectively. $\mathcal {H} = 1$ implies kH(k) = 2M(k). The nonlinear terms allow coupling only among modes k and p satisfying p + q = k, where q is an arbitrary mediating wavenumber in the three-wave interaction, helicity is conserved for each mode, so pH(p) = 2M(p) and qH(q) = 2M(q). Since energy is conserved in the interaction, M(k) = M(p) + M(q), and thus pH(p) + qH(q) = kH(k). Because helicity must also be conserved, H(k) = H(p) + H(q), and

$$\begin{equation} k = \frac{p H(p) + q H(q)}{H(p) + H(q)}. \end{equation} \tag{ 8 }$$

If k is the final wavenumber, one of p or q is the starting wavenumber, and Equation (8) demands that k is less than or equal to the maximum of p or q. Since k is smaller than or equal to the parent wavenumber, this corresponds to an inverse cascade of magnetic energy. Pouquet et al. (1976) explicitly relate this to the current helicity, and construct a magnetic α_M term that back-reacts on the flow as this cascade proceeds.

This formalism gives a heuristic explanation for the existence of a mean-field dynamo driven by small-scale, helical flows, and incorporating the magnetic helicity constraint into mean-field dynamo models allows them to correctly predict the saturation behavior of α² dynamos (Brandenburg 2001). Here, we want to understand the MRI dynamo action in our simulations in light of the constraints imposed by helicity conservation. We explore both boundary terms and scale transfer of H.

We begin by testing the Pencil Code's numerical conservation of magnetic helicity. For a system with periodic boundary conditions, the second term in Equation (6) is zero, and the only contribution to changes in helicity can come from resistively limited current helicity fluctuations. Therefore, if our simulations numerically conserve helicity, a run with periodic boundary conditions should show helicity variations only on a resistive timescale. In order to establish that helicity conservation is robust in our simulations, we ran a simulation with Re = 3200 and Pm = 2, with periodic boundary conditions, and tracked the time evolution of magnetic and current helicities (upper panel of Figure 7). We used a simple forward Euler scheme to integrate dH/dt = −2ηC with volume average data from the simulations. The H(t) that results from this integration is shown in Figure 7 as the green triangles, while the H calculated directly during the run is given by the blue solid line. The agreement is quite good, considering the crudeness of the integration method and the sparseness of the data (it is only sampled at intervals of 100 time steps). As a result, we are reasonably confident that our simulations conserve magnetic helicity.

Zoom In Zoom Out Reset image size

Moving to the case with VF boundary conditions, the lower panel of Figure 7 again shows the current helicity and the now-gauge-dependent quantity H = 〈A · B〉. Because of the presence of the gauge in the second term of Equation (6), this quantity is not physically meaningful itself. However, it is well defined and provides a clue as to the behavior of the system. Once again, we integrate −2ηC using C in the same way. We expect that if magnetic helicity ejection occurs as a result of MRI turbulence, the value of H resulting from this integration will not track the actual value from the simulation, and indeed the lower panel of Figure 7 shows that it does not. Because of our gauge choice, we know that the flux of magnetic helicity flux density must be vertical, and this figure confirms that significant amounts of helicity escape. Furthermore, H fluctuates considerably more frequently in the VF case than in the run with periodic boundary conditions, showing that the timescale for variation of the global magnetic helicity decreases considerably when a helicity flux is allowed.

How does transport work in the periodic case, when there is no boundary flux? The only way for net helicity to change in this case is via resistive effects. Typically, resistively limited systems lead to catastrophic quenching of the dynamo effect, where the saturated magnetic energy is ∝Rm^−1/2 (Vainshtein & Cattaneo 1992). We have shown that this is not the case in the stratified MRI: boundary conditions do not make significant differences in the transport (Figure 8) and for VF boundaries, Figure 4 demonstrates that the saturated large-scale field strength does not decline strongly with increasing Rm. Vishniac & Cho (2001) suggested a potential solution to this situation: helicity may not need to be ejected entirely across a system boundary (as in, say, a coronal mass ejection in the solar dynamo). It could instead be transported spectrally, transferring from small to large scales. Indeed, Brandenburg (2001) shows that exactly this occurs for a helically driven turbulence simulation with no shear or rotation—the prototypical α² dynamo. However, in this case, while equipartition fields are built even for periodic boundary conditions, the timescale required to do so is t_sat∝η⁻¹: instead of catastrophic quenching for low resistivity, there is a catastrophic timescale problem instead (Brandenburg 2001). While our data are not conclusive on the relationship between Rm, Pm, and the cycle period of large-scale magnetic fields, it certainly does not suggest an inverse relationship between t_sat and Rm.

Zoom In Zoom Out Reset image size

The two runs in Figure 7, taken together, tell an intriguing tale: the stratified MRI can transport angular momentum and build large-scale magnetic energy even without a global helicity flux, but if one is allowed, the system does transport helicity across the domain boundary. Ultimately, as any boundary condition is in some sense an approximation of reality, we must understand the details of the accretion disk dynamo independent of the choice of such conditions.

The MRI is difficult to analyze in terms of a simple, two-scale mean-field model: it does not present a single energy injection scale at which we could expect small-scale helicity to be generated (Davis et al. 2010). This lack of clear scale separation makes the dynamo coefficients extracted via the test field method rather noisy and somewhat difficult to interpret, though it has been performed for the MRI by several groups (Brandenburg 2005; Gressel 2010). We forgo the test field method here, since we should be able to see a crude difference between the current helicity at the small and large scales if this is in fact what allows dynamo action, and thus transport. It is also clear from comparing Figures 7 and 4 that the timescale for building mean fields for the VF boundary condition case is not related in any obvious way to the flux of magnetic helicity. Furthermore, we do not see any significant changes in angular momentum transport with differing boundary conditions that do and do not allow a helicity flux. We address this in the next section.

We want to understand if scale transfer is responsible for dynamo action in closed, stratified MRI simulations. To do so, we compute the spectrum for the relative helicity, $\mathcal {H}$, on spherical shells in k-space, in runs with closed boundary conditions (i.e., periodic or perfect conductor). This ensures that magnetic helicity is gauge-independent by removing the surface terms in Equation (6).

In Figure 9, we show the helicity spectrum for Re = 3200, Pm = 2. The field is not strongly helical at any scale, reaching only $\mathcal {H} \sim 0.2$ at the smallest scales and remaining much lower on the largest scales. By contrast, the α² dynamo driven by a fully helical fluid forcing function has $\mathcal {H} \sim 1$ at all scales (Brandenburg 2001). In that case, something similar to the classic inverse cascade of magnetic energy due to helicity conservation (Frisch et al. 1975; Pouquet et al. 1976, and Section 4.2) occurs. In our case, however, we do not see a significant difference in the properties of the dynamo when the boundary conditions change between those that do not allow a flux of magnetic helicity and those that do, despite the fact that our results suggest that the MRI does in fact eject helicity when given the chance (see Figure 7). This resolves that observation: the MRI-generated field is not strongly helical, even when both kinetic and magnetic α effects occur. Thus, the constraints placed on the field by the conservation of magnetic helicity do not dominate its formation.

Zoom In Zoom Out Reset image size

Finally, the fact that the relative helicity is peaked at small scales suggests that helicity constraints might be more important in the unstratified MRI, as suggested by Käpylä & Korpi (2011). In that case, the generation of field is not affected by dynamo waves and magnetic buoyancy, and so the helical turbulence might wind the field into a helically limited state. The small scales in Figure 9, where coherent dynamo waves and magnetic buoyancy effects are less important, show a rise in helicity, consistent with this idea.

There remains some incongruity between our results for α_M(z) and those of Gressel (2010). The only significant difference between his simulations and ours is his use of outflow boundary conditions on the fluid (the magnetic field boundary conditions are identical) and his use of a slightly larger domain. By running a model with a larger vertical extent, we have shown that above the ±2H z-boundaries of our standard domain, the profile reverses and matches Gressel's. However, our results show a much stronger reversal near the midplane than his does, and this appears to be robust regardless of the domain size. The MRI coupled with outflow fluid boundary conditions leads to a magnetized wind, and this, together with the increased efficiency of the Parker instability with height, likely explains the remaining discrepancy. The role of fluid boundary conditions (which we have not varied here) should also be examined, in order to better understand how winds interact with a disk dynamo (e.g., Vishniac & Cho 2001).

The overall pattern of dynamo-generated fields in the z–t plane for the MRI is determined by the boundary conditions (Brandenburg et al. 1995), and this can be understood in terms of analytic mean-field theory (Brandenburg & Subramanian 2005). Nonetheless, Figure 8 shows that these various dynamo modes do not affect angular momentum transport. Thus, all of the boundary conditions studied here for the stratified MRI do an equally good job of generating magnetic flux for the MRI to continue to grow on. We have identified the fact that the magnetic fields are not fully helical at any scale as a likely reason why the boundary conditions do not determine the total angular momentum transport. Future work should address this point in more detail. Measuring the helicity spectrum in unstratified shearing boxes would determine if the dramatic differences between periodic and VF boundary conditions in that case are indeed related to a helicity-limited field evolution. Further theoretical work on the relationship between the mean-field dynamo mechanisms in unstratified domains (likely the incoherent-α effect of Vishniac & Brandenburg 1997) and stratified domains (the α − Ω effect) and their expected degree of helicity is also required in order to complete this picture of dynamo-mediated MRI transport in the zero-net-flux case.

We have demonstrated that understanding angular momentum transport via the MRI is strongly tied to the details of MHD turbulence and dynamo action. It is worth commenting briefly on some results for isotropic (non-shearing), non-rotating MHD turbulence, a much better studied system. Recently, several groups have demonstrated that non-local interactions in k-space are important at large scales in MHD turbulence, cross-coupling large-scale velocity fields with small-scale magnetic fields (e.g., Mininni et al. 2005; Lessinnes et al. 2009; Cho 2010). If such an analysis holds for the stratified MRI dynamo, it could explain the existence of a critical magnetic Reynolds number Rm_crit: if the large-scale velocity fields are coupled to the magnetic dissipation range, a much more efficient energy sink appears than if they are coupled to a magnetic inertial range. In the latter case, the large-scale velocities could contribute to small-scale helicity production, continuing to drive large-scale dynamo action. This could be verified by a shell-transfer analysis of the stratified MRI, which we will pursue in a future publication.

It is possible that our hyperdiffusive terms might be biasing our results. However, in a series of small-scale dynamo simulations, groups led by Schekochihin and Brandenburg have found no evidence that hyperviscosity plays a significant role (Schekochihin et al. 2005, 2007). There is a well-understood mechanism by which hyperresistivity causes an increase in saturated mean-field strengths for α² dynamos (Brandenburg & Sarson 2002). The main effect of hyperresistivity is simply to modify the timescales and length scales over which resistive effects occur, leading to differences from the magnetic helicity theory developed for Laplacian resistivity (Brandenburg 2001). Once the hyperresistive scales are properly accounted for, magnetic helicity conservation indeed correctly predicts the saturation field strength. Thus, our results are unlikely to be significantly affected by this mechanism, as we do not offer a detailed theory for the saturation amplitude for the MRI. Furthermore, the Brandenburg & Sarson (2002) simulations did not include a standard, Laplacian resistivity term, while ours do.

More directly, our resolution study also suggests that hyperdiffusion is not a dominant effect: for integrated quantities, our 64 zones/H and 128 zones/H runs are converged (see Figure 2 and Table 1). Because the hyperdiffusive coefficients are smaller for larger resolution, this suggests that the values of the coefficients make little difference. Nevertheless, we cannot rule out any systematic hyperdiffusive influence on the dynamo behavior without running simulations at high enough resolution to simultaneously resolve the MRI and dissipate energy fast enough at the grid scale using only regular Laplacian diffusion operators at all times and through the entire spatial domain. While this is practical for unstratified MRI simulations, because the distribution of u is relatively narrow, in stratified simulations where the distribution of velocities is much broader, it is hard to maintain stability without hyperdiffusivity. Direct numerical simulations of stratified MRI thus require very high resolution and are therefore extremely expensive and are not practical for the parameter study considered here.