SYMMETRIES, SCALING LAWS, AND CONVERGENCE IN SHEARING-BOX SIMULATIONS OF MAGNETO-ROTATIONAL INSTABILITY DRIVEN TURBULENCE

The magneto-rotational instability (MRI) is commonly invoked to explain the origin of turbulence in electrically conducting disks (Balbus & Hawley 1991). Although the instability itself is amenable to an analytic treatment, much of what is currently known about the nonlinear development of this instability is based on numerical simulations. In particular, the efficiency of MRI-driven turbulence at transporting angular momentum, which is the crucial quantity that ultimately controls the accretion rate, is known almost exclusively through numerical simulations. Because of the difficulties inherent in simulating an entire disk the majority of the numerical work has been based on the local model known as the shearing-box approximation (SBA; Hawley et al. 1995). The advantages of the SBA are numerous: the cylindrical geometry of the full disk is replaced with the Cartesian geometry of a rectangle, simple periodic and shear-periodic boundary conditions can be applied, and most importantly, whatever numerical resolution is available can be deployed to resolving the ensuing turbulence. However, a number of issues arise that have led several authors to question the applicability of the SBA to the study of disk turbulence. Regev & Umurhan (2008) have noted some inconsistencies in the formulations of the SBA with uniform background states; they have also questioned the assumption of locality that underlies the derivation of the SBA. Related concerns about locality were expressed by Bodo et al. (2008). However, the biggest issue about the SBA remains the so-called problem of convergence.

There are two configurations in which MRI-driven turbulence is studied numerically: one in which there is a net magnetic flux threading the layer and the other in which there is none. In the former case, if the uniform component of the magnetic field is vertical, say, there is a linearly unstable mode with a well-defined vertical wavenumber of maximum growth rate that sets the scale of the instability; in the latter case no such state exists, and the MRI must set in as a subcritical instability. Similar considerations apply if the initial field is azimuthal. What exactly determines the characteristic scale of the turbulence in the no-net-flux case is, at the moment, an open question. The conceptual appeal of the no-flux case is that it offers the possibility of a universal state of MRI turbulence in which the disk becomes self-magnetized through dynamo action, so that the angular momentum transport depends on the disk properties but not on the amount of flux threading the disk. Clearly much effort has been devoted to determining if such a universal state exists, and the value of the associated turbulent transport. Simply stated, the problem of convergence is that the angular momentum transport measured in numerical simulations based on the SBA with homogeneous background state and zero mean flux depends on numerical resolution, and decreases as the resolution increases.

This effect was first noted by Fromang & Papaloizou (2007) in a series of numerical experiments with the ZEUS code and subsequently confirmed by other authors using different codes (e.g., Pessah et al. 2007; Simon et al. 2009; Guan et al. 2009). There are a number of important points that should be made. Most of the evidence for the convergence problem is based on codes with no explicit viscosity or magnetic diffusivity in which the dissipation arises solely from truncation errors. An increase in resolution is interpreted as a decrease in effective dissipation, so that the convergence problem can equivalently be stated as a decrease in the effective transport with decreasing dissipation. Recently, however, Fromang (2010) indicated that the problem does not arise when explicit viscosity and magnetic diffusivity are included. Also, the problem is seen in simulations with homogeneous background state and periodic boundary conditions in the vertical. Simulations by Käpylä & Korpi (2011) suggest that the convergence problem is absent if different boundary conditions are applied in which the magnetic field is vertical on the upper and lower boundaries. Also, the problem seems to disappear when stratification is included (Davis et al. 2010), or the aspect ratio is such that the computational domain consists of a sufficiently tall box (J. M. Stone 2010, private communication). Finally, even when the effect is clearly observed there is no universal agreement about the rate at which the transport decreases with resolution. A number of questions naturally arise: what is the origin of the convergence problem? Is it a numerical artifact associated with the indiscriminate use of ideal solvers? Is it a real physical phenomenon related to subcritial dynamo instabilities? Is it that the SBA in its simplest form is too idealized to describe MRI turbulence in a disk? What is the asymptotic scaling of the angular momentum transport with resolution/dissipation, and what can we learn from it?

In this paper, we address some of these issues. We revisit numerically the problem of MRI-driven turbulence in shearing boxes with uniform background state and periodic boundary conditions in the vertical. We formulate the problem in a slightly unconventional way that is designed to emphasize the role of symmetries in the SBA and to make it easier to distinguish between numerical and physical effects.

In the ideal limit, the shearing-box equations are scale invariant, in the sense that doubling the number of grid points in a simulation or doubling the size of the computational domain are entirely equivalent. This gives rise to two equivalent limiting procedures: letting the grid spacing vanish for finite domain size, and letting the domain size become infinite for finite grid spacing. Although the two limits are equivalent, the latter makes it easier to interpret the results in terms of some classical results from dynamo theory. Using this formulation we consider the results of high-resolution numerical experiments and reach the following conclusions. The convergence problem is real in the sense that the transport by MRI-driven turbulence does decrease with decreasing dissipation and vanishes for vanishing dissipation. The problem is related to two aspects of the SBA with periodic boundary conditions: one is the inherent symmetry of the SBA and the other is the absence of an effective inverse cascade in the dynamo solutions. We find that there exists an asymptotic regime in which the effects of the symmetries are apparent but it is realized at very high resolution or, as we shall see presently, at very large system size, suggesting that most of the simulations described in the current literature are not in that regime. Finally, we argue that even though the SBA in its simplest form is probably unable to capture the physics of MRI-driven turbulence in disks, it nevertheless gives us very useful clues about which effects are likely to be important, and how to proceed to improve our models.

The paper is organized as follows. In the next section, we discuss the zero-flux case; we review the standard formulations and propose a slightly different approach that better separates physical and numerical quantities. We present numerical evidence that suggests the existence of an asymptotic state in which the transport becomes independent of system size. We then discuss cases with small but finite flux and show that they match smoothly with the zero-flux cases. Finally, we discuss the implications for modeling astrophysical disks.

