1 - INAF, Osservatorio Astronomico di Torino, Strada Osservatorio 20, Pino Torinese, Italy
2 - Dipartimento di Fisica Generale, Universitá di Torino, via Pietro Giuria 1, 10125 Torino, Italy
3 - Department of Astronomy and Astrophysics, The University of Chicago, 5640 S. Ellis ave., Chicago IL 60637, USA

Abstract
Aims. We study the changes in the properties of turbulence driven by the magnetorotational instability in a shearing box, as the computational domain size in the radial direction is varied relative to the height.
Methods. We perform 3D simulations in the shearing box approximation, with a net magnetic flux, and we consider computational domains with different aspect ratios.
Results. We find that in boxes of aspect ratio unity the transport of angular momentum is strongly intermittent and dominated by channel solutions in agreement with previous work. In contrast, in boxes with larger aspect ratios, the channel solutions and the associated intermittent behavior disappear.
Conclusions. There is strong evidence that, as the aspect ratio becomes larger, the characteristics of the solution become aspect ratio independent. We conclude that shearing box calculations with an aspect ratio of unity or near unity may introduce spurious effects.

Key words: accretion, accretion disks - instabilities - turbulence - magnetic fields - magnetohydrodynamics (MHD)

Understanding the process of angular momentum transport is one of the fundamental issues in the physics of accretion disks. Turbulence was recognized as the main source of the required enhanced transport (Shakura & Sunyaev 1973), but its driving mechanism remained unclear until Balbus & Hawley (1991) proposed the magnerotational instability (MRI). Much of what is presently known about MRI driven turbulence has been obtained in the shearing box approximation in which only a small periodic patch of the disk is simulated (see e.g. Sano et al. 2004; Hawley et al. 1995; Lesur & Longaretti 2007; Sano & Stone 2002; Sano & Inutsuka 2001; Balbus 2003; Turner et al. 2003). The advantage of this local approach, as opposed to a global disk simulation, is the possibility of reaching much higher resolutions at the same computational cost. The shearing box approach is of course meaningful only if it captures the characteristics of MRI driven turbulence, and of the related angular momentum transport in a full disk. A critical discussion of the validity of the shearing box approximation can be found in Regev & Umurhan (2007).

Ideally one should compare local and global results; however, available global disk simulations (see e.g. Hawley 2001; Hawley et al. 2001) have relatively low resolution while simulations with adequately high resolution are still beyond present capabilities. Given the present resources, one should at least check the self consistency of the shearing box results. The shearing box equations are obtained as a formal local expansion in the limit of large radii and small gap (Hill 1878), but there is no guarantee that the solutions thus obtained satisfy the same locality conditions. This should be verified a posteriori, by checking, for example, that the properties of the solutions do not depend on the size of the computational domain.

$\begin{figure} \par\includegraphics[width=16cm,clip]{9730fig1.eps} \end{figure}$	Figure 1: Time histories of the Maxwell stresses averaged over the computational box and normalized to the pressure for the low resolution simulations. The three panels correspond to cases with L=1, 4, and 8. Note that time is in units of the rotation time.
Open with DEXTER

It is useful to distinguish two related issues. One concerns the size of the computational domain relative to some characteristic length scale of the problem determined by the physical parameters. For instance, in the MRI case, the vertical wavelength of the mode of maximum growth rate depends on the strength of the applied uniform magnetic field. For a given field strength then, the vertical size of the computational domain is chosen to contain some multiple of that wavelength. In this case, varying the field strength is in some sense equivalent to varying the (vertical) box size. A related issue concerns the aspect ratio, i.e. the size of the computational domain in the radial and azimuthal directions relative to the vertical size. The modes of maximum growth rate do not offer much guidance in this respect since they only vary in the vertical. Clearly a very thin box may introduce spurious effects, a very wide box rapidly becomes computationally expensive. A typical compromise is to adopt boxes with aspect ratios between the radial and vertical direction of unity. Presumably, the basis for this choice is that once nonlinear processes saturate the exponential growth, the typical resulting structures will have a characteristic size in the radial direction comparable to or not much larger than that in the vertical direction. This being the case, a computational domain of unit aspect ratio is adequate to capture the properties of the solution. Strangely, this point has not received careful consideration and there has been no systematic study of the dependence of the solutions on the aspect ratio. Furthermore, many studies in computational domains of unit aspect ratio have shown intermittent behavior associated with the formation and subsequent disruption of ``channel'' solutions (Sano & Inutsuka 2001; Hawley & Balbus 1992). These are the nonlinear analogues of the exponentially growing linear modes. Recently a number of authors have addressed the effect of the presence or absence of channel solutions in different contexts, for example, by introducing vertical gravity and stratification (Coppi & Keyes 2003) or changing the radial boundary conditions (Liu et al. 2006; Umurhan et al. 2007a; Umurhan et al. 2007b). Here we note that, interestingly, in the nonlinear regime, the observed channel solution is not necessarily the one that corresponds to the linear mode of maximum growth rate. In fact, often, the observed channel solution is the one whose vertical extent fills the computational domain once. In this case, it is not a priori obvious what aspect ratio should be used. In this paper we address these issues by carrying out a systematic study of shearing box results as a function of aspect ratio. In particular, we want to see if the results observed in boxes of unit aspect ratio are representative of more extended systems, or if they display peculiarities induced by an overly constrained geometry.

