LETTER TO THE EDITOR
MHD simulations of the magnetorotational instability in a shearing box with zero net flux: the case Pm = 4
S. Fromang^{1,2}
1  CEA, Irfu, SAp, Centre de Saclay, 91191 GifsurYvette,
France
2  UMR AIM, CEACNRSUniv. Paris VII, Centre de Saclay, 91191 GifsurYvette, France
Received 18 February 2010 / Accepted 12 April 2010
Abstract
Aims. This letter investigates the transport properties of
MHD turbulence induced by the magnetorotational instability at large
Reynolds numbers Re when the magnetic Prandtl number Pm 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 Re is gradually increased from 3125 to 12 500 while simultaneously keeping Pm
= 4. To ensure that the small scale features of the flow are well
resolved, the resolution varies from 128 cells per disk
scaleheight to 512 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 15% in all simulations. There is no systematic evolution with Re.
For the best resolved model, the kinetic energy power spectrum
tentatively displays a powerlaw range with an exponent 3/2, 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 welldefined 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 Pm=4 in unstratified zero net flux shearing boxes.
Key words: accretion, accretion disks  magnetohydrodynamics (MHD)  methods: numerical
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 1998,1991). 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 MRIpowered 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 (Guan et al. 2009; Simon et al. 2009). This result, although it raised the concern that MRIinduced 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 MRIinduced 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 Pm, 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 Pm. If unchecked, this decreasing would mean that MRIinduced MHD turbulence is ineffective at transporting angular momentum without a mean flux, even in systems that have Pm values higher than unity. Here, high resolution numerical simulations in which Re and Rm are simultaneously increased while keeping their ratio Pm constant are used to examine if this is indeed the case.
2 Numerical setup
In the simulations described below, the nonideal MHD equations (i.e. including viscosity and resistivity ) are solved in the unstratified shearing box (Goldreich & LyndenBell 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 c_{0}, 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 Re=3125 (hereafter labeled model Re3125) to Re=6250 (model Re6250) and finally Re=12 500 (model Re12500). 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 Re3125 is identical to model 128Re3125Pm4 of Fromang et al. (2007), for which different diagnostics have shown that 128 cells per scaleheight are sufficient when using Zeus. Thus the resolutions (N_{x},N_{y},N_{z})=(128,192,128), (256,384,256) and (512,768,512) are adopted respectively for model Re3125, Re6250 and Re12500^{}.
Table 1: Properties of the simulations and time averaged value of .
For model Re12500, 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 Re3125 was run from t=0 to t=150 orbits, starting from the initial state described above and identical to that used by Fromang et al. (2007). At t=60, 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 t=60 and t=150 orbits. This constitutes model Re6250. This procedure was repeated at time t=90 orbits to produce model Re12500. The latter was run between t=90 and t=135 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 (N_{x},N_{y},N_{z}) and the third the Reynolds number Re for that run. All models share the same value Pm=4. Finally, the last column in Table 1 gives timeaveraged 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 nonaxisymmetric density waves propagate radially in the box (Heinemann & Papaloizou 2009b,a) on top of smaller scales velocity and magnetic field turbulent fluctuations, the latter exhibiting a tangled structure typical of Pm 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 Re3125 ( blue dotted line), Re6250 ( green dashed line) and Re12500 ( red solid line). The three curves are consistent with the same time averaged value for , independently of the Reynolds number. 

Open with DEXTER 
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 Re3125, Re6250 and Re12500 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 Re 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 , timeaveraged between 90 and 130 orbits, are reported for the different models. The rate of angular momentum transport appears to be somewhat smaller in model Re6250 than in model Re3125 and Re12500. Nevertheless, the difference between the three simulations remains less than 25%. 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 Re 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 Re12500, time averaged over twenty snapshots between t=90 and t=120. The dotted line shows a power law line with index 3/2 for the purpose of comparison. Bottom panel: kinetic energy power spectra compensated by 1 ( dotted line), 3/2 ( solid line) and 5/3 ( dashed line). Both panels are suggestive of a k^{3/2} spectrum in the range 30<k<100. 

Open with DEXTER 
Figure 3: Top panel: plot of for model Re3125 ( dotted line), Re6250 ( dashed line) and Re12500 ( solid line). The dotted line shows a powerlaw line with the index 3/2 for comparison. As Re and Rm increase, the kinetic energy display an increasing region wellfitted by a power law with the index 3/2, while the viscous cutoff region moves to higher k 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). 

Open with DEXTER 
The top panel of Fig. 2 shows the shellaveraged kinetic energy power spectrum (solid line) and magnetic energy power spectrum (dashed line) for model Re12500. 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 powerlaw behavior for wavenumber 20<k<100. For the purpose of comparison, the dotted line shows a pure powerlaw with the index 3/2 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 powerlaw part of the spectrum with an injection length . Similarly, the small scale end can be associated with the viscous cutoff length and is found to be . This is about 32 cells at that resolution and is thus wellresolved 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 16 cells and also well resolved, which shows that numerical dissipation is most likely negligible in this simulation. Given the still limited resolution of model Re12500, the reliability of the powerlaw 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 powerlaw extends over less than a decade in wavenumber. Nevertheless, the dashed line, which unambiguously rises over the interval of 10<k<100, excludes a k^{5/3} spectrum and rather suggests an exponent larger than 5/3. The dotted line on the other hand suggests 1 as an upper limit. Finally, the solid line suggests k^{3/2} as a tentative fit for the powerlaw range of the spectrum (30<k<100). Finally, Fig. 3 (top panel) compares the shape of in model Re3125 (dotted line), Re6250 (dashed line) and Re12500 (solid line). For all models, the kinetic energy power spectrum peaks at 20. For larger wavenumbers, the k^{3/2} powerlaw 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 40, 5060 and 7080 respectively when Re=3125, 6250 and 12500. In other words, the scale at which most of the magnetic energy is located moves toward smaller and smaller scales as Rm 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 Pm(Schekochihin et al. 2002b,a). 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 powerspectra 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 Re3125 ( left panel), Re6250 ( middle panel) and Re12500 ( right panel). Its structure is only weakly dependent on the Reynolds number, suggesting a welldefined injection length. 