We begin by discussing the case with no mean flux. Conceptually, this is the simplest case since for an ideal incompressible fluid, the SBA is completely scale invariant. As we shall see presently, this symmetry plays an important role in determining the asymptotic form of the solutions.

2.1. Formulation

The first step is to cast the SBA equations in dimensionless form. Often this is accomplished by selecting units that depend on the sound speed; although this choice is convenient from an astrophysical point of view it is not the most suitable to bring out the natural symmetries of the equations. Instead, we begin by considering an incompressible fluid, i.e., a fluid with infinite sound speed, with finite viscosity and magnetic diffusivity, and we shall return to the compressible case presently. A detailed presentation of the SBA can be found in Hawley et al. (1995); (see also Lesur & Longaretti 2007). The equations in dimensional form can be written as

$$\begin{eqnarray} \frac{\partial {\mathbf v}}{\partial t} + {\mathbf v}\cdot \nabla {\mathbf v}+ 2 \Omega \times {\mathbf v}&=& \frac{{\mathbf B}\cdot \nabla {\mathbf B}}{4 \pi \rho } - \frac{1}{\rho } \nabla \left(\frac{{\mathbf B}^2}{8 \pi } + P \right)\nonumber\\ && -\, \nabla (2 A \Omega x^2) + \nu \nabla ^2 {\mathbf v}, \end{eqnarray} \tag{ 1 }$$

$$\begin{equation} \frac{\partial {\mathbf B}}{\partial t} + {\mathbf v}\cdot \nabla {\mathbf B}- {\mathbf B}\cdot \nabla {\mathbf v}= \eta \nabla ^2 {\mathbf B}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \nabla \cdot {\mathbf v}=\nabla \cdot {\mathbf B}= 0, \end{equation} \tag{ 3 }$$

where B, v, ρ, and P denote, respectively, magnetic field, velocity, density, and pressure. The local angular velocity Ω = Ωe_z and shear rate

$$\begin{equation} A \equiv \frac{R}{2} \frac{\partial \Omega }{\partial R} \end{equation} \tag{ 4 }$$

are considered constants. For a Keplerian disk Ω∝R^−3/2 and A = −(3/4)Ω. The velocity can be decomposed as the sum of the base Keplerian flow and the fluctuations:

$$\begin{equation} {\mathbf v}= 2 A x {\mathbf e_y}+ {\mathbf u}; \end{equation} \tag{ 5 }$$

likewise, the pressure can be decomposed as the sum of the average and the fluctuations:

$$\begin{equation} P = P_0 + p. \end{equation} \tag{ 6 }$$

With these decompositions, the dimensional shearing-box equations become

$$\begin{eqnarray} \frac{D {\mathbf u}}{D t} + A x \frac{\partial {\mathbf u}}{\partial y} + A u_x {\mathbf e_y}+ 2 \Omega \times {\mathbf u}&=& \frac{{\mathbf B}\cdot \nabla {\mathbf B}}{4 \pi \rho } - \frac{1}{\rho } \nabla \left(\frac{{\mathbf B}^2}{8 \pi } + p \right)\nonumber\\ && +\, \nu \nabla ^2 {\mathbf u}, \end{eqnarray} \tag{ 7 }$$

$$\begin{equation} \frac{D {\mathbf B}}{D t} + A x \frac{\partial {\mathbf B}}{\partial y} - A B_x {\mathbf e_y}- {\mathbf B}\cdot \nabla {\mathbf u}= \eta \nabla ^2 {\mathbf B}, \end{equation} \tag{ 8 }$$

$$\begin{equation} \nabla \cdot {\mathbf u}= \nabla \cdot {\mathbf B}= 0, \end{equation} \tag{ 9 }$$

with D/Dt = ∂/∂t + u · ∇. The above system of equations contains only three dimensional parameters: the rotation frequency Ω, the viscosity ν, and the magnetic diffusivity η.

By adopting the following units of time, length, velocity, and magnetic field intensity

$$\begin{equation} \tau = \frac{1}{\Omega }, \quad {\mathcal L} =l_D \equiv \sqrt{\frac{\nu }{\Omega }}, \quad u^* = \sqrt{\nu \Omega }, \quad B^* = \sqrt{\rho \nu \Omega } \end{equation} \tag{ 10 }$$

the equations can be written in dimensionless form as

$$\begin{eqnarray} \frac{D \hat{{\mathbf u}}}{D \hat{t}} {-} \frac{3}{4} \hat{x} \frac{\partial \hat{{\mathbf u}}}{\partial \hat{y}} -\frac{3}{4} \hat{u}_x {\mathbf e_y}+ 2 {\mathbf e_z}\times \hat{{\mathbf u}}&=& \frac{1}{4 \pi } \hat{{\mathbf B}}\cdot \hat{\nabla }\hat{{\mathbf B}} \nonumber\\ && - \hat{\nabla }\left(\frac{\hat{{\mathbf B}}^2}{8 \pi } + \hat{p} \right) \,{+}\, \hat{\nabla }^2 \hat{{\mathbf u}},\nonumber\\ \end{eqnarray} \tag{ 11 }$$

$$\begin{equation} \frac{D \hat{{\mathbf B}}}{D \hat{t}} -\frac{3}{4} \hat{x} \frac{\partial \hat{{\mathbf B}}}{\partial \hat{y}} + \frac{3}{4} \hat{B}_x {\mathbf e_y}- \hat{{\mathbf B}}\cdot \nabla \hat{{\mathbf u}}+ \frac{1}{P_m} \hat{\nabla }^2 \hat{{\mathbf B}}= 0, \end{equation} \tag{ 12 }$$