2 Formulation

We perform a series of 3D, compressible, isothermal numerical simulations in the shearing box approximation (for a detailed description of the shearing box model, see Hawley et al. 1995). The Cartesian coordinates x, y, z are defined around a fiducial point comoving with angular velocity $\Omega$ and corresponding, respectively, to the radial, azimuthal and vertical directions. An equilibrium solution is given by a uniform shear flow, $\vec{v} = -q\Omega x\hat{\vec{e}}_y$ , with q = 3/2 for a Keplerian disk, threaded by a uniform vertical magnetic field, $\vec{B} = B_0\hat{\vec{e}}_z$ , in which the density $\rho$ and the pressure p are constant ( $\rho = 1$ , $p = c_{\rm s}^2\rho$ ). We have assumed, as done in most of the literature on the subject, the simplest case with no vertical stratification and gravity. This may limit the applicability to real disks (Regev & Umurhan 2007), but it is adequate for our purpose of studying the aspect ratio dependence. This equilibrium is initially perturbed by small amplitude, spatially uncorrelated azimuthal velocity perturbations. The box has size $L~\times~~4~\times~~1$ in the x, y and z directions, respectively, with L varying from 1 to 8.

We take $1/\Omega$ as the unit of time (note however that in the plots time is in units of the rotation time $2 \pi / \Omega$ ), and in all the simulations the sound speed $c_{\rm s}$ and the plasma $\beta=p_{\rm gas}/p_{\rm mag}$ are fixed to values of 4.56 and 10⁴, respectively. With these choices the fastest growing MRI mode has a vertical wavelength close to 1/3.

Computations are carried out at both low and high resolutions (32 and 128 zones per unit length, respectively) on equally spaced grids. The MHD equations are solved in conservative form using the isothermal MHD module available in the PLUTO code (Mignone et al. 2007; Mignone 2007). The latter is a finite volume, Riemann solver based code in which the evolution of the magnetic field is carried out using the constrained transport method of Balsara & Spicer (1999). In the present formulation we do not include explicit viscosity or magnetic diffusivity, and rely solely on numerical dissipation to limit the potential unbounded growth of the solutions (for a theoretical discussion of numerical dissipation in Riemann solver based schemes see, for example, Toro (1999), for a discussion of the role of numerical dissipation in turbulent simulations see Grinstein et al. 2007).

3 Results

Our main objective is to study how the properties of the turbulent solutions vary as the (radial) aspect ratio L is varied from 1 to 4 to 8. We recall that in simulations with an aspect ratio of unity, the angular momentum transport shows an intermittent behavior with episodes of high transport. These correspond to states, loosely referred to as ``channel'' solutions, in which the velocity and magnetic perturbations are highly spatially correlated. Strictly speaking, the channel solutions are exponentially growing exact solutions of the incompressible shearing box equations in which the velocity and magnetic field perturbations depend sinusoidally on z, have directions at right angle to each other, and a growth rate determined by the angle they make to the azimuthal direction (Goodman & Xu 1994). More commonly the term ``channel solutions'' is used to describe a behavior with similar properties to that observed in the exact solutions. It was shown by Goodman & Xu (1994) that the channel solutions, as their amplitudes grow, become unstable to parasitic instabilities. It is now believed that the formation of near channel solutions and their subsequent disruption by the parasitic instabilities are the basis for the observed intermittent behavior (Sano & Inutsuka 2001). The vertical wavenumber determines the angle $\gamma_{\rm c}$ between the direction of the magnetic perturbations and the azimuthal direction, and hence the growth rate. Although, in principle, channel solutions exist with any wavenumber, the most commonly observed ones have a vertical extent equal to the vertical box size, for which the angle $\gamma_{\rm c} \approx 23^\circ$ . We shall use these results below.

