MHD simulations of the magnetorotational instability in a shearing box with zero net flux: the case Pm=4

MHD simulations of the magnetorotational instability in a shearing box with zero net flux: the case

Sébastien Fromang CEA, Irfu, SAp, Centre de Saclay, F-91191 Gif-sur-Yvette, France UMR AIM, CEA-CNRS-Univ. Paris VII, Centre de Saclay, F-91191 Gif-sur-Yvette, France.
Accepted; Received; in original form;
Key Words.:
Accretion, accretion disks - MHD - Methods: numerical
offprints: S.Fromang


Aims:This letter investigates the transport properties of MHD turbulence induced by the magnetorotational instability at large Reynolds numbers when the magnetic Prandtl number is larger than unity.

Methods:Three MHD simulations of the magnetorotational instability (MRI) in the unstratified shearing box with zero net flux are presented. These simulations are performed with the code Zeus and consider the evolution of the rate of angular momentum transport as is gradually increased from to while simultaneously keeping . To ensure that the small scale features of the flow are well resolved, the resolution varies from cells per disk scaleheight to cells per scaleheight. The latter constitutes the highest resolution of an MRI turbulence simulation to date.

Results:The rate of angular momentum transport, measured using the parameter, depends only very weakly on the Reynolds number: is found to be about with variations around this mean value bounded by in all simulations. There is no systematic evolution with . For the best resolved model, the kinetic energy power spectrum tentatively displays a power-law range with an exponent , while the magnetic energy is found to shift to smaller and smaller scales as the magnetic Reynolds number increases. A couple of different diagnostics both suggest a well–defined injection length of a fraction of a scaleheight.

Conclusions:The results presented in this letter are consistent with the MRI being able to transport angular momentum efficiently at large Reynolds numbers when in unstratified zero net flux shearing boxes.

1 Introduction

Angular momentum transport in accretion disks has been an outstanding issue in theoretical astrophysics for decades. To date the most likely mechanism appears to be MHD turbulence driven by the magnetorotational instability (MRI, Balbus & Hawley, 1991, 1998). Several numerical simulations have been performed to study its properties. The most popular approach is to work in the local approximation, using the shearing box model, as pioneered by Hawley et al. (1995), Hawley et al. (1996) or Brandenburg et al. (1995). These early simulations have shown that MRI–powered MHD turbulence is a robust mechanism that transports angular momentum outward. The rate of transport, measured by the famous parameter (Shakura & Sunyaev, 1973) depends on the field geometry but is always positive, indicating outward flux of angular momentum. The results obtained in the 1990’s however obviously suffered from the limited computational resources available at that time. With no mean vertical magnetic field threading the shearing box (a field geometry referred to as the zero net flux case), Fromang & Papaloizou (2007) recently demonstrated with the code Zeus (Hawley & Stone, 1995) that it is indeed a problem: decreases by a factor of two each time the resolution is doubled. This behavior has since been shown to be very robust as it has been confirmed by simulations performed with codes using different algorithms (Simon et al., 2009; Guan et al., 2009). This result, although it raised the concern that MRI–induced transport could vanish at infinite resolution, was interpreted as an indication that the small scale behavior of the flow is an important ingredient to determine the rate of MRI–induced angular momentum transport: small scale explicit dissipation coefficients, namely viscosity and resistivity, need to be included in the simulations. With such calculations Lesur & Longaretti (2007) showed that, for a nonzero vertical mean magnetic field, rises with the magnetic Prandtl number , the ratio of viscosity over resistivity. This result is actually very general: it is independent of the field geometry and was also found for a mean toroidal magnetic field (Simon & Hawley, 2009) and in the zero net flux case of interest here (Fromang et al., 2007). Recently Simon et al. (2009) measured the numerical dissipation properties of the code Athena (Gardiner & Stone, 2008; Stone et al., 2008). They found that an increase in resolution amounts to an increase of the numerical Reynolds numbers, while keeping the effective magnetic Prandtl number (i.e. the ratio between the numerical viscosity and the numerical resistivity) roughly constant and equal to about two. In light of these results a possible interpretation of the findings of Fromang & Papaloizou (2007) is that is decreasing when the physical Reynolds number increases at fixed . If unchecked, this decreasing would mean that MRI–induced MHD turbulence is ineffective at transporting angular momentum without a mean flux, even in systems that have values higher than unity. Here, high resolution numerical simulations in which and are simultaneously increased while keeping their ratio constant are used to examine if this is indeed the case.

2 Numerical setup