where the hatted quantities are dimensionless and P_m = ν/η is the magnetic Prandtl number. It is important to note that the induction equation (12) depends only on P_m, while the momentum equation (11) has no adjustable parameters at all. This is because in the SBA with no mean flux, the diffusive scale l_d is the only intrinsic length that can be defined using the dimensional parameters of the problem. Clearly, in order to solve Equations (11) and (12) one also needs to define a computational domain and suitable boundary conditions, thereby introducing an external length scale L that characterizes the size of the computational box; the latter, however, has no direct physical meaning. In particular, in the case of periodic boundary conditions (shearing periodic in the radial direction), formally one is considering an infinite domain, with L representing a long-wavelength cutoff. The domain size can itself be non-dimensionalized to define the additional parameter $R \equiv L/{\mathcal L}$. For fixed magnetic Prandtl number, R—the domain size in units of the diffusive scale—is the only adjustable parameter.

We note that the non-dimensional form of the equations above is different from what is typically found in the literature. In general, the box size is used as the unit of length and the square of the parameter R is then identified with the Reynolds number. However, in a turbulent flow, the Reynolds number should be more properly defined in terms of the integral scale of the turbulence. The question then becomes whether asymptotically, for R → ∞, in MRI turbulence the integral scale is linked to the box size, or to the diffusive scale, or to a combination of both. We note that if the solutions become asymptotically independent of the box size, the system of Equations (11) and (12) does not depend on any non-dimensional parameters, apart from the Prandtl number, i.e., for fixed Prandtl number there is a universal solution.

In the zero net flux case, dynamo processes must effectively maintain the magnetic field against dissipation; however the magnetic field is itself responsible for generating, through the MRI, the turbulence that sustains it. It is well known that, in a (small-scale) turbulent dynamo, the magnetic energy is concentrated near the resistive scale. In this case, one expects that the MRI-driven velocity may also be characterized by the diffusive scale. In order for the velocity to have a significant component comparable to the box size an efficient mechanism is required that can transfer energy backward, from the resistive scale to the largest available scale, i.e., an inverse cascade. Whether such an inverse cascade exists is at the heart of the problem.

As an aside, we note that the condition for the validity of the SBA is that the characteristic scales of the solution be much smaller than the box size, lest the use of periodic boundary conditions is not justified. Therefore, if there were cases in which the shearing-box results would show a dependence of the solutions on the box size R as R → ∞, one should actually reformulate these cases by introducing some global scale characteristic of the whole disk.

One of the main objectives of MRI studies is to determine the efficiency of the transport of angular momentum by Maxwell and Reynolds stresses. It is customary to define the quantity

$$\begin{equation} \Sigma \equiv \overline{\left\langle u_x u_y - \frac{B_x B_y}{4 \pi \rho }\right\rangle } \end{equation} \tag{ 13 }$$

that represents the box and time-averaged value of the total stresses (hereinafter overbars denote time averages, while angle brackets denote a box average). This quantity has dimensions of square velocity; thus with our choice of non-dimensionalization it is measured in units of Ων. For fixed magnetic Prandtl number we can write

$$\begin{equation} \Sigma \sim f(R) \Omega ^2 l_D^2 = f(R) \Omega \nu. \end{equation} \tag{ 14 }$$

Here, f(R) is a dimensionless function that accounts for the possible dependence of the solutions on the externally imposed box size. When R = O(1) the box size is comparable to the diffusive scale, and, on general grounds, one expect a strong dependence of f(R) on R. It is not clear, a priori, what dependence one should expect as R → ∞. It is important, however, to be clear about the relationship between the asymptotic form of f(R) and the transport efficiency. In turbulence theory, it is customary to refer to "turbulent transport" to describe a transport process that becomes independent of diffusion for small diffusion, and to "collisional transport" to describe one such process that remains proportional to the diffusivities when these become small. Because, here, Σ is measured in diffusive units (Ων), if f(R) ∼ constant as R → ∞ the angular momentum transport is entirely collisional. Conversely, in order for the transport to be "turbulent" in the sense defined above, f(R) ∼ R² as R → ∞. In fact, any asymptotic dependence of f(R) weaker than R² will ultimately lead to a transport that vanishes for vanishing viscosity.

Typically, in the literature, the transport efficiency is measured in terms of the parameter α. In the incompressible case this is defined as the total stresses measured in units of L²Ω² (Lesur & Longaretti 2007); in other words

$$\begin{equation} \alpha \equiv \frac{\Sigma }{\Omega ^2 L^2} = f(R) \frac{1}{R^2}, \end{equation} \tag{ 15 }$$

since R² = L²Ω/ν. Again, we recover the result that in order for α to have asymptotically finite (constant) value the function f(R) must diverge quadratically as R becomes large.

When compressibility is taken into account, an additional physical parameter is introduced, namely the sound speed c_s, that can be used to define a new length scale H = Ωc_s. When vertical gravity is accounted for, H represents the scale height of the disk. In studies such as the present one, in which vertical gravity is neglected, this length has no direct meaning; nevertheless it introduces a new non-dimensional parameter $R_{H} = H/{\mathcal L}$, that must be taken into account in the scaling arguments discussed above. As before we can write

$$\begin{equation} \Sigma = F(R,R_H) \Omega \nu, \end{equation} \tag{ 16 }$$

where F(R, R_H) is a dimensionless function. For R_H = O(1), again we expect a strong dependence of F on R_H. However for R_H large, F must become asymptotically independent of R_H, since this corresponds to the incompressible limit in which F must approach f.

Physically, this corresponds to situations in which H is much larger than the characteristic scales of the solutions. Again, if there is no dependence of these scales on the externally imposed box size L, the only available length is ${\mathcal L}$ and the characteristic scales of the solutions will be proportional to ${\mathcal L}$. For constant sound speed and therefore constant ratio H/L, the limit R → ∞ corresponds to R_H → ∞, and therefore the stress behavior should be the same as in the incompressible case. In this case α is defined as in Shakura & Sunyaev (1973),

$$\begin{equation} \alpha \equiv \frac{\Sigma }{\Omega H c_s} = \frac{\Sigma }{ c_s^2}, \end{equation} \tag{ 17 }$$

and we therefore expect again the scaling

$$\begin{equation} \alpha \sim \frac{1}{R^2}. \end{equation} \tag{ 18 }$$