We start our discussion with the low resolution results and compare three cases with different aspect ratios. Figure 1 shows the the volume averaged Maxwell stresses $\langle w_{xy}\rangle=-\langle B_x B_y\rangle$ as a function of time. In what follows we concentrate on the Maxwell stresses because they give the largest contribution to the total stresses exceeding the Reynolds stresses on average by a factor of about 5. We have normalized the stresses to the pressure, thus the values shown in the figure correspond to the $\alpha$ parameter (Shakura & Sunyaev 1973). The left panel, corresponding to L = 1, shows the characteristic intermittent behavior described above. A striking difference appears if we compare this case with the other two having larger aspect ratios in which the spikes are absent. A more quantitative view of the differences between the three cases is afforded by their respective probability distribution functions shown in Fig. 2. The case L = 1 (solid curve) is quite distinct from the other two with a long tail corresponding to the episodes of enhanced transport. Careful examination shows that a small tail is still present in the case with L=4 (dashed curve), and absent in the case with L=8 (dot-dashed curve). However, we see that the differences between the last two cases are much smaller than the differences to the L = 1 case, indicating that we have almost reached a converged behavior. Because of the presence of the peaks, the time-averaged value of $\alpha$ for L=1 exceeds the other two cases by approximately a factor of 2. It is natural to assume that the distinctive behavior of the case with L=1 is due to the presence of channel solutions that are otherwise absent in the simulations with larger aspect ratios. This can be verified by looking for typical imprints of the channel solutions such as a high correlation coefficient between directional components of the magnetic field, or a specific relationship between the vertical wavenumber and the angle between the magnetic field and the azimuthal direction. In order to define these quantities it is useful to introduce the idea of a scatter diagram whereby for each gridpoint the value of B_y is plotted as a function of B_x, Fig. 4. For a pure channel solution this procedure would yield a straight line through the origin with slope $\tan \gamma_{\rm c}$ . In the general case, the points will not lie on a straight line; still, a least-square straight line fit through the points can be used to define the average angle, and the spread about this line would give a measure of the correlation coefficient R. Equivalently, these quantities can be defined by

$\begin{displaymath}\tan \gamma = -\frac{\langle B_x B_y\rangle}{\langle B_x^2\ra... ... B_y\rangle^2}{\langle B_x^2\rangle \langle B_y^2\rangle}\cdot \end{displaymath}$

(1)

$\begin{figure} \par\includegraphics[width=8cm,clip]{9730fig2.eps}\end{figure}$	Figure 2: Probability distribution functions of the Maxwell stresses for the same cases as in Fig. 1. The three curves correspond to the three values of the aspect ratio: L = 1 - solid curve, L = 4 - dashed curve, and L = 8 - dot-dashed.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=14.5cm,clip]{9730fig3.eps}\end{figure}$	Figure 3: Time histories of the average slope $\tan \gamma$ ( upper row) and correlation coefficient R ( lower row) for the same cases as in Fig. 1. The straight horizontal lines indicate the theoretical value of $\tan \gamma_{\rm c}$ for the exact channel solution with wavelength unity. Note that time is in units of the rotation time.
Open with DEXTER

Differences between the various cases can be further illustrated by scatter plots of the distribution of B_y as a function of B_x like those presented in Fig. 4. The two panels on the left show two examples of distributions for the case L = 1, one corresponding to a maximum of $\langle w_{xy}\rangle$ and the other corresponding to a minimum. This confirms that in going from the minimum to the maximum the magnetic field fluctuations increase their intensity and the x and y components tend to become more correlated. The velocity fluctuations show a similar behavior, with almost no correlation in the minimum state. In contrast, for L = 4, there are no significant variations in the distributions during the evolution, see Fig. 4 where we display only one representative example that shows a similarity to the minimum state of the case with L=1. The figure also makes apparent why the Reynolds stresses are significantly lower than the Maxwell stresses. The difference partly arises because the intensity of the fluctuations is smaller, but to a somewhat larger extent because the velocity fluctuations, except during a channel solution phase, remain uncorrelated.

$\begin{figure} \par\includegraphics[width=14.2cm,clip]{9730fig4.eps}\end{figure}$	Figure 4: Scatter diagrams of horizontal magnetic fluctuations (dark tones) and velocities (lighter tones). The first panel corresponds to the case with L=1 at an instant near a maximum in the stress. The second panel corresponds to the same case near a minimum. The third panel corresponds to a representative instant for the case with L=4.
Open with DEXTER

Another property of the distributions not entirely evident from the figure is that, while for L = 4 and at the minimum of the L = 1 case the distributions peak at zero value, at the maximum of the L = 1 case there are two peaks around the largest values of the fluctuations. This is clear in Fig. 5 where we show the probability distribution function of the azimuthal magnetic field component B_y. The panels show the behavior at an instant near a maximum in the stress (solid line) and a minimum (dashed line) for two different aspect ratios. In the left panel, we also show (dot-dashed line) the corresponding distribution for the channel solution. The presence of the double peak distribution in the case with L = 1 (left panel) can be considered again as an indication of the presence of the channel solution around the maxima in the stress, as it can be seen by comparing the solid and the dot-dashed lines.

