EDP Sciences
Free Access
Issue
A&A
Volume 514, May 2010
Article Number L5
Number of page(s) 4
Section Letters
DOI https://doi.org/10.1051/0004-6361/201014284
Published online 19 May 2010
A&A 514, L5 (2010)

LETTER TO THE EDITOR

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

S. Fromang1,2

1 - CEA, Irfu, SAp, Centre de Saclay, 91191 Gif-sur-Yvette, France
2 - UMR AIM, CEA-CNRS-Univ. Paris VII, Centre de Saclay, 91191 Gif-sur-Yvette, 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 $\alpha $ parameter, depends only very weakly on the Reynolds number: $\alpha $ is found to be about $7 \times 10^{-3}$ 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 power-law 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 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 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 MRI-powered MHD turbulence is a robust mechanism that transports angular momentum outward. The rate of transport, measured by the famous $\alpha $ 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: $\alpha $ 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 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, $\alpha $ 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 $\alpha $ is decreasing when the physical Reynolds number increases at fixed Pm. If unchecked, this decreasing $\alpha $ would mean that MRI-induced 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 non-ideal MHD equations (i.e. including viscosity $\nu$ and resistivity $\eta$) 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 $\Omega$ (thus defining the orbital time $T_{\rm orb} = 2 \pi / \Omega$), the equation of state is isothermal with the sound speed c0, and the size of the box is fixed to $(L_x,L_y,L_z)=(H,\pi H,H)$, where $H=c_0/\Omega$ 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 $Pm=\nu/\eta=4$. The Reynolds number $Re=c_0H/\nu$ 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 (Nx,Ny,Nz)=(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 $\alpha $.

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 (Nx,Ny,Nz) 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 time-averaged values of $\alpha $ 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 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

\begin{figure}
\par\includegraphics[scale=0.38]{14284f1.eps}
\end{figure} 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 $\alpha _{\rm Max}$, 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 $\alpha $ parameter, the sum of the Reynolds stress tensor $\alpha_{\rm Rey}$ and the Maxwell stress tensor $\alpha _{\rm Max}$. All three coefficients are calculated as in Fromang & Papaloizou (2007). The time history of $\alpha _{\rm Max}$ 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 $\alpha _{\rm Max}$ as the resolution was increased. Here, no such systematic evolution is found as Re goes up. Indeed $\alpha _{\rm Max}$ 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 $\alpha $, time-averaged 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 $\alpha $ is of the order of $7 \times 10^{-3}$ and is fairly independent of the Reynolds number. At the very least, a systematic evolution of $\alpha $ with Re is ruled out by the simulations.

3.2 Power spectrum

\begin{figure}
\par\includegraphics[scale=0.4]{14284f2.ps}\par\includegraphics[scale=0.4]{14284f3.ps}
\end{figure} 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

\begin{figure}
\par\includegraphics[scale=0.4]{14284f4.ps}\par\includegraphics[scale=0.4]{14284f5.ps}
\end{figure} Figure 3:

Top panel: plot of $E_{\rm K}$ for model Re3125 ( dotted line), Re6250 ( dashed line) and Re12500 ( solid line). The dotted line shows a power-law line with the index -3/2 for comparison. As Re and Rm increase, the kinetic energy display an increasing region well-fitted by a power law with the index -3/2, while the viscous cut-off region moves to higher k values. Bottom panel: same as the top panel, but for the quantity $k E_{\rm M}$. 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 shell-averaged kinetic energy power spectrum $E_{\rm K}$ (solid line) and magnetic energy power spectrum $E_{\rm M}$ (dashed line) for model Re12500. The latter is rather flat over about a decade in wavenumber (from $k \sim 10$ to $k
\sim 100$) and is larger than the kinetic energy over that range. By contrast, the former displays a clear power-law behavior for wavenumber 20<k<100. For the purpose of comparison, the dotted line shows a pure power-law 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 power-law part of the spectrum with an injection length $l_{\rm inj} \sim 2\pi / k_{\rm min} \sim
0.3~H$. Similarly, the small scale end can be associated with the viscous cut-off length and is found to be $l_{\rm visc} \sim 2\pi / k_{\rm visc}
\sim 0.06~H$. This is about 32 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 $l_{\rm res} \sim
Pm^{-1/2}l_{\rm visc}$ (Schekochihin et al. 2004). Thus, $l_{\rm res}$ 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 power-law exponent mentioned above can however be questioned: for that purpose, the bottom panel of Fig. 2 displays three compensated spectra of $E_{\rm K}$, $kE_{\rm K}$ (dotted line), $k^{3/2}E_{\rm K}$ (solid line) and $k^{5/3}E_{\rm K}$ (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 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 power-law range of the spectrum (30<k<100). Finally, Fig. 3 (top panel) compares the shape of $E_{\rm K}$ in model Re3125 (dotted line), Re6250 (dashed line) and Re12500 (solid line). For all models, the kinetic energy power spectrum peaks at $k \sim 10$-20. For larger wavenumbers, the k-3/2 power-law becomes more and more apparent as the Reynolds number increases. The bottom panel of Fig. 3 plots the quantity $k E_{\rm M}$ 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 $k_{\rm peak} \sim 30$-40, 50-60 and 70-80 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 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

\begin{figure}
\par\includegraphics[scale=0.22]{14284f6.eps}\includegraphics[sca...
...284f7.eps}\includegraphics[scale=0.22]{14284f8.eps}
\vspace*{2.5mm}
\end{figure} Figure 4:

Structure of the correlation function $\xi _v$ in the $(\Delta _x,\Delta _y)$ 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 well-defined injection length.

Open with DEXTER

The shape of the kinetic energy power-spectrum described above suggests an injection length $l_{\rm inj}$ 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 two-points-correlation function (Davis et al. 2010; Guan et al. 2009):

\begin{displaymath}\xi_v(\mbox{\boldmath {$\Delta x$ }})= \langle \Sigma_i \delt...
... {$\Delta x$ }})\rangle/ \langle\Sigma_i \delta
v_i^2 \rangle.
\end{displaymath} (1)