2.2. Numerical Considerations

In most of the MRI studies, as well as in this paper, the analysis of MRI turbulence is performed using ideal codes, where dissipative effects only arise due to numerical truncation. In general, in this case it is not known how to write the dissipation term explicitly; however, it can be said that the amount of dissipation depends on the ratio between the local scale of variation of the physical quantities and the cell size. Structures comparable to the cell size are subject to strong dissipation that decreases when the scale of the structures increases. Different numerical schemes, however, will behave in different ways; in particular the reduction of dissipation when the local scale is decreased will be sharper for high-order schemes and gentler for lower-order ones. In (10), we defined the dissipation scale in terms of the viscosity and the shear as $l_D = {\mathcal L} = \sqrt{\nu /\Omega }$. In the numerical approach one is not able to get an explicit expression for the dissipation scale; nevertheless it is known that l_D is related to the cell size δ in a way that, however, depends on the scheme. In addition, it is known that l_D (in units of δ) will be larger for lower-order schemes and smaller for the high-order ones. Thus, in the non-dimensionalization of the equations, in principle one should use l_D, but in practice one can only use δ. With the two dimensional parameters Ω and δ, it is possible to define the units of time τ, length ${\mathcal L}$, velocity u*, and magnetic field B* as

$$\begin{equation} \tau \,{=}\, \frac{1}{\Omega }, \quad {\mathcal L} \,{=}\, \delta, \quad u^* \,{=}\, \Omega \delta, \quad B^* \,{=}\, \sqrt{\rho } \Omega \delta. \end{equation} \tag{ 19 }$$

By the same considerations discussed above, the scaling of stresses can be written as

$$\begin{equation} \Sigma \sim f_\delta (N) \Omega ^2 \delta ^2, \end{equation} \tag{ 20 }$$

where N = L/δ represents the number of grid points and the function f_δ(N) depends on the numerical scheme and accounts for a possible dependence of the solutions on the externally imposed box size. The scaling for α is correspondingly given by

$$\begin{equation} \alpha \sim f_\delta (N) \frac{1}{N^2}. \end{equation} \tag{ 21 }$$

If there is no asymptotic dependence on the box size, i.e., f_δ(N) ∼ constant as N → ∞, in this limit we obtain αN² = constant where the constant will be different from scheme to scheme in view of the fact that different numerical schemes, as discussed above, have different relationships between the dissipation length l_D and the cell size δ. The proper unit for Σ is given by Ω²l²_D, and consequently, we expect

$$\begin{equation} \alpha N^2 \sim \left(\frac{l_D}{\delta }\right)^2. \end{equation} \tag{ 22 }$$

Since l_D is expected to be smaller for high-order codes the same should be true for the value of αN². It should be noted that this discussion ignores possible differences in the effective Prandtl number for the different schemes, which can account for additional residual differences between codes.

For the compressible case, we can repeat the same reasoning as above: we have an additional parameter N_H = H/δ that represents the number of grid points over H; however the dependence on this parameter should disappear, for constant sound speed, in the limit N → ∞.

2.3. Numerical Results

We have performed a series of compressible isothermal shearing-box simulations for the zero net flux case, with different resolutions and with three different numerical schemes available in the PLUTO code (Mignone et al. 2007). The main difference between the three schemes is the order of accuracy of the reconstruction; more precisely we employed a scheme with total variation diminishing linear reconstruction (second-order accurate for smooth flows), a parabolic reconstruction (PPM), and a monotonicity-preserving fifth-order reconstruction (MP5) (for a description of this scheme, see Suresh & Huynh 1997; Mignone & Tzeferacos 2010; Mignone et al. 2010). Our goal is to study the asymptotic behavior of the solutions as the separation of scales between the box size and the dissipation length becomes large. This can be achieved by increasing the resolution for a fixed scheme or, for a fixed resolution, by using a scheme with higher order of accuracy.

In the shearing-box approach, the Cartesian coordinates x, y, and z refer, respectively, to the radial, azimuthal, and vertical directions. Our computational box has aspect ratio L_x:L_y:L_z = 1:π:1. As in Fromang & Papaloizou (2007) we use cells elongated in the azimuthal direction, with aspect ratio 1:2:1. For each of the schemes we used four different resolutions: 32, 64, 128, and 256 points in the vertical direction, the corresponding grid sizes are therefore 32 × 50 × 32, 64 × 100 × 64, 128 × 200 × 128, and 256 × 400 × 256. The sound speed is chosen such that H = L, the initial magnetic field is (0, 0, B₀sin (2πx)), with B₀ corresponding to β = 1500, and random noise in the y-component of the velocity is introduced initially to start the growth of the instability.

As is well known, the time histories of quantities like the box-averaged magnetic energy or the box-averaged Maxwell or Reynolds stresses have an initial transient followed by the establishment of a stationary state in which the quantities fluctuate around a well-defined mean value. The amplitude of the fluctuations decreases with increasing resolution. To estimate the mean value we average over the simulation time excluding the initial transient phase. The simulation time varies from case to case and it is always larger than 100 revolutions. We estimate the error of the mean value by subdividing the averaging interval into 10 subintervals and then computing the standard deviation of the subinterval averages. The resulting estimate of the error is always less than 10%.

2.3.1. Transport

We start our analysis by considering the behavior of α, as defined in Equation (17), as a function of resolution. Figure 1 shows αN² as a function of N for all the cases considered. The quantity displayed represents the function f_δ(N) defined in Equation (21), which accounts for a possible dependence of the solutions on the externally imposed box size. Consonant with our objectives, we are interested in the asymptotic behavior of f_δ(N) as N → ∞.

Zoom In Zoom Out Reset image size

We recall that in order for α to be resolution independent, f_δ(N) should diverge as N². The results of Fromang & Papaloizou (2007) suggested that f_δ(N)∝N. A constant value of f_δ(N), on the other hand, means that the solution is asymptotically independent of the box size, that the scaling of α is determined solely by the equations, and that α decreases with resolution.