Recently, the work of Fromang & Papaloizou (2007) has raised some concerns about the dependence of numerical studies of MRI turbulence on resolution. Accordingly, we have repeated the simulations with L=1, 4 at a higher resolution of 128 gridpoints per unit length. We find the results to be qualitatively the same. For example, Fig. 6 shows the time histories of $\langle w_{xy}\rangle$ for the higher resolution simulations. The behavior is strikingly similar to that of Fig. 1, with a strong presence of the channel solutions in the L=1 case and a complete absence of such solutions in the runs with larger aspect ratio. We have also checked other indicators, like the probability distribution functions and the scatter diagrams and found them to be qualitatively the same as in the cases described above. We should note however that quantitatively we do observe a dependence on resolution; in particular the average value of the Maxwell stresses increases with increasing resolution (cf. the zero flux case Fromang & Papaloizou 2007). This follows mainly from an increase in the peak values of magnetic fluctuations, which is not surprising since the effective magnetic Reynolds number associated with numerical dissipation increases with increasing resolution.

4 Conclusion

Channel solutions represent an highly correlated state for which the transport is very efficient and are a characteristic feature that strongly affects the behavior of the angular momentum transport in most of the shearing box simulations presented so far in the literature. There is however the question of whether the dominance of this state could not be due to the constraints imposed by the typical shearing box approach. Recently, several authors have addressed this question considering the effects of vertical gravity and stratification (Coppi & Keyes 2003), and of different boundary conditions (Liu et al. 2006; Umurhan et al. 2007a; Umurhan et al. 2007b). The general result is that the channel solution is not robust and disappears whenever one introduces some modification to the typical shearing box approach and, as a consequence, the angular momentum transport is less efficient.

In this paper we have focused our analysis on the effects of a change of the box aspect ratio and the results above show two different behaviors depending on it. For aspect ratios close to unity the solutions are strongly affected by the channel solutions, with frequent episodes of high correlation and efficient transport. For larger aspect ratios the channel solutions disappear, the system remains in the state with lower correlation and the average transport of angular momentum is correspondingly reduced. Further increase in aspect ratio does not lead to any significant changes. Thus we conclude that the less correlated state is more likely to be representative of the extended system.

$\begin{figure} \par\includegraphics[width=8cm,clip]{9730fig5.eps} \end{figure}$	Figure 5: Probability distribution functions for the azimuthal component of magnetic field B_y. The left panel refers to the case with L=1. Solid and dashed lines correspond to the first two panels in Fig. 4 while the dot-dashed line corresponds to the exact channel solution. The right panel corresponds to the third panel in Fig. 4.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=7cm,clip]{9730fig6.eps} \end{figure}$	Figure 6: Time histories of the Maxwell stresses averaged over the computational box and normalized to the pressure for the high resolution simulations. The two panels correspond to cases with L=1, and 4. Note that time is in units of the rotation time.
Open with DEXTER

It is interesting to speculate why the channel solutions disappear from boxes with larger aspect ratios. One possibility is to assume that nontrivial correlations are necessary to form the channel solutions and that these cannot develop faster than the sound crossing time. Accordingly, in a large box it takes longer to form the channel solutions; if the rate at which they are destroyed is fixed in a sufficiently large box they may never form. We do not believe this is the correct explanation, since we have found the same behavior in simulations with a sound speed ten times larger. Another, more likely explanation is that the secondary instabilities responsible for the destruction of the channel solutions are to some extent suppressed at small aspect ratios, but can achieve their maximal growth rate once the aspect ratio is large enough. In support of this idea, we note that the parasitic instabilities studied by Goodman & Xu (1994) in general have a horizontal size larger than the vertical size of the channel solution on which they grow. In a box of unit aspect ratio containing a single channel they may simply be stable.

Whatever the reason for the differences, clearly, the obvious conclusion is that studies of turbulence driven by the MRI with net flux should be conducted in shearing boxes sufficiently large to allow the solution to develop naturally. The exact domain size depends on the particular problem and the quantities under consideration. However, it seems clear that boxes with aspect ratios close to unity over-emphasize the role of the channel solutions and may lead to misleading results.

On the other hand our discussion does not put in doubt the fact that MRI may be adequate to produce the turbulence necessary to support the required ``viscosity'' to transport angular momentum. However, the uncertainties discussed by many authors about the applicability of the shearing box approximation for the astrophysical problem of accretion disks suggest that the extension of the shearing box results to the full disk has to be tested on global simulations with sufficient resolution, that may soon be feasible.

Aspect ratio dependence in magnetorotational instability shearing box simulations

1 Introduction

2 Formulation

3 Results

4 Conclusion

References