Here $\langle$.$ \rangle$ stands for a volume average, $\delta v_i(\mbox{\boldmath {$x$ }})$corresponds to the velocity fluctuations in the direction i and the sum is over spatial coordinates. Isocontours of $\xi _v$ in the plane $\Delta z=0$ 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), $\xi _v$ has an ellipsoidal shape for all models. The tilt angle $\theta_v$ of the major axis is $\theta_v \sim 8$ 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 $(\lambda_{\min},\lambda_{\max},\lambda_{z})=(0.08,0.45,0.08)H$, 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 $\lambda_{\rm min}$ or $\lambda_{z}$ and $l_{\rm visc}$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 MRI-driven MHD turbulence are used to demonstrate this result: when Pm =4, the dependence of $\alpha $ on the Reynolds number is very weak. In all models, $\alpha \sim 7 \times 10^{-3}$ to within about 15%. This result unambiguously shows that the decrease of $\alpha $ 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 $\alpha $ at infinite Reynolds numbers for a magnetic Prandtl number higher than unity. Note that this weak dependence of $\alpha $ 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 two-points-correlation function of the velocity both suggest a well-defined injection length $l_{\rm inj}$ of a few tens of a scaleheight. For the range of the Reynolds numbers Re that can be probed with current resources, $l_{\rm inj}$ seems to be independent of Re. 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 -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.

Acknowledgements
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).

References

Footnotes

...12500[*]
Model Re12500, with 512 cells per scaleheight, constitutes the highest resolution published so far of MRI induced turbulence. With about $2 \times 10^8$ 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 $\alpha $.

All Figures

  \begin{figure}
\par\includegraphics[scale=0.38]{14284f1.eps}
\end{figure} 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 $\alpha _{\rm Max}$, independently of the Reynolds number.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[scale=0.4]{14284f2.ps}\par\includegraphics[scale=0.4]{14284f3.ps}
\end{figure} 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

  \begin{figure}
\par\includegraphics[scale=0.4]{14284f4.ps}\par\includegraphics[scale=0.4]{14284f5.ps}
\end{figure} Figure 3:

Top panel: plot of $E_{\rm K}$ for model Re3125 ( dotted line), Re6250 ( dashed line) and Re12500 ( solid line). The dotted line shows a power-law line with the index -3/2 for comparison. As Re and Rm increase, the kinetic energy display an increasing region well-fitted by a power law with the index -3/2, while the viscous cut-off region moves to higher k values. Bottom panel: same as the top panel, but for the quantity $k E_{\rm M}$. 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

  \begin{figure}
\par\includegraphics[scale=0.22]{14284f6.eps}\includegraphics[sca...
...284f7.eps}\includegraphics[scale=0.22]{14284f8.eps}
\vspace*{2.5mm}
\end{figure} Figure 4:

Structure of the correlation function $\xi _v$ in the $(\Delta _x,\Delta _y)$ 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 well-defined injection length.

Open with DEXTER
In the text


Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.