In Figure 1, the squares refer to the linear scheme, the stars to the PPM, and the triangles to the MP5. It is obvious that, as the order of accuracy of the scheme is increased, the curves tend to be flatter for increasing N. However, while the linear and PPM schemes show a continually decreasing slope and do not seem to have yet reached an asymptotic behavior, the MP5 scheme, by contrast, appears to be leveling out. The results of MP5 suggest that f_δ(N) approaches a constant value as N → ∞ and therefore that α has the asymptotic natural scaling ∼1/N² arising from the equations alone. Our expectation is that, given sufficient resolution, a similar behavior would be observed with the other schemes.

It is important to note that even if all the curves corresponding to different schemes eventually level off, they will not in general have identical asymptotic values. The reason for this is that we do not know the explicit form of the numerical dissipation for each scheme; accordingly we have used the cell size instead of the dissipation length to non-dimensionalize the equations. Indeed, we note that in agreement with Equation (22) the higher-order scheme MP5 has a lower value of the non-dimensional stresses than the other schemes.

2.3.2. Structures

The behavior of α as a function of resolution suggests that the only scale that determines the properties of the solutions is the dissipative scale, which, in numerical studies, is related to the cell size. It is therefore important to examine the behavior of the size of the typical observed structures as a function of resolution. In order to characterize the scales of magnetic structures, we measure the average scales of variation of the field in the directions parallel and perpendicular to itself. We define the quantities

$$\begin{equation} l_\parallel = \overline{\left(\frac{\langle {|} {\mathbf B}|^4\rangle }{\langle |{\mathbf B}\cdot \nabla {\mathbf B}|^2\rangle } \right)^{1/2}}, \qquad l_\perp = \overline{ \left(\frac{\langle {|} {\mathbf B}|^4\rangle }{\langle | {\mathbf J}\times {\mathbf B}|^2\rangle } \right)^{1/2} }, \end{equation} \tag{ 23 }$$

which represent measures of the characteristic scales of the magnetic field in the parallel direction and along the direction of maximum gradient, respectively (see Schekochihin et al. 2004). In Figures 2 and 3, we plot, respectively, Nl_⊥ and Nl_∥ as functions of N. As before, the three curves refer to the three different numerical schemes and the symbols have the same meaning as in the previous figure. The lengths are normalized with respect to the cell size and the quantities plotted therefore represent the average number of grid points on which the observed structures extend. As expected, the magnetic structures are highly anisotropic with l_∥/l_⊥ ∼ 10. Considering the dependence on resolution, we see again that the linear and PPM schemes show an increase of the number of grid points over which magnetic structures extend, while the MP5 scheme appears to tend to a constant number of grid points independent of resolution, in both the longitudinal and transverse directions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This analysis suggests that the magnetic field is characterized by elongated structures extending transversally over a few grid points. In order to get further insight into the three-dimensional structure one can define a third scale in the direction orthogonal to both B and B × J and that points approximately in the vertical direction:

$$\begin{equation} l_{{\mathbf B}\cdot {\mathbf J}} = \overline{ \left(\frac{\langle | {\mathbf B}|^4\rangle }{\langle |{\mathbf B}\,{\bm \cdot}\, {\mathbf J}|^2\rangle } \right)^{1/2} }. \end{equation} \tag{ 24 }$$

The ratios of the scales along the three directions are typically l_⊥:l_∥:l_{B · J} ∼ 1:10:5 indicating that the magnetic structures can be thought of as thin sheets. This result is also discussed by Guan et al. (2009), who base their analysis on the correlation function. They find that the correlation lengths scale as N^−2/3 and that therefore, measured in units of the cell size, show an increase with resolution, suggesting that one could expect some kind of transition. We have also repeated this convergence analysis on the correlation lengths and again we find that while the linear and PPM schemes show a behavior similar to that discussed by Guan et al. (2009), the MP5 scheme at a resolution of 256 grid points shows a convergent behavior, i.e., the correlation lengths scale as the cell size. A better impression of the magnetic field structure can be obtained by looking at Figure 4, where we show a two-dimensional cut of the Maxwell stress distribution in the x–y-plane at z = 0. The three-dimensional sheets appear in the figure as filaments of high Maxwell stress and high magnetic field intensity, separated by regions of low magnetic field.

Zoom In Zoom Out Reset image size

We can now discuss this filamentary structure from a more physical point of view. From the y-component of the induction equation it can be seen that the field structure is determined by the competition between three terms: the stretching by the background shear, the stretching due to velocity fluctuations (B · ∇v term) and dissipation, and the transverse size is determined by the balance between these terms. The background shear tends to stretch the field along the azimuthal direction, simultaneously decreasing its perpendicular scale of variation until dissipation sets in. These magnetic filaments, however, give rise to a local transport of momentum and, therefore, tend to reduce the local shear. One therefore expects that the effect of the nonlinear term B · ∇v is to reduce the stretching by the background shear. The observed increase of the normalized transverse scale with increasing resolution may be related to the increase of this effect, which however ultimately has to be limited since it cannot completely suppress the background shear. Thus, asymptotically, the transverse scale has to remain constant and be determined by the balance between shear stretching and dissipation. The MP5 scheme, which is closer to the asymptotic behavior, accordingly shows a tendency of the transverse scale to become constant. This leads to the conclusion that a possible dependence on the scale of the box may appear only in the parallel scale. However, the ratio l_∥/l_⊥ seems to be independent of resolution, i.e., the magnetic structures do not appear to become more elongated as the resolution increases and all the characteristic lengths appear to become asymptotically constant when measured in units of the dissipative length.

2.3.3. Probability Distribution Functions

Further details on the solutions can be examined by considering the behavior of the probability distribution function (pdf) of magnetic field intensity and of the second-order structure functions. Our purpose is twofold: on the one hand we want to show that, indeed, if we measure all quantities in their proper units, there is a convergence to a well-defined solution when we increase the number of grid points and, on the other hand, we want to better characterize this solution. We start by examining the pdf of magnetic field intensity. Figure 5 compares the two MP5 cases with the highest resolution; the pdfs are normalized to unity and the field values are measured in units of B*.