Model Resolution
Table 1: Properties of the simulations and time averaged value of . The errors on are computed following Longaretti & Lesur (2010): the time history of is divided in N bins of size . is varied between and orbits. For each N, the standard deviation is computed according to , where is the mean value in bin i. For large N, scales like . The errors reported on use that scaling to estimate when orbits, the time duration over which the mean values of are calculated.

In the simulations described below, the non–ideal MHD equations (i.e. including viscosity and resistivity ) are solved in the unstratified shearing box (Goldreich & Lynden-Bell, 1965) by the code Zeus (Hawley & Stone, 1995). The setup is identical to that used by Fromang et al. (2007): the shearing box rotates around the central point mass with angular velocity (thus defining the orbital time ), the equation of state is isothermal with the sound speed , and the size of the box is fixed to , where is the disk scaleheight. As mentioned in the introduction, the magnetic flux threading the disk vanishes in all directions. Three simulations are presented here. They share the same value for the magnetic Prandtl number . The Reynolds number is gradually increased from (hereafter labeled model ) to (model ) and finally (model ). The resolution is increased at the same time as the Reynolds number to ensure that the smallest scale features of the flow are always resolved. Model is identical to model of Fromang et al. (2007), for which different diagnostics have shown that cells per scaleheight are sufficient when using Zeus. Thus the resolutions , and are adopted respectively for model , and 111Model , with cells per scaleheight, constitutes the highest resolution published so far of MRI induced turbulence. With about cells, the simulation required over million timesteps to be completed and a total of about CPU hours on the CEA supercomputer BULL Novascale 3045 hosted in France by CCRT. . For model , it was found that early transients associated with the linear instability kept affecting the flow for long times, resulting in prohibitively long simulations. For the computational cost of that simulation to remain acceptable, the following procedure was used: model was run from to orbits, starting from the initial state described above and identical to that used by Fromang et al. (2007). At , the flow was interpolated on a grid twice finer. The dissipation coefficients were reduced by a factor of two and the model was restarted between and orbits. This constitutes model . This procedure was repeated at time orbits to produce model . The latter was run between and orbits. The properties of the three models are summarized in Table 1: the first column gives the label of the model, the second column reports its resolution and the third the Reynolds number for that run. All models share the same value . Finally, the last column in Table 1 gives time–averaged values of that are discussed in the subsections below.

3 Flow properties

In the three models flow features typical of unstratified shearing boxes simulations are recovered: weakly non–axisymmetric density waves propagate radially in the box (Heinemann & Papaloizou, 2009a, b) on top of smaller scales velocity and magnetic field turbulent fluctuations, the latter exhibiting a tangled structure typical of values higher than unity (Schekochihin et al., 2004). Below we concentrate on the transport properties of the turbulence, the shape of the kinetic and magnetic energy power spectra and the the two points correlation function.

3.1 Angular momentum transport

Figure 1: Time history of the Maxwell stress tensor for model (blue dotted line), (green dashed line) and (red solid line). The three curves are consistent with the same time averaged value for , independently of the Reynolds number.

The angular momentum transport properties of the turbulence in the three models are assessed by calculating the parameter, the sum of the Reynolds stress tensor and the Maxwell stress tensor . All three coefficients are calculated as in Fromang & Papaloizou (2007). The time history of is shown in Fig. 1 for models , and respectively, using a dotted, a dashed and a solid line. The result is dramatically different from the results of Fromang & Papaloizou (2007) who found without explicit dissipation a monotonic decrease of as the resolution was increased. Here, no such systematic evolution is found as goes up. Indeed appears to vary only very weakly with the Reynolds number. This is confirmed by the last column of Table 1 in which the values of , time–averaged between and orbits, are reported for the different models. The rate of angular momentum transport appears to be somewhat smaller in model than in model and . Nevertheless, the difference between the three simulations remains less than . Taken together, the three measurements suggest that is of the order of and is fairly independent of the Reynolds number. At the very least, a systematic evolution of with is ruled out by the simulations.

3.2 Power spectrum

Figure 2: Top panel: kinetic (solid line) and magnetic (dashed line) energy power spectrum for model , time averaged over twenty snapshots between and . The dotted line shows a power law line with index for the purpose of comparison. Bottom panel: kinetic energy power spectra compensated by (dotted line), (solid line) and (dashed line). Both panels are suggestive of a spectrum in the range .
Figure 3: Top panel: plot of for model (dotted line), (dashed line) and (solid line). The dotted line shows a power–law line with the index for comparison. As and increase, the kinetic energy display an increasing region well–fitted by a power law with the index , while the viscous cut–off region moves to higher values. Bottom panel: Same as the top panel, but for the quantity . On both panels the insets reproduce the results of Fromang & Papaloizou (2007) for model STD64 (dotted line), STD128 (dashed line) and STD256 (solid line).