Open with DEXTER 
The shape of the kinetic energy powerspectrum described above
suggests an injection length
that appears to be independent
of Re 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 twopointscorrelation function (Davis et al. 2010; Guan et al. 2009):
(1) 
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 Re3125, Re6250 and Re12500 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 Sect. 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 MRIdriven MHD turbulence are used to demonstrate this result: when Pm =4, the dependence of on the Reynolds number is very weak. In all models, to within about 15%. 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 Re for Pm>1 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 twopointscorrelation function of the velocity both suggest a welldefined injection length of a few tens of a scaleheight. For the range of the Reynolds numbers Re that can be probed with current resources, seems to be independent of Re. At the highest resolution achieved here, the kinetic energy power spectrum displays a powerlaw scaling over almost a decade in wavenumber. However, given the limited extent of the powerlaw range, the precise exponent of this powerlaw cannot be accurately determined: an exponent of 3/2 appears to be consistent with the data, while a 5/3 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 (Goldreich & Sridhar 1995; Iroshnikov 1963; Kraichnan 1965) 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 Rm 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.
AcknowledgementsThe 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).
References
 Balbus, S., & Hawley, J. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S., & Hawley, J. 1998, Rev. Mod. Phys., 70, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Nordlund, A., Stein, R. F., et al. 1995, ApJ, 446, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., & Papaloizou, J. 2007, A&A, 476, 1113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., Papaloizou, J., Lesur, G., et al. 2007, A&A, 476, 1123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gardiner, T. A., & Stone, J. M. 2008, J. Comput. Phys., 227, 4123 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Guan, X., Gammie, C. F., Simon, J. B., et al. 2009, ApJ, 694, 1010 [NASA ADS] [CrossRef] [Google Scholar]
 Haugen, N. E. L., Brandenburg, A., & Dobler, W. 2003, ApJ, 597, L141 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J., & Stone, J. 1995, Comput. Phys. Commun., 89, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690 [NASA ADS] [CrossRef] [Google Scholar]
 Heinemann, T., & Papaloizou, J. C. B. 2009a, MNRAS, 397, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Heinemann, T., & Papaloizou, J. C. B. 2009b, MNRAS, 397, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Iroshnikov, P. S. 1963, AZh, 40, 742 [NASA ADS] [Google Scholar]
 Kraichnan, R. H. 1965, Phys. Fluids, 8, 1385 [NASA ADS] [CrossRef] [Google Scholar]
 Lesur, G., & Longaretti, P.Y. 2007, MNRAS, 378, 1471 [NASA ADS] [CrossRef] [Google Scholar]
 Longaretti, P., & Lesur, G. 2010, A&A, in press [Google Scholar]
 Rincon, F., Ogilvie, G. I., Proctor, M. R. E., et al. 2008, Astron. Nachr., 329, 750 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Boldyrev, S. A., & Kulsrud, R. M. 2002a, ApJ, 567, 828 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Maron, J. L., Cowley, S. C., et al. 2002b, ApJ, 576, 806 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Simon, J. B., & Hawley, J. F. 2009, ApJ, 707, 833 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Hawley, J. F., & Beckwith, K. 2009, ApJ, 690, 974 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
Footnotes
 ...12500^{}
 Model Re12500, with 512 cells per scaleheight, constitutes the highest resolution published so far of MRI induced turbulence. With about cells, the simulation required over 1.4 million timesteps to be completed and a total of about 350 000 CPU hours on the CEA supercomputer BULL Novascale 3045 hosted in France by CCRT.
All Tables
Table 1: Properties of the simulations and time averaged value of .
All Figures
Figure 1: Time history of the Maxwell stress tensor for model Re3125 ( blue dotted line), Re6250 ( green dashed line) and Re12500 ( red solid line). The three curves are consistent with the same time averaged value for , independently of the Reynolds number. 

Open with DEXTER  
In the text 
Figure 2: Top panel: kinetic ( solid line) and magnetic ( dashed line) energy power spectrum for model Re12500, time averaged over twenty snapshots between t=90 and t=120. The dotted line shows a power law line with index 3/2 for the purpose of comparison. Bottom panel: kinetic energy power spectra compensated by 1 ( dotted line), 3/2 ( solid line) and 5/3 ( dashed line). Both panels are suggestive of a k^{3/2} spectrum in the range 30<k<100. 

Open with DEXTER  
In the text 
Figure 3: Top panel: plot of for model Re3125 ( dotted line), Re6250 ( dashed line) and Re12500 ( solid line). The dotted line shows a powerlaw line with the index 3/2 for comparison. As Re and Rm increase, the kinetic energy display an increasing region wellfitted by a power law with the index 3/2, while the viscous cutoff region moves to higher k 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). 

Open with DEXTER  
In the text 
Figure 4: Structure of the correlation function in the plane for model Re3125 ( left panel), Re6250 ( middle panel) and Re12500 ( right panel). Its structure is only weakly dependent on the Reynolds number, suggesting a welldefined injection length. 

Open with DEXTER  
In the text 
Copyright ESO 2010