Zoom In Zoom Out Reset image size

The pdfs have an exponential character; the mean value of |B|/B* is about 5.5 with a difference of less than 1% between the two cases at different resolution. The contribution to the magnetic energy and to the Maxwell stresses comes mainly from the high field values; values to the right of the vertical line plotted in the figure contribute to about 70% of the total magnetic energy and to about 80% of the total Maxwell stresses, while occupying only about 20% of the volume. The higher contribution to the Maxwell stresses with respect to the magnetic energy is explained by the fact that these higher field values correspond to the magnetic structures discussed above. For these, there is a strong correlation between the radial and azimuthal components of the field, while for lower field values the correlation is much lower. The fact that the volume fraction supporting the transport is independent of resolution explains the decrease in the amplitude of the fluctuations when resolution is increased (Fromang & Papaloizou 2007). In fact, as resolution is increased and the volume fraction remains the same, the number of points that contributes to the average gets larger and therefore the fluctuations of the average become smaller.

In Figure 6, we plot the time-averaged second-order longitudinal and transverse structure functions for the y-component of magnetic field, for the two MP5 cases at the highest resolution. They are explicitly defined as

$$\begin{equation} \begin{array}{l} \dsty \overline{S_{2l}} = \overline{\langle (B_y (x, y+h, z) - B_y(x,y,z))^2 \rangle }, \\ \ms \ms \dsty \overline{S_{2t}} = \overline{\langle (B_y (x+h, y, z) - B_y(x,y,z))^2 \rangle }.\end{array} \end{equation} \tag{ 25 }$$

Following our discussion, we compute S_2l and S_2t in units of B*² and h in units of δ. Again, it can be seen that, using these units, the difference between the two cases at different resolution is very small, indicating a convergent behavior. The structure functions in both directions become flat at scales larger than few grid points. This plateau indicates that the magnetic field values become uncorrelated at these scales. A similar conclusion could be reached from the behavior of the averaged spectra, like those shown in Fromang & Papaloizou (2007), that become flat at small wavenumbers.

Zoom In Zoom Out Reset image size

From these two results we conclude that the solution is characterized by the superposition of uncorrelated magnetic sheets. Those with high field intensity also have a high degree of correlation between magnetic field components and, consequently, give the highest contribution to both the total magnetic energy and the Maxwell stresses. An increase of the domain size does not change the properties of these sheets, but simply increases the available volume and therefore the number of sheets, leaving unchanged all the average properties. In a similar way there should not be any modification of the solution by changing the aspect ratio of the computational box, and we have indeed tested this by examining two further cases with respective aspect ratios of 1:2π:1 and 2:π:1, with MP5 and 128 points in the vertical direction, finding no changes in all the averaged properties.

We now turn to the case in which a net magnetic flux threads the computational domain. In this case, the magnetic field has a uniform component of strength B₀, say, that can be used to define a new dimensional length scale λ_A given by

$$\begin{equation} \lambda _A = \frac{1}{\Omega } \frac{B_0}{\sqrt{4 \pi \rho }}. \end{equation} \tag{ 26 }$$

Accordingly, the averaged stresses will show a dependence also on this scale and keeping the same units as in the previous section, in general we can write

$$\begin{equation} \Sigma = g \left(\frac{\lambda _A}{l_D}, \frac{L}{l_D} \right) \Omega ^2 l_D^2. \end{equation} \tag{ 27 }$$

Here, we are mostly interested in the form of the function g when L/l_D = R → ∞. On general grounds, we expect it will take different forms for different relationships between the three scales L, λ_A, and l_D. We start our discussion by considering what is known from numerical simulations in the available literature. The case that has received the greatest attention is that of a mean vertical field; the results have been summarized in Pessah et al. (2007). For this case the system is linearly unstable to the MRI; the (vertical) wavelength of maximum growth rate is given by $\lambda _M = 2 \pi \sqrt{16/15} \lambda _A$; the instability exists only for wavelengths $\lambda > 1/\sqrt{3} \lambda _A$. If λ_A is larger than $\sqrt{3} L$, no unstable wavelengths fit in the box, the system becomes stable, and the stresses drop to zero. Pessah et al. (2007) show that as λ_A increases, the stresses scale linearly with λ_A, reach a maximum, and eventually, drop to zero when the above condition is met. Furthermore, it has been found that, at fixed λ_A/L and at low resolution, the stresses increase with resolution, i.e., with R (see e.g., Silvers 2008; Bodo et al. 2008); thus, for high enough resolution, since they cannot grow indefinitely, they should converge to a value independent of R. Also, we have to note that Bodo et al. (2008) have found a dependence of the solutions on the aspect ratio L_x/L_z of the computational box, with a convergent behavior for high enough values of L_x/L_z. Putting these considerations together, we conclude that for l_D ≪ λ_A ≲ L,

$$\begin{equation} \Sigma \sim g \left(\frac{\lambda _A}{L} \right) \Omega ^2 L^2, \end{equation} \tag{ 28 }$$

where the function g is initially proportional to λ_A/L, then reaches a maximum (at a fixed value of λ_A/L), and finally drops to zero for $\lambda _A > \sqrt{3} L$. Thus, the value of Σ at the maximum scales as Ω²L².

In the opposite regime, when l_D ∼ λ_A ≪ L, the averaged stresses tend toward the value obtained in the zero-flux case (Pessah et al. 2007). The results of the previous section suggest that, for large R, Σ is independent of L and scales as Ω²l²_D.