The top panel of Fig. 2 shows the shell–averaged kinetic energy power spectrum (solid line) and magnetic energy power spectrum (dashed line) for model . The latter is rather flat over about a decade in wavenumber (from to ) and is larger than the kinetic energy over that range. By contrast, the former displays a clear power–law behavior for wavenumber . For the purpose of comparison, the dotted line shows a pure power–law with the index that nicely fits the solid line of the plot. By analogy with hydrodynamic turbulence it is tempting to associate the large scale end of the power–law part of the spectrum with an injection length . Similarly, the small scale end can be associated with the viscous cut–off length and is found to be . This is about cells at that resolution and is thus well-resolved by the code. Furthermore, results obtained in the kinematic regime of incompressible and homogeneous MHD turbulence suggest that the resistive length (Schekochihin et al., 2004). Thus, is of order cells and also well resolved, which shows that numerical dissipation is most likely negligible in this simulation. Given the still limited resolution of model , the reliability of the power–law exponent mentioned above can however be questioned: for that purpose, the bottom panel of Fig. 2 displays three compensated spectra of , (dotted line), (solid line) and (dashed line) respectively. First, the figure illustrates the difficulty of a reliable determination of the exponent. Indeed, the power–law extends over less than a decade in wavenumber. Nevertheless, the dashed line, which unambiguously rises over the interval of , excludes a spectrum and rather suggests an exponent larger than . The dotted line on the other hand suggests as an upper limit. Finally, the solid line suggests as a tentative fit for the power–law range of the spectrum (). Finally, Fig. 3 (top panel) compares the shape of in model (dotted line), (dashed line) and (solid line). For all models, the kinetic energy power spectrum peaks at . For larger wavenumbers, the power–law becomes more and more apparent as the Reynolds number increases. The bottom panel of Fig. 3 plots the quantity for the three simulations. The peak of each curve thus provides an estimate of the scale at which magnetic energy is located. It is found to lie at , and respectively when , and . In other words, the scale at which most of the magnetic energy is located moves toward smaller and smaller scales as is increased. This is different from the results reported by Haugen et al. (2003), but not unexpected given existing theories of small scale dynamos with large (Schekochihin et al., 2002a, b). On both panels, the small insets plot the spectra obtained by Fromang & Papaloizou (2007) without explicit dissipation. Aside from the decrease of their amplitude with resolution, the most noticable differences with the results presented here are twofold: first, the kinetic energy power–spectra appear flatter at intermediate wavenumbers. In addition, there is more energy (both kinetic and magnetic) at the smallest scales of the box.

3.3 Correlation length

Figure 4: Structure of the correlation function in the plane for model (left panel), (middle panel) and (right panel). Its structure is only weakly dependent on the Reynolds number, suggesting a well–defined injection length.

The shape of the kinetic energy power–spectrum described above suggests an injection length that appears to be independent of for the range of the Reynolds numbers investigated here. However, the shell average involved in its derivation washes out all information about the anisotropy of the turbulence. This can be investigated using the two–points–correlation function (Guan et al., 2009; Davis et al., 2010):


Here . stands for a volume average, corresponds to the velocity fluctuations in the direction i and the sum is over spatial coordinates. Isocontours of in the plane are shown in Fig. 4. From left to right, the different panels correspond to models , and respectively. As found by Guan et al. (2009) and Davis et al. (2010), has an ellipsoidal shape for all models. The tilt angle of the major axis is degrees, which agrees with the results of Guan et al. (2009). Following the procedure outlined by these authors (i.e. by fitting the shape of the correlation function in a given direction by an exponential), the correlation lengths along the major, minor and vertical axis of the ellipsoid are found to be , independently of the Reynolds number. This suggests once more that the injection length of the turbulence is only weakly dependent on the Reynolds number. At the same time, the small difference between or and as quoted in section 3.2 is a warning that the injection range and the dissipative range might overlap in these simulations.

4 Conclusion

Here zero net flux high resolution numerical simulations of MRI–driven MHD turbulence are used to demonstrate this result: when , the dependence of on the Reynolds number is very weak. In all models, to within about . This result unambiguously shows that the decrease of with resolution reported by Fromang & Papaloizou (2007) is a numerical artifact that contains no physical information about the nature of the MHD turbulence in accretion disks. Quite differently, the present simulations are consistent with a nonzero value of at infinite Reynolds numbers for a magnetic Prandtl number higher than unity. Note that this weak dependence of with for is also suggested by the data recently reported by Simon & Hawley (2009) and Longaretti & Lesur (2010) respectively for a mean azimuthal and vertical magnetic field.