In order to connect smoothly the two portions with, respectively, l_D ∼ λ_A ≪ L and l_D ≪ λ_A ≲ L, there must be an intermediate range of values of λ_A in which the stresses scale as Ω²λ²_A. In fact, when λ_A is sufficiently far both from the dissipation scale and from the size of the box, the transport should depend only on λ_A and, therefore, by dimensional arguments it should have a quadratic scaling. Increasing R, the range over which λ_A is sufficiently far both from the dissipation scale and from the size of the box increases and, therefore, the region with quadratic scaling should increase. Conversely, decreasing R should cause the interval to shrink and eventually to disappear for low values of R. In Figure 7, we show the results of calculations with a net vertical flux performed at a resolution of 128 points in the vertical direction, with the MP5 scheme. The horizontal dashed line corresponds to the averaged stress value for the zero-flux case, the dashed line corresponds to a quadratic scaling, the solid line corresponds to a linear scaling, and the vertical solid line marks the stability boundary; to the right of it, the stresses vanish.

Zoom In Zoom Out Reset image size

It is apparent from the figure that indeed for small enough λ_A, a range of values exists with quadratic scaling, while for larger λ_A the curve reverts to the linear scaling discussed above. We have deliberately restricted our investigation to values of λ_A for which we have quadratic scaling and the transition from quadratic to linear.

Larger values of λ_A, for which the stresses increase linearly up to a maximum and drop to zero to the right of the stability boundary, have already been discussed extensively in the literature (see e.g., Pessah et al. 2007). We note here that the portion with quadratic scaling has never before been reported in the literature (however, Longaretti & Lesur 2010 find evidence of a scaling steeper than linear); this is most likely due to limited resolution since, as discussed above, the quadratic portion disappears for low values of R. We believe that the reason why we were able to detect it is because of the high order of accuracy of the MP5 scheme, which for the same resolution produces a much higher effective value of R.

In this paper, we have addressed the problem of convergence in numerical studies of MRI-driven turbulence in the homogeneous SBA. Any simulation that is physically useful should have some meaningful property that becomes independent of the number of grid points as the latter becomes large. This does not seem to be the case here. As noted by several authors, in the weak flux regime, i.e., when the total magnetic flux threading the computational box is small or zero, the dimensional angular momentum transport decreases as the resolution increases, eventually becoming arbitrarily small (Fromang & Papaloizou 2007; Pessah et al. 2007; Simon et al. 2009; Guan et al. 2009). There are at least two possible frameworks in which to address this issue: one in terms of Reynolds number, the other in terms of computational domain size. In both cases one is interested in the asymptotic behavior of the solutions as either the Reynolds number or the domain size become large. Clearly, both frameworks are related and to some extent it is a matter of preference which one is picked. In this paper, we have chosen to discuss the convergence problem in terms of domain size. The reason is that, because of the symmetries of the homogeneous SBA, there is no a priori characteristic outer scale for the velocity, which leads to an ambiguity in the definition of the Reynolds number. On the other hand there is a natural inner scale—i.e., the dissipation length—that can be used to non-dimensionalize the equations. In this case, the Reynolds number is fixed, and of order unity, and the only remaining free parameter is the domain size. With this choice, increasing the resolution and increasing the domain size are entirely equivalent operations. It is then possible to define the characteristic scale of the solutions, once they are computed, in terms of, say, the (inverse) scale at which the velocity or magnetic spectrum peaks; the latter quantity is conceptually useful since it is related to the effective angular momentum transport. The problem of convergence then relates to the behavior of the characteristic scale of the solutions as the computational domain size increases.

As we discussed in Section 2, there are two limiting possibilities: one in which the solution scale diverges with the domain size—measured in units of the grid size or the dissipation length—the other in which it becomes asymptotically independent of it. The former can lead to a net turbulent transport that is independent of diffusivity, the latter to a quasi-collisional transport that scales as some fixed multiple of the diffusivity. The evidence from the numerical simulations, in particular Figure 1, is that it is the latter and not the former that gives the correct asymptotic behavior. In other words, the characteristic scale of the solutions becomes independent of the size of the computational domain. If we accept this result as correct, there are a number of issues that naturally arise.

The first concerns the physical interpretation of the result. In the weak-field regime, MRI-driven turbulence should be considered within the framework of a subcritical dynamo instability since the magnetic field necessary for the development of the MRI must be generated self-consistently by the MRI turbulence itself. In dynamo theory, it is customary to distinguish between two types of dynamo action, small-scale and large-scale. Small-scale dynamos generate magnetic field with characteristic scale comparable to or smaller than the characteristic scale of the velocity field. The only requirements for the operation of such dynamos are that the magnetic Reynolds number be sufficiently high, and that the underlying velocity not be too symmetric. Clearly, both conditions can be satisfied here provided the system size is large. Large-scale dynamos, on the other hand, generate magnetic field with characteristic scale larger than that of the velocity. The latter types of dynamo are associated with flows lacking reflectional symmetry, i.e., helical flows and an inverse cascade of magnetic helicity. Recently, Tobias et al. (2011) have argued that in the nonlinear regime it is often difficult to distinguish unambiguously between large- and small-scale dynamo action and have introduced instead the idea of a system-scale dynamo, as one that produces magnetic structures comparable in size to the system scale irrespective of the scale of the velocity. The fact that, here, the scale at which the magnetic spectrum peaks is comparable with the dissipation scale and it is asymptotically independent of the system size indicates that the dynamo operating here is of the small-scale type (Vainshtein & Kichatinov 1986; Boldyrev & Cattaneo 2004). Dynamos of this type have no manifest inverse cascade of (magnetic) helicity; they can lead to the production of magnetic energy but not of substantial magnetic flux. In the specific context of MRI-driven turbulence a small-scale dynamo will not give rise to a "turbulent" transport of angular momentum. When regarded in the general scheme of dynamo theory, this result is actually not that surprising since the underlying flow is not strongly helical, and therefore there is no a priori reason to expect an efficient inverse cascade.

This brings us to the second issue, which relates to the astrophysical significance of these types of MRI solutions. Clearly, it is to be expected that this type of dynamo excitation will be prevalent in most (electrically conducting) disks. Its role in the transport of angular momentum, and therefore in regulating the accretion rate, however will be negligible. With things as they stand, the only dynamo solutions that might have an impact on the angular momentum transport are those with an efficient inverse cascade. Interestingly, this implies that by their very nature such solutions should not be treated in the context of a local approximation. We find ourselves in a paradoxical situation in which if a solution is astrophysically relevant it cannot properly be described by a local model, and if it can be described by a local model then it is astrophysically irrelevant. The third issue concerns the role of numerics in influencing what we understand, or think we understand, in computational astrophysics. It is interesting that the asymptotic behavior that eventually emerges according to the present study is in many ways the most natural that is consonant with a subcritical small-scale dynamo instability developing in a system where there are not many reasons to expect otherwise. Likewise, the scaling behavior for the cases with weak imposed fluxes is asymptotically the one that makes most sense given the symmetries of the equations. Yet, their numerical realizations were somewhat long in the making, and had not been reported in the literature to date. The reason is that these asymptotic regimes require a substantial dynamic range to manifest themselves. In this particular case, this was achieved by a combination of high-resolution and high-accuracy numerics. If we had tried to obtain the correct limiting behaviors with the less accurate codes, it would have been prohibitive.⁴ Finally, it behooves us to compare the conclusions of the present study with those of related studies in which different scaling behaviors have been observed. Of particular relevance are: the work of Käpylä & Korpi (2011) who also consider unstratified shearing boxes but with non-periodic boundary conditions; the work of Fromang (2010) who has the same set-up as ours but with explicit diffusivities; and the work of Davis et al. (2010) who include stratification. In all these cases the authors report that the problem of convergence as stated above is absent. In the present context, these conclusions can be attributed to one or more of the following possibilities: the simulations are not in the asymptotic regime because of limited resolution and/or accuracy; the extra physics has engendered an inverse cascade, i.e., the type of dynamo action has changed from small-scale to large-scale; the extra physics has changed the nature of the subcritical dynamo instability and has caused it to become of the system-size type. The last two possibilities are related but distinct. They both refer to circumstances that lead to the production of large-scale magnetic fields, but in a large-scale dynamo the velocity correlation length remains small and it is only the magnetic structures that become large; in a system-scale dynamo both the velocity and magnetic structures grow in size until they occupy the largest available scale—the system size—hence the name (Tobias et al. 2011).

We begin by discussing the work of Käpylä & Korpi (2011). They consider unstratified shearing-box simulations with boundary conditions at the upper and lower boundaries corresponding to impenetrable, stress-free velocities, and vertical magnetic fields. The distinctive feature of these conditions is that they allow a flux of magnetic helicity in and out of the box. This has been advocated by several authors as a feature that might be beneficial to large-scale dynamo action (see, e.g., Vishniac & Cho 2001). The calculations are carried out with explicit diffusivities and the authors report that the effective angular momentum transport is substantial and insensitive to variations in the values of both the magnetic Reynolds and Prandtl numbers. The morphology of the solutions is dramatically different from the periodic cases, with two large current boundary layers forming in the vicinity of the upper boundaries and the generation of a coherent nearly uniform azimuthal field. All these findings are consistent with the activation of a large-scale-type dynamo. It remains to be seen if these solutions persist at higher resolutions, as well as how they map to a global geometry. These issues notwithstanding, the Käpylä & Korpi (2011) solutions are quite remarkable and could potentially have considerable astrophysical significance. Next, we discuss stratification (Davis et al. 2010). Its presence drives two new effects: it introduces a new spatial scale in the problem, namely the vertical scale height that breaks the symmetry of the homogeneous, periodic models, and it also introduces buoyancy forces associated with thermal, pressure, and magnetic pressure fluctuations. It is not unlikely that these new effects could modify the type of dynamo action and drive flows and field coherent on scales comparable with the scale height. If this is indeed the case, it would be interesting to see how the solutions behave as the simulation size is increased for fixed scale height. Such a study would be numerically very demanding but astrophysically very useful.

Finally, we consider the recent results of Fromang (2010) who reports that when explicit diffusivities are included, the angular momentum transport becomes independent of resolution, i.e., the convergence problem goes away. If this result has a physical basis, it implies that the presence or absence of an inverse cascade and the type of associated dynamo action (large-scale as opposed to small-scale) depends on the detailed form of the diffusivities. The dynamo is small-scale for all three types of numerical diffusivities considered here and in other similar studies, and large-scale for "physical" diffusivities. Although this possibility is not inconceivable, it is somewhat peculiar. It implies that the type of dynamo action is determined not by the symmetries of the problem, or large-scale features, but by the micro-physics. If correct it is a remarkable result. Another possibility is that the simulations are not in the asymptotic regime. To explore this possibility we note that the study by Fromang (2010) was carried out by incorporating explicit diffusivities into the ZEUS code. In order for the explicit diffusivities to be correctly represented they must give rise to boundary layers whose size is not smaller than those arising from truncation errors. Thus, the dynamic range of a code with numerical diffusivities gives an upper bound on the dynamic range achievable with the same code, the same resolution, and explicit diffusivities. Now, ZEUS has an order of accuracy comparable to the code with quadratic reconstruction: thus we can ask the following question. What is the resolution that we need with the quadratic code to see the same asymptotic behavior we observe with the MP5 code? This can be estimated by comparing the curves in Figure 1 and taking the ratio in resolutions for which the curves corresponding to the PPM and MP5 codes have similar slopes. This exercise gives a factor close to four. A more careful analysis comparing boundary layers in a channel flow problem yields a similar estimate. Thus, our conclusion is that our quadratic code requires a resolution exceeding 1000 grid points to reproduce the asymptotic behavior observed by the MP5 code with 256 grid points. Since the highest resolution in the work by Fromang (2010) was 512 grid points it is possible that his results may not yet be in the asymptotic regime.

This work was supported in part by the National Science Foundation-sponsored Center for Magnetic Self Organization at the University of Chicago. Calculations were performed at CINECA (Bologna, Italy) thanks to support by INAF and at CASPUR (Rome, Italy).