In addition, a number of statistical properties of the turbulence are reported. The kinetic energy power spectrum of the turbulence and the two–points–correlation function of the velocity both suggest a well–defined injection length of a few tens of a scaleheight. For the range of the Reynolds numbers that can be probed with current resources, seems to be independent of . At the highest resolution achieved here, the kinetic energy power spectrum displays a power–law scaling over almost a decade in wavenumber. However, given the limited extent of the power–law range, the precise exponent of this power–law cannot be accurately determined: an exponent of appears to be consistent with the data, while a exponent seems too steep. Nevertheless, as suggested in Sect. 3.3, the separation between the forcing and the dissipative scales might still be marginal. This is why a detailed comparison of these exponents with existing MHD turbulence theories (Iroshnikov, 1963; Kraichnan, 1965; Goldreich & Sridhar, 1995) is probably premature at this stage. Higher resolution simulations are definitively needed. Finally, the shape of the magnetic energy power spectrum shows that magnetic energy is mostly located at small scales and shifts to smaller and smaller scales as increases, as expected from small scale dynamo theory (Schekochihin et al., 2002a). This is consistent with the scenario postulated by Rincon et al. (2008) of a large scale MRI forcing that generates and coexists with a small scale dynamo.


The author acknowledges insightful discussions with F. Rincon, G. Lesur and P.–Y. Longaretti and is indebted to S. Pires for her help in analyzing the data presented here. These simulations were granted access to the HPC resources of CCRT under the allocation x2008042231 made by GENCI (Grand Equipement National de Calcul Intensif).


  • Balbus & Hawley (1991) Balbus, S. & Hawley, J. 1991, ApJ, 376, 214
  • Balbus & Hawley (1998) Balbus, S. & Hawley, J. 1998, Rev.Mod.Phys., 70, 1
  • Brandenburg et al. (1995) Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741
  • Davis et al. (2010) Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52
  • Fromang & Papaloizou (2007) Fromang, S. & Papaloizou, J. 2007, A&A, 476, 1113
  • Fromang et al. (2007) Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123
  • Gardiner & Stone (2008) Gardiner, T. A. & Stone, J. M. 2008, Journal of Computational Physics, 227, 4123
  • Goldreich & Lynden-Bell (1965) Goldreich, P. & Lynden-Bell, D. 1965, MNRAS, 130, 125
  • Goldreich & Sridhar (1995) Goldreich, P. & Sridhar, S. 1995, ApJ, 438, 763
  • Guan et al. (2009) Guan, X., Gammie, C. F., Simon, J. B., & Johnson, B. M. 2009, ApJ, 694, 1010
  • Haugen et al. (2003) Haugen, N. E. L., Brandenburg, A., & Dobler, W. 2003, ApJ, 597, L141
  • Hawley & Stone (1995) Hawley, J. & Stone, J. 1995, Comput. Phys. Commun., 89, 127
  • Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
  • Hawley et al. (1996) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690
  • Heinemann & Papaloizou (2009a) Heinemann, T. & Papaloizou, J. C. B. 2009a, MNRAS, 397, 52
  • Heinemann & Papaloizou (2009b) Heinemann, T. & Papaloizou, J. C. B. 2009b, MNRAS, 397, 64
  • Iroshnikov (1963) Iroshnikov, P. S. 1963, AZh, 40, 742
  • Kraichnan (1965) Kraichnan, R. H. 1965, Physics of Fluids, 8, 1385
  • Lesur & Longaretti (2007) Lesur, G. & Longaretti, P.-Y. 2007, MNRAS, 378, 1471
  • Longaretti & Lesur (2010) Longaretti, P. & Lesur, G. 2010, A&A, accepted
  • Rincon et al. (2008) Rincon, F., Ogilvie, G. I., Proctor, M. R. E., & Cossu, C. 2008, Astronomische Nachrichten, 329, 750
  • Schekochihin et al. (2002a) Schekochihin, A. A., Boldyrev, S. A., & Kulsrud, R. M. 2002a, ApJ, 567, 828
  • Schekochihin et al. (2004) Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276
  • Schekochihin et al. (2002b) Schekochihin, A. A., Maron, J. L., Cowley, S. C., & McWilliams, J. C. 2002b, ApJ, 576, 806
  • Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  • Simon & Hawley (2009) Simon, J. B. & Hawley, J. F. 2009, ApJ, 707, 833
  • Simon et al. (2009) Simon, J. B., Hawley, J. F., & Beckwith, K. 2009, ApJ, 690, 974
  